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

name = "lrslib"
version = v"7.0"
versionstr = "070"

# Collection of sources required to build lrslibBuilder
sources = [
    "http://cgm.cs.mcgill.ca/~avis/C/lrslib/archive/lrslib-$versionstr.tar.gz" =>
    "e9f12b29be89b3ad6984f3a9cf83f5776ac06edc57b0716649e63395e5ac4dfe",
]

# Bash recipe for building across all platforms
script = raw"""
cd lrslib-070
extraargs=""
cflags="-O3 -Wall"

# 32bit linux, arm and windows:
if [[ $target == i686* ]] || [[ $target == arm* ]]; then
# no 128 bit ... patch makefile
sed -i -e 's#lrs-shared: ${SHLINK} lrs.o#lrs-shared: ${SHLINK} lrs64.o#' makefile
sed -i -e 's#$(CC) lrs.o -o $@ -L . -llrs#$(CC) lrs64.o -o $@ -L . -llrs#' makefile
sed -i -e 's#redund-shared: ${SHLINK}  redund.o#redund-shared: ${SHLINK}  redund64.o#' makefile
sed -i -e 's#$(CC) redund.o -o $@ -L . -llrs#$(CC) redund64.o -o $@ -L . -llrs#' makefile
sed -i -e 's#SHLIBOBJ=lrslong1-shr.o lrslong2-shr.o lrslib1-shr.o lrslib2-shr.o#SHLIBOBJ=lrslong1-shr.o lrslib1-shr.o#' makefile
fi

if [[ $target == *apple* ]]; then
sed -i -e 's#-Wl,-soname=#-install_name #' makefile
extraargs="SONAME=liblrs.0.dylib SHLINK=liblrs.dylib SHLIB=liblrs.0.0.0.dylib"
elif [[ $target == *freebsd* ]]; then
export CC="$CC $LDFLAGS"
elif [[ $target == *mingw* ]]; then
extraargs="SONAME=liblrs.0.dll SHLINK=liblrs.dll SHLIB=liblrs.0.0.0.dll"
cflags="$cflags -DSIGNALS -DTIMES"
fi

make prefix=$prefix INCLUDEDIR=$prefix/include LIBDIR=$prefix/lib CFLAGS="$cflags" $extraargs install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "liblrs", :liblrs)
    ExecutableProduct(prefix, "lrs", :lrs)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: benlorenz/lrslibBuilder

max_issues_repo_name: benlorenz/lrslibBuilder

max_forks_repo_name: benlorenz/lrslibBuilder


text:
using Oceananigans:
    RegularCartesianGrid,
    Field, CellField, FaceField, FaceFieldX, FaceFieldY, FaceFieldZ, EdgeField,
    VelocityFields, TracerFields, PressureFields, SourceTerms

# Increment and decrement integer a with periodic wrapping. So if n == 10 then
# incmod1(11, n) = 1 and decmod1(0, n) = 10.
@inline incmod1(a, n) = ifelse(a==n, 1, a + 1)
@inline decmod1(a, n) = ifelse(a==1, n, a - 1)

@inline δx_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[decmod1(i, g.Nx), j, k]
@inline δx_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[incmod1(i, g.Nx), j, k] - f[i, j, k]
@inline δx_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[incmod1(i, g.Nx), j, k] - f[i, j, k]
@inline δx_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[decmod1(i, g.Nx), j, k]

@inline δy_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, decmod1(j, g.Ny), k]
@inline δy_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, incmod1(j, g.Ny), k] - f[i, j, k]
@inline δy_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, incmod1(j, g.Ny), k] - f[i, j, k]
@inline δy_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k] - f[i, decmod1(j, g.Ny), k]

@inline function δz_c2f(g::RegularCartesianGrid, f, i, j, k)
    if k == 1
        return 0
    else
        @inbounds return f[i, j, k-1] - f[i, j, k]
    end
end

@inline function δz_f2c(g::RegularCartesianGrid, f, i, j, k)
    if k == g.Nz
        @inbounds return f[i, j, g.Nz]
    else
        @inbounds return f[i, j, k] - f[i, j, k+1]
    end
end

@inline function δz_e2f(g::RegularCartesianGrid, f, i, j, k)
    if k == g.Nz
        @inbounds return f[i, j, g.Nz]
    else
        @inbounds return f[i, j, k] - f[i, j, k+1]
    end
end

@inline function δz_f2e(g::RegularCartesianGrid, f, i, j, k)
    if k == 1
        return 0
    else
        @inbounds return f[i, j, k-1] - f[i, j, k]
    end
end

@inline avgx_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[i, j, k] + f[decmod1(i, g.Nx), j, k])
@inline avgx_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[incmod1(i, g.Nx), j, k] + f[i, j, k])
@inline avgx_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[i, j, k] + f[decmod1(i, g.Nx), j, k])

@inline avgy_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[i, j, k] + f[i, decmod1(j, g.Ny), k])
@inline avgy_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[i, incmod1(j, g.Ny), k] + f[i, j, k])
@inline avgy_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds 0.5 * (f[i, j, k] + f[i, decmod1(j, g.Ny), k])

@inline avg_xy(g::RegularCartesianGrid, u, i, j, k) = 0.5 * (avgy_f2c(g, u, i, j, k) + avgy_f2c(g, u, incmod1(i, g.Nx), j, k))

@inline function avgz_c2f(g::RegularCartesianGrid, f, i, j, k)
    if k == 1
        @inbounds return f[i, j, k]
    else
        @inbounds return  0.5 * (f[i, j, k] + f[i, j, k-1])
    end
end

@inline function avgz_f2c(g::RegularCartesianGrid, f, i, j, k)
    if k == g.Nz
        @inbounds return 0.5 * f[i, j, k]
    else
        @inbounds return 0.5 * (f[i, j, incmod1(k, g.Nz)] + f[i, j, k])
    end
end

@inline function avgz_f2e(g::RegularCartesianGrid, f, i, j, k)
    if k == 1
        @inbounds return f[i, j, k]
    else
        @inbounds return 0.5 * (f[i, j, k] + f[i, j, k-1])
    end
end

function avgx_4(f, Nx, i, j, k)
    @inbounds (f[i, j, k] + f[decmod1(i, Nx), j, k] -
		                   (f[incmod1(i, Nx), j, k] - f[i, j, k] -
                            f[decmod1(i, Nx), j, k] + f[decmod2(i, Nx), j, k]) / 6.0) * 0.5
end

function avgy_4(f, Ny, i, j, k)
    @inbounds (f[i, j, k] + f[i, decmod1(j, Ny), k] -
		                   (f[i, incmod1(j, Ny), k] - f[i, j, k] -
                            f[i, decmod1(j, Ny), k] + f[i, decmod2(j, Ny), k]) / 6.0) * 0.5
end

function avgz_4(f, Nz, i, j, k)
	if k == 1
		@inbounds return f[i, j, 1]
	else
		@inbounds return (f[i, j, k] + f[i, j, max(1, k-1)] -
		                              (f[i, j, min(Nz, k+1)] - f[i, j, k] -
							           f[i, j, max(1, k-1)] + f[i, j, max(1, k-2)]) / 6.0 ) * 0.5
    end
    nothing
end

@inline function div_f2c(g::RegularCartesianGrid, fx, fy, fz, i, j, k)
    (δx_f2c(g, fx, i, j, k) / g.Δx) + (δy_f2c(g, fy, i, j, k) / g.Δy) + (δz_f2c(g, fz, i, j, k) / g.Δz)
end

@inline function div_c2f(g::RegularCartesianGrid, fx, fy, fz, i, j, k)
    (δx_c2f(g, fx, i, j, k) / g.Δx) + (δy_c2f(g, fy, i, j, k) / g.Δy) + (δz_c2f(g, fz, i, j, k) / g.Δz)
end

@inline function δx_f2c_ab̄ˣ(g::RegularCartesianGrid, a, b, i, j, k)
    @inbounds (a[incmod1(i, g.Nx), j, k] * avgx_c2f(g, b, incmod1(i, g.Nx), j, k) -
               a[i,                j, k] * avgx_c2f(g, b, i,                j, k))
end

@inline function δy_f2c_ab̄ʸ(g::RegularCartesianGrid, a, b, i, j, k)
    @inbounds (a[i, incmod1(j, g.Ny), k] * avgy_c2f(g, b, i, incmod1(j, g.Ny), k) -
               a[i,                j, k] * avgy_c2f(g, b, i, j,                k))
end

@inline function δz_f2c_ab̄ᶻ(g::RegularCartesianGrid, a, b, i, j, k)
    if k == g.Nz
        @inbounds return a[i, j, k] * avgz_c2f(g, b, i, j, k)
    else
        @inbounds return (a[i, j,   k] * avgz_c2f(g, b, i, j,   k) -
                          a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1))
    end
end

@inline function div_flux(g::RegularCartesianGrid, u, v, w, Q, i, j, k)
    if k == 1
        @inbounds return (δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) - ((w[i, j, 2] * avgz_c2f(g, Q, i, j, 2)) / g.Δz)
    else
        return (δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)
    end
end

@inline function δx_c2f_ūˣūˣ(g::RegularCartesianGrid, u, i, j, k)
    avgx_f2c(g, u, i, j, k)^2 - avgx_f2c(g, u, decmod1(i, g.Nx), j, k)^2
end

@inline function δy_e2f_v̄ˣūʸ(g::RegularCartesianGrid, u, v, i, j, k)
    avgx_f2e(g, v, i, incmod1(j, g.Ny), k) * avgy_f2e(g, u, i, incmod1(j, g.Ny), k) -
    avgx_f2e(g, v, i,                j, k) * avgy_f2e(g, u, i,                j, k)
end

@inline function δz_e2f_w̄ˣūᶻ(g::RegularCartesianGrid, u, w, i, j, k)
    if k == g.Nz
        @inbounds return avgx_f2e(g, w, i, j, k) * avgz_f2e(g, u, i, j, k)
    else
        @inbounds return avgx_f2e(g, w, i, j,   k) * avgz_f2e(g, u, i, j,   k) -
                         avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1)
    end
end

@inline function u∇u(g::RegularCartesianGrid, u, v, w, i, j, k)
    (δx_c2f_ūˣūˣ(g, u, i, j, k) / g.Δx) + (δy_e2f_v̄ˣūʸ(g, u, v, i, j, k) / g.Δy) + (δz_e2f_w̄ˣūᶻ(g, u, w, i, j, k) / g.Δz)
end

@inline function δx_e2f_ūʸv̄ˣ(g::RegularCartesianGrid, u, v, i, j, k)
    avgy_f2e(g, u, incmod1(i, g.Nx), j, k) * avgx_f2e(g, v, incmod1(i, g.Nx), j, k) -
    avgy_f2e(g, u, i,                j, k) * avgx_f2e(g, v, i,                j, k)
end

@inline function δy_c2f_v̄ʸv̄ʸ(g::RegularCartesianGrid, v, i, j, k)
    avgy_f2c(g, v, i, j, k)^2 - avgy_f2c(g, v, i, decmod1(j, g.Ny), k)^2
end

@inline function δz_e2f_w̄ʸv̄ᶻ(g::RegularCartesianGrid, v, w, i, j, k)
    if k == g.Nz
        @inbounds return avgy_f2e(g, w, i, j, k) * avgz_f2e(g, v, i, j, k)
    else
        @inbounds return avgy_f2e(g, w, i, j,   k) * avgz_f2e(g, v, i, j,   k) -
                         avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1)
    end
end

@inline function u∇v(g::RegularCartesianGrid, u, v, w, i, j, k)
    (δx_e2f_ūʸv̄ˣ(g, u, v, i, j, k) / g.Δx) + (δy_c2f_v̄ʸv̄ʸ(g, v, i, j, k) / g.Δy) + (δz_e2f_w̄ʸv̄ᶻ(g, v, w, i, j, k) / g.Δz)
end

@inline function δx_e2f_ūᶻw̄ˣ(g::RegularCartesianGrid, u, w, i, j, k)
    avgz_f2e(g, u, incmod1(i, g.Nx), j, k) * avgx_f2e(g, w, incmod1(i, g.Nx), j, k) -
    avgz_f2e(g, u, i,                j, k) * avgx_f2e(g, w, i,                j, k)
end

@inline function δy_e2f_v̄ᶻw̄ʸ(g::RegularCartesianGrid, v, w, i, j, k)
    avgz_f2e(g, v, i, incmod1(j, g.Ny), k) * avgy_f2e(g, w, i, incmod1(j, g.Ny), k) -
    avgz_f2e(g, v, i,                j, k) * avgy_f2e(g, w, i,                j, k)
end

@inline function δz_c2f_w̄ᶻw̄ᶻ(g::RegularCartesianGrid, w, i, j, k)
    if k == 1
        return 0
    else
        return avgz_f2c(g, w, i, j, k-1)^2 - avgz_f2c(g, w, i, j, k)^2
    end
end

@inline function u∇w(g::RegularCartesianGrid, u, v, w, i, j, k)
    (δx_e2f_ūᶻw̄ˣ(g, u, w, i, j, k) / g.Δx) + (δy_e2f_v̄ᶻw̄ʸ(g, v, w, i, j, k) / g.Δy) + (δz_c2f_w̄ᶻw̄ᶻ(g, w, i, j, k) / g.Δz)
end

@inline δx²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δx_c2f(g, f, incmod1(i, g.Nx), j, k) - δx_c2f(g, f, i, j, k)
@inline δy²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δy_c2f(g, f, i, incmod1(j, g.Ny), k) - δy_c2f(g, f, i, j, k)

@inline function δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k)
    if k == g.Nz
        return δz_c2f(g, f, i, j, k)
    else
        return δz_c2f(g, f, i, j, k) - δz_c2f(g, f, i, j, k+1)
    end
end

@inline function κ∇²(g::RegularCartesianGrid, Q, κh, κv, i, j, k)
    ((κh/g.Δx^2) * δx²_c2f2c(g, Q, i, j, k)) + ((κh/g.Δy^2) * δy²_c2f2c(g, Q, i, j, k)) + ((κv/g.Δz^2) * δz²_c2f2c(g, Q, i, j, k))
end

@inline δx²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2c(g, f, i, j, k) - δx_f2c(g, f, decmod1(i, g.Nx), j, k)
@inline δy²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2c(g, f, i, j, k) - δy_f2c(g, f, i, decmod1(j, g.Ny), k)

@inline δx²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2e(g, f, incmod1(i, g.Nx), j, k) - δx_f2e(g, f, i, j, k)
@inline δy²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2e(g, f, i, incmod1(j, g.Ny), k) - δy_f2e(g, f, i, j, k)

@inline function δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k)
    if k == g.Nz
        return δz_f2e(g, f, i, j, k)
    else
        return δz_f2e(g, f, i, j, k) - δz_f2e(g, f, i, j, k+1)
    end
end

@inline function δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k)
    if k == 1
        return 0
    else
        return δz_f2c(g, f, i, j, k-1) - δz_f2c(g, f, i, j, k)
    end
end

@inline function 𝜈∇²u(g::RegularCartesianGrid, u, 𝜈h, 𝜈v, i, j, k)
    ((𝜈h/g.Δx^2) * δx²_f2c2f(g, u, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, u, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, u, i, j, k))
end

@inline function 𝜈∇²v(g::RegularCartesianGrid, v, 𝜈h, 𝜈v, i, j, k)
    ((𝜈h/g.Δx^2) * δx²_f2e2f(g, v, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2c2f(g, v, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, v, i, j, k))
end

@inline function 𝜈∇²w(g::RegularCartesianGrid, w, 𝜈h, 𝜈v, i, j, k)
    ((𝜈h/g.Δx^2) * δx²_f2e2f(g, w, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, w, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2c2f(g, w, i, j, k))
end

@inline function ∇²_ppn(g::RegularCartesianGrid, f, i, j, k)
	(δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)
end

ext: jl


text:
export Microphone

struct Microphone end

function show(io::IO, ::MIME"text/html", microphone::Microphone)
    mic_id = randstring('a':'z')
    mic_btn_id = randstring('a':'z')
    microphone
    withtag(() -> (), io, :audio, :id => mic_id)
    print(io, """<input type="button" id="$mic_btn_id" class="mic-button" value="Stop">""")
    withtag(io, :script) do
        print(io, """
            const player = document.getElementById('$mic_id');
            const stop = document.getElementById('$mic_btn_id');

            const handleSuccess = function(stream) {
            const context = new AudioContext({ sampleRate: 44100 });
            const analyser = context.createAnalyser();
            const source = context.createMediaStreamSource(stream);

            source.connect(analyser);

            const bufferLength = analyser.frequencyBinCount;

            let dataArray = new Float32Array(bufferLength);
            let animFrame;

            const streamAudio = () => {
                animFrame = requestAnimationFrame(streamAudio);
                analyser.getFloatTimeDomainData(dataArray);
                player.value = dataArray;
                player.dispatchEvent(new CustomEvent("input"));
            }

            streamAudio();

            stop.onclick = e => {
                source.disconnect(analyser);
                cancelAnimationFrame(animFrame);
            }
            }

            navigator.mediaDevices.getUserMedia({ audio: true, video: false })
            .then(handleSuccess)
        """
        )
    end
    withtag(io, :style) do
        print(io, """
            .mic-button {
                background-color: darkred;
                border: none;
                border-radius: 6px;
                color: white;
                padding: 15px 32px;
                text-align: center;
                text-decoration: none;
                display: inline-block;
                font-size: 16px;
                font-family: "Alegreya Sans", sans-serif;
                margin: 4px 2px;
                cursor: pointer;
            }
        """
        )
    end
end

get(microphone::Microphone) = microphone


text:
"""
    print_tree(tree; kwargs...)
    print_tree(io::IO, tree; kwargs...)
    print_tree(f::Function, io::IO, tree; kwargs...)

Print a text representation of `tree` to the given `io` object.

# Arguments

* `f::Function` - custom implementation of [`printnode`](@ref) to use. Should have the
  signature `f(io::IO, node)`.
* `io::IO` - IO stream to write to.
* `tree` - tree to print.
* `maxdepth::Integer = 5` - truncate printing of subtrees at this depth.
* `indicate_truncation::Bool = true` - print a vertical ellipsis character beneath
  truncated nodes.
* `charset::TreeCharSet` - [`TreeCharSet`](@ref) to use to print branches.
* `printkeys::Union{Bool, Nothing}` - Whether to print keys of child nodes (using
  `pairs(children(node))`). A value of `nothing` uses [`printkeys_default`](@ref) do decide the
  behavior on a node-by-node basis.

# Examples

```julia
julia> tree = [1:3, "foo", [[[4, 5], 6, 7], 8]];

julia> print_tree(tree)
Vector{Any}
├─ UnitRange{Int64}
│  ├─ 1
│  ├─ 2
│  └─ 3
├─ "foo"
└─ Vector{Any}
   ├─ Vector{Any}
   │  ├─ Vector{Int64}
   │  │  ├─ 4
   │  │  └─ 5
   │  ├─ 6
   │  └─ 7
   └─ 8

julia> print_tree(tree, maxdepth=2)
Vector{Any}
├─ UnitRange{Int64}
│  ├─ 1
│  ├─ 2
│  └─ 3
├─ "foo"
└─ Vector{Any}
   ├─ Vector{Any}
   │  ⋮
   │
   └─ 8

julia> print_tree(tree, charset=AbstractTrees.ASCII_CHARSET)
Vector{Any}
+-- UnitRange{Int64}
|   +-- 1
|   +-- 2
|   \\-- 3
+-- "foo"
\\-- Vector{Any}
    +-- Vector{Any}
    |   +-- Vector{Int64}
    |   |   +-- 4
    |   |   \\-- 5
    |   +-- 6
    |   \\-- 7
    \\-- 8
```

"""
function print_tree end

"""
    printnode(io::IO, node)

Print a compact representation of a single node.  By default, this prints `nodevalue(node)`.

**OPTIONAL**: This can be extended for custom types and controls how nodes are shown
in [`print_tree`](@ref).
"""
printnode(io::IO, node) = show(IOContext(io, :compact => true, :limit => true), nodevalue(node))

"""
    repr_node(node; context=nothing)

Get the string representation of a node using [`printnode`](@ref). This works
analagously to `Base.repr`.

`context` is an `IO` or `IOContext` object whose attributes are used for the
I/O stream passed to `printnode`.
"""
function repr_node(node; context=nothing)
    buf = IOBuffer()
    io = context === nothing ? buf : IOContext(buf, context)
    printnode(io, node)
    return String(take!(buf))
end

const CharArg = Union{AbstractString, Char}

"""
    TreeCharSet(mid, terminator, skip, dash, trunc, pair)

Set of characters (or strings) used to pretty-print tree branches in [`print_tree`](@ref).

## Fields
- `mid`: "Forked" branch segment connecting to middle children.
- `terminator`: Final branch segment connecting to last child.
- `skip`: Vertical branch segment.
- `dash`: Horizontal branch segmentt printed to the right of `mid` and `terminator`.
- `trunc`: Used to indicate the subtree has been truncated at the maximum depth.
- `pair`: Printed between a child node and its key.
"""
struct TreeCharSet
    mid::String
    terminator::String
    skip::String
    dash::String
    trunc::String
    pair::String

    function TreeCharSet(mid::CharArg, terminator::CharArg, skip::CharArg, dash::CharArg, trunc::CharArg, pair::CharArg)
        return new(String(mid), String(terminator), String(skip), String(dash), String(trunc), String(pair))
    end
end

"""
    TreeCharSet(base::TreeCharSet; fields...)

Create a new `TreeCharSet` by modifying select fields of an existing instance.
"""
function TreeCharSet(base::TreeCharSet;
                     mid = base.mid,
                     terminator = base.terminator,
                     skip = base.skip,
                     dash = base.dash,
                     trunc = base.trunc,
                     pair = base.pair,
                    )
    return TreeCharSet(mid, terminator, skip, dash, trunc, pair)
end

"""
    TreeCharSet(name=:unicode)

Generate one of the default tree character sets.  Valid options are `:unicode` (default) and `:ascii`.
"""
function TreeCharSet(name::Symbol=:unicode)
    if name == :unicode
        TreeCharSet("├", "└", "│", "─", "⋮", " ⇒ ")
    elseif name == :ascii
        TreeCharSet("+", "\\", "|", "--", "...", " => ")
    else
        throw(ArgumentError("unrecognized dfeault TreeCharSet name: $name"))
    end
end

"""
    shouldprintkeys(children)::Bool

Whether a collection of children should be printed with its keys by default.

The base behavior is to print keys for all collections for which `keys()` is defined, with the
exception of `AbstractVector`s and tuples.
"""
shouldprintkeys(ch) = applicable(keys, ch)
shouldprintkeys(::AbstractVector) = false
shouldprintkeys(::Tuple) = false
shouldprintkeys(::Base.Generator) = false

"""
    print_child_key(io::IO, key)

Print the key for a child node.
"""
print_child_key(io::IO, key) = show(io, key)
print_child_key(io::IO, key::CartesianIndex) = show(io, Tuple(key))

branchwidth(cs::TreeCharSet) = sum(textwidth.((cs.mid, cs.dash)))

function print_tree(printnode::Function, io::IO, node;
                    maxdepth::Integer=5,
                    indicate_truncation::Bool=true,
                    charset::TreeCharSet=TreeCharSet(),
                    printkeys::Union{Bool,Nothing}=nothing,
                    depth::Integer=0,
                    prefix::AbstractString="",
                   )
    # Get node representation as string
    str = repr_node(node, context=io)

    # Copy buffer to output, prepending prefix to each line
    for (i, line) in enumerate(split(str, '\n'))
        i ≠ 1 && print(io, prefix)
        println(io, line)
    end

    # Node children
    c = children(node)

    # No children?
    isempty(c) && return

    # Reached max depth?
    if depth ≥ maxdepth
        # Print truncation char(s)
        if indicate_truncation
            println(io, prefix, charset.trunc)
            println(io, prefix)
        end
        return
    end

    # Print keys?
    this_printkeys = applicable(keys, c) && (isnothing(printkeys) ? shouldprintkeys(c) : printkeys)

    # Print children
    s = Iterators.Stateful(this_printkeys ? pairs(c) : c)

    while !isempty(s)
        child_prefix = prefix

        if this_printkeys
            child_key, child = popfirst!(s)
        else
            child = popfirst!(s)
            child_key = nothing
        end

        print(io, prefix)

        # Last child?
        if isempty(s)
            print(io, charset.terminator)
            child_prefix *= " " ^ (textwidth(charset.skip) + textwidth(charset.dash) + 1)
        else
            print(io, charset.mid)
            child_prefix *= charset.skip * " " ^ (textwidth(charset.dash) + 1)
        end

        print(io, charset.dash, ' ')

        # Print key
        if this_printkeys
            buf = IOBuffer()
            print_child_key(IOContext(buf, io), child_key)
            key_str = String(take!(buf))

            print(io, key_str, charset.pair)

            child_prefix *= " " ^ (textwidth(key_str) + textwidth(charset.pair))
        end

        print_tree(printnode, io, child;
                   maxdepth, indicate_truncation, charset, printkeys,
                   depth=depth+1, prefix=child_prefix
                  )
    end
end

print_tree(io::IO, node; kw...) = print_tree(printnode, io, node; kw...)
print_tree(node; kw...) = print_tree(stdout, node; kw...)

"""
    repr_tree(tree; context=nothing, kw...)

Get the string result of calling [`print_tree`](@ref) with the supplied arguments.

The `context` argument works as it does in `Base.repr`.
"""
function repr_tree(tree; context=nothing, kw...)
    buf = IOBuffer()
    io = context === nothing ? buf : IOContext(buf, context)
    print_tree(io, tree; kw...)
    return String(take!(buf))
end


text:
##Univariate kerneldensity and regression
using Distributions
x=rand(Normal(10), 500)
xeval=linspace(minimum(x), maximum(x), 100)
h = bwlscv(x, gaussiankernel)
@test h>0
denvalues=kerneldensity(x, xeval=xeval)
@test all(denvalues .>= 0)

y=2 .* x.^2 .+ rand(Normal(0,5),500)
b0, b1 = linreg(x, y)
regfit = b0 .+ b1 .* x
yfit0=localconstant(x, y, xeval=x)
@test sumabs2(regfit .- y) > sumabs2(yfit0.-y)

yfit1=locallinear(x, y, xeval=x)
@test sumabs2(regfit .- y) > sumabs2(yfit1.-y)

yfit0=npr(x, y, xeval=x, reg=localconstant)
@test sumabs2(regfit .- y) > sumabs2(yfit0.-y)

yfit1=npr(x, y, xeval=x, reg=locallinear)
@test sumabs2(regfit .- y) > sumabs2(yfit1.-y)

#confidence band
yfit1=npr(x, y, xeval=xeval, reg=locallinear)
cb=bootstrapCB(x, y, xeval=xeval)
@test mean(vec(cb[1,:]) .<= yfit1 .<= vec(cb[2,:])) > .8

###Bounded gamma kernel density and regression
x = rand(Gamma(4,2), 500)
xeval = linspace(0.01,20, 100)
h = bwlscv(x, gammakernel)
@test h>0
denvalues = kerneldensity(x, xeval=xeval, kernel=gammakernel, lb=0.0)
@test all(denvalues .> 0)

y=2 .* x.^2 + x.*rand(Normal(0, 5), 500)
b0, b1 = linreg(x, y)
regfit = b0 .+ b1 .* x
yfit0=npr(x, y, xeval=x, reg=localconstant, kernel=gammakernel,lb=0.0)
@test sumabs2(regfit .- y) > sumabs2(yfit0.-y)

yfit1=npr(x, y, xeval=x, reg=locallinear, kernel=gammakernel, lb=0.0)
@test sumabs2(regfit .- y) > sumabs2(yfit1.-y)

#bounded beta kernel density and regression
x = rand(Beta(4,2), 500) * 10
xeval = linspace(0, 10, 100)
h = bwlscv(x./10, betakernel)
@test h>0
denvalues=kerneldensity(x, xeval=xeval, kernel=betakernel,h=h, lb=0.0,ub=10.0)
@test all(denvalues .> 0)

##Bounded local regression
y=2 .* x.^2 + x.*rand(Normal(0, 5), 500)
b0, b1 = linreg(x, y)
regfit = b0 .+ b1 .* x
yfit0=npr(x, y, xeval=x, reg=localconstant, kernel=betakernel,lb=0.0, ub=10.0)
@test sumabs2(regfit .- y) > sumabs2(yfit0.-y)

yfit1=npr(x, y, xeval=x, reg=locallinear, kernel=betakernel, lb=0.0,ub=10.0)
@test sumabs2(regfit .- y) > sumabs2(yfit1.-y)

yfit1=npr(x, y, xeval=xeval, reg=locallinear, kernel=betakernel, lb=0.0,ub=10.0)
cb=bootstrapCB(x, y, xeval=xeval,reg=locallinear, kernel=betakernel, lb=0.0,ub=10.0)
@test mean(vec(cb[1,:]) .<= yfit1 .<= vec(cb[2,:])) > .8


text:
#
# Optional time-consuming Benchmarks
#

using DataFrames
using Benchmark

benchmarks = [ "io.jl"]

# TODO: Print summary to stdout_stream, while printing results
#       to file with appends.
#println("Running benchmarks:")

for benchmark in benchmarks
#    println(" * $(benchmark)")
    include(benchmark)
end


text:
run(`julia julia/julia.jl`) # logs to file

# move to rust directory and bench
cd("rust/src")
run(`cargo +nightly bench`) # logs to console

# back to home directory
cd("...")

cd("python")

max_stars_repo_name: JeannotJeannot/Learn

max_issues_repo_name: JeannotJeannot/Learn

max_forks_repo_name: JeannotJeannot/Learn


text:
using Test

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

# Iterations:
# linear5 : > 1000, < 2000
# linear9 : > 3000, < 4000
# linear15: > 20000, Don't know if ever converges so we exclude it
import CDCS
const OPTIMIZER_CONSTRUCTOR = MOI.OptimizerWithAttributes(CDCS.Optimizer, MOI.Silent() => true, "maxIter" => 4000)
const OPTIMIZER = MOI.instantiate(OPTIMIZER_CONSTRUCTOR)

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "CDCS"
end

@testset "supports_allocate_load" begin
    @test MOIU.supports_allocate_load(OPTIMIZER, false)
    @test !MOIU.supports_allocate_load(OPTIMIZER, true)
end

const BRIDGED = MOI.instantiate(OPTIMIZER_CONSTRUCTOR, with_bridge_type=Float64)
const CONFIG = MOIT.TestConfig(atol=3e-2, rtol=3e-2)

@testset "Unit" begin
    MOIT.unittest(BRIDGED, CONFIG, [
        # `NumberOfThreads` not supported.
        "number_threads",
        # `TimeLimitSec` not supported.
        "time_limit_sec",
        # Need to investigate...
        "solve_with_lowerbound", "solve_affine_deletion_edge_cases", "solve_blank_obj", "solve_single_variable_dual_min", "solve_single_variable_dual_max",
        # Need https://github.com/JuliaOpt/MathOptInterface.jl/issues/529
        "solve_qp_edge_cases",
        # Error using cdcs_hsde.preprocess (line 14)
        # No variables in your problem?
        "solve_unbounded_model",
        # Integer and ZeroOne sets are not supported
        "solve_integer_edge_cases", "solve_objbound_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end

@testset "Continuous linear problems" begin
    MOIT.contlineartest(BRIDGED, CONFIG, [
        # Need to investigate...
        "linear12", "linear15"])
end

@testset "Continuous conic problems" begin
    MOIT.contconictest(BRIDGED, CONFIG, [
        # ITERATION_LIMIT
        "normone1v", "normone1f",
        # rotatedsoc2: Returns Inf and -Inf instead of infeasibility certificate
        "rotatedsoc2",
        # Need to investigate...
        "psdt3", "psds3",
        # Unsupported cones
        "pow", "dualpow", "rootdets", "exp", "dualexp", "logdet", "normspec", "normnuc", "relentr"])
end


text:
# -*- coding: utf-8 -*-

# modularize so I can 'reload' in repl
using Base.Test
# reload("Nodes.jl")  # for repl debugging reload things in the order they would call each other
using Parsimonious
import Parsimonious.DEFUNCT_visit  # for overloading
# import Nodes: errstring # for testing internal stuff

type HtmlFormatter <: NodeVisitor end

#
# class HtmlFormatter(NodeVisitor):
#     """Visitor that turns a parse tree into HTML fragments"""
#
#     def visit_bold_open(self, node, visited_children):
#         return ''

DEFUNCT_visit(::HtmlFormatter, ::LeafNode{:bold_open}) = ""

#     def visit_bold_close(self, node, visited_children):
#         return ''
#

DEFUNCT_visit(::HtmlFormatter, ::LeafNode{:bold_close}) = ""

#     def visit_text(self, node, visited_children):
#         """Return the text verbatim."""
#         return node.text

DEFUNCT_visit(::HtmlFormatter, node::LeafNode{:text}) = nodetext(node)

#
#     def visit_bold_text(self, node, visited_children):
#         return ''.join(visited_children)

DEFUNCT_visit(::HtmlFormatter, ::ParentNode{:bold_text}, visited_children) = join(visited_children, "")

type ExplosiveFormatter <: NodeVisitor end

#
#
# class ExplosiveFormatter(NodeVisitor):
#     """Visitor which raises exceptions"""
#
#     def visit_boom(self, node, visited_children):
#         raise ValueError

# BoundsError is supposed to get wrapped with VisitationError
DEFUNCT_visit(::ExplosiveFormatter, ::LeafNode{:boom}) = throw(BoundsError)

#
#
# def test_visitor():
#     """Assert a tree gets visited correctly.
#
#     We start with a tree from applying this grammar... ::
#
#         bold_text  = bold_open text bold_close
#         text       = ~'[a-zA-Z 0-9]*'
#         bold_open  = '(('
#         bold_close = '))'
#
#     ...to this text::
#
#         ((o hai))
#
#     """
#     text = '((o hai))'
#     tree = Node('bold_text', text, 0, 9,
#                 [Node('bold_open', text, 0, 2),
#                  Node('text', text, 2, 7),
#                  Node('bold_close', text, 7, 9)])
#     result = HtmlFormatter().visit(tree)
#     eq_(result, 'o hai')

text = "((o hai))"
tree = Node("bold_text", text, 1, 9,
            (Node("bold_open", text, 1, 2),
              Node("text", text, 3, 7),
              Node("bold_close", text, 8, 9)))
result = DEFUNCT_visit(HtmlFormatter(), tree)
@test result == "o hai"

# def test_visitation_exception():
#     assert_raises(VisitationError,
#                   ExplosiveFormatter().visit,
#                   Node('boom', '', 0, 0))

n = Node("boom", "", 1, 0)
@test_throws DEFUNCT_visit(ExplosiveFormatter(), n)
try
    DEFUNCT_visit(ExplosiveFormatter(), n)
catch ex
    if isa(ex, VisitationError)
        @test true
    else
        println("Rethrowing in test_nodes.jl which is bad")
# TODO: Hard to fix under the new node Union paradigm
#        rethrow(ex)
    end
end

#
#
# def test_str():

txt = "o hai"
# stringexample = "<Node called 'text' matching 'o hai'>\n"
stringexample = "1 LeafNode{:text}                                                       'o hai'\n"
n = Node("text", txt, 1, 5)
@show stringexample
@show string(n)
@test string(n) == stringexample

#
# def test_repr():
#     """Test repr of ``Node``."""
#     s = u'hai ö'
#     boogie = u'böogie'
#     n = Node(boogie, s, 0, 3, children=[
#             Node('', s, 3, 4), Node('', s, 4, 5)])
#     eq_(repr(n), """s = {hai_o}\nNode({boogie}, s, 0, 3, children=[Node('', s, 3, 4), Node('', s, 4, 5)])""".format(hai_o=repr(s), boogie=repr(boogie)))

s = "hai ö"
boogie = "böogie"
n = Node(boogie, s, 1, 3, (
        Node("", s, 4, 3), Node("", s, 5, 4)))
shouldbe = "\"s = \\\"hai ö\\\" ; ParentNode{:böogie}(1, 3, [Node{:}(4, 3), Node{:}(5, 4)])\""
# TODO: who gives a fuck?
#@test repr(n) == shouldbe

## More test I wrote ##

# start and end must be within the text
# TODO: put back assertions in constructors
# @test_throws Node("myexpr", "123456", 0, 2)
# @test_throws Node("myexpr", "123456", 2, 7)
@test isa(Node("myexpr", "123456", 1, 6), AnyNode)

# _end may be less than start to indicate a 0-length capture
@test isa(Node("myexpr", "123456", 1, 0), AnyNode)

# Some nodes for testing

#         123456789012345
mytext = "this is my text"
copytext = string(mytext)
n = Node("myexpr", mytext, 5, 8)
n2 = Node("myexpr", mytext, 5, 8)
n3 = Node("myexpr2", mytext, 5, 8)
nct = Node("myexpr", copytext, 5, 8)

# length(n) is the number of children n has
@test length(n) == 0

@test isa(n, AnyNode)
@test nodetext(n) == " is "

# n == nct even though text is a copy, it is an exact copy
@test !is(mytext, copytext)
@test n == nct

@test isempty(Node())
@test !isempty(n)
@test Node() == Node()  # Because immutable

@test isequal(n, n)
@test isequal(n, n2)
@test !isequal(n, n3)
@test isequal(Node(), Node())

# add some children to n2
# Thinking of making children immutable ... yeah, that should be fine
n2 = Node("myexpr", mytext, 5, 8, (n3, n3, n3, n3))

# no longer equal
@test !isequal(n, n2)

# length(n) is number of children
@test length(n2.children) == 4
@test length(n2) == 4

# Print an error string when pretty printing tree traceback
@test Parsimonious.errstring(n, n2) == ""
@test match(r"Error", Parsimonious.errstring(n2, n2)) != nothing
@test match(r"Error", Parsimonious.errstring(n, n2)) == nothing

# Indent nodes when pretty printing
@test Parsimonious.indent(Parsimonious.indent("foo","|"),"|") == "||foo"
@test Parsimonious.indent("foo\nbar","  ") == "  foo\n  bar"

# TODO: test RegexNode

type SomeVisitor <: NodeVisitor end
@test_throws DEFUNCT_visit(SomeVisitor(), n2)
@test lift_child(SomeVisitor(), n2, ["foo"]) == "foo"

# keyword syntax
# many tests because I thought for a long time this can't work
# I strongly suspect it is slow
# keyword syntax is totally out
#@time begin for i in 1:1000
#    nwc = Node(T="withchildren", fulltext=mytext, start=1, _end=4, children=[n n2 n3])
#    @test length(nwc) == 3
#    nwd = Node(T="nochildren", fulltext=mytext, start=1, _end=4)
#    @test length(nwd) == 0
#    m=match(r"\S+", mytext)
#    nwe = Node(match=m, _end=4, start=1, fulltext=mytext, T="withchildren")
#    @test isa(nwe.match, RegexMatch)
#    nwf = Node("withchildren", mytext, 1, 0, children=[n n2 n3])
#    nwg = Node("withchildren", mytext, 1, 0, match=m, children=[n n2 n3])
#    nwh = Node("withchildren", mytext, 1, 0, match=m)
#end
#end

@test textlength(Node("foo","",1,0)) == 0
@test textlength(Node("foo","f", 1,1)) == 1
@test textlength(Node("foo","foo", 1,3)) == 3

txt = "asdf;lkj"
n1 = Node("string", txt, 1, 8)
n2 = Node("string", txt, 1, 8, (Node("alpha", txt, 1, 4), Node("sym", txt, 5, 5), Node("not", txt, 6, 5), Node("alpha", txt, 6, 8)))
for n in (n1, n2)
    print(n)
end

####################### TEST
@test name(n2) == "string"
####################### TEST

########### TEST
@test nodetext(n2) == "asdf;lkj"
########### TEST

########### TEST
@test isempty(n2) == false
@test isempty(EmptyNode()) == true
# @test isempty(nothing) == true
@test isempty(Node("asdf", "asdf", 1, 0)) == false
########### TEST

########### TEST
@test n2 == n2
@test Node("asdf", "texty", 1, 2) == Node("asdf", "texty", 1, 2)
@test Node("asdf", "texty", 1, 3) != Node("asdf", "texty", 1, 2)
@test Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 1),)) != Node("asdf", "texty", 1, 2)
@test Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 1),)) != Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 2),))
@test Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 2),)) == Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 2),))
@test Node("asdf", "texty", 1, 2, (Node("bsdf", "texty", 1, 2),)) != EmptyNode()
@test EmptyNode() == EmptyNode()
########### TEST

####################### TEST
showerror(STDOUT, VisitationError(n2, BoundsError()))
[println(n) for n in n2]  # iteration
####################### TEST

####################### TEST
type TestVisitor <: NodeVisitor end
DEFUNCT_visit(v::TestVisitor, n::ParentNode{:lit}, visited_children) = lift_child(v,n,visited_children)
DEFUNCT_visit(v::TestVisitor, n::LeafNode{:generic}, _) = nodetext(n)
@test Parsimonious.visit_on_the_way_down(TestVisitor(), Node("lit", "texty", 1, 2)) == []
@test DEFUNCT_visit(TestVisitor(), Node("generic", "asdf", 1, 2)) == "as"
@test nodetext(Node("generic", "asdf", 1, 2)) == "as"
@test DEFUNCT_visit(TestVisitor(), Node("generic", "asdf", 1, 2)) == "as"
@test DEFUNCT_visit(TestVisitor(), Node("lit", "asdf", 1, 2, (Node("generic", "asdf", 1, 2),))) == "as"
@test_throws DEFUNCT_visit(TestVisitor(), Node("notimplemented", "asdf", 1, 2)) == "as"
####################### TEST

s = "ελληνική"
nu = Node("", s, chr2ind(s, 1), chr2ind(s, 2))
@test pos(nu) == chr2ind(s, 1)
@test _end(nu) == chr2ind(s, 2)
@test endof(nu.match) == chr2ind(s, 2)
@test length(nu.match) == 2

nu2 = Node("",s, chr2ind(s, 2), chr2ind(s, 1))
@test pos(nu2) == chr2ind(s, 2)
@show _end(nu2)
@show pos(nu2)
@show chr2ind(s, 1)
xdump(nu2)
@show s[_end(nu2)]
@test _end(nu2) == chr2ind(s, 1)
@test length(nu2.match) == 0
@test endof(nu2.match) == 0
@show chr2ind(s, 2)
nu3 = Node("",s, chr2ind(s, 2), chr2ind(s, 2))
@test pos(nu2) == chr2ind(s, 2)
@test _end(nu3) == chr2ind(s, 2)
@test length(nu3.match) == 1
@test endof(nu3.match) == 1


text:
using SwitchTimeOpt

# Plotting settings
using PyPlot, PyCall
close("all")
# Import seaborn for nice plots
@pyimport seaborn as sns
sns.set_palette("hls", 8)
sns.set_context("paper", font_scale=1.5)
sns.set_style("whitegrid")
# Use Latex Labels in Plots
plt[:rc]("text", usetex=true)
plt[:rc]("font", family="serif")

maxiter = 20
using Ipopt
solver = IpoptSolver(
  print_level=0,
  max_iter=maxiter)

### Define system parameters
# Time vector
t0 = 0.0
tf = 12.0

# Size of the system
nx = 4

# Integer input
uvec = [repmat([0.0; 1.0], 4, 1); 0.0]'

# Number of switching times
N = size(uvec,2) - 1

# Cost funcction matrix
C = [1.0 0.0 -1.0 0.0;
     0.0 1.0 0.0  -1.0]
Q = C'*C

# Define initial state
x0 = [0.5; 0.7; 1; 1]

### Define system dynamics
function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
  n = length(x)
  f = zeros(n)
  f[1] = x[1] - x[1]*x[2] - 0.4*x[1]*u[1]
  f[2] = -x[2] + x[1]*x[2] - 0.2*x[2]*u[1]
  return f

end

function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
  df = [1.0-x[2]-0.4*u[1]       -x[1]                   0    0;
        x[2]                     -1+x[1]-0.2*u[1]       0    0;
        0                        0                      0    0;
        0                        0                      0    0]
end

### Generate and solve problems with different grid points
ngrid = [100; 150; 200; 250]

idx200 = (find(ngrid .== 200))[1]

# Preallocate vectors for results
objode45 = Array{Float64}(length(ngrid))
objlin = Array{Float64}(length(ngrid))
objiterates = Array{Float64}(maxiter+1, length(ngrid))
cputime = Array{Float64}(length(ngrid))
nobjeval = Array{Int}(length(ngrid))
ngradeval = Array{Int}(length(ngrid))
nhesseval = Array{Int}(length(ngrid))
tauopt = Array{Float64}(N, length(ngrid))
xsim = Array{Float64}(4, 10^4, length(ngrid))
xlinsim = Array{Float64}(4, 10^4, length(ngrid))
usim = Array{Float64}(1, 10^4, length(ngrid))

# Initialize the problem first
m = stoproblem(x0, nldyn, nldyn_deriv, uvec)

for i = 1:length(ngrid)  # Iterate over all grid points numbers

  m = stoproblem(
    x0,                 # Initial state
    nldyn,              # Nonlinear dynamics
    nldyn_deriv,        # Nonlinear dynamics derivative
    uvec,               # Vector of integer inputs
    ngrid = ngrid[i],   # Number of grid points
    t0 = t0,            # Initial time
    tf = tf,            # Final time
    Q = Q,              # Cost matrix
    solver = solver)

  # Solve optimization
  solve!(m)

  # Obtain values
  tauopt[:, i] = gettau(m)
  objlin[i] = getobjval(m)
  cputime[i] = getsoltime(m)
  nobjeval[i] = getnobjeval(m)
  ngradeval[i] = getngradeval(m)
  nhesseval[i] = getnhesseval(m)

  # Simulate system
  # Nonlinear simulation
  xsim[:, :, i], xpts, objode45[i], t = simulate(m, tauopt[:, i])
  usim[:, :, i], _ = simulateinput(m, t)

  # Linearized simulation
  xlinsim[:, :, i], _, Jlinsim, _ = simulatelinearized(m, tauopt[:, i], t)

  # Save objective function iterates
  objiterates[:, i] = m.STOev.obj[2:end]

end

###  Print results
@printf("RESULTS\n")
@printf("-------\n\n")

@printf("+==================================================================================+\n")
@printf("|  ngrid   |   objode45  |   objlin    | deltaobj [%%] |  nobjeval  |  cputime [s]  |\n")
@printf("+==================================================================================+\n")

for i = 1:length(ngrid)
  @printf("|  %7.i | %9.4f   | %9.4f   |   %6.3f     | %9.i  | %12.4f  | \n", ngrid[i], objode45[i], objlin[i], 100*norm(objode45[i]- objlin[i])/objode45[i], nobjeval[i], cputime[i])
  @printf("+----------------------------------------------------------------------------------+\n")
end

@printf("\nOptimal Switching Times for ngrid = 200\n")
@printf("--------------------------------------\n")
@printf("tauopt = "); show(round.(tauopt[:, idx200],3)); @printf("\n")

# Generate plots for ngrid = 25
t = linspace(t0, tf, 10000)

figure()
subplot(3,1,1)
plot(t, xlinsim[1,:, idx200], sns.xkcd_rgb["grass green"], linestyle = "dashdot")
plot(t, xsim[1,:, idx200], sns.xkcd_rgb["denim blue"])
ylim(0, 1.75)
yticks([0; 1; ])
ylabel(L"x_1")

subplot(3,1,2)
plot(t, xlinsim[2,:, idx200], sns.xkcd_rgb["grass green"], linestyle = "dashdot")
plot(t, xsim[2, :, idx200], sns.xkcd_rgb["denim blue"])
ylabel(L"x_2")
ylim(0, 1.75)
yticks([0; 1; ])

subplot(3,1,3)
plot(t, usim[1,:, idx200], sns.xkcd_rgb["denim blue"])
ylim(-0.2, 1.2)
yticks([0; 1])
ylabel(L"u")
xlabel(L"$\mathrm{Time}\; [s]$")

# Save figure
# tight_layout()
# savefig("fishing_problem.pdf")

# Do not return anything
nothing


text:
using ErrorBudget
using Test

@testset "ErrorBudget.jl" begin
    # Write your tests here.
end


text:
function simulatetrajectories(mpc::MPCSystem1D,
                              osc::OfflineSystemControl1D,
                              ωmodel::UnivariateDistribution,
                              x0,
                              numsimulations::Int = 1000,
                              includemean::Bool = true;
                              ωtrue::UnivariateDistribution = ωmodel,
                              verbose::Bool = false,
                              optimizer::Optim.AbstractOptimizer = LBFGS(),
                              aguessinit::Vector = fill(0.5*(mpc.system.amin+mpc.system.amax), mpc.system.T))
    @assert mpc.system === osc.system
    system = mpc.system

    osctrajectories = Vector{DynamicSystemTrajectory1D}(numsimulations)
    mpctrajectories = Vector{DynamicSystemTrajectory1D}(numsimulations)
    atmpc = copy(aguessinit)

    for sim = 1:numsimulations
        trajosc = DynamicSystemTrajectory1D(system)
        trajmpc = DynamicSystemTrajectory1D(system)
        initializestate!(trajosc, x0)
        initializestate!(trajmpc, x0)
        aguess = copy(aguessinit)

        for t = 0:system.T-1
            numscenarios = mpc.numscenarios[t+1]
            includemean = includemean || (numscenarios == 1) # Always predict the future using the mean when only one sample
            if t > 0
                aguess = zeros(system.T-t)
                numsteps = length(aguess)

                for step = 1:numsteps
                    idxs = [(step+1)+numsteps*(scenario-1) for scenario=1:numscenarios]
                    aguess[step] = mean(atmpc[idxs])
                end
            end

            atosc = osc.policy(t, trajosc.state[t+1])
            atmpc = onlinedecision(mpc, trajmpc.state[t+1], t, aguess, ωmodel, includemean;
                                   verbose = verbose, optimizer = optimizer)

            # Update guess for next time
            w = rand(ωtrue)
            step!(trajosc, atosc, t, w)
            step!(trajmpc, atmpc[1], t, w)
        end
        osctrajectories[sim] = trajosc
        mpctrajectories[sim] = trajmpc
    end

    return osctrajectories, mpctrajectories
end

function simulatetrajectories{T<:Real}(oscvec::Vector{OfflineSystemControl1D{T}},
                                       ωmodel::UnivariateDistribution,
                                       x0,
                                       numsimulations::Int = 1000;
                                       ωtrue::UnivariateDistribution = ωmodel)
    numosc = length(oscvec)
    systems = [osc.system for osc in oscvec]
    @assert length(unique(systems)) == 1
    system = systems[1]

    osctrajectories = tuple([Vector{DynamicSystemTrajectory1D}(numsimulations)
                             for k = 1:numosc]...)

    for sim = 1:numsimulations
        trajosc = Vector{DynamicSystemTrajectory1D}(numosc)
        for k = 1:numosc
            trajosc[k] = DynamicSystemTrajectory1D(system)
            initializestate!(trajosc[k], x0)
        end

        for t = 0:system.T-1
            atosc = [oscvec[k].policy(t, trajosc[k].state[t+1]) for k = 1:numosc]
            # Update guess for next time
            w = rand(ωtrue)
            for k = 1:numosc
                step!(trajosc[k], atosc[k], t, w)
            end

        end
        for k = 1:numosc
            osctrajectories[k][sim] = trajosc[k]
        end
    end

    return osctrajectories
end

function simulatetrajectories(olfc::OLFCSystem1D,
                              osc::OfflineSystemControl1D,
                              ωmodel::UnivariateDistribution,
                              x0,
                              numsimulations::Int = 1000,
                              includemean::Bool = true;
                              ωtrue::UnivariateDistribution = ωmodel,
                              verbose::Bool = false,
                              optimizer::Optim.AbstractOptimizer = LBFGS(),
                              aguessinit::Vector = fill(0.5*(olfc.system.amin+olfc.system.amax), olfc.system.T))
    @assert olfc.system === osc.system
    system = olfc.system

    osctrajectories = Vector{DynamicSystemTrajectory1D}(numsimulations)
    olfctrajectories = Vector{DynamicSystemTrajectory1D}(numsimulations)
    atolfc = copy(aguessinit)

    for sim = 1:numsimulations
        trajosc = DynamicSystemTrajectory1D(system)
        trajolfc = DynamicSystemTrajectory1D(system)
        initializestate!(trajosc, x0)
        initializestate!(trajolfc, x0)
        aguess = copy(aguessinit)

        for t = 0:system.T-1
            numscenarios = olfc.numscenarios[t+1]
            includemean = includemean || (numscenarios == 1) # Always predict the future using the mean when only one sample
            if t > 0
                aguess = atolfc[2:end]
            end

            atosc = osc.policy(t, trajosc.state[t+1])
            atolfc = onlinedecision(olfc, trajolfc.state[t+1], t, aguess, ωmodel, includemean;
                                    verbose = verbose, optimizer = optimizer)

            # Update guess for next time
            w = rand(ωtrue)
            step!(trajosc, atosc, t, w)
            step!(trajolfc, atolfc[1], t, w)
        end
        osctrajectories[sim] = trajosc
        olfctrajectories[sim] = trajolfc
    end

    return osctrajectories, olfctrajectories
end


text:
using DataInterpolations, FileIO, CSV, DataFrames, Dates, Statistics
using StatsPlots, JLD2, Flux, ForwardDiff, GalacticOptim, LinearAlgebra
using Parameters, Distributions, OrdinaryDiffEq, DiffEqCallbacks, LogExpFunctions
using RecursiveArrayTools, Interpolations, MCMCChains, Suppressor, Interpolations
import KenyaCoVaccines

##Load methods and group data into correct age groups

include("generation_methods.jl");
include("inference_methods.jl");
@load("data/p_ID.jld2")
@load("data/p_IH.jld2")
@load("data/cleaned_linelist20211109_deaths_c_age__date_of_lab_confirmation.jld2")

##Gather infection predictions by county

fitfiles = readdir("modelfits", join = true)
inf_predictions = Vector{NamedTuple}(undef, 47)
prev_inf_predictions = Vector{Vector}(undef, 47)

for (i, filename) in enumerate(fitfiles)
    model = @suppress_err load(filename)["model"]
    KenyaCoVaccines.change_prob!(model, zeros(500); startdate = Date(2020, 12, 1), enddate = Date(2021, 9, 27))
    println("Making infection predictions for county $(model.areaname)")
    inf_predictions[i] = gather_post_mean_infection_rates_by_type_novac(model, Date(2021, 9, 27))
    prev_inf_predictions[i] = model.init_sero .* mean(model.MCMC_results.chain[:init_scale]) .* model.N
end

##

@load("data/fitted_mort_by_county_age_variant.jld2", mort_scales, mort_age, mort_var)

##unscaled hosp pred infections

unscaled_hosp_pred_inf = zeros(300)
unscaled_hosp_pred_inf_wt = zeros(300)
unscaled_hosp_pred_inf_ab = zeros(300)
unscaled_hosp_pred_inf_delta = zeros(300)

temp_hosp_pred = zeros(300)
for i = 1:47
    for a = 1:6
        temp_hosp_pred .= inf_predictions[i].infections_by_age_first_pred[a] .* mort_scales[i] .* mort_age[a]
        temp_hosp_pred .+= 0.15 .* inf_predictions[i].infections_by_age_reinf_pred[a] .* mort_scales[i] .* mort_age[a]
        unscaled_hosp_pred_inf_wt .+= temp_hosp_pred .* inf_predictions[i].wt_frq_pred
        unscaled_hosp_pred_inf_ab .+= (mort_var[1]) .* temp_hosp_pred .* inf_predictions[i].alpha_beta_frq_pred
        unscaled_hosp_pred_inf_delta .+= (mort_var[2]) .* temp_hosp_pred .* inf_predictions[i].delta_frq_pred

        temp_hosp_pred .*= (inf_predictions[i].wt_frq_pred .+ (mort_var[1]) .* inf_predictions[i].alpha_beta_frq_pred .+ (mort_var[2]) .* inf_predictions[i].delta_frq_pred)
        unscaled_hosp_pred_inf .+= temp_hosp_pred

    end
end

##Plots

scatter(sum(deaths_data.deaths[deaths_data.dates.>=Date(2020, 12, 25), :, :], dims = [2, 3])[:], lab = "")
plot!(unscaled_hosp_pred_inf, lab = "", lw = 3)

##Gather hospitalisation data

kenya_MoH = CSV.File("data/2022-01-10_MOH_Kenya_user_tweets__EXTRACTED_cleaned_dates.csv") |> DataFrame
kenya_MoH.dates

kenya_MoH.rel_time = [(d - Date(2020, 12, 1)).value for d in kenya_MoH.dates]
idxs_hosp_valid = kenya_MoH.in_hospital .> 0
idxs_ICU_valid = (kenya_MoH.intensive_care_unit .> 0) .& (kenya_MoH.intensive_care_unit .< 4000)
idxs_suppoxy_valid = kenya_MoH.supplementary_oxygen .> 0

##

###
### CONSTRUCT DELAY DISTRIBUTIONS
###

d_incubation = LogNormal(1.644, 0.363)#Lauer estimate
# d_duration_in_hosp = Weibull(1.572, 17.819)#Haw et al
d_duration_in_hosp = Gamma(8.0472378052287681, 0.0017022772867586302 * 365)
mean(d_duration_in_hosp)
std(d_duration_in_hosp)
# shape = 53.273947475158401, scale = 0.0003599898843749761

# d_ICUstay = Gamma(2.66, 3.42)#Fit to IQR
d_ICUstay = Gamma(53.273947475158401, 0.0003599898843749761 * 365) #From Cameline
mean(d_ICUstay)
std(d_ICUstay)
d_recovery = Exponential(2.4)
#lag distributions
p_IS = [cdf(d_incubation, t) - cdf(d_incubation, t - 1) for t = 1:100]#infection to symptoms
p_ICU = [cdf(d_ICUstay, t) - cdf(d_ICUstay, t - 1) for t = 1:100]#ICU to leave ICU
p_HR = [cdf(d_duration_in_hosp, t) - cdf(d_duration_in_hosp, t - 1) for t = 1:100]#Hospital to discharge
p_ICUR = KenyaCoVaccines.simple_conv(p_ICU, p_HR)#ICU to discharge assuming its the sum of the two
p_SH = [0.2 for t = 1:5] #Moghadas et al uniform 1-5 day estimate
p_R = [cdf(d_recovery, t) - cdf(d_recovery, t - 1) for t = 1:1000]
#Upper tail functions
Q_HR = vcat([1.0], [1 - cumsum(p_HR)[t] for t = 1:100])
Q_ICUH = vcat([1.0], [1 - cumsum(p_ICU)[t] for t = 1:100])
Q_ICUR = vcat([1.0], [1 - cumsum(p_ICUR)[t] for t = 1:100])
Q_R = vcat([1.0], [1 - cumsum(p_R)[t] for t = 1:1000])
F_R = 1 .- Q_R

##Functions for creating fit
"""
    function pred_ICU_incidence(H_ICU, η_hosp, unscaled_pred_inf, p_IH)

Predict daily incidence at ICU
"""
function pred_ICU_incidence(H_ICU, η_hosp, unscaled_pred_inf, p_IH)
    T = eltype(H_ICU)
    _p_IH = normalize([exp(η_hosp * t) * p_IH[t] for t = 1:50], 1)
    pred_ICU_inc = zeros(T, size(unscaled_pred_inf))
    for t = 1:length(pred_ICU_inc)
        for s = 1:min(t - 1, 50)
            pred_ICU_inc[t] += H_ICU * unscaled_pred_inf[t-s] * _p_IH[s]
        end
    end
    return pred_ICU_inc
end

pred_ICU_inc = pred_ICU_incidence(3.0, 0.0, unscaled_hosp_pred_inf, p_IH)

plot(pred_ICU_inc)
scatter!(kenya_MoH.rel_time[idxs_ICU_valid], kenya_MoH.intensive_care_unit[idxs_ICU_valid])

"""
    function pred_ICU_occupancy(pred_ICU_inc, Q_ICUH)

Predict the number of people in ICU on any day.
"""
function pred_ICU_occupancy(pred_ICU_inc, Q_ICUH)
    pred_ICU_occupancy = zeros(eltype(pred_ICU_inc), length(pred_ICU_inc))
    for t = 1:length(pred_ICU_occupancy)
        for s = 1:min(t - 1, 100)
            pred_ICU_occupancy[t] += pred_ICU_inc[t-s] * Q_ICUH[s]
        end
    end
    return pred_ICU_occupancy
end
pred_ICU_inc = pred_ICU_incidence(1.5, 0.5, unscaled_hosp_pred_inf, p_IH)
pred_ICU_occ = pred_ICU_occupancy(pred_ICU_incidence(1.5, 0.5, unscaled_hosp_pred_inf, p_IH), Q_ICUH)

plot(pred_ICU_occ)
scatter!(kenya_MoH.rel_time[idxs_ICU_valid], kenya_MoH.intensive_care_unit[idxs_ICU_valid])

function regularization_ICU(H_ICU, η_hosp, α)
    reg_loss = -logpdf(Normal(log(1.5), 0.2), H_ICU)
    reg_loss += -logpdf(Normal(0.0, 0.2), η_hosp)
    reg_loss += -logpdf(Normal(log(0.25), 0.5), α)

    return reg_loss
end

function logpdf_mean_alpha_param_NB(data, mean, alpha)
    σ² = clamp(mean + (alpha * mean^2), 0.00001, 1e6)
    p_negbin = clamp(1 - (alpha * mean^2 / σ²), 1e-10, 0.9999999999)
    r_negbin = clamp(1 / alpha, 0, 1000000)
    return logpdf(NegativeBinomial(r_negbin, p_negbin), data)
end

function loss_ICU_occupancy(H_ICU, η_hosp, α, unscaled_pred_inf, p_IH, Q_ICUH, ts, ICU_occ_data)
    T = eltype(H_ICU)
    loss = T(0.0)
    loss += regularization_ICU(H_ICU, η_hosp, α)
    pred_ICU_occ = pred_ICU_occupancy(pred_ICU_incidence(H_ICU, η_hosp, unscaled_pred_inf, p_IH), Q_ICUH)
    for (k, t) in enumerate(ts)
        if t >= 90 && t <= 300
            loss -= logpdf_mean_alpha_param_NB(ICU_occ_data[k], pred_ICU_occ[t], α)
        end
    end
    return loss, pred_ICU_occ
end

@time loss, pred = loss_ICU_occupancy(1.5, 0.3, 0.25, unscaled_hosp_pred_inf, p_IH, Q_ICUH, kenya_MoH.rel_time[idxs_ICU_valid], kenya_MoH.intensive_care_unit[idxs_ICU_valid])

##Use GalacticOptim

ICU_fixed_data = (unscaled_pred_inf = unscaled_hosp_pred_inf, p_IH = p_IH, Q_ICUH = Q_ICUH, ts = kenya_MoH.rel_time[idxs_ICU_valid], ICU_occ_data = kenya_MoH.intensive_care_unit[idxs_ICU_valid])

function ICU_loss(x, p)
    H_ICU = exp(x[1])
    η_hosp = x[2]
    α = exp(x[3])
    @unpack unscaled_pred_inf, p_IH, Q_ICUH, ts, ICU_occ_data = p
    loss, pred = loss_ICU_occupancy(H_ICU, η_hosp, α, unscaled_pred_inf, p_IH, Q_ICUH, ts, ICU_occ_data)
    return loss, pred
end

_loss, _pred = ICU_loss([log(1.5), 0.3, log(0.25)], ICU_fixed_data)

recorded_losses = []
function plot_ICU_occ(x, loss, pred)
    println("loss = $(loss), alpha = $(exp(x[3]))")
    push!(recorded_losses, loss)
    plt = scatter(kenya_MoH.rel_time[idxs_ICU_valid], kenya_MoH.intensive_care_unit[idxs_ICU_valid], lab = "ICU occ data")
    plot!(plt, pred, lw = 3, lab = "pred")
    display(plt)
    return false
end

##Fit
x0 = [log(1.5), 0.3, log(0.25)]

f = OptimizationFunction(ICU_loss, GalacticOptim.AutoForwardDiff())
prob = OptimizationProblem(f, x0, ICU_fixed_data)
sol = solve(prob, Flux.ADAM(0.001); cb = plot_ICU_occ, maxiters = 3000)

##
x = sol.u
H_ICU_to_death_fit = exp(x[1])
η_ICU = x[2]
alpha_ICU = exp(x[3])
# @save("data/fitted_rel_risk_ICU.jld2", H_ICU_to_death_fit, η_ICU, alpha_ICU)

## HOSPITALISATION
best_loss_ICU, best_pred_ICU = ICU_loss(sol.u, ICU_fixed_data)
best_pred_ICU_inc = pred_ICU_incidence(H_ICU_to_death_fit, η_ICU, ICU_fixed_data.unscaled_pred_inf, ICU_fixed_data.p_IH)
best_pred_H_inc_from_ICU = KenyaCoVaccines.simple_conv(best_pred_ICU_inc, p_ICU)

function regularization_hosp(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, α)
    reg_loss = -logpdf(Normal(log(15.0), 5.0), H_hosp)
    reg_loss = -logpdf(Normal(log(15.0), 5.0), H_var_ab)
    reg_loss = -logpdf(Normal(log(17.0), 5.0), H_var_delta)

    reg_loss += -logpdf(Normal(0.0, 0.2), η_hosp)
    reg_loss += -logpdf(Normal(0.0, 0.2), η_hosp_delta)

    reg_loss += -logpdf(Normal(log(0.25), 0.5), α)

    return reg_loss
end

function pred_hosp_incidence(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, unscaled_pred_inf_wt, unscaled_pred_inf_ab, unscaled_pred_inf_delta, p_IH)
    T = eltype(H_hosp)
    _p_IH = normalize([exp(η_hosp * t) * p_IH[t] for t = 1:50], 1)
    _p_IH_delta = normalize([exp(η_hosp_delta * t) * p_IH[t] for t = 1:50], 1)

    pred_hosp_inc = zeros(T, size(unscaled_pred_inf_wt))
    for t = 1:length(pred_ICU_inc)
        for s = 1:min(t - 1, 50)
            pred_hosp_inc[t] += (H_hosp * unscaled_pred_inf_wt[t-s] + H_var_ab * unscaled_pred_inf_ab[t-s]) * _p_IH[s]
            pred_hosp_inc[t] += (H_var_delta * unscaled_pred_inf_delta[t-s]) * _p_IH_delta[s]

        end
    end
    return pred_hosp_inc
end

pred_hosp_inc = pred_hosp_incidence(10.0, 10.0, 15.0, 0.3, 0.6, unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, ICU_fixed_data.p_IH)
plot(pred_hosp_inc)

function pred_hosp_occupancy(pred_hosp_inc, best_pred_H_inc_from_ICU, Q_HR)
    pred_hosp_occupancy = zeros(eltype(pred_hosp_inc), length(pred_hosp_inc))
    for t = 1:length(pred_hosp_occupancy)
        for s = 1:min(t - 1, 100)
            # pred_hosp_occupancy[t] += (pred_hosp_inc[t-s] + best_pred_H_inc_from_ICU[t-s]) * Q_ICUH[s]
            pred_hosp_occupancy[t] += (pred_hosp_inc[t-s] + best_pred_H_inc_from_ICU[t-s]) * Q_HR[s]
        end
    end
    return pred_hosp_occupancy
end

pred_hosp_occ = pred_hosp_occupancy(pred_hosp_inc, best_pred_H_inc_from_ICU, Q_HR)
plot(pred_hosp_occ)
scatter!(kenya_MoH.rel_time[idxs_hosp_valid], kenya_MoH.in_hospital[idxs_hosp_valid])

function loss_hosp_occupancy(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, α, unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, p_IH, Q_HR, ts, hosp_occ_data)
    T = eltype(H_hosp)
    loss = T(0.0)
    loss += regularization_hosp(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, α)
    pred_hosp_occ = pred_hosp_occupancy(pred_hosp_incidence(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, p_IH), best_pred_H_inc_from_ICU, Q_HR)
    for (k, t) in enumerate(ts)
        if t >= 90 && t <= 300
            loss -= logpdf_mean_alpha_param_NB(hosp_occ_data[k], pred_hosp_occ[t], α)
        end
    end
    return loss, pred_hosp_occ
end

loss, pred = loss_hosp_occupancy(10.0, 10.0, 15.0, 0.4, 0.25, 0.4, unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, ICU_fixed_data.p_IH, Q_HR, kenya_MoH.rel_time[idxs_hosp_valid], kenya_MoH.in_hospital[idxs_hosp_valid])
plot!(pred)

hosp_fixed_data = (unscaled_hosp_pred_inf_wt = unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab = unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta = unscaled_hosp_pred_inf_delta,
    p_IH = p_IH, Q_HR = Q_HR, ts = kenya_MoH.rel_time[idxs_hosp_valid], hosp_occ_data = kenya_MoH.in_hospital[idxs_hosp_valid])

function hosp_loss(x, p)
    H_hosp = exp(x[1])
    H_var_ab = exp(x[2])
    H_var_delta = exp(x[3])
    η_hosp = x[4]
    η_hosp_delta = x[5]
    α = exp(x[6])
    @unpack unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, p_IH, Q_HR, ts, hosp_occ_data = p
    loss, pred = loss_hosp_occupancy(H_hosp, H_var_ab, H_var_delta, η_hosp, η_hosp_delta, α, unscaled_hosp_pred_inf_wt, unscaled_hosp_pred_inf_ab, unscaled_hosp_pred_inf_delta, p_IH, Q_HR, ts, hosp_occ_data)
    return loss, pred
end
_loss, _pred = hosp_loss([log(10.0), log(10.0), log(15.0), 0.4, 0.4, log(0.25)], hosp_fixed_data)
plot(_pred)
scatter!(kenya_MoH.rel_time[idxs_hosp_valid], kenya_MoH.in_hospital[idxs_hosp_valid])

recorded_losses_hosp = []
function plot_hosp_occ(x, loss, pred)
    println("loss = $(loss), alpha = $(exp(x[6]))")
    push!(recorded_losses_hosp, loss)
    plt = scatter(kenya_MoH.rel_time[idxs_hosp_valid], kenya_MoH.in_hospital[idxs_hosp_valid], lab = "Hosp occ data")
    plot!(plt, pred, lw = 3, lab = "pred")
    display(plt)
    return false
end

##
x0_h = [log(10.0), log(10.0), log(15.0), 0.4, 0.4, log(0.25)]

f_h = OptimizationFunction(hosp_loss, GalacticOptim.AutoForwardDiff())
prob_h = OptimizationProblem(f_h, x0_h, hosp_fixed_data)
sol_h = solve(prob_h, Flux.ADAM(0.001); cb = plot_hosp_occ, maxiters = 3000)
scatter(recorded_losses)

x = sol_h.u

H_hosp = exp(x[1])
H_var_ab = exp(x[2])
H_var_delta = exp(x[3])
η_hosp = x[4]
η_hosp_delta = x[5]
α = exp(x[6])
alpha_hosp = α
@save("data/fitted_rel_risk_hosp.jld2", H_hosp, H_var_ab, H_var_delta, η_hosp, alpha_hosp)


text:
g = EvolvingGraph{Node{Int}, Int}()
add_edge!(g, 1, 2, 2001)
add_edge!(g, 1, 3, 2001)
add_edge!(g, 4, 5, 2002)
add_edge!(g, 2, 3, 2002)
add_edge!(g, 5, 6, 2003)

ml = evolving_graph_to_matrices(g)
for (i,m) in enumerate(ml)
    println("Matrix $i")
    println(m)
end

A = sparse_adjacency_matrix(ml, 1)

@test A == sparse_adjacency_matrix(g, 2001)


text:
module Conv2DM
    include("../tools/glorotuniform.jl")
    using .GlorotUniform
    using LoopVectorization

    mutable struct Conv2D
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple
        kernel_size::Tuple{Int64, Int64}
        activation_function::Module

        filters::Array{Float32}
        biases::Array{Float32}

        padding::Int64
        filter::Int64
        strides::Tuple{Int64, Int64}

        ∇weights::Array{Float32}
        Vdw::Array{Float32}
        Sdw::Array{Float32}
        ∇biases_c::Array{Float32}
        Vdb::Array{Float32}
        Sdb::Array{Float32}

        padding_input::Array{Float32}
        value::Array{Float32}
        output::Array{Float32}

        ∇biases::Array{Float32}
        pre_δ::Array{Float32}
        δ::Array{Float32}

        function Conv2D(;input_shape::Tuple, filter::Int64, padding::Int64=0, kernel_size::Tuple{Int64, Int64}, strides::Tuple{Int64, Int64}=(1,1), activation_function::Module, randomization::Bool=true)
            conv_num_per_col = (input_shape[1]+2*padding-kernel_size[1])÷strides[1]+1
            conv_num_per_row = (input_shape[2]+2*padding-kernel_size[2])÷strides[2]+1
            output_shape = (conv_num_per_col, conv_num_per_row, filter)

            input_size = 1
            for i in input_shape
                input_size *= i
            end

            filter_size = (filter, input_shape[3], kernel_size[1], kernel_size[2])
            filters = !randomization ? zeros(Float32, filter_size) : Array{Float32}(GUMatrix(input_size, conv_num_per_row*conv_num_per_col*filter, filter_size))
            biases = !randomization ? zeros(Float32, filter) : Array{Float32}(GUMatrix(input_size, conv_num_per_row*conv_num_per_col*filter, filter))

            ∇weights = zeros(Float32, filter_size)
            Vdw = zeros(Float32, filter_size)
            Sdw = zeros(Float32, filter_size)
            ∇biases_c = zeros(Float32, filter)
            Vdb = zeros(Float32, filter)
            Sdb = zeros(Float32, filter)
            new(save_Conv2D, load_Conv2D, activate_Conv2D, init_Conv2D, update_Conv2D, input_shape, output_shape, kernel_size, activation_function, filters, biases, padding, filter, strides, ∇weights, Vdw, Sdw, ∇biases_c, Vdb, Sdb)
        end
    end

    function init_Conv2D(layer::Conv2D, mini_batch::Int64)
        layer.padding_input = zeros(Float32, (layer.input_shape[1]+2*layer.padding, layer.input_shape[2]+2*layer.padding, layer.input_shape[3], mini_batch))
        layer.value = zeros(Float32, layer.output_shape..., mini_batch)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)

        layer.∇biases = zeros(Float32, layer.output_shape..., mini_batch)
        layer.pre_δ = zeros(Float32, (layer.input_shape[1]+2*layer.padding, layer.input_shape[2]+2*layer.padding, layer.input_shape[3], mini_batch))
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_Conv2D(layer::Conv2D, input::Array{Float32})
        @avxt for i in axes(input, 1), j in axes(input, 2), c in axes(input, 3), b in axes(input, 4)
            layer.padding_input[i+layer.padding, j+layer.padding, c, b] = input[i, j, c, b]
        end

        x, y = layer.strides
        @avxt for i in axes(layer.value, 1), j in axes(layer.value, 2), f in axes(layer.value, 3), b in axes(layer.value, 4)
            s = 0.0f0
            for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                s += layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, b] * layer.filters[f, c, k₁, k₂]
            end
            layer.value[i, j, f, b] = s + layer.biases[f]
        end

        layer.activation_function.func!(layer.output, layer.value)
    end

    function update_Conv2D(layer::Conv2D, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=1)
        x, y = layer.strides
        layer.activation_function.get_∇biases!(layer.∇biases, layer.value, δₗ₊₁)
        @avxt layer.pre_δ .= 0.0f0

        @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3), b in axes(layer.∇biases, 4)
            for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                layer.pre_δ[(i-1)*x+k₁, (j-1)*y+k₂, c, b] += layer.filters[f, c, k₁, k₂]*layer.∇biases[i, j, f, b]
            end
        end

        @avxt for i in axes(layer.δ, 1), j in axes(layer.δ, 2), c in axes(layer.δ, 3), b in axes(layer.δ, 4)
            layer.δ[i, j, c, b] = layer.pre_δ[layer.padding+i, layer.padding+j, c, b]
        end

        if optimizer=="SGD"
            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.filters[f, c, k₁, k₂] -= α*layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]*direction
                end
                layer.biases[f] -= α*layer.∇biases[i, j, f, 1]*direction
            end

        elseif optimizer=="Minibatch_GD"
            batch_size = size(layer.∇biases, 4)
            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3), b in axes(layer.∇biases, 4)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.filters[f, c, k₁, k₂] -= α*layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, b]*layer.∇biases[i, j, f, b]/batch_size*direction
                end
                layer.biases[f] -= α*layer.∇biases[i, j, f, b]/batch_size*direction
            end

        elseif optimizer=="Adam"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt layer.∇weights .= 0.0f0
            @avxt layer.∇biases_c .= 0.0f0

            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.∇weights[f, c, k₁, k₂] += layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]
                end
                layer.∇biases_c[f] += layer.∇biases[i, j, f, 1]
            end

            @avxt for i in eachindex(layer.filters)
                layer.Vdw[i] = β₁*layer.Vdw[i] + (1-β₁)*layer.∇weights[i]
                layer.Sdw[i] = β₂*layer.Sdw[i] + (1-β₂)*layer.∇weights[i]^2
                layer.filters[i] -= α*(layer.Vdw[i]/(1-β₁^t))/(sqrt(layer.Sdw[i]/(1-β₂^t))+ϵ)*direction
            end
            @avxt for i in eachindex(layer.biases)
                layer.Vdb[i] = β₁*layer.Vdb[i] + (1-β₁)*layer.∇biases_c[i]
                layer.Sdb[i] = β₂*layer.Sdb[i] + (1-β₂)*layer.∇biases_c[i]^2
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end

        elseif optimizer=="AdaBelief"
            t = parameters[1]
            β₁ = parameters[2]
            β₂ = parameters[3]
            ϵ = parameters[4]

            @avxt layer.∇weights .= 0.0f0
            @avxt layer.∇biases_c .= 0.0f0

            @avxt for i in axes(layer.∇biases, 1), j in axes(layer.∇biases, 2), f in axes(layer.∇biases, 3)
                for k₁ in axes(layer.filters, 3), k₂ in axes(layer.filters, 4), c in axes(layer.filters, 2)
                    layer.∇weights[f, c, k₁, k₂] += layer.padding_input[(i-1)*x+k₁, (j-1)*y+k₂, c, 1]*layer.∇biases[i, j, f, 1]
                end
                layer.∇biases_c[f] += layer.∇biases[i, j, f, 1]
            end

            @avxt for i in eachindex(layer.filters)
                layer.Vdw[i] = β₁*layer.Vdw[i] + (1-β₁)*layer.∇weights[i]
                layer.Sdw[i] = β₂*layer.Sdw[i] + (1-β₂)*(layer.∇weights[i]-layer.Vdw[i])^2
                layer.filters[i] -= α*(layer.Vdw[i]/(1-β₁^t))/(sqrt(layer.Sdw[i]/(1-β₂^t))+ϵ)*direction
            end
            @avxt for i in eachindex(layer.biases)
                layer.Vdb[i] = β₁*layer.Vdb[i] + (1-β₁)*layer.∇biases_c[i]
                layer.Sdb[i] = β₂*layer.Sdb[i] + (1-β₂)*(layer.∇biases_c[i]-layer.Vdb[i])^2
                layer.biases[i] -= α*(layer.Vdb[i]/(1-β₁^t))/(sqrt(layer.Sdb[i]/(1-β₂^t))+ϵ)*direction
            end
        end
    end

    function save_Conv2D(layer::Conv2D, file::Any, id::String)
        write(file, id, "Conv2D")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"filter", layer.filter)
        write(file, id*"padding", layer.padding)
        write(file, id*"kernel_size", collect(layer.kernel_size))
        write(file, id*"strides", collect(layer.strides))
        write(file, id*"activation_function", layer.activation_function.get_name())

        write(file, id*"filters", layer.filters)
        write(file, id*"biases", layer.biases)
    end

    function load_Conv2D(layer::Conv2D, file::Any, id::String)
        layer.filters = read(file, id*"filters")
        layer.biases = read(file, id*"biases")
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        filter = read(file, id*"filter")
        padding = read(file, id*"padding")
        kernel_size = tuple(read(file, id*"kernel_size")...)
        strides = tuple(read(file, id*"strides")...)
        return (input_shape=input_shape, filter=filter, padding=padding, kernel_size=kernel_size, strides=strides)
    end
end


text:
import Base.Random: rand, srand
import RNG: AbstractRNG, gen_seed

"""
```julia
AbstractXorshift1024 <: AbstractRNG{UInt64}
```

The base abstract type for `Xorshift1024`, `Xorshift1024Star` and `Xorshift1024Plus`.
"""
abstract AbstractXorshift1024 <: AbstractRNG{UInt64}

for (star, plus) in (
        (false, false),
        (false, true),
        (true, false),
    )
    rng_name = Symbol(string("Xorshift1024", star ? "Star" : plus ? "Plus" :""))
    @eval begin
        type $rng_name <: AbstractXorshift1024
            s0::UInt64
            s1::UInt64
            s2::UInt64
            s3::UInt64
            s4::UInt64
            s5::UInt64
            s6::UInt64
            s7::UInt64
            s8::UInt64
            s9::UInt64
            s10::UInt64
            s11::UInt64
            s12::UInt64
            s13::UInt64
            s14::UInt64
            s15::UInt64
            p::Int
            function $rng_name(seed::NTuple{16, UInt64}=gen_seed(UInt64, 16))
                r = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
                srand(r, seed)
                r
            end
        end
    end # If no another @eval: type definition not allowed inside a local scope. (It's a bug of julia)
    @eval begin
        function $rng_name(seed::Integer...)
            l = length(seed)
            @assert 0 ≤ l ≤ 16
            if l == 0
                return $rng_name(gen_seed(UInt64, 16))
            end
            $rng_name(map(x -> x % UInt64, (seed..., [i for i in l+1:16]...)))
        end

        @inline function xorshift_next(r::$rng_name)
            p = r.p
            s0 = unsafe_load(Ptr{UInt64}(pointer_from_objref(r)), p + 1)
            p = (p + 1) % 16
            s1 = unsafe_load(Ptr{UInt64}(pointer_from_objref(r)), p + 1)
            s1 $= s1 << 31
            s1 $= s1 >> 11
            s1 $= s0 >> 30
            s1 $= s0
            unsafe_store!(Ptr{UInt64}(pointer_from_objref(r)), s1, p + 1)
            r.p = p
            $(star ? :(s1 * 2685821657736338717) :
              plus ? :(s1 + s0) : :s1)
        end
    end
end

function srand(r::AbstractXorshift1024, seed::NTuple{16, UInt64})
    @inbounds for i in 1:16
        unsafe_store!(Ptr{UInt64}(pointer_from_objref(r)), seed[i], i)
    end
    r.p = 0
    for i in 1:16
        xorshift_next(r)
    end
    r
end

function srand(r::AbstractXorshift1024, seed::Integer...)
    l = length(seed)
    @assert 0 ≤ l ≤ 16
    if l == 0
        srand(r, gen_seed(UInt64, 16))
    end
    srand(r, map(x -> x % UInt64, (seed..., [i for i in l+1:16]...)))
end

@inline rand(r::AbstractXorshift1024, ::Type{UInt64}) = xorshift_next(r)


text:
# run with julia --track-allocation=user

using Profile
using GeometricIntegrators

#set_config(:nls_solver, NewtonSolver)

const Δt    = 0.1
const ntime = 10000

function f(t, x, v)
    v[1] = x[2]
    v[2] = sin(x[1])
    nothing
end

function profile_pendulum(tableau, file, prof=false)
    x0 = [acos(0.4), 0.0]
    ode = ODE(f, x0)
    int = Integrator(ode, tableau, Δt)
    sol = Solution(ode, Δt, ntime; filename=file)

    integrate!(int, sol)
    write_to_hdf5(sol)

    @time set_initial_conditions!(sol, ode)
    integrate!(int, sol)
    @time write_to_hdf5(sol)

    if prof
        Profile.clear()
        Profile.clear_malloc_data()

        @profile set_initial_conditions!(sol, ode)
        integrate!(int, sol)
        @profile write_to_hdf5(sol)
    end

    close(sol)
    rm(file)
end

profile_pendulum(getTableauERK4(), "test.hdf5", true)


text:
filter_ext(path, ext::AbstractString) = filter(fn -> splitext(fn)[2] == ext, readdir(path))
filter_ext(path, exts::Union{Tuple,AbstractArray}) =
    filter(fn -> splitext(fn)[2] in exts, readdir(path))
filter_ext(path, ext::Nothing) = readdir(path)

cleankeys(name) = (_cleankey(name),)
function cleankeys(keys::Union{NamedTuple,Tuple,AbstractArray})
    Tuple(map(_cleankey, keys, ntuple(i -> i, length(keys))))
end

function _cleankey(name::Union{Symbol,AbstractString,Name,NoName}, i=1)
    if name in (NoName(), Symbol(""), Name(Symbol("")))
        Symbol("layer$i")
    else
        Symbol(name)
    end
end

noindex_to_sampled(A) = rebuild(A; dims=noindex_to_sampled(dims(A)))
function noindex_to_sampled(dims::DimTuple)
    map(dims) do d
        lookup(d) isa NoLookup ? set(d, Sampled) : d
    end
end

function maybe_typemin_as_missingval(filename::String, A::AbstractRaster{T}) where T
    if ismissing(missingval(A))
        newmissingval = typemin(Missings.nonmissingtype(T))
        base, ext = splitext(filename)
        A1 = replace_missing(A, newmissingval)
        if missing isa eltype(A1)
            A1 = replace_missing(A, missing)
        end
        @warn "`missing` cant be written to $ext, typemin for `$(eltype(A1))` of `$newmissingval` used instead"
        return A1
    elseif missing isa eltype(A)
        A1 = replace_missing(A, missingval)
    else
        return A
    end
end

# We often need to convert the locus and the lookup in the same step,
# as doing it in the wrong order can give errors.
# function convert_locus_lookup(M1::Type{<:LookupArray}, L1::Type{<:Locus}, dim::Dimension)
#     _convert(S1, L1, sampling(dim), locus(dim), span(dim), dim)
# end

# _convert(::Type{M1}, ::Type{L1}, lookup::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2<:L1} = dim
# _convert(::Type{M1}, ::Type{L1}, lookup::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2} = shiftlocus(L1(), dim)
# _convert(::Type{M1}, ::Type{L1}, lookup::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2} =
#     _convert_by_locus(M1, L1, lookup, l2, span, dim)

# _convert_by_locus(M1, ::Type{Center}, lookup, l2::Union{Start,End}, span, dim) =
#     _convert_by_lookup(M1, dim)
# _convert_by_locus(M1, L1::Type{Union{Start,End}}, lookup, l2::Center, span, dim) =
#     _convert_by_lookup(M1, L1, dim)
# _convert_by_locus(M1, L1::Type{Start}, lookup, l2::End, span, dim) =
#     convertlookup(M1, shiftlocus(L1, dim))
# _convert_by_locus(M1, L1::Type{End}, lookup, l2::Start, span, dim) =
#     convertlookup(M1, shiftlocus(L1, dim))
# _convert_by_locus(M1, ::Type{L1}, lookup, l2::L2, span, dim) where {L1,L2<:L1} =
#     convertlookup(M1, dim)

# # Projected will have an accurate center point on an equal scale, but Mapped may not.
# # So we always shift the locus while in Projected lookup to avoid errors.
# _convert_by_lookup(::Type{Mapped}, dim) = convertlookup(Mapped, shiftlocus(Center(), dim))
# _convert_by_lookup(::Type{Projected}, dim) = shiftlocus(Center(), convertlookup(Projected, dim))

_missingval_or_missing(x) = missingval(x) isa Nothing ? missing : missingval(x)

maybe_eps(dims::DimTuple) = map(maybe_eps, dims)
maybe_eps(dim::Dimension) = maybe_eps(eltype(dim))
maybe_eps(::Type) = nothing
maybe_eps(T::Type{<:AbstractFloat}) = _default_atol(T)

_writeable_missing(filename::Nothing, T) = missing
_writeable_missing(filename::AbstractString, T) = _writeable_missing(T)
function _writeable_missing(T)
    missingval = typemin(Missings.nonmissingtype(T))
    @info "`missingval` set to typemin of $missingval"
    return missingval
end

# Map filename suffix over a stack
function mapargs(f, st::AbstractRasterStack, args...)
    layers = map(values(st), args...) do A, mappedargs...
        f(A, mappedargs...)
    end
    return DD.rebuild_from_arrays(st, Tuple(layers))
end

_without_mapped_crs(f, x) = _without_mapped_crs(f, x, mappedcrs(x))
_without_mapped_crs(f, x, ::Nothing) = f(x)
function _without_mapped_crs(f, dims::DimTuple, mappedcrs::GeoFormat)
    dims1 = setmappedcrs(dims, nothing)
    x = f(dims1)
    if x isa DimTuple
        x = setmappedcrs(x, mappedcrs)
    end
    return x
end
function _without_mapped_crs(f, A::AbstractRaster, mappedcrs::GeoFormat)
    A = setmappedcrs(A, nothing)
    x = f(A)
    if x isa AbstractRaster
        x = setmappedcrs(x, mappedcrs)
    end
    return x
end
function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat)
    st1 = map(A -> setmappedcrs(A, nothing), st)
    x = f(st1)
    if x isa AbstractRasterStack
        x = map(A -> setmappedcrs(A, mappedcrs(st)), x)
    end
    return x
end

_warn_disk() = @warn "Disk-based objects may be very slow here. User `read` first."


text:
invoke_map!{F, N}(f::F, dest, As::Vararg{NullableArray, N}) =
    invoke(map!, Tuple{F, AbstractArray, Vararg{AbstractArray, N}}, f, dest, As...)

"""
    map(f, As::NullableArray...)

Call `map` with nullable lifting semantics and return a `NullableArray`.
Lifting means calling function `f` on the the values wrapped inside `Nullable` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the source `As` arrays as all
`NullableArray`s. Thus, calling `map` on arguments consisting
of both `Array`s and `NullableArray`s will fall back to the standard implementation
of `map` (i.e. without lifting).
"""
function Base.map{F}(f::F, As::NullableArray...)
    # These definitions are needed to avoid allocation due to splatting
    @inline f2(x1) = lift(f, (x1,))
    @inline f2(x1, x2) = lift(f, (x1, x2))
    @inline f2(x1, x2, x3) = lift(f, (x1, x2, x3))
    @inline f2(x1, x2, x3, x4) = lift(f, (x1, x2, x3, x4))
    @inline f2(x1, x2, x3, x4, x5) = lift(f, (x1, x2, x3, x4, x5))
    @inline f2(x1, x2, x3, x4, x5, x6) = lift(f, (x1, x2, x3, x4, x5, x6))
    @inline f2(x1, x2, x3, x4, x5, x6, x7) = lift(f, (x1, x2, x3, x4, x5, x6, x7))
    @inline f2(x...) = lift(f, x)

    T = nullable_broadcast_eltype(f, As...)
    dest = similar(NullableArray{T}, size(As[1]))
    invoke_map!(f2, dest, As...)
end

"""
    map!(f, dest::NullableArray, As::NullableArray...)

Call `map!` with nullable lifting semantics.
Lifting means calling function `f` on the the values wrapped inside `Nullable` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the destination `dest` array as well as the
source `As` arrays as all `NullableArray`s. Thus, calling `map!` on a arguments
consisting of both `Array`s and `NullableArray`s will fall back to the standard implementation
of `map!` (i.e. without lifting).
"""
function Base.map!{F}(f::F, dest::NullableArray, As::NullableArray...)
    # These definitions are needed to avoid allocation due to splatting
    @inline f2(x1) = lift(f, (x1,))
    @inline f2(x1, x2) = lift(f, (x1, x2))
    @inline f2(x1, x2, x3) = lift(f, (x1, x2, x3))
    @inline f2(x1, x2, x3, x4) = lift(f, (x1, x2, x3, x4))
    @inline f2(x1, x2, x3, x4, x5) = lift(f, (x1, x2, x3, x4, x5))
    @inline f2(x1, x2, x3, x4, x5, x6) = lift(f, (x1, x2, x3, x4, x5, x6))
    @inline f2(x1, x2, x3, x4, x5, x6, x7) = lift(f, (x1, x2, x3, x4, x5, x6, x7))
    @inline f2(x...) = lift(f, x)

    invoke_map!(f2, dest, As...)
end

# This definition is needed to avoid dispatch loops going back to the above one
function Base.map!{F}(f::F, dest::NullableArray)
    f2(x1) = lift(f, (x1,))
    invoke_map!(f2, dest, dest)
end


text:
################## ic updates from store for emulation problems simulation #################

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{InitialTimeDurationOn, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_aux_variable_value(store, TimeDurationOn(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{InitialTimeDurationOff, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_aux_variable_value(store, TimeDurationOff(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {T <: InitialCondition{DevicePower, S}} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_variable_value(store, ActivePowerVariable(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{DeviceStatus, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_variable_value(store, OnVariable(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{DeviceAboveMinPower, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val =
            get_variable_value(store, PowerAboveMinimumVariable(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{InitialEnergyLevel, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_variable_value(store, EnergyVariable(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{InitialEnergyLevelUp, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_variable_value(store, EnergyVariableUp(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end

function update_initial_conditions!(
    ics::Vector{T},
    store::EmulationModelStore,
    ::Dates.Millisecond,
) where {
    T <: InitialCondition{InitialEnergyLevelDown, S},
} where {S <: Union{Float64, PJ.ParameterRef}}
    for ic in ics
        var_val = get_variable_value(store, EnergyVariableDown(), get_component_type(ic))
        set_ic_quantity!(ic, get_last_recorded_value(var_val)[get_component_name(ic)])
    end
    return
end


text:
  module utils

  using Reexport

  export isinside,intersection,clipSH,tuple2intpoint,isoverlap

  import Luxor: Point
using Clipper

isinside(p, a, b) = (b.x - a.x) * (p.y - a.y) > (b.y - a.y) * (p.x - a.x)

function intersection(a, b, s, f)
    dc = [a.x - b.x, a.y - b.y]
    dp = [s.x - f.x, s.y - f.y]
    n1 = a.x * b.y - a.y * b.x
    n2 = s.x * f.y - s.y * f.x
    n3 = 1.0 / (dc[1] * dp[2] - dc[2] * dp[1])
    Point((n1 * dp[1] - n2 * dc[1]) * n3, (n1 * dp[2] - n2 * dc[2]) * n3)
end

function clipSH(spoly, cpoly)
    outarr = spoly
    q = cpoly[end]
    for p in cpoly
        # @show p
        inarr = outarr
        outarr = Point[]
        # DEBUG if inarr is empty, I assume the two poly does not overlap
        # if isempty(inarr) return Point[] end
        s = inarr[end]
        for vtx in inarr
            if isinside(vtx, q, p)
                if !isinside(s, q, p)
                    push!(outarr, intersection(q, p, s, vtx))
                end
                push!(outarr, vtx)
            elseif isinside(s, q, p)
                push!(outarr, intersection(q, p, s, vtx))
            end
            s = vtx
        end
        q = p
    end
    outarr
end

function isoverlap_old(spoly, cpoly)
    if isempty(clipSH(spoly, cpoly)) return false
    else return true
    end
end

function isoverlap_clip(poly0, polys...)
    path1 = Vector{IntPoint}()
    for p in poly0
        push!(path1, p)
    end
    c = Clip()
    add_path!(c, path1, PolyTypeSubject, true)

    for poly in polys
        path2 = Vector{IntPoint}()
        for p in poly
            push!(path2, p)
        end
        add_path!(c, path2, PolyTypeClip, true)
    end

    result, polys = execute(c,
                            # ClipTypeUnion,
                            # ClipTypeDifference,
                            ClipTypeIntersection,
                            PolyFillTypeEvenOdd,
                            PolyFillTypeEvenOdd)
    # @show polys
    # FIXME isempty(polys) does not work
    return length(polys) > 0
end

function tuple2intpoint(t)
    map(t) do (x,y)
        IntPoint(round(x*1e6),round(y*1e6))
    end
end

function isoverlap(poly0, polys...)
    isoverlap_clip(tuple2intpoint(poly0), map(tuple2intpoint, polys)...)
end

  end

max_stars_repo_name: bhdl/bhdl

max_issues_repo_name: bhdl/bhdl

max_forks_repo_name: bhdl/bhdl


text:
using LinearAlgebra

include("gaussnum.jl")
include("ode-solver.jl")

# Define & Solve ODE

function lotkavolterra(x,θ)
    α, β, γ, δ = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

struct GaussODEProblem{F,T<:Tuple{Number,Number},U<:AbstractMatrix,P<:AbstractVector,I} <: AbstractODEProblem
    f::F
    tspan::T
    u0::U
    θ::P
    θ_idx::I
    size::Int
end
Base.size(prob::GaussODEProblem) = size(prob.u0,1)
statesize(prob::GaussODEProblem) = prob.size

function GaussODEProblem(f,tspan,u0::V,θ::V) where V<:AbstractVector{<:GaussNum}
    θ_idx = findall(isuncertain, θ)

    μ = vcat(mean.(u0), mean.(θ[θ_idx]))
    d = length(μ)
    Σ = vcat(var.(u0), var.(θ[θ_idx])) |> Diagonal |> Matrix
    Σ½ = cholesky(Σ).L
    Qp = sqrt(d)*[I(d) -I(d)]
    Xp = μ .+ Σ½*Qp

    GaussODEProblem(f,tspan,Xp,mean.(θ),θ_idx,length(u0))
end

function setmean!(prob::GaussODEProblem,xi)
    prob.θ[prob.θ_idx] .= xi[(statesize(prob)+1):end]
    xi[1:statesize(prob)]
end

struct GaussODESolver{S<:ODESolver} <: ODESolver
    solver::S
end

function (s::GaussODESolver)(prob::GaussODEProblem, Xp, t)
    Xp = reduce(hcat, map(eachcol(Xp)) do xi
        ui = setmean!(prob, xi)
        ui = s.solver(prob, ui, t)[1]
        xi[1:statesize(prob)] .= ui
        xi
    end)
    Xp, t+s.solver.dt
end

function filter(prob::GaussODEProblem, solver::GaussODESolver, data, σy)
    t = prob.tspan[1]; u = prob.u0
    us = [u]; ts = [t]

    for (ty, y) in data
        while t <= ty
            (u,t) = solver(prob, u, t)
            push!(us,u)
            push!(ts,t)
        end

        yₚ = u[1,:] # prediction. we measure only the first variable
        μᵧ = mean(yₚ)
        μᵤ = mean(u, dims=2)

        Σᵧᵧ = cov(yₚ, corrected=false) + σy
        Σᵤᵤ = cov(u, dims=2, corrected=false)

        Σᵤᵧ = cov(u, yₚ, dims=2, corrected=false)
        K   = Σᵤᵧ*inv(Σᵧᵧ)

        μ = vec(μᵤ + K*(y .- μᵧ))
        Σ = Hermitian(Σᵤᵤ - K*Σᵤᵧ')

        d = size(prob)
        Qp = sqrt(d)*[I(d) -I(d)]
        u = μ .+ cholesky(Σ).L * Qp
    end
    ts, us
end

function lotkavolterra(x,θ)
    α, β, γ, δ = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

θ = [0.1±0.01, 0.2, 0.3, 0.2]
u0 = [1.0±0.2, 1.0±0.2]
tspan = (0., 100.)
prob = GaussODEProblem(lotkavolterra,tspan,u0,θ)
solver = GaussODESolver(RK2(0.1))

using JLD2
raw_data = load("../../lecture_12/lotkadata.jld2")
ts = Float32.(raw_data["t"][1:100:end])
ys = Float32.(raw_data["u"][:,1:100:end])
data = [(t,y[1]) for (t,y) in zip(ts,eachcol(ys))]
t, us = solve(prob,solver)
t, us = filter(prob, solver, data, 0.01)

gus = map(us) do u
    u = u[1:statesize(prob),:]
    μ = mean(u,dims=2) |> vec
    σ = std(u,dims=2) |> vec
    GaussNum.(μ,σ)
end
gus = reduce(hcat,gus)
p1 = plot(t, gus[1,:], lw=3)
plot!(p1, t, gus[2,:], lw=3)
plot!(p1, raw_data["t"], raw_data["u"][1,:], lw=3, alpha=0.5) |> display

solver = RK2(0.1)
prob = ODEProblem(lotkavolterra,tspan,u0,θ)

function Base.rand(x::AbstractVector{<:GaussNum{T}}) where T
    mean.(x) .+ std.(x) .* randn(T,length(x))
end
Base.rand(prob::ODEProblem) = ODEProblem(prob.f, prob.tspan, rand(prob.u0), rand(prob.θ))

p2 = plot()
Us = map(1:200) do i
    t, us = solve(rand(prob),solver.solver)
    us = reduce(hcat,us)
end

plot!(p2, t, GaussNum.(mean(Us)[1,:],std(Us)[1,:]), lw=3)
plot!(p2, t, GaussNum.(mean(Us)[2,:],std(Us)[2,:]), lw=3)
display(plot(p1,p2))


text:
function wavenumbers_mesh(ω::T, k_effs::Vector{Complex{T}}, medium::PhysicalMedium{Dim}, species::Species{T};
        symmetry::AbstractSymmetry = PlanarSymmetry{Dim}(),
        tol::T = 1e-5,
        mesh_refine::T = T(0.4),
        inner_optimizer = LBFGS(),
        verbose::Bool = false,
        max_Rek::T = maximum(real(k1) for k1 in k_effs),
        min_Rek::T = minimum(real(k1) for k1 in k_effs),
        max_Imk::T = maximum(imag(k1) for k1 in k_effs),
        min_Imk::T = minimum(imag(k1) for k1 in k_effs),
        kws...) where {T<:Number,Dim}

    eff_medium = effective_medium(medium, species)
    k0 = ω/eff_medium.c
    if isnan(k0) k0 = ω + T(0)*im end
    # abs(k0) can be used to non-dimensionlise k_vec
    kscale = abs(k0)

    low_tol = 1e-2*minimum( (k1 == k2) ? Inf : abs(k1-k2) for k1 in k_effs, k2 in k_effs) / kscale
    low_tol = max(low_tol, 1e-4, tol) # a tolerance used for a first pass with time_limit

    # the dispersion equation is given by: `dispersion(k1,k2) = 0` where k_eff = k1 + im*k2.
    dispersion_dim = dispersion_equation(ω, medium, species, symmetry; tol = low_tol, kws...)
    dispersion(vec::Vector{T}) = dispersion_dim((vec[1] + vec[2] * im) .* kscale)

    # non-dimensionalise
    min_Rek = min_Rek / kscale
    min_Imk = min_Imk / kscale
    max_Rek = max_Rek / kscale
    max_Imk = max_Imk / kscale

    kx = LinRange(min_Rek, max_Rek, Int(round(length(k_effs)/(2*mesh_refine))))# tree shape makes this division by 2 natural
    ky = LinRange(min_Imk, max_Imk, Int(round(length(k_effs)/(2*mesh_refine))))
    k_mesh = [[x,y] for x in kx, y in ky][:]

    # Include k_effs in the mesh, we do not assume these are solutions
    k_vecs = [[real(keff),imag(keff)] for keff in k_effs] ./ kscale
    k_mesh = [k_mesh; k_vecs]

    # make slightly bigger box to be constrained within
    min_Rek = (min_Rek < zero(T)) ? min_Rek*T(1.1) : min_Rek*T(0.9)
    min_Imk = (min_Imk < zero(T)) ? min_Imk*T(1.1) : min_Imk*T(0.9)
    max_Rek = (max_Rek > zero(T)) ? max_Rek*T(1.001) : max_Rek*T(0.999)
    max_Imk = (max_Imk > zero(T)) ? max_Imk*T(1.001) : max_Imk*T(0.999)

    filter!(keff -> min_Rek <= keff[1] <= max_Rek && min_Imk <= keff[2] <= max_Imk, k_mesh)

    # Find all wavenumbers
    lower = [min_Rek, min_Imk]
    upper = [max_Rek, max_Imk]

    new_ks = [
        optimize(dispersion, lower, upper, kvec,
            Fminbox(inner_optimizer), Optim.Options(x_tol=low_tol, g_tol = low_tol^3)).minimizer
    for kvec in k_mesh]
    deleteat!(new_ks, findall(dispersion.(new_ks) .> low_tol))
    new_ks = reduce_kvecs(new_ks, low_tol)

    filter!(keff -> lower[1] <= keff[1] <= upper[1] && lower[2] <= keff[2] <= upper[2], new_ks)

    # Here we refine the new roots
    new_ks = map(new_ks) do k_vec
        # res = optimize(dispersion, k_vec; g_tol = tol^2.0, f_tol = tol^4.0, x_tol=tol)
        res = optimize(dispersion, lower, upper, k_vec, Fminbox(inner_optimizer),
            Optim.Options(g_tol = tol^3.0, x_tol=tol))
        if res.minimum < T(100)*tol || (Optim.converged(res) && res.minimum < T(10)*low_tol)
            res.minimizer
        else
            [zero(T),-one(T)]
        end
    end

    if verbose println("New roots from mesh refiner:",new_ks) end

    # group together wavenumbers which are closer than tol
    deleteat!(k_vecs, findall(v-> v == [zero(T),-one(T)], k_vecs) )
    k_vecs = reduce_kvecs(new_ks, T(10)*tol)

    # Finally delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
    deleteat!(k_vecs, findall([k_vec[2] < -low_tol for k_vec in k_vecs]))
    deleteat!(k_vecs, findall([-low_tol < abs(k_vec[2])/k_vec[1] < zero(T) for k_vec in k_vecs]))
    # deleteat!(k_vecs, find(k_vec[2] < tol && k_vec[1] < tol for k_vec in k_vecs))

    k_effs = kscale .* [k_vec[1] + k_vec[2]*im for k_vec in k_vecs]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

# Alternative solvers
# res = nlsolve(detMM!,initial_k_eff; iterations = 10000, factor=2.0)
# k_eff_nl = res.zero[1] + im*res.zero[2]
# lower = [0.,0.]; upper = [T(2)*k0,k0]
# result = optimize(dispersion, initial_k0, lower, upper; g_tol = tol^2.0, f_tol = tol^4.0)

max_stars_repo_name: arturgower/EffectiveWaves


text:
module CartesianGridsTests

using FillArrays
using Test

using Gridap.Helpers
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Geometry
using Gridap.ReferenceFEs

domain = (0,1)
partition = (1,)
grid = CartesianGrid(domain,partition)
test_grid(grid)
test_grid(UnstructuredGrid(grid))

pmin = Point(0.0,-1.0)
pmax = Point(1.0,2.0)
partition = (3,4)
grid1 = CartesianGrid(pmin,pmax,partition)

domain = (0.0,1.0,-1.0,2.0)
partition = (3,4)
grid2 = CartesianGrid(domain,partition)

domain = [0.0,1.0,-1.0,2.0]
partition = [3,4]
grid3 = CartesianGrid(domain,partition)

origin = Point(0,-1)
sizes = (1/3,3/4)
partition = (3,4)
grid4 = CartesianGrid(origin,sizes,partition)

for grid in (grid1,grid2,grid3,grid4)
  test_grid(grid)
  @test is_oriented(grid) == true

  desc = get_cartesian_descriptor(grid)
  @test desc.origin == origin
  @test desc.sizes == sizes
  @test desc.partition == partition

  @test num_cell_dims(grid) == 2
  @test num_point_dims(grid) == 2
  @test num_cells(grid) == 3*4
  @test num_nodes(grid) == (3+1)*(4+1)

  x = get_node_coordinates(grid)
  @test length(x) == 20
  @test x[13] == Point(0.0,1.25)
  @test x[4] == Point(1.0,-1.0)
  @test x[1] == Point(0.0,-1.0)
  @test x[end] == Point(1.0,2.0)

  t = get_cell_node_ids(grid)
  @test t[1] == [1,2,5,6]
  @test t[11] == [14,15,18,19]

  cell_type = get_cell_type(grid)
  @test isa(cell_type,Fill)

  reffes = get_reffes(grid)
  @test length(reffes) == 1
  @test reffes[1] == QUAD4

  ugrid = UnstructuredGrid(grid)
  test_grid(ugrid)
  @test is_oriented(ugrid) == true
end

domain = (0,1,0,1)
partition = (2,2)

grid = CartesianGrid(domain,partition)

map = get_cell_map(grid)
∇map = lazy_map(∇,map)

x = [Point(0.5,0.5),]
ax = Fill(x,partition)
r = Vector{Point{2,Float64}}[[(0.25, 0.25)], [(0.75, 0.25)], [(0.25, 0.75)], [(0.75, 0.75)]]
∇r = Vector{TensorValue{2,2,Float64,4}}[
  [(0.5, 0.0, 0.0, 0.5)], [(0.5, 0.0, 0.0, 0.5)],
  [(0.5, 0.0, 0.0, 0.5)], [(0.5, 0.0, 0.0, 0.5)]]
#test_array_of_fields(map,ax,r,grad=∇r)
@test all(lazy_map(test_field,map,ax,r))
@test all(lazy_map(test_field,∇map,ax,∇r))

mx = lazy_map(evaluate,map,ax)
∇mx = lazy_map(evaluate,∇map,ax)

test_array(mx,reshape(r,size(map)))
test_array(∇mx,reshape(∇r,size(map)))
@test isa(∇mx,Fill)

domain = (0,1,0,1,0,1)
partition = (2,2,2)
grid = CartesianGrid(domain,partition)
tgrid = simplexify(grid)
test_grid(tgrid)

grid = compute_reference_grid(HEX8,4)
test_grid(grid)

grid = compute_linear_grid(HEX8)
test_grid(grid)

end # module


text:
using TensorFlow
include("mnist_loader.jl")

loader = DataLoader()

sess = Session(Graph())

x = placeholder(Float32)
y_ = placeholder(Float32)

W = Variable(zeros(Float32, 784, 10))
b = Variable(zeros(Float32, 10))

run(sess, global_variables_initializer())

y = nn.softmax(x*W + b)

cross_entropy = reduce_mean(-reduce_sum(y_ .* log(y), axis=[2]))
train_step = train.minimize(train.GradientDescentOptimizer(.00001), cross_entropy)

correct_prediction = indmax(y, 2) .== indmax(y_, 2)
accuracy=reduce_mean(cast(correct_prediction, Float32))

for i in 1:1000
    batch = next_batch(loader, 100)
    run(sess, train_step, Dict(x=>batch[1], y_=>batch[2]))
end

testx, testy = load_test_set()

println(run(sess, accuracy, Dict(x=>testx, y_=>testy)))

max_stars_repo_name: octonion/tensorflow

max_issues_repo_name: octonion/tensorflow

max_forks_repo_name: octonion/tensorflow