text:
module HeightOrderedGrid

using StaticArrays

export HOGrid

"""
    immutable HOGrid{T, M <: AbstractMatrix}

A spatial index suitable for 3D data clustered about a plane (so-called
2.5D data) such as large-scale point clouds of terrain. Allows for fast querying
of points withing a radius of a specified location.
"""
immutable HOGrid{T}
    grid::Matrix{Vector{SVector{3,T}}}
    x0::T
    y0::T
    spacing::T
end

"""
    HOGrid(points::Vector{SVector{3,T}}, gridspacing) -> grid

Instantiate a height-ordered grid of `points` with grid spacing `gridspacing`.
During the process, the `points` will be copied and reordered to maximise cache
coherency on spatial querying.
"""
function HOGrid{T}(points::Vector{SVector{3,T}}, gridspacing)
    if length(points) == 0
        error("Cannot create HOGrid for zero points")
    end

    (xmin, xmax, ymin, ymax) = bounds(points)
    xmin -= eps(xmin) # make sure we never fall into grid 0
    ymin -= eps(ymin)

    n_x = Int(cld(xmax - xmin, gridspacing))
    n_y = Int(cld(ymax - ymin, gridspacing))

    grid = Matrix{Vector{SVector{3,T}}}((n_x, n_y))
    for i = 1:n_x
        for j = 1:n_y
            grid[i,j] = Vector{SVector{3,T}}() # Should this be done sparsely?
        end
    end

    # Stash each point into the correct grid cell
    for p ∈ points
        @inbounds x = p[1]
        @inbounds y = p[2]

        grid_x = Int(cld(x - xmin, gridspacing))
        grid_y = Int(cld(y - ymin, gridspacing))

        @inbounds push!(grid[grid_x, grid_y], p)
    end

    # Sort all the grid cells on height
    for g ∈ grid
        sort!(g; lt = height_less)
    end

    return HOGrid(grid, xmin, ymin, gridspacing)
end

immutable SpatialQuery{T}
    grid::HOGrid{T}
    query_point::SVector{3,T}
    radius::T
    gx_min::Int
    gx_max::Int
    gy_min::Int
    gy_max::Int
end

function SpatialQuery{T}(grid::HOGrid{T}, query_point::SVector{3,T}, radius::T)
    r² = radius * radius

    @inbounds qx = query_point[1]
    @inbounds qy = query_point[2]

    # Make a square cutoff
    xmin = qx - radius
    xmax = qx + radius
    ymin = qy - radius
    ymax = qy + radius

    # Make a square grid
    gx_min = Int(cld(xmin - grid.x0, grid.spacing))
    if gx_min < 1
        gx_min = 1
    end
    gx_max = Int(cld(xmax - grid.x0, grid.spacing))
    if gx_max > size(grid.grid, 1)
        gx_max = size(grid.grid, 1)
    end

    gy_min = Int(cld(ymin - grid.y0, grid.spacing))
    if gy_min < 1
        gy_min = 1
    end
    gy_max = Int(cld(ymax - grid.y0, grid.spacing))
    if gy_max > size(grid.grid, 2)
        gy_max = size(grid.grid, 2)
    end

    return SpatialQuery(grid, query_point, radius, gx_min, gx_max, gy_min, gy_max)
end

@inline height{T <: Number}(p::SVector{3,T}) = @inbounds return p[3]
@inline height{T <: Number}(h::T) = h

function Base.start{T}(sq::SpatialQuery{T})
    gx = sq.gx_min
    gy = sq.gy_min
    r² = sq.radius * sq.radius

    @inbounds qx = sq.query_point[1]
    @inbounds qy = sq.query_point[2]
    @inbounds qz = sq.query_point[3]

    while true
        if gy > sq.gy_max
            break
        end

        # Fit a sphere and find the vertical range to query, if any
        if qy > sq.grid.y0 + gy * sq.grid.spacing
            dy = qy - (sq.grid.y0 + gy * sq.grid.spacing)
        elseif qy >= sq.grid.y0 + (gy - 1) * sq.grid.spacing
            dy = zero(T)
        else
            dy = (sq.grid.y0 + (gy - 1) * sq.grid.spacing) - qy
        end

        while true
            if gx > sq.gx_max
                gx = sq.gx_min
                break
            end

            # Fit a sphere and find the vertical range to query, if any
            if qx > sq.grid.x0 + gx * sq.grid.spacing
                dx = qx - (sq.grid.x0 + gx * sq.grid.spacing)
            elseif qx >= sq.grid.x0 + (gx - 1) * sq.grid.spacing
                dx = zero(T)
            else
                dx = (sq.grid.x0 + (gx - 1) * sq.grid.spacing) - qx
            end

            h² = dx*dx + dy*dy
            if h² > r²
                gx += 1
                continue
            end
            dz = sqrt(r² - h²)

            # Now find the first element in this range
            @inbounds cell = sq.grid.grid[gx, gy]
            i = searchsortedfirst(cell, qz - dz, Base.By(height))

            # Now loop until we find one
            while true
                if i > size(cell, 1)
                    break
                end

                @inbounds p = cell[i]
                diff = p - sq.query_point
                if dot(diff, diff) <= r²
                    return (gx, gy, i, dz)
                end

                if height(p) - qz > dz
                    break
                end

                i += 1
            end

            gx += 1
        end

        gy += 1
    end

    return (gx, gy, i, dz)
end

function Base.next{T}(sq::SpatialQuery{T}, state::Tuple{Int,Int,Int,T})
    (gx, gy, i, dz) = state
    r² = sq.radius * sq.radius
    @inbounds qx = sq.query_point[1]
    @inbounds qy = sq.query_point[2]
    @inbounds qz = sq.query_point[3]

    # If we are calling next, then gx, gy, i is a valid state
    @inbounds cell = sq.grid.grid[gx, gy]
    @inbounds p_now = cell[i]

    # increment i and check the same cell before continuing
    while true
        i += 1

        if i > size(cell, 1)
            break
        end

        @inbounds p = cell[i]
        diff = p - sq.query_point
        if dot(diff, diff) <= r²
            return (p_now, (gx, gy, i, dz))
        end

        if height(p) - qz > dz
            break
        end
    end

    # increment gx and check other cells
    gx += 1

    while true
        if gy > sq.gy_max
            break
        end

        # Fit a sphere and find the vertical range to query, if any
        if qy > sq.grid.y0 + gy * sq.grid.spacing
            dy = qy - (sq.grid.y0 + gy * sq.grid.spacing)
        elseif qy >= sq.grid.y0 + (gy - 1) * sq.grid.spacing
            dy = zero(T)
        else
            dy = (sq.grid.y0 + (gy - 1) * sq.grid.spacing) - qy
        end

        while true
            if gx > sq.gx_max
                gx = sq.gx_min
                break
            end

            # Fit a sphere and find the vertical range to query, if any
            if qx > sq.grid.x0 + gx * sq.grid.spacing
                dx = qx - (sq.grid.x0 + gx * sq.grid.spacing)
            elseif qx >= sq.grid.x0 + (gx - 1) * sq.grid.spacing
                dx = zero(T)
            else
                dx = (sq.grid.x0 + (gx - 1) * sq.grid.spacing) - qx
            end

            h² = dx*dx + dy*dy
            if h² > r²
                gx += 1
                continue
            end
            dz = sqrt(r² - h²)

            # Now find the first element in this range
            @inbounds cell = sq.grid.grid[gx, gy]
            i = searchsortedfirst(cell, qz - dz, Base.By(height))

            # Now loop until we find one
            while true
                if i > size(cell, 1)
                    break
                end

                @inbounds p = cell[i]
                diff = p - sq.query_point
                if dot(diff, diff) <= r²
                    return (p_now, (gx, gy, i, dz))
                end

                if height(p) - qz > dz
                    break
                end

                i += 1
            end

            gx += 1
        end

        gy += 1
    end

    return (p_now, (gx, gy, i, dz))
end

Base.done{T}(sq::SpatialQuery{T}, state::Tuple{Int,Int,Int,T}) = state[2] > sq.gy_max

immutable SpatialIterator{T, ZMin, ZMax}
    grid_cells::Vector{Vector{SVector{3,T}}}
    dz::Vector{T}
    query_point::SVector{3,T}
    radius::T
    zmin::ZMin
    zmax::ZMax
end

# Returns an spatial iterator for getting points within a radius of a point
@inline function inrange{T}(grid::HOGrid{T}, query_point, radius; zmin = nothing, zmax = nothing)
    grid_cells = Vector{Vector{SVector{3,T}}}()
    dz = Vector{T}()
    r² = radius * radius

    @inbounds qx = query_point[1]
    @inbounds qy = query_point[2]

    # Make a square cutoff
    xmin = qx - radius
    xmax = qx + radius
    ymin = qy - radius
    ymax = qy + radius

    # Make a square grid
    g_xmin = Int(cld(xmin - grid.x0, grid.spacing))
    if g_xmin < 1
        g_xmin = 1
    end
    g_xmax = Int(cld(xmax - grid.x0, grid.spacing))
    if g_xmax > size(grid.grid, 1)
        g_xmax = size(grid.grid, 1)
    end

    g_ymin = Int(cld(ymin - grid.y0, grid.spacing))
    if g_ymin < 1
        g_ymin = 1
    end
    g_ymax = Int(cld(ymax - grid.y0, grid.spacing))
    if g_ymax > size(grid.grid, 2)
        g_ymax = size(grid.grid, 2)
    end

    # Refine the square grid into a cicular one and compose the cells
    for g_x ∈ g_xmin:g_xmax
        # Check if qx is in this cell, or to the left or right
        if qx > grid.x0 + g_x * grid.spacing
            dx = qx - (grid.x0 + g_x * grid.spacing)
        elseif qx >= grid.x0 + (g_x - 1) * grid.spacing
            dx = zero(T)
        else
            dx = (grid.x0 + (g_x - 1) * grid.spacing) - qx
        end

        for g_y ∈ g_ymin:g_ymax
            # Check if qy is in this cell, or to above or below
            if qy > grid.y0 + g_y * grid.spacing
                dy = qy - (grid.y0 + g_y * grid.spacing)
            elseif qy >= grid.y0 + (g_y - 1) * grid.spacing
                dy = zero(T)
            else
                dy = (grid.y0 + (g_y - 1) * grid.spacing) - qy
            end

            hdist² = dx*dx + dy*dy
            if hdist² <= r²
                @inbounds push!(grid_cells, grid.grid[g_x, g_y])
                #push!(horizontal_distance², hdist²)
                push!(dz, sqrt(r² - hdist²))
            end
        end
    end

    SpatialIterator(grid_cells, dz, query_point, radius, zmin, zmax)
end

@inline function Base.start(si::SpatialIterator)
    if length(si.grid_cells) == 0
        return (1,0)
    end

    (g, i) = getfirst(si)
    return (g, i)
end

@inline function Base.next(si::SpatialIterator, gi::Tuple{Int, Int})
    (g, i) = getnext(gi[1], gi[2], si)
    @inbounds return (si.grid_cells[gi[1]][gi[2]], (g, i))
end

@inline function Base.done(si::SpatialIterator, gi::Tuple{Int, Int})
    return gi[1] > size(si.grid_cells, 1)
end

@inline function getfirst{T}(si::SpatialIterator{T, Void, Void})
    g = 1
    @inbounds cell = si.grid_cells[g]
    @inbounds dz = si.dz[g]
    # Interesting new (undocumented) interface to sorting in Julia Base:
    # Note that the old method is currently not type stable in Julia 0.5 and
    # results in a noticeable slowdown of the entire algorithm!
    #i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz); by = p -> p[3]) # by is applied to comparitor as well... (???)
    i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz), Base.By(p -> p[3]))
    r² = si.radius * si.radius

    while true # Loop over g
        while true # Loop over i
            if i > size(cell, 1)
                break
            end

            @inbounds p = cell[i]
            @inbounds if p[3] > si.query_point[3] + dz
                break
            end

            diff = si.query_point - p
            if dot(diff, diff) <= r²
                return (g, i)
            end

            i += 1
        end

        g += 1
        if g > size(si.grid_cells, 1)
            return (g, i)
        else
            @inbounds cell = si.grid_cells[g]
            @inbounds dz = si.dz[g]
            #i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz); by = p -> p[3]) # by is applied to comparitor as well... (???)
            i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz), Base.By(p -> p[3]))
        end
    end
end

@inline function getnext{T}(g::Int, i::Int, si::SpatialIterator{T, Void, Void})
    @inbounds cell = si.grid_cells[g]
    i += 1
    r² = si.radius * si.radius
    @inbounds dz = si.dz[g]

    while true # Loop over g

        while true # Loop over i
            if i > size(cell, 1)
                break
            end

            @inbounds p = cell[i]
            @inbounds if p[3] > si.query_point[3] + dz
                break
            end

            diff = si.query_point - p
            if dot(diff, diff) <= r²
                return (g, i)
            end

            i += 1
        end

        g += 1
        if g > size(si.grid_cells, 1)
            return (g, i)
        else
            @inbounds cell = si.grid_cells[g]
            @inbounds dz = si.dz[g]
            #i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz); by = p -> p[3]) # by is applied to comparitor as well... (???)
            i = searchsortedfirst(cell, si.query_point - SVector(0.0, 0.0, dz), Base.By(p -> p[3]))
        end
    end
end

@inline function height_less(p1::SVector{3}, p2::SVector{3})
    @inbounds z1 = p1[3]
    @inbounds z2 = p2[3]
    return z1 < z2
end

function bounds{T}(points::Vector{SVector{3,T}})
    xmin = typemax(T)
    xmax = typemin(T)
    ymin = typemax(T)
    ymax = typemin(T)

    for p ∈ points
        @inbounds x = p[1]
        if x < xmin
            xmin = x
        end
        if x > xmax
            xmax = x
        end

        @inbounds y = p[2]
        if y < ymin
            ymin = y
        end
        if y > ymax
            ymax = y
        end
    end

    return (xmin, xmax, ymin, ymax)
end

end # module


text:
abstract type AbstractExecution end
# Execution API #
# ------------------------------------------------------------
nthetas(lshaped::AbstractLShapedSolver) = nthetas(lshaped, lshaped.execution)
initialize_subproblems!(lshaped::AbstractLShapedSolver, scenarioproblems::AbstractScenarioProblems, x::AbstractVector, subsolver::MPB.AbstractMathProgSolver) = initialize_subproblems!(lshaped.execution, scenarioproblems, x, subsolver)
finish_initilization!(lshaped::AbstractLShapedSolver) = finish_initilization!(lshaped, lshaped.execution)
solve_master!(lshaped::AbstractLShapedSolver) = solve_master!(lshaped, lshaped.execution)
iterate!(lshaped::AbstractLShapedSolver) = iterate!(lshaped, lshaped.execution)
start_workers!(lshaped::AbstractLShapedSolver) = start_workers!(lshaped, lshaped.execution)
close_workers!(lshaped::AbstractLShapedSolver) = close_workers!(lshaped, lshaped.execution)
resolve_subproblems!(lshaped::AbstractLShapedSolver) = resolve_subproblems!(lshaped, lshaped.execution)
current_decision(lshaped::AbstractLShapedSolver) = current_decision(lshaped, lshaped.execution)
calculate_objective_value(lshaped::AbstractLShapedSolver) = calculate_objective_value(lshaped, lshaped.execution)
timestamp(lshaped::AbstractLShapedSolver) = timestamp(lshaped, lshaped.execution)
incumbent_decision(lshaped::AbstractLShapedSolver, t::Integer, regularizer::AbstractRegularization) = incumbent_decision(lshaped, t, regularizer, lshaped.execution)
incumbent_objective(lshaped::AbstractLShapedSolver, t::Integer, regularizer::AbstractRegularization) = incumbent_objective(lshaped, t, regularizer, lshaped.execution)
incumbent_trustregion(lshaped::AbstractLShapedSolver, t::Integer, regularizer::AbstractRegularization) = incumbent_trustregion(lshaped, t, regularizer, lshaped.execution)
readd_cuts!(lshaped::AbstractLShapedSolver, consolidation::Consolidation) = readd_cuts!(lshaped, consolidation, lshaped.execution)
subobjectives(lshaped::AbstractLShapedSolver) = subobjectives(lshaped, lshaped.execution)
set_subobjectives(lshaped::AbstractLShapedSolver, Qs::AbstractVector) = set_subobjectives(lshaped, Qs, lshaped.execution)
model_objectives(lshaped::AbstractLShapedSolver) = model_objectives(lshaped, lshaped.execution)
set_model_objectives(lshaped::AbstractLShapedSolver, Qs::AbstractVector) = set_model_objectives(lshaped, Qs, lshaped.execution)
fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::AbstractScenarioProblems) = fill_submodels!(lshaped, scenarioproblems, lshaped.execution)

# ------------------------------------------------------------
include("common.jl")
include("serial.jl")
include("distributed.jl")
include("synchronous.jl")
include("asynchronous.jl")

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.3

using Markdown
using InteractiveUtils

# ╔═╡ 504df9a8-0aff-11eb-13b1-6d8a03e123aa
using Pkg, DrWatson

# ╔═╡ 515f176e-0aff-11eb-342c-9f6e5c550932
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 425b24f6-0aff-11eb-124e-8f7e7dc68c84
md"## Clip-05-50-54t.jl"

# ╔═╡ 2998907c-0aff-11eb-279d-67723fa39fb3
md"### snippet 5.50"

# ╔═╡ a86c956c-0b00-11eb-0e54-a7fa8640895e
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; missingstring = "NA")
	df.clade |> unique |> sort
end;

# ╔═╡ a86d54ca-0b00-11eb-217f-798f08670430
md"### snippet 5.51"

# ╔═╡ a87df190-0b00-11eb-3747-c35abd52dda3
md"##### You can't just turn Strings into Integers in Julia but hashing them should give the same result"

# ╔═╡ a87e7c9e-0b00-11eb-12f2-53d3c65e745c
df.clade_id = Int.(indexin(df.clade, unique(df.clade)))

# ╔═╡ a88a6312-0b00-11eb-1238-0d168fc40850
md"### snippet 5.52"

# ╔═╡ a88af58e-0b00-11eb-3729-4d96fd03ce11
df.K = zscore(df.kcal_per_g)

# ╔═╡ a896bdae-0b00-11eb-09c8-4185f6e08bea
md"##### Clade DynamicPPL model"

# ╔═╡ a8989c8c-0b00-11eb-2964-8bb8b153bf44
@model function ppl5_9(clade_id, K)
    σ ~ Exponential(1)
    α ~ filldist(Normal(0, 0.5), length(unique(clade_id)))
    μ = α[clade_id]
    K ~ MvNormal(μ, σ)
end

# ╔═╡ a8a04290-0b00-11eb-3ea7-afbefe340644
begin
	m5_9t = ppl5_9(df.clade_id, df.K)
	quap5_9t = quap(m5_9t)
	post5_9t = DataFrame(rand(quap5_9t.distr, 1000)', quap5_9t.params)
	Text(precis(post5_9t; io=String))
end

# ╔═╡ a8a809bc-0b00-11eb-3d71-9fdca69bcc10
# TODO: plot

md"### snippet 5.54"

# ╔═╡ a8b74526-0b00-11eb-2e7e-43231317823c
df.house = rand(1:4, nrow(df));

# ╔═╡ a8ba5926-0b00-11eb-17d3-ed90f3311578
@model function ppl5_10(clade_id, house, K)
    σ ~ Exponential(1)
    a ~ filldist(Normal(0, 0.5), length(unique(clade_id)))
    h ~ filldist(Normal(0, 0.5), length(unique(house)))
    μ = a[clade_id] + h[house]
    K ~ MvNormal(μ, σ)
end

# ╔═╡ a8c94a96-0b00-11eb-0290-ef5f4f51848e
begin
	m5_10t = ppl5_10(df.clade_id, df.house, df.K)
	quap5_10t = quap(m5_10t)
	Text(precis(DataFrame(rand(quap5_10t.distr, 1000)', quap5_10t.params); io=String))
end

# ╔═╡ a8cc6bb8-0b00-11eb-1f59-099af38de865
md"## End of clip-05-45-54.jl"

# ╔═╡ Cell order:
# ╟─425b24f6-0aff-11eb-124e-8f7e7dc68c84
# ╠═504df9a8-0aff-11eb-13b1-6d8a03e123aa
# ╠═515f176e-0aff-11eb-342c-9f6e5c550932
# ╟─2998907c-0aff-11eb-279d-67723fa39fb3
# ╠═a86c956c-0b00-11eb-0e54-a7fa8640895e
# ╟─a86d54ca-0b00-11eb-217f-798f08670430
# ╟─a87df190-0b00-11eb-3747-c35abd52dda3
# ╠═a87e7c9e-0b00-11eb-12f2-53d3c65e745c
# ╟─a88a6312-0b00-11eb-1238-0d168fc40850
# ╠═a88af58e-0b00-11eb-3729-4d96fd03ce11
# ╟─a896bdae-0b00-11eb-09c8-4185f6e08bea
# ╠═a8989c8c-0b00-11eb-2964-8bb8b153bf44
# ╠═a8a04290-0b00-11eb-3ea7-afbefe340644
# ╠═a8a809bc-0b00-11eb-3d71-9fdca69bcc10
# ╠═a8b74526-0b00-11eb-2e7e-43231317823c
# ╠═a8ba5926-0b00-11eb-17d3-ed90f3311578
# ╠═a8c94a96-0b00-11eb-0290-ef5f4f51848e
# ╟─a8cc6bb8-0b00-11eb-1f59-099af38de865


text:
module plutonotebooks

end # module

max_stars_repo_name: matutenun/projectbinder

max_issues_repo_name: matutenun/projectbinder

max_forks_repo_name: matutenun/projectbinder


text:
module FrameworkPacote

using LinearAlgebra
using NLPModels
using SolverTools

include("gradiente.jl")

end # module


text:
using Statistics: mean
using Oceananigans.Grids: Center, Face
using Oceananigans.Fields: Field, location
using Oceananigans.OutputReaders: FieldTimeSeries, FieldDataset

"""
    coarse_grain(Φ, n, ::Center)

Average or coarse grain a `Center`-centered field `Φ` down to size `n`. `Φ` is required to have evenly spaced points and `n` needs to evenly divide `length(Φ)`.
"""
function coarse_grain(Φ, n, ::Center)
    N = length(Φ)
    Δ = Int(N / n)
    Φ̅ = similar(Φ, n)
    for i in 1:n
        Φ̅[i] = mean(Φ[Δ*(i-1)+1:Δ*i])
    end
    return Φ̅
end

"""
    coarse_grain(Φ, n, ::Face)

Average or coarse grain a `Face`-centered field `Φ` down to size `n`. `Φ` is required to have evenly spaced points. The values at the left and right endpoints of `Φ` will be preserved in the output.
"""
function coarse_grain(Φ, n, ::Face)
    N = length(Φ)
    Φ̅ = similar(Φ, n)
    Δ = (N-2) / (n-2)

    Φ̅[1] = Φ[1]
    Φ̅[n] = Φ[N]

    if isinteger(Δ)
        Φ̅[2:n-1] .= coarse_grain(Φ[2:N-1], n-2, Center())
    else
        for i in 2:n-1
            i1 = 2 + (i-2)*Δ
            i2 = 2 + (i-1)*Δ

            # Like modf but with ::Int integer part.
            f1, i1 = rem(i1, 1), trunc(Int, i1)
            f2, i2 = rem(i2, 1), trunc(Int, i2)

            left_contrib = (1 - f1) * Φ[i1]
            right_contrib = f2 * Φ[i2]
            middle_contrib = sum(Φ[i1+1:i2-1])

            Φ̅[i] = (left_contrib + middle_contrib + right_contrib) / Δ
        end
    end

    return Φ̅
end

@deprecate coarse_grain(Φ, n, loc::Type{Center}) coarse_grain(Φ, n, loc::Center)
@deprecate coarse_grain(Φ, n, loc::Type{Face}) coarse_grain(Φ, n, loc::Face)

"""
    coarse_grain(field::Field, new_grid; dims=3)

Coarse grain a `field` onto a `new_grid` along `dims`. Returns a new `Field`.
"""
function coarse_grain(field::Field{X, Y, Center}, new_grid; dims=3) where {X, Y}

    # TODO: Generalize to x and y.
    @assert dims == 3
    @assert new_grid.Nx == 1
    @assert new_grid.Ny == 1

    # TODO: Generalize `coarse_grain` to non-integer ratios.
    r = field.grid.Nz / new_grid.Nz
    @assert isinteger(r)
    r = Int(r)

    coarse_field = Field(location(field)..., field.architecture, new_grid, field.boundary_conditions)

    coarse_data = zeros(size(coarse_field))
    field_interior = interior(field)[1, 1, :]
    coarse_data[1, 1, :] .= OceanParameterizations.coarse_grain(field_interior, new_grid.Nz, Center())

    set!(coarse_field, coarse_data)

    return coarse_field
end

function coarse_grain(field::Field{X, Y, Face}, new_grid; dims=3) where {X, Y}

    # TODO: Generalize to x and y.
    @assert dims == 3
    @assert new_grid.Nx == 1
    @assert new_grid.Ny == 1

    r = field.grid.Nz / new_grid.Nz

    coarse_field = Field(location(field)..., field.architecture, new_grid, field.boundary_conditions)

    coarse_data = zeros(size(location(coarse_field), new_grid))
    field_interior = interior(field)[1, 1, :]
    coarse_data[1, 1, :] .= OceanParameterizations.coarse_grain(field_interior, new_grid.Nz+1, Face())

    set!(coarse_field, coarse_data)

    return coarse_field
end

function coarse_grain(fts::FieldTimeSeries, new_grid)
    fts_new = FieldTimeSeries(new_grid, location(fts), fts.times; ArrayType=Array{Float32})

    Nt = size(fts, 4)
    for n in 1:Nt
        fts_new.data[:, :, :, n] .= coarse_grain(fts[n], new_grid).data
    end

    return fts_new
end

function coarse_grain(fds::FieldDataset, new_grid)
    coarse_fields = Dict{String, FieldTimeSeries}(
        name => coarse_grain(fts, new_grid)
        for (name, fts) in fds.fields
    )
    return FieldDataset(coarse_fields, fds.metadata, fds.filepath)
end


text:
#==============================================================================#
# SCFT Simulation
#==============================================================================#

"""
	struct SCFTOptions

	SCFTOptions(steps = 1000, nsout = 0, nsave = 0, ftol = 1e-5, stol = 1e-5, domain = false, savepath = "")

Input options for an SCFT calculation.
"""
@with_kw struct SCFTOptions
	nsteps    :: Int = 100
	nsout     :: Int = 0
	nsave	  :: Int = 0
	ftol      :: Float64 = 1e-5
	stol      :: Float64 = 1e-4
	savepath  :: String = ""
end

"""
	struct SCFTResults

Container for output summary from an SCFT calculation.
"""
@with_kw struct SCFTResults
	time      :: Float64
	nsteps    :: Int
	converged :: Bool
	ferr      :: Float64
	serr      :: Float64
	fhelm     :: Float64
end

#==============================================================================#

function isconverged(ferr, ftol, serr, stol, vcell = false)
	fconverged = false
	sconverged = true
	if ferr < ftol
		fconverged = true
		if vcell
			if serr > stol
				sconverged = false
			end
		end
	end
	return fconverged && sconverged
end

"""
	scft!(sys, updater; domain = nothing, opts = SCFTOptions())

Run an SCFT calculation to relax the system configuration with the provided field updater.
The system is checked with the `validate` method prior to the calculation
to ensure all fields are properly initialized.

Returns a summary in the form of an `SCFTResults` object.
"""
function scft!(sys::FieldSystem, updater::AbstractFieldUpdater;
	domain::Option{DomainUpdater} = nothing, opts::SCFTOptions = SCFTOptions())
	simtime = @elapsed begin
		# Validate and initial setup
		validate(sys)
		residuals!(sys)
		setup!(updater, sys)

		# Iterate SCFT
		step = 0
		ferr = field_error(sys)
		fhelm = free_energy(sys)

		serr = 0.0
		vcell = !isnothing(domain)
		if vcell
			setup!(domain, sys)
			serr = stress_error(domain)
		end

		converged = isconverged(ferr, opts.ftol, serr, opts.stol, vcell)
		if !converged
			# Print header of output
			println("--------------------------------------------------------------")
		end

		while step < opts.nsteps && !converged
			if opts.nsout > 0 && step % opts.nsout == 0
				@printf("Step (%-5d / %5d): ferr = %8e, fhelm = %.5e\n", step, opts.nsteps, ferr, fhelm)
			end

			# Update the fields
			step!(updater)

			# Update the cell if necessary
			if vcell && step % domain.nskip == 0
				scfstress!(domain)
				serr = stress_error(domain)
				if serr > opts.stol
					step!(domain; nsteps = domain.nper, tol = opts.stol)
					serr = stress_error(domain)
					@printf("Updating cell dimensions: serr = %.5e\n", serr)
				end
			end

			# Track error
			ferr = field_error(sys)
			fhelm = free_energy(sys)

			if isnan(ferr)
				error("Diverging field updater -- breaking.") # Better way to check for this?
			end

			# Check for convergence of fields and cell
			converged = isconverged(ferr, opts.ftol, serr, opts.stol, vcell)
			if converged
				@printf("Converged @ step = %d: ferr = %8e, fhelm = %.5e\n", step, ferr, fhelm)
				if vcell
					@printf("Converged cell parameters: %s\n", sys.cell.params)
				end
				break
			end

			# Check if output should be saved
			if opts.nsave > 0 && step % opts.nsave == 0
				if !isempty(opts.savepath)
					savefields(sys, opts.savepath)
				end
			end

			step += 1
		end

		# Make sure system state is current
		residuals!(sys)
		ferr = field_error(sys)
		fhelm = free_energy(sys)

		if vcell
			serr = stress_error(domain)
		end

		if !isempty(opts.savepath)
			savefields(sys, opts.field_out)
		end
	end

	return SCFTResults(simtime, step, converged, ferr, serr, fhelm)
end


text:
function perform_haar_transf(M, Jmin, dir)
    n = size(M)[1]
    Mf = copy(M)
    if dir == 1
        MW = copy(M)

        for j in Jmin : Int(log2(n))
            p = Int(n/2^(j-1))
            sel = 1 : p
            even = 1 : 2 : p
            odd = 2 : 2 : p
            # average/ difference along X
            MW[sel, sel, sel] = cat(1, (MW[even, sel, sel] + MW[odd, sel, sel])./sqrt(2), (MW[even, sel, sel] - MW[odd, sel, sel])./sqrt(2))
            # average/ difference along Y
            MW[sel, sel, sel] = cat(2, (MW[sel, even, sel] + MW[sel, odd, sel])./sqrt(2), (MW[sel, even, sel] - MW[sel, odd, sel])./sqrt(2))
            # average/ difference along Z
            MW[sel, sel, sel] = cat(3, (MW[sel, sel, even] + MW[sel, sel, odd])./sqrt(2), (MW[sel, sel, even] - MW[sel, sel, odd])./sqrt(2))
            Mf = copy(MW)
        end
    else

        M1 = copy(M)

        for j in Int(log2(n)) : -1 : Jmin
            p = Int(n/2^j)
            sel = 1 : p
            sel1 = 1:2*p
            selw = p+1:2*p
            even = 1:2:2*p
            odd = 2:2:2*p
            # average/ difference along X
            A = M1[sel, sel1, sel1]
            D = M1[selw, sel1, sel1]
            M1[even, sel1, sel1] = (A + D)./sqrt(2)
            M1[odd, sel1, sel1] = (A - D)./sqrt(2)
            # average/ difference along Y
            A = M1[sel1, sel, sel1]
            D = M1[sel1, selw, sel1]
            M1[sel1, even, sel1] = (A + D)./sqrt(2)
            M1[sel1, odd, sel1] = (A - D)./sqrt(2)
            # average/ difference along Z
            A = M1[sel1, sel1, sel]
            D = M1[sel1, sel1, selw]
            M1[sel1, sel1, even] = (A + D)./sqrt(2)
            M1[sel1, sel1, odd] = (A - D)./sqrt(2)
            Mf = copy(M1)
        end
    end

    return Mf
end


text:
# Edge
# ====

struct Edge
    labels::ByteSet
    precond::Precondition
    actions::ActionList
end

function Edge(labels::ByteSet)
    return Edge(labels, Precondition(), ActionList())
end

function Edge(labels::ByteSet, actions::ActionList)
    return Edge(labels, Precondition(), actions)
end


text:
using DataStructures

function check(x, as)
    N = length(as)
    tmps = [(a < x) ? -1 : 1 for a in as]
    ss = vcat(0, cumsum(tmps))
    ps = sortperm(ss)
    bs = FenwickTree{Int}(length(ps))

    q_inv = 0
    for p in ps
        q_inv += prefixsum(bs, p)
        inc!(bs, p, 1)
    end

    return 2q_inv ≥ N * (N +1) ÷ 2
end

function main()
    N = parse(Int, readline())
    as = parse.(Int, split(readline()))

    x̲, x̄ = minimum(as), maximum(as)
    while x̄ - x̲ > 1
        x = (x̲ + x̄) ÷ 2

        if check(x, as)
            x̲ = x
        else
            x̄ = x
        end
    end

    if check(x̄, as)
        println(x̄)
    else
        println(x̲)
    end

end

main()


text:
# A SearchSpace determines which candidate points can be
# considered in a search/optimization. This base class has very few restrictions
# and can allow varying number of dimensions etc.
abstract SearchSpace

# SearchSpace with a fixed number of dimensions. The vast majority of problems.
abstract FixedDimensionSearchSpace <: SearchSpace

# In a ContinuousSearchSpace each dimension has a continuous range of
# valid values.
abstract ContinuousSearchSpace <: FixedDimensionSearchSpace

ranges(css::ContinuousSearchSpace) = [range_for_dim(css, i) for i in 1:numdims(css)]

# Access individual range for a dim.
range_for_dim(css::ContinuousSearchSpace, i::Int) = (mins(css)[i], maxs(css)[i])

# Generate a number of individuals by random sampling in the search space.
function rand_individuals(css::ContinuousSearchSpace, numIndividuals)
  # Basically min + delta * rand(), but broadcast over the columns...
  broadcast(+, mins(css), broadcast(*, deltas(css), rand(numdims(css), numIndividuals)))
end

# Generate a number of individuals via latin hypercube sampling. This should
# be the default when creating the initial population.
function rand_individuals_lhs(css::ContinuousSearchSpace, numIndividuals)
  BlackBoxOptim.Utils.latin_hypercube_sampling(mins(css), maxs(css), numIndividuals)
end

# Generate one random candidate.
function rand_individual(css::ContinuousSearchSpace)
  rand_individuals(css, 1)[:,1:1]
end

# True iff ind is within the search space.
function isinspace(ind, css::ContinuousSearchSpace)
  all(mins(css) .<= ind .<= maxs(css))
end

# In a RangePerDimSearchSpace each dimension is specified as a range of valid
# values.
type RangePerDimSearchSpace <: ContinuousSearchSpace
  # We save the ranges as individual mins, maxs and deltas for faster access later.
  mins::Array{Float64,2}
  maxs::Array{Float64,2}
  deltas::Array{Float64,2}

  function RangePerDimSearchSpace(ranges)
    mins = map(t -> t[1], ranges)''
    maxs = map(t -> t[2], ranges)''
    new(mins, maxs, (maxs - mins))
  end
end
mins(rss::RangePerDimSearchSpace) = rss.mins
maxs(rss::RangePerDimSearchSpace) = rss.maxs
deltas(rss::RangePerDimSearchSpace) = rss.deltas
numdims(rss::RangePerDimSearchSpace) = size(mins(rss), 1)

diameters(rss::RangePerDimSearchSpace) = deltas(rss)

# Convenience function to create symmetric search spaces.
symmetric_search_space(numdims, range = (0.0, 1.0)) = RangePerDimSearchSpace([range for i in 1:numdims])

# Create a feasible point (i.e. within the search space) given one which is
# outside.
function feasible(v, ss::RangePerDimSearchSpace)
  minimum(hcat(maxs(ss), maximum(hcat(mins(ss), v), 2)), 2)
end


text:
"""
    depth_first_impl(g, i, j, verbose = true)
    depth_first_impl(g, i_key, i_t, j_key, j_t, verbose = true)

Find the shortest temporal path between TimeNode `i` and `j` on an evolving graph `g` using DFS. Alternatively, we could inpute the node keys `i_key`, `j_key` and node timestamps `i_t`, `j_t` respectively.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = EvolvingGraph{Node{String}, Int}()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_bunch_of_edges!(g, [("A", "B", 1), ("B", "F", 1), ("B", "G", 1), ("C", "E", 2), ("E", "G", 2), ("A", "B", 2), ("A", "B", 3), ("C", "F", 3), ("E","G", 3)])
Directed EvolvingGraph 6 nodes, 9 static edges, 3 timestamps

julia> depth_first_impl(g, "A", 1, "F", 3)
Current path: TimeNode(A, 1)
Current path: TimeNode(A, 1)->TimeNode(B, 1)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(F, 1)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(F, 1)->TimeNode(F, 3)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(G, 1)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(G, 1)->TimeNode(G, 2)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(G, 1)->TimeNode(G, 3)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(B, 2)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(B, 2)->TimeNode(B, 3)
Current path: TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(B, 3)
Current path: TimeNode(A, 1)->TimeNode(A, 2)
Current path: TimeNode(A, 1)->TimeNode(A, 2)->TimeNode(B, 2)
Current path: TimeNode(A, 1)->TimeNode(A, 2)->TimeNode(B, 2)->TimeNode(B, 3)
Current path: TimeNode(A, 1)->TimeNode(A, 2)->TimeNode(A, 3)
Current path: TimeNode(A, 1)->TimeNode(A, 2)->TimeNode(A, 3)->TimeNode(B, 3)
Current path: TimeNode(A, 1)->TimeNode(A, 3)
Current path: TimeNode(A, 1)->TimeNode(A, 3)->TimeNode(B, 3)
TimeNode(A, 1)->TimeNode(B, 1)->TimeNode(F, 1)->TimeNode(F, 3)
```
"""
function depth_first_impl(g::Union{AbstractEvolvingGraph, IntAdjacencyList}, i::Union{TimeNode, Tuple{Int,Int}}, j::Union{TimeNode, Tuple{Int,Int}}, p = TemporalPath(), shortest = TemporalPath(); verbose = true)

    p = deepcopy(p)
    push!(p, i)
    if verbose
        println("Current path: $(p)")
    end

    if i == j
        return p
    end

    for n in forward_neighbors(g, i)
        if !(n in p)
            n_t = typeof(n) <: TimeNode ? node_timestamp(n) : n[2]
            j_t = typeof(j) <: TimeNode ? node_timestamp(j) : j[2]
            if n_t <= j_t
                if shortest == TemporalPath() || (length(p) < length(shortest))
                    new_p = depth_first_impl(g, n, j, p, shortest, verbose = verbose)
                    if new_p != TemporalPath()
                        shortest = new_p
                    end
                end
            end
        end
    end
    return shortest

end
function depth_first_impl{V,E,T,KV}(g::EvolvingGraph{V,E,T,KV}, key_i::KV, t_i::T, key_j::KV, t_j::T)
    if !((key_i, t_i) in keys(g.active_node_indexof)) || !((key_j, t_j) in keys(g.active_node_indexof))
        return p
    end

    index_i = g.active_node_indexof[(key_i, t_i)]
    index_j = g.active_node_indexof[(key_j, t_j)]

    i = TimeNode(index_i, key_i, t_i)
    j = TimeNode(index_j, key_j, t_j)
    return depth_first_impl(g, i, j)
end

function depth_first_impl(g::IntAdjacencyList, key_i::Int, t_i::Int, key_j::Int, t_j::Int)
    return depth_first_impl(g, (key_i, t_i), (key_j, t_j))
end


text:
const MultivectorArray = MVector{2^N}
const COMBINATIONS = [[], ga.combinations(1:N)...]

KVector{K}(coefs::MVector) where {K} = KVector{K,typeof(coefs)}(coefs, linear_index(K))
KVector{K,V}(coefs::V) where {K,V<:AbstractVector} = KVector{K,V}(coefs, linear_index(K))

zero(::Type{<:KVector{K}}, T) where {K} = KVector{K}(zeros(MVector{binomial(N, K),T}))
zero(::Type{<:Multivector}, T) = Multivector(zeros(MultivectorArray))

fill(v, T::Type{<:KVector{K}}) where {K} = T(fill(v, MVector{binomial(N, K)}))
fill(v, ::Type{<:Multivector}) = Multivector(fill(v, MultivectorArray))

linear_index(grade::Integer) = linear_index(COMBINATIONS, grade)
linear_index(indices) = linear_index(COMBINATIONS, indices)

indices_from_linear_index(index::Integer) = indices_from_linear_index(COMBINATIONS, index)

grade_projection(x::Blade, ::Val{K}) where {K} = grade(x) == K ? x : K > N ? scalar(zero(eltype(x))) : scalar(zero(eltype(x)))

grade(b::Blade) = length(indices_from_linear_index(b.index))

function grade(mv::Multivector)
    higher_nonzero_index = findlast(≉(0), mv.coefs)
    isnothing(higher_nonzero_index) && return N
    findfirst(≤(higher_nonzero_index), linear_index.(0:N) .+ 1) - 1
end

dual(x) = x / I

(∨)(xs...) = ∧(dual.(xs)...) ⋅ I

reverse(x::Multivector) = Multivector(map(x -> reverse_sign(length(indices_from_linear_index(x[1]))) * x[2], enumerate(x.coefs)))

"""
Get `Blade` element with index `(i, j)`.
"""
getindex(mv::Multivector, i::CartesianIndex) = getindex(mv.coefs, linear_index(SVector(i.I)))

show(io::IO, @nospecialize b::Blade) = print(io, b.coef, "v$(join(map(subscript, indices_from_linear_index(b.index))))")

@generated op_result_type(::typeof(*), T1, T2) = Multivector{MultivectorArray{promote_type(eltype(T1), eltype(T2))}}

kvector(mv::Multivector, ::Val{0}) = scalar(first(mv.coefs))

function kvector(mv::Multivector, ::Val{K}) where {K}
    start = linear_index(K)
    KVector{K}(mv.coefs[start:start + binomial(N, K) - 1])
end

kvectors(mv::Multivector) = (kvector(mv, Val(k)) for k ∈ 0:N)


text:
import StatisticalNonlocality: chebyshev, fourier
import StatisticalNonlocality: fourier_nodes, fourier_wavenumbers

@testset "Fourier Grid and Wavenumbers" begin
    for N in 1:5
        x = fourier_nodes(N, a = 0, b = N)
        @test all(x - collect(0:(N - 1)) .≤ eps(N * 1.0))
    end

    N = 4
    k = fourier_wavenumbers(N, L = 2π)
    @test all(k - [0, 1, -2, -1] .≤ eps(N * 1.0))
end

@testset "Fourier Derivative" begin
    N = 8
    D, x = fourier(N)
    # only wavenumbers 0, 1, 2, ..., div(N,2)-1 are well represented
    # differentiation matrix is exact for n = 0, sin(x), cos(x), sin(2x), cos(2x)
    # all the wave up to div(N,2)-1
    for n in 1:(div(N, 2) - 1)
        @test all(abs.(D * sin.(n * x) - n * cos.(n * x)) .≤ eps(10.0 * N))
        @test all(abs.(D * cos.(n * x) + n * sin.(n * x)) .≤ eps(10.0 * N))
    end
end

@testset "1D Helmholtz Equation" begin
    N = 8
    a, b = 0, 2π
    k = fourier_wavenumbers(N, L = b - a)
    Dx, x = fourier(N, a = a, b = b)
    H = Dx * Dx - I
    rhs = sin.(x) + sin.(k[3] * x)
    numerical = H \ rhs
    analytic = sin.(x) / (-1 - k[2]^2) + sin.(k[3] * x) / (-1 - k[3]^2)
    tolerance = 1e2 * eps(maximum(analytic))
    @test norm(numerical - analytic, Inf) ≤ tolerance
end

@testset "Chebyshev Grid and Wavenumbers" begin
    N = 2
    D, x = chebyshev(N)
    xᴬ = [1, 0, -1]
    Dᴬ = [3/2 -2 1/2; 1/2 0 1/2; -1/2 2.0 -3/2]
    @test all(x - xᴬ .≤ eps(N * 1.0))
    @test all(D .- Dᴬ .≤ eps(N * 1.0))
end

@testset "Chebyshev Derivative" begin
    N = 8
    D, x = chebyshev(N)
    for n in 3:(N - 1)
        @test all(D * (x .^ n) - n * x .^ (n - 1) .≤ eps(N * 3.0))
    end
end

@testset "Chebyshev Poisson" begin
    N = 8
    D, x = chebyshev(N)

    Δ = D * D
    Δ[1, :] .= 0.0
    Δ[1, 1] = 1.0
    Δ[end, :] .= 0.0
    Δ[end, end] = 1.0

    rhs = ones(N + 1) * 2
    rhs[1] = 0.0
    rhs[end] = 0.0
    numerical = Δ \ rhs
    analytic = x .^ 2 .- 1
    tolerance = 1e2 * eps(norm(analytic, Inf))
    @test norm(numerical - analytic, Inf) ≤ tolerance
end

@testset "Fourier-Chebyshev Heat Equation" begin
    N = 8
    M = 8
    Dz, z = chebyshev(M)
    a, b = 0, 2π
    k = fourier_wavenumbers(N, L = b - a)
    Dx, x = fourier(N, a = a, b = b)
    x = reshape(x, (N, 1))
    z = reshape(z, (1, M + 1))

    ∂z = kron(Dz, I + zeros(N, N))
    ∂x = kron(I + zeros(M + 1, M + 1), Dx)

    Δ = ∂x * ∂x + ∂z * ∂z
    bΔ = copy(Δ)
    # Figure out boundary indices without thinking
    zlifted = sparse(kron(Diagonal(z[:]), I + zeros(N, N)))
    ∂Ω¹ = zlifted.rowval[zlifted.nzval .== z[1]]
    ∂Ω² = zlifted.rowval[zlifted.nzval .== z[end]]
    ∂Ω = vcat(∂Ω¹, ∂Ω²)

    # apply Dirichlet boundary conditions
    for ∂i in ∂Ω
        bΔ[∂i, :] .= 0.0
        bΔ[∂i, ∂i] = 1.0
    end

    # Simple Test: Δϕ = 1 & ϕ(boundary) = 1 => z^2 /2 + 1/2
    rhs = ones(size(Δ)[2])
    numerical = bΔ \ rhs
    analytic = @. x * 0 + z^2 / 2 + 1 / 2
    numerical = reshape(numerical, size(analytic))
    e1 = norm(analytic[:] - numerical[:], Inf)
    tolerance = 1e3 * eps(norm(analytic[:], Inf))
    @test e1 < tolerance

    # Simple Test 2: Δϕ = 2 * sin(2 * x) + (z^2 - 1) * 4 * sin(2x)
    # & ϕ(boundary) = 0
    # => (z^2 - 1) * sin(2 x)
    rhs = @. 2 * sin(k[3] * x) - (z^2 - 1) * k[3]^2 * sin(k[3] * x)
    rhs = rhs[:] # flatten rhs
    # homogenous boundary conditions
    for ∂i in ∂Ω
        rhs[∂i] = 0.0
    end
    numerical = bΔ \ rhs
    analytic = @. (z^2 - 1) * sin(k[3] * x)
    numerical = reshape(numerical, size(analytic))
    e2 = norm(analytic[:] - numerical[:], Inf)
    tolerance = 1e2 * eps(norm(analytic[:], Inf))
    @test e2 < tolerance

    # Simple Test 3: A block diagonal equation of 2: [Δ 0 ; 0 Δ]
    Δ = ∂x^2 + ∂z^2
    bΔ = [Δ 0*Δ; 0*Δ Δ]

    zlifted = sparse(kron(Diagonal(z[:]), I + zeros(N, N)))
    zlifted = [zlifted 0*zlifted; 0*zlifted zlifted]
    ∂Ω¹ = zlifted.rowval[zlifted.nzval .== z[1]]
    ∂Ω² = zlifted.rowval[zlifted.nzval .== z[end]]
    ∂Ω = vcat(∂Ω¹, ∂Ω²)

    rhs1 = @. 2 * sin(k[3] * x) - (z^2 - 1) * k[3]^2 * sin(k[3] * x)
    rhs2 = @. 2 * sin(k[4] * x) - (z^2 - 1) * k[4]^2 * sin(k[4] * x)
    rhs = vcat(rhs1[:], rhs2[:]) # flatten rhs
    # homogenous boundary conditions
    for ∂i in ∂Ω
        bΔ[∂i, :] .= 0.0
        bΔ[∂i, ∂i] = 1.0
        rhs[∂i] = 0.0
    end
    numerical = qr(bΔ) \ rhs
    analytic1 = @. (z^2 - 1) * sin(k[3] * x)
    analytic2 = @. (z^2 - 1) * sin(k[4] * x)
    analytic = vcat(analytic1[:], analytic2[:])
    e3 = norm(analytic - numerical, Inf)
    tolerance = 1e3 * eps(norm(analytic[:], Inf))
    @test e3 < tolerance
end

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
#=##############################################################################
# DESCRIPTION
    Special grid types
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2018
  * License   : MIT License
=###############################################################################

################################################################################
# TRIANGULAR SURFACE GRID TYPE
################################################################################

"""
  `GridTriangleSurface(orggrid, dimsplit)`

Receives a 3D surface grid (like the one in `Paneled Wing Example` in the
documentation), which by construction is made of nonplanar quadrilateral panels,
and creates a surface grid of planar triangular panels by splitting every
original quadrilateral panel into triangles.

  **Arguments**
  * `orggrid`         : Original quadrilateral surface grid.
  * `dimsplit`        : Dimension along which to split the quadrilaterals.
"""
mutable struct GridTriangleSurface <: AbstractGrid

  # User inputs
  orggrid::Grid                       # Original quadrilateral surface grid
  dimsplit::Int64                     # Dimension to split

  # Properties
  dims::Int64                         # Number of dimensions
  nnodes::Int64                       # Number of nodes
  ncells::Int64                       # Number of cells
  field::Dict{String, Dict{String, Any}}  # Calculated fields

  # Internal data
  _ndivsnodes::Tuple                  # Number of nodes in each coordinate
  _ndivscells::Tuple                  # Number of cells in each coordinate
  _override_vtkcelltype::Int64        # Option for overriding vtk outputs

  GridTriangleSurface(orggrid, dimsplit,
                      dims=orggrid.dims,
                        nnodes=orggrid.nnodes,
                        ncells=2*orggrid.ncells,
                        field=Dict{String, Dict{String, Any}}(),
                      _ndivsnodes=orggrid._ndivsnodes,
                        _ndivscells=_ndivscells(orggrid, dimsplit),
                        _override_vtkcelltype=5
      ) = _checkGridTriangleSurface(orggrid, dimsplit) ? new(orggrid, dimsplit,
                      dims,
                        nnodes,
                        ncells,
                        field,
                      _ndivsnodes,
                        _ndivscells,
                        _override_vtkcelltype
      ) : error("Logic error!")
end

get_node(self::GridTriangleSurface, i::Int64) = get_node(self.orggrid, i)
get_node(self::GridTriangleSurface, coor::Array{Int64,1}) = get_node(self.orggrid, coor)

function get_cell(self::GridTriangleSurface, i::Int64)
  if i>self.ncells
    error("Requested invalid cell index $i; max is $(self.ncells).")
  end
  dims = Tuple(d != 0 ? d : 1 for d in self._ndivscells)
  return get_cell(self, collect(Tuple(CartesianIndices(dims)[i])))
end

function get_cell(self::GridTriangleSurface, coor::Array{Int64,1})
  # ERROR CASES
  if length(coor)!=self.dims
    error("$(self.dims)-dimensional grid requires $(self.dims) coordinates,"*
            " got $(length(coor)).")
  end
  for (dim, i) in enumerate(coor)
    if i>self._ndivscells[dim]
      if i==1 && self._ndivscells[dim]==0
        nothing
      else
        error("Requested cell $coor but max cell in"*
              " $dim-dimension is $(self._ndivscells[dim])")
      end
    end
  end

  # Converts this coordinates to the coordinates of the quadrilateral panel
  quadcoor = Int64[ceil(ind/2^(i==self.dimsplit)) for (i,ind) in enumerate(coor)]

  # Gets the nodes of the quadrilateral panel
  quadnodes = get_cell(self.orggrid, quadcoor)

  # Splits the quadrilateral into a triangle
  if coor[self.dimsplit]%2!=0
    return [quadnodes[1], quadnodes[2], quadnodes[3]]
  else
    return [quadnodes[3], quadnodes[4], quadnodes[1]]
  end
end

"""
    `get_area(self::GridTriangleSurface, i_or_coor::Union{Int, Array{Int,1}})`
Returns the area of the i-th cell.
"""
function get_area(self::GridTriangleSurface, i)
    cell = get_cell(self, i)
    A = get_node(self, cell[1])
    B = get_node(self, cell[2])
    C = get_node(self, cell[3])
    return 0.5*norm(cross(B-A, C-A))
end

"""
    `get_area(self::GridTriangleSurface)`
Returns the area of the entire grid.
"""
function get_area(self::GridTriangleSurface)
    A = 0
    for i in 1:self.ncells
        A += get_area(self, i)
    end
    return A
end

"""
    `get_volume(self::GridTriangleSurface)`
Returns the volume encloused by the grid using Green's theorem. See
https://en.wikipedia.org/wiki/Green%27s_theorem#Area_Calculation
"""
function get_volume(self::GridTriangleSurface)
    V = 0
    for i in 1:self.ncells
        A = get_area(self, i)
        n = get_normal(self, i)
        r = get_cellcenter(self, i)
        V += 1/3 * dot(r, A*n)
    end
    return V
end

"""
    `get_centroid(self::GridTriangleSurface)`
Returns the centroid of the volume encloused by the grid using Green's theorem.
See https://en.wikipedia.org/wiki/Green%27s_theorem#Area_Calculation
"""
function get_centroid(self::GridTriangleSurface)
    R = zeros(self.dims)
    for i in 1:self.ncells
        A = get_area(self, i)
        n = get_normal(self, i)
        r = get_cellcenter(self, i)
        R += 3/4*r * (1/3 * dot(r, A*n))
    end
    return R/get_volume(self)
end

"""
  `get_normal(self::GridTriangleSurface, i::Int64)`

Returns the normal vector of the i-th panel.
"""
function get_normal(self::GridTriangleSurface, i::Int64)
  return _calc_normal(get_cellnodes(self, i))
end
function get_normal(self::GridTriangleSurface, coor::Array{Int64,1})
  return get_normal(self, sub2ind(self._ndivsnodes, coor...))
end

"""
  `get_tangent(self::GridTriangleSurface, i::Int64)`

Returns the tangential vector of the i-th panel.
"""
function get_tangent(self::GridTriangleSurface, i::Int64)
  return _calc_tangent(get_cellnodes(self, i))
end
function get_tangent(self::GridTriangleSurface, coor::Array{Int64,1})
  return get_tangent(self, sub2ind(self._ndivsnodes, coor...))
end

"""
  `get_unitvectors(self::GridTriangleSurface, i::Int64)`

Returns the orthogonal system `(t, o, n)` of the i-th panel, with `t` the
tangent vector, `o` the oblique vector (which is also tangent), and `n` the
normal vector.
"""
function get_unitvectors(self::GridTriangleSurface, i::Int64)
  return _calc_unitvectors(get_cellnodes(self, i))
end
function get_unitvectors(self::GridTriangleSurface, coor::Array{Int64,1})
  return get_unitvectors(self, sub2ind(self._ndivsnodes, coor...))
end

function lintransform!(self::GridTriangleSurface, args...; optargs...)
  lintransform!(self.orggrid, args...; optargs...)
end

function transform!(self::GridTriangleSurface, args...; optargs...)
  transform!(self.orggrid, args...; optargs...)
end
function transform2!(self::GridTriangleSurface, args...; optargs...)
  transform!(self.orggrid, args...; optargs...)
end
function transform3!(self::GridTriangleSurface, args...; optargs...)
  transform!(self.orggrid, args...; optargs...)
end

##### INTERNAL FUNCTIONS  ######################################################
function _checkGridTriangleSurface(orggrid::Grid, dimsplit::Int64)
  # Number of quasi-dimensions
  qdims = length( [1 for ndiv in orggrid._ndivscells if ndiv==0] )

  if qdims!=1
    error("Expected one quasi-dimension, found $qdims.")

  elseif orggrid.dims!=3
    error("Expected a three-dimensional grid, got $(orggrid.dims)-dimensional.")

  elseif !(dimsplit in [1,2])
    error("Invalid split dimension $dimsplit")

  end

  return true
end

function _ndivscells(orggrid::Grid, dimsplit::Int64)
  return Tuple([
            2^(i==dimsplit)*divs for (i,divs) in enumerate(orggrid._ndivscells)
              ])
end

function _calc_tangent(nodes::Array{Array{T,1},1}) where{T<:Real}
  t = nodes[2] - nodes[1]
  return t/norm(t)
end

function _calc_normal(nodes::Array{Array{T,1},1}) where{T<:Real}
  n = cross(nodes[2]-nodes[1], nodes[3]-nodes[1])
  return n/norm(n)
end

function _calc_unitvectors(nodes::Array{Array{T,1},1}) where{T<:Real}
  t = _calc_tangent(nodes)
  n = _calc_normal(nodes)
  o = cross(n,t)
  return t, o, n
end
##### END OF TRIANGULAR GRID ###################################################


text:
module JuDocTemplates

export newsite

const JUDOC_PATH = splitdir(pathof(JuDocTemplates))[1] # .../JuDocTemplates/src
const TEMPL_PATH = joinpath(JUDOC_PATH, "templates")
const LIST_OF_TEMPLATES = ("sandbox", "basic", "jemdoc", "just-the-docs", "hyde",
                           "hypertext", "lanyon", "pure-sm", "tufte", "vela")

"""
    newsite(topdir; template="basic", cd=true)

Generate a new folder (an error is thrown if it already exists) that contains the skeleton of a
website that can be processed by JuDoc. The user can specify a `template` out of the list of
available templates.

* `template="basic"` specify the name of the desired template,
* `cd=true` specify whether to change the current directory to the newly created folder or not.
* `verbose=true` specify whether to display information or not.

### Example
```julia
newsite("MyNewWebsite", template="pure-sm")
```
"""
function newsite(topdir::String="TestWebsite";
                 template::String="basic",
                 changedir::Bool=true,
                 verbose::Bool=true)

    template = lowercase(template)
    template ∈ LIST_OF_TEMPLATES || throw(ArgumentError("Template $template doesn't exist."))

    # create the top-directory
    topdir = mkdir(topdir)

    # create the sub-directories
    for foldername ∈ ("libs", "assets", "src")
        subdir = joinpath(topdir, foldername)
        # common/foldername
        cp(joinpath(TEMPL_PATH, "common", foldername), subdir)
        # template/foldername
        template_subdir = joinpath(TEMPL_PATH, template, foldername)
        isdir(template_subdir) && mergefolders(template_subdir, subdir)
    end

    # Pkg.jl does something odd with file permissions, restoring to 644
    # otherwise there may be files that are read-only which is annoying
    for (root, _, files) ∈ walkdir(topdir)
        for file in files
            run(`chmod 644 $(joinpath(root, file))`)
        end
    end

    # move to the directory if relevant
    changedir && cd(topdir)

    # display information as adequate
    verbose && begin
        print("✓ Website folder generated at \"$(topdir)\"")
        println(ifelse(changedir, " (now the current directory)." , "."))
        println("→ Use serve() from JuDoc to see the website in your browser.")
    end

    return nothing
end

#
# utils
#

"""
    mergefolders(src, dst)

Internal function to looks at what's inside `src/` and put it in `dst/`. If there are paths that
match, the files are merged. It is assumed that files will not clash, if they clash files in `dst`
are replaced. See also [`newsite`](@ref).
"""
function mergefolders(src, dst)
    for (root, _, files) ∈ walkdir(src)
        for file ∈ files
            newpath = replace(root, Regex("^$(escape_string(src))")=>"$dst")
            isdir(newpath) || mkpath(newpath)
            newpathfile = joinpath(newpath, file)
            cp(joinpath(root, file), newpathfile; force=true)
        end
    end
end

end # module


text:
## MNIST data using Non-negative matrix Factorization
using MLDatasets
train_x, train_y = MNIST.traindata()
test_x,  test_y  = MNIST.testdata()

## Create the input for NMF
m = size(train_x,1)*size(train_x,2)
n = size(train_x,3)
X = Float64.(reshape(train_x, m, n))
A = copy(X') # input

## Show the average image
using Images
heatmap(Gray.(reshape(sum(A;dims=1)/size(A,1),28,28)'),
  framestyle=:none,yflip=true,colorbar=true)

## Show a random set of images
imgs = Gray.(reshape(A[1:25,:]',28,28,25))
p = plot([heatmap(imgs[:,:,i]') for i=1:25]...,framestyle=:none,yflip=true)
## Turn this into a function
myheatmap(X) = heatmap(Gray.(X'./maximum(X)),framestyle=:none,yflip=true)
plotimgs(X) = plot([myheatmap(X[:,:,i]) for i=1:size(X,3)]...)
plotimgs(reshape(A[1:25,:]',28,28,25))
## Use NMF to find structure
#using ImageMagick
using NMF
## Make a rank 5 factorization
k = 5
W, H = NMF.nndsvd(A, k; variant= :ar)
# alginst = NMF.MultUpdate{Float64}(obj=:mse, maxiter=1000, verbose=false)
alginst = NMF.ALSPGrad{Float64}(maxiter=100)
r = NMF.solve!(alginst, A, W, H)
## show the images
plotimgs(reshape(r.H',28,28,k))
## Turn this into a method
function nmf(A,k)
  W, H = NMF.nndsvd(A, k; variant= :ar)
  alginst = NMF.ALSPGrad{Float64}(maxiter=100)
  r = NMF.solve!(alginst, A, W, H)
  p = plotimgs(reshape(r.H',28,28,k))
  savefig("nmf-mnist-$k.pdf")
  return p, r
end
p5, r5 = nmf(A,5)
p5

## Run a set of these
p8, r8 = nmf(A,8)
p10, r10 = nmf(A,10)
p15, r15 = nmf(A,15)
p25, r25 = nmf(A,25)
p35, r35 = nmf(A,35)
p50, r50 = nmf(A,50)
p100, r100 = nmf(A,100)

## Use these to re-assemble digits
ndig = 60
Xs = reshape((r10.W*r10.H)[1:ndig,:]',28,28,ndig)
As = reshape(A[1:ndig,:]',28,28,ndig)
anim = @animate for i=1:60
  plot(myheatmap(Xs[:,:,i]),myheatmap(As[:,:,i]))
end
gif(anim, "nmf-mnist-10-approx.gif",fps=3)

##
ndig = 60
Xs = reshape((r50.W*r50.H)[1:ndig,:]',28,28,ndig)
As = reshape(A[1:ndig,:]',28,28,ndig)
anim = @animate for i=1:60
  plot(myheatmap(Xs[:,:,i]),myheatmap(As[:,:,i]))
end
gif(anim, "nmf-mnist-50-approx.gif",fps=3)


text:
module BridgeSDEInference

using Bridge, StaticArrays, Distributions
using Statistics, Random, LinearAlgebra
using ForwardDiff
using ForwardDiff: value

# fitzHughNagumo.jl
export FitzhughDiffusion, FitzhughDiffusionAux, ℝ
export regularToAlter, alterToRegular, regularToConjug, conjugToRegular, display

# types.jl
export ImproperPrior, NoChangePt, SimpleChangePt

# mcmc.jl
export mcmc, PartObs, FPT, FPTInfo, ConjugateUpdt, MetropolisHastingsUpdt

# ODE solvers:
export Ralston3, RK4, Tsit5, Vern7

# euler_maruyama_dom_restr.jl
export forcedSolve!, forcedSolve

# random_walk.jl
export RandomWalk

# priors.jl
export Priors

# blocking_schedule.jl
export NoBlocking, ChequeredBlocking

# starting_pt.jl
export KnownStartingPt, GsnStartingPt

# radial_ornstein_uhlenbeck.jl
export RadialOU, RadialOUAux

# lorenz_system.jl
export Lorenz, LorenzAux

# lorenz_system_const_vola.jl
export LorenzCV, LorenzCVAux

# adaptation.jl
export Adaptation, NoAdaptation

# setup.jl
export MCMCSetup, set_observations!, set_imputation_grid!, set_transition_kernels!
export set_priors!, set_mcmc_params!, set_blocking!, set_solver!, initialise!

export Workspace

include("types.jl")

_DIR = "stochastic_process"
include(joinpath(_DIR, "bounded_diffusion_domain.jl"))
include(joinpath(_DIR, "guid_prop_bridge.jl"))
include(joinpath(_DIR, "path_to_wiener.jl"))

_DIR = "solvers"
include(joinpath(_DIR, "vern7.jl"))
include(joinpath(_DIR, "tsit5.jl"))
include(joinpath(_DIR, "rk4.jl"))
include(joinpath(_DIR, "ralston3.jl"))
include(joinpath(_DIR, "euler_maruyama_dom_restr.jl"))

_DIR = "transition_kernels"
include(joinpath(_DIR, "random_walk.jl"))

_DIR = "mcmc_extras"
include(joinpath(_DIR, "adaptation.jl"))
include(joinpath(_DIR, "blocking_schedule.jl"))
include(joinpath(_DIR, "first_passage_times.jl"))
include(joinpath(_DIR, "starting_pt.jl"))

_DIR = "examples"
include(joinpath(_DIR, "fitzhugh_nagumo.jl"))
include(joinpath(_DIR, "radial_ornstein_uhlenbeck.jl"))
include(joinpath(_DIR, "lorenz_system.jl"))
include(joinpath(_DIR, "lorenz_system_const_vola.jl"))

_DIR = "mcmc"
include(joinpath(_DIR, "priors.jl"))
include(joinpath(_DIR, "setup.jl"))
include(joinpath(_DIR, "workspace.jl"))
include(joinpath(_DIR, "conjugate_updt.jl"))
include(joinpath(_DIR, "mcmc.jl"))
include(joinpath(_DIR, "repeated.jl"))

end


text:
using LinearAlgebra
using StaticArrays
using Plots; pyplot()
using FastClosures
using NPZ
using LaTeXStrings

include("MagnusExpansion.jl")
include("AuxFunctions.jl")
include("HarmonicCorrectionSnap.jl")
include("snap.jl")

function run(tf)
	σ = pauli_matrices()
	χ_GHz = (-8281.0, 0.0, 0.0)
	χ = χ_GHz./abs(χ_GHz[1])
	N = 10
	levels = SA[0, 4]
	α = (pi, pi)./2
	Ncorr = 10

	system_params = Dict("chi"=>χ, "N"=>N, "levels"=>levels, "alpha"=>α,
						 "tf"=>tf, "Ncorr"=>Ncorr)

	Identity = Matrix{Complex{Float64}}(I, 2*N, 2*N)
	Uideal = U(system_params, σ; ideal_dynamics=true, matrix=true)
	Ur = U(system_params, σ; ideal_dynamics=false, matrix=true)
	F0 = fidelity(Ur, Uideal, 2*N)
	# "First" order correction
	magnus_order = 1
	max_order = 2
	coeffs = get_coeffs(magnus_order, max_order, Ncorr, system_params)

	Ucorr = U(system_params, σ, coeffs, max_order; frame="qubit", matrix=true)
	F1 = fidelity(Ucorr, Uideal, 2*N)
	# "Second" order correction
	magnus_order = 2
	coeffs2 = get_coeffs(coeffs, magnus_order, max_order, Ncorr, system_params)

	coeffs2 = coeffs + coeffs2
	Ucorr = U(system_params, σ, coeffs2, max_order; frame="qubit", matrix=true)
	F2 = fidelity(Ucorr, Uideal, 2*N)

	return F0, F1, F2
end

npoints = 100
tf0, tf1 = 15, 200
r = (tf1/tf0)^(1/(npoints-1))
tf_array = tf0*r.^collect(0: npoints-1)
data = zeros(npoints, 4)
for i = 1: npoints
	tf = tf_array[i]
	println(tf)
	data[i, 1] = tf
	data[i, 2: end] .= run(tf)
	println(data[i, 2: end])
end

x = data[:,1]
y = log10.(1 .- data[:, [2, 3, 4]])
plot(x, y,
	 linewidth = 2,
	 labels = ["No\ncorrection" "2nd order\nMagnus" "4th order\nMagnus"])
plot!(tickfont = font(12, "Serif"),
      xlabel = "gate time "*L"|\chi| t_\mathrm{f}",
      ylabel = "Fidelity error "*L"\varepsilon",
      guidefont = font(14, "Serif"),
	  legendfont = font(12, "Serif"),
	  background_color_legend = false,
      margin = 5Plots.mm)
savefig("fidelity_error.pdf")
npzwrite("data.npy", data)

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
function elementary_rules(rule::UInt8)
    rule_map = zeros(Bool, 2^3)
    for i = 1:length(rule_map)
        rule_map[i] = rule >> (i - 1) & 1
    end
    rule_map
end
function elementary_reversible_rules(rule::UInt8)
    rule_map = zeros(Bool, 2^4)
    for p = 0:1
        for i = 1:2^3
            rule_map[p * 8 + i] = rule >> (i - 1) & 1
            if p == 1
                rule_map[p * 8 + i] = !rule_map[p * 8 + i]
            end
        end
    end
    rule_map
end

export elementary
function elementary(rule::Int, steps::Int, initial_state::Vector{Bool}=Bool[ true ])
    return elementary(UInt8(rule), Int64(steps), initial_state)
end
function elementary(rule::UInt8, steps::Int64, initial_state::Vector{Bool}=Bool[ true ])
    rule_map = elementary_rules(rule)

    current_state = initial_state
    uns = false
    states = AutomataBinaryHistory()

    push!(states, current_state, uns)

    for i = 2:steps
        next_state = zeros(Bool, length(current_state) + 2)
        for j = 1:length(next_state)
            next_state[j] = rule_map[state_at(current_state, uns, j-2) << 2 | state_at(current_state, uns, j-1) << 1 | state_at(current_state, uns, j) + 1]
        end
        uns = rule_map[uns ? 8 : 1]
        current_state = next_state
        push!(states, current_state, uns)
    end

    return states
end
function elementary_reversible(rule::UInt8, steps::UInt64, initial_state::AbstractArray{Bool, 1}=Bool[ true ], initial_state_above::AbstractArray{Bool, 1}=Bool[ false, false, false ])
    rule_map = elementary_reversible_rules(rule)

    current_state = initial_state
    uns = false
    pre_uns = false
    states = AutomataBinaryHistory()

    push!(states, current_state, uns)

    for i = 2:steps
        next_state = zeros(Bool, length(current_state) + 2)
        for j = 1:length(next_state)
            next_state[j] = rule_map[(i == 2 ? state_at(initial_state_above, pre_uns, j) : state_at(states[Int(i-2)][1], pre_uns, j-2)) << 3 | state_at(current_state, uns, j-2) << 2 | state_at(current_state, uns, j-1) << 1 | state_at(current_state, uns, j) + 1]
        end
        pre_uns = uns
        uns = rule_map[pre_uns << 3 | uns << 2 | uns << 1 | uns + 1]
        current_state = next_state
        push!(states, current_state, uns)
    end

    return states
end


text:
#InspectDR: Compute grid/tick locations
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#Allowed mantissa values for grid step size:
const GRID_MANTSTEPS = DReal[1, 2, 2.5, 5, 10]
const GRID_MANTSTEP_EXPOFFSET = [0, 0, -1, 0, 0] #Ensure all significant digits show up.
const GRID_MANTTHRESH = (GRID_MANTSTEPS[1:end-1] + GRID_MANTSTEPS[2:end]) / 2

#Maps desired # of minor grid lines to optimal divisor for each possibility in GRID_MANTSTEPS
#TODO: Figure out algorithmically?
#GRID_MINORDIV[index(GRID_MANTSTEPS), DESIRED_NUMBER_OF_MINOR_GRIDLINES]
#NOTE: dividing major step size by 5 gives 4 grid lines (off-by-one)
const GRID_MINORDIV = DReal[
#	1  2  3  4  5  6  7  8  9  - desired # of grid lines
	2  2  4  5  5  5  10 10 10; #1
	2  2  4  4  4  8  8  8  10; #2
	1  1  5  5  5  5  10 10 10; #2.5
	2  2  5  5  5  5  10 10 10; #5
	2  2  4  5  5  5  10 10 10; #10 (matches 1)
]

#==Types
===============================================================================#

abstract type AbstractGridLines end

#Identifies where to place ticks/grid lines
mutable struct GridLines <: AbstractGridLines
	displaymajor::Bool
	displayminor::Bool
	major::Vector{DReal}
	minor::Vector{DReal}
	rnginfo::RangeDisplayInfo
end
GridLines(displaymajor=true, displayminor=true) =
	GridLines(displaymajor, displayminor, [], [], NoRangeDisplayInfo())

mutable struct GridLinesAuto <: AbstractGridLines
	displaymajor::Bool
	displayminor::Bool
end

mutable struct UndefinedGridLines <: AbstractGridLines
	minline::DReal
	maxline::DReal
end

abstract type PlotGrid end

#Rectilinear grid (ex: normal cartesian +logarithmic, ...):
mutable struct GridRect <: PlotGrid
	xlines::AbstractGridLines
	ylines::AbstractGridLines
#TODO: Make it possible for user to specify grids/labels as well.
end
GridRect(;vmajor=false, vminor=false, hmajor=false, hminor=false) =
	GridRect(GridLinesAuto(vmajor, vminor), GridLinesAuto(hmajor, hminor))

mutable struct GridSmith <: PlotGrid
	zgrid::Bool #Set to false for Y (admittance)-grid (x -> -x)
	ref::Float64 #Scaling factor (Zref/Yref)
	majorR::Vector{DReal} #Major constant resistance/admittance circles
	minorR::Vector{DReal} #Minor constant resistance/admittance circles
	minorX::Vector{DReal} #Minor constant reactance circles
	labelR::Vector{DReal}
	labelX::Vector{DReal}
end

#TODO: Curvilinear grid (ex: polar plots):
mutable struct GridCurv <: PlotGrid
end

#==Constructor-like functions
===============================================================================#
InputXfrm2D(xs::AxisScale, ys::AxisScale, g::PlotGrid) = InputXfrm2D(xs, ys) #Grid usually does not matter.

#Input data must be split into real/imag components:
InputXfrm2D(xs::AxisScale, ys::AxisScale, g::GridSmith) = InputXfrm2D(InputXfrm1DSpec{:real}(), InputXfrm1DSpec{:imag}())

function GridSmith(zgrid::Bool, ref::Float64)
	labelbase = DReal[0.2, 0.4, 0.6, 2, 4, 10] #WANTCONST
	minorextra = DReal[0.8, 1.5, 3, 6, 20] #WANTCONST
	#TODO: change with ext??
	majorR = DReal[0, 1]
	minorR = vcat(labelbase, minorextra)
	minorX = vcat(labelbase, minorextra, DReal[1])
	labelR = vcat(labelbase, DReal[0, 1])
	labelX = vcat(labelbase, DReal[1, 0.8, 1.5])
	return GridSmith(zgrid, ref, majorR, minorR, minorX, labelR, labelX)
end

function GridSmith(t::Symbol, ref::Float64)
	if !in(t, [:Z, :Y])
		MethodError(GridSmith, "Unrecognized Smith chart: :$t.")
	end
	return GridSmith(:Z == t, ref)
end
GridSmith(t::Symbol; ref::Number=1.0) = GridSmith(t, Float64(ref))

#==Main algorithm
===============================================================================#

#Generate grid step size resulting in pretty tick labels
#tgtminor: Targeted number of minor grid lines/ticks
function linearstep_pretty(tgtstep::DReal, tgtminor::Int)
	#TODO: Rounding errors in calculation of major cause display issues with tick labels

	lstep = log10(tgtstep)
	ilstep = floor(Int, lstep)
	stepexp = 10.0^(ilstep)
	stepmant = tgtstep / stepexp

	istep = length(GRID_MANTSTEPS)
	for i in 1:length(GRID_MANTTHRESH)
		thresh = GRID_MANTTHRESH[i]
		if stepmant < thresh
			istep = i
			break
		end
	end

	mstep = GRID_MANTSTEPS[istep]
	expoffset = GRID_MANTSTEP_EXPOFFSET[istep]
	tgtminor = clamp(tgtminor, 0, 9)
	major = mstep*stepexp
	minor = major
	if tgtminor > 0
		minor /= GRID_MINORDIV[istep, tgtminor]
	end

	return (major, minor, expoffset)
end

#Compute grid line configuration for a linear scale (dB scale is basically linear as well)
#tgtmajor: Targeted number of major grid lines
#tgtminor: Targeted number of minor grid lines
function gridlines(scale::AxisScale, min::DReal, max::DReal, displaymajor::Bool, displayminor::Bool, tgtmajor::DReal, tgtminor::Int)
	if !isfinite(min) || !isfinite(max)
		return UndefinedGridLines(min, max)
	elseif abs(max - min) == 0 #TODO: specify Δmin?
		return UndefinedGridLines(min, max)
	end

	grd = GridLines(displaymajor, displayminor)
	span = max - min

	#Compute major grid lines:
	tgtstep = span / (tgtmajor+1)
	(stepmajor, stepminor, expoffset) = linearstep_pretty(tgtstep, tgtminor)
	major1 = ceil(Int, min/stepmajor)*stepmajor
	grd.major = collect(major1:stepmajor:max)
	grd.rnginfo = LinearRangeDisplayInfo(
		base10exp(major1), base10exp(grd.major[end]),
		base10exp(stepmajor) + expoffset
	)

	#Compute minor grid lines:
	minor1 = ceil(Int, min/stepminor)*stepminor
	minorall = minor1:stepminor:max

	#Throw away major points...
	#Yuk... is there not a more elegant solution?
	major_i1 = round(Int, abs(major1-minor1) / stepminor)+1 #Index of first major grid line
	major_ipitch = round(Int, abs(stepmajor/stepminor)) #Index pitch for major grid lines
	npts = length(minorall)
	if major_i1 <= npts
		nredundant = div(npts-major_i1, major_ipitch)+1
		npts = npts-nredundant
	else #Just in case
		npts = 0
	end
	grd.minor = Array{DReal}(undef, npts)

	nextmajor = major_i1
	j = 1
	for i in 1:length(minorall)
		if nextmajor == i
			nextmajor += major_ipitch
		else
			grd.minor[j] = minorall[i]
			j+=1
		end
	end

	return grd
end

#Compute grid line configuration for a log10 scale
#Ignore keywords
function gridlines(scale::LogScale{10}, logmin::DReal, logmax::DReal, displaymajor::Bool, displayminor::Bool; kwargs...)
	if !isfinite(logmin) || !isfinite(logmax)
		return UndefinedGridLines(logmin, logmax)
	elseif abs(logmax - logmin) == 0 #TODO: specify Δmin?
		return UndefinedGridLines(logmin, logmax)
	end
	grd = GridLines(displaymajor, displayminor)

	major1 = ceil(Int, logmin)
	majorend = floor(Int, logmax)
	grd.major = collect(major1:1:majorend)

	offsets = collect(log10.(2:1:9))

	if length(grd.major) > 0
		minor = (grd.major[1]-1) .+ offsets

		for curline in grd.major
			append!(minor, curline .+ offsets)
		end

		#Strip out values beyond specified range:
		istart = 1; iend = 0
		for i in 1:length(minor)
			if minor[i] < logmin
				istart = i+1
			end
			if minor[i] < logmax
				iend = i
			end
		end

		grd.minor = minor[istart:iend]
	end

	return grd
end

#Compute grid line configuration for a ln scale
function gridlines(scale::LogScale{:e}, min::DReal, max::DReal, displaymajor::Bool, displayminor::Bool)
	#TODO: Limit major grid lines to step sizes >= e^1?
	return GridLines(displaymajor, displayminor)
end

#Compute grid line configuration for a log2 scale
function gridlines(scale::LogScale{2}, min::DReal, max::DReal, displaymajor::Bool, displayminor::Bool)
	#TODO: Limit major grid lines to step sizes >= 2^1?
	return GridLines(displaymajor, displayminor)
end

#==High-level interface
===============================================================================#

#Evaluate grid lines
#-------------------------------------------------------------------------------
function _eval(alines::GridLinesAuto, s::AxisScale, emin::DReal, emax::DReal)
	return gridlines(s, emin, emax, alines.displaymajor, alines.displayminor)
end
function _eval(alines::GridLinesAuto, s::LinScale, emin::DReal, emax::DReal)
	lines = gridlines(s, emin, emax,
		alines.displaymajor, alines.displayminor, s.tgtmajor, s.tgtminor
	)
	return lines
end

#Default: use specified grid lines:
_eval(lines::AbstractGridLines, s::AxisScale, emin::DReal, emax::DReal) = lines

#Default: Nothing to evaluate (use grid as-specified).
_eval(grid::PlotGrid, xs::AxisScale, ys::AxisScale, ext::PExtents2D) = grid

function _eval(grid::GridRect, xs::AxisScale, ys::AxisScale, ext::PExtents2D)
	xlines = _eval(grid.xlines, xs, ext.xmin, ext.xmax)
	ylines = _eval(grid.ylines, ys, ext.ymin, ext.ymax)
	return GridRect(xlines, ylines)
end

#Get grid used to extract coordinates, etc on a graph:
coord_grid(grid::GridRect, xs::AxisScale, ys::AxisScale, ext::PExtents2D) = _eval(grid, xs, ys, ext)
function coord_grid(grid::GridSmith, xs::AxisScale, ys::AxisScale, ext::PExtents2D)
	autogrid = GridRect(vmajor=false, vminor=false, hmajor=false, hminor=false)
	return _eval(autogrid, xs, ys, ext)
end

#Last line


text:
"""
    cll()

Clears the current legend.
"""
cll() = (reset!(gca().legend); PreviewFigure(gcf()))

"""
    $SIGNATURES

Update the properties of an existing legend object present on `axes`. If none
exist then a new one is created with the given properties.
"""
function legend(vd::Option{Vector{DrawingHandle{T}}}=nothing,
                labels::Option{Vector{String}}=nothing;
                axes=nothing, opts...) where {T}
    axes isa Axes3D && throw(NotImplementedError("Legend for Axes3D"))
    axes=check_axes(axes)
    # create a new legend object
    axes.legend = Legend()
    # check if handles are given, if none, try to apply a default legend
    if isnothing(vd)
        isnothing(labels) || throw(ArgumentError("Cannot pass labels without handles."))
        axes.legend.handles = [DrawingHandle(d) for d ∈ axes.drawings]
        axes.legend.labels  = fill("", length(axes.drawings))

        # check if any object has a label, otherwise fill with default
        onehaslabel = false
        for (k,d) ∈ enumerate(axes.drawings)
            if isa(d, Union{Scatter2D,Bar2D})
                isempty(d.labels) || (axes.legend.labels[k] = d.labels; onehaslabel = true)
            elseif isa(d, Union{Fill2D, Hist2D})
                isempty(d.label)  || (axes.legend.labels[k] = d.label; onehaslabel = true)
            end
        end
        if onehaslabel
            # remove the entries that don't have a label
            mask = @. !(isempty(axes.legend.labels))
            axes.legend.handles = axes.legend.handles[mask]
            axes.legend.labels  = axes.legend.labels[mask]
        else
            s_ctr = 1
            f_ctr = 1
            h_ctr = 1
            b_ctr = 1
            for (k,d) ∈ enumerate(axes.drawings)
                if d isa Scatter2D
                    axes.legend.labels[k] = ["plot $(s_ctr+e-1)" for e ∈ 1:d.nobj]
                    s_ctr += d.nobj
                elseif d isa Fill2D
                    axes.legend.labels[k] = "fill $f_ctr"
                    f_ctr += 1
                elseif d isa Hist2D
                    axes.legend.labels[k] = "hist $h_ctr"
                    h_ctr += 1
                elseif d isa Bar2D
                    axes.legend.labels[k] = ["bar $(b_ctr+e-1)" for e ∈ 1:d.nobj]
                    b_ctr += d.nobj
                end
            end
        end
    # handles are given, check the labels
    else
        isnothing(labels) && throw(ArgumentError("Labels must be provided along with handles"))
        length(vd) == length(labels) || throw(ArgumentError("There must be as many labels given "*
                                                            "as drawing handles"))
        axes.legend.handles = vd
        axes.legend.labels  = labels
    end
    set_properties!(axes.legend; defer_preview=true, opts...)
    return PreviewFigure(gcf())
end
legend(d::DrawingHandle{T}, s::String; o...) where T = legend([d],[s]; o...)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Furlongs

export Furlong

# Here we implement a minimal dimensionful type Furlong, which is used
# to test dimensional correctness of various functions in Base.

# represents a quantity in furlongs^p
struct Furlong{p,T<:Number} <: Number
    val::T
    Furlong{p,T}(v::Number) where {p,T} = new(v)
end
Furlong(x::T) where {T<:Number} = Furlong{1,T}(x)
Furlong(x::Furlong) = x
(::Type{T})(x::Furlong) where {T<:Number} = T(x.val)::T
Furlong{p}(v::Number) where {p} = Furlong{p,typeof(v)}(v)
Furlong{p}(x::Furlong{q}) where {p,q} = (@assert(p==q); Furlong{p,typeof(x.val)}(x.val))
Furlong{p,T}(x::Furlong{q}) where {T,p,q} = (@assert(p==q); Furlong{p,T}(T(x.val)))

Base.promote_type(::Type{Furlong{p,T}}, ::Type{Furlong{p,S}}) where {p,T,S} =
    Furlong{p,promote_type(T,S)}

Base.one(x::Furlong{p,T}) where {p,T} = one(T)
Base.one(::Type{Furlong{p,T}}) where {p,T} = one(T)
Base.oneunit(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(one(T))
Base.oneunit(x::Type{Furlong{p,T}}) where {p,T} = Furlong{p,T}(one(T))
Base.zero(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(zero(T))
Base.zero(::Type{Furlong{p,T}}) where {p,T} = Furlong{p,T}(zero(T))
Base.iszero(x::Furlong) = iszero(x.val)
Base.float(x::Furlong{p}) where {p} = Furlong{p}(float(x.val))
Base.eps(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(eps(T))
Base.eps(::Furlong{p,T}) where {p,T<:AbstractFloat} = eps(Furlong{p,T})
Base.floatmin(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmin(T))
Base.floatmin(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmin(Furlong{p,T})
Base.floatmax(::Type{Furlong{p,T}}) where {p,T<:AbstractFloat} = Furlong{p}(floatmax(T))
Base.floatmax(::Furlong{p,T}) where {p,T<:AbstractFloat} = floatmax(Furlong{p,T})
Base.conj(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(conj(x.val))

# convert Furlong exponent p to a canonical form
canonical_p(p) = isinteger(p) ? Int(p) : Rational{Int}(p)

Base.abs(x::Furlong{p}) where {p} = Furlong{p}(abs(x.val))
Base.abs2(x::Furlong{p}) where {p} = Furlong{canonical_p(2p)}(abs2(x.val))
Base.inv(x::Furlong{p}) where {p} = Furlong{canonical_p(-p)}(inv(x.val))

for f in (:isfinite, :isnan, :isreal, :isinf)
    @eval Base.$f(x::Furlong) = $f(x.val)
end
for f in (:real,:imag,:complex,:+,:-)
    @eval Base.$f(x::Furlong{p}) where {p} = Furlong{p}($f(x.val))
end

import Base: +, -, ==, !=, <, <=, isless, isequal, *, /, //, div, rem, mod, ^
for op in (:+, :-)
    @eval function $op(x::Furlong{p}, y::Furlong{p}) where {p}
        v = $op(x.val, y.val)
        Furlong{p}(v)
    end
end
for op in (:(==), :(!=), :<, :<=, :isless, :isequal)
    @eval $op(x::Furlong{p}, y::Furlong{p}) where {p} = $op(x.val, y.val)
end
for (f,op) in ((:_plus,:+),(:_minus,:-),(:_times,:*),(:_div,://))
    @eval function $f(v::T, ::Furlong{p}, ::Union{Furlong{q},Val{q}}) where {T,p,q}
        s = $op(p, q)
        Furlong{canonical_p(s),T}(v)
    end
end
for (op,eop) in ((:*, :_plus), (:/, :_minus), (://, :_minus), (:div, :_minus))
    @eval begin
        $op(x::Furlong{p}, y::Furlong{q}) where {p,q} =
            $eop($op(x.val, y.val),x,y)
        $op(x::Furlong{p}, y::S) where {p,S<:Number} = $op(x,Furlong{0,S}(y))
        $op(x::S, y::Furlong{p}) where {p,S<:Number} = $op(Furlong{0,S}(x),y)
    end
end
# to fix an ambiguity
//(x::Furlong, y::Complex) = x // Furlong{0,typeof(y)}(y)
for op in (:rem, :mod)
    @eval begin
        $op(x::Furlong{p}, y::Furlong) where {p} = Furlong{p}($op(x.val, y.val))
        $op(x::Furlong{p}, y::Number) where {p} = Furlong{p}($op(x.val, y))
    end
end
Base.sqrt(x::Furlong) = _div(sqrt(x.val), x, Val(2))

end

max_stars_repo_name: Marlin-Na/julia

max_issues_repo_name: Marlin-Na/julia

max_forks_repo_name: Marlin-Na/julia


text:
using LTC
using Plots
gr()
using BenchmarkTools
using DiffEqSensitivity
using OrdinaryDiffEq
using DiffEqFlux
using GalacticOptim

include("half_cheetah_data_loader.jl")

function traintest(n, solver=VCABM(), sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
	function loss(x,y,m)
	  Flux.reset!(m)
	  #m = re(θ)
	  #ŷ = m.(x)
	  ŷ = map(xi -> m(xi), x)
	  #ŷ = [m(xi)[end-m.cell.wiring.n_motor+1:end, :] for xi in x]
	  #ŷ = m.(x)
	  #length(ŷ) < length(y) && return Inf
	  #sum(Flux.Losses.mse.(ŷ, y; agg=mean))
	  sum(sum([(ŷ[i] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y)#, ŷ
	  #sum([sum((ŷ[i] .- y[i]) .^ 2) for i in 1:length(y)])/length(y)
	end

  function cb(x,y,l,m)
    println(l)
    # pred = m.(x)
    # # isnan(l) && return false
    # fig = plot([ŷ[1,1] for ŷ in pred])
	#   plot!(fig, [ŷ[end,1] for ŷ in pred])
    # plot!(fig, [yi[1,1] for yi in y])
	#   plot!(fig, [yi[end,1] for yi in y])
    # display(fig)
    return false
  end

  function lg(p,x,y)
	m = model
	Flux.reset!(m)
	ŷ = map(xi -> m(xi,p), x)
    sum(sum([(ŷ[i] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y), ŷ
  end
  cbg = function (p,l,pred;doplot=false) #callback function to observe training
   display(l)
    # plot current prediction against data
    if doplot
  	  pl = scatter(first(train_dl)[1],label="data")
	  scatter!(pl,pred[1,:],label="prediction")
	  display(plot(pl))
	end
	return false
  end

  train_dl, test_dl, valid_dl = get_dl(batchsize=32, seq_len=32)
  ncp = LTC.LTCNet(NCPWiring(17,2,
          n_sensory=4, n_inter=4, n_command=7, n_motor=2,
          sensory_in=-1, rec_sensory=0, sensory_inter=2, sensory_command=0, sensory_motor=0,
          inter_in=2, rec_inter=2, inter_command=3, inter_motor=1,                       # inter_in = sensory_out
          command_in=0, rec_command=4, command_motor=2,                   # command_in = inter_out
          motor_in=0, rec_motor=3), solver, sensealg)
  model1 = ncp
  model = ncp

  model2 = DiffEqFlux.FastDense(2,17)
  #model = (x,p) -> model1(x,p)[end-16:end, :]
  model = DiffEqFlux.FastChain(model1, (x,p) -> x[end-1:end,:], model2)
  # model2 = Dense(2,17)
  # model = Chain(model1, x -> x[end-1:end, :], model2)
  # model = model1
  # model = DiffEqFlux.FastChain(ncp,(x,p) -> x[end-1:end, :],DiffEqFlux.FastDense(2,17))
  #θ = Flux.params(model1,model2)
  # θ = Flux.params(model)
  # @show sum(length.(θ))
  pp = initial_params(model)
  @show length(pp)
  # pp = DiffEqFlux.initial_params(model)
  #lower,upper = get_bounds(model)
  lower,upper = [],[]
  @show length(lower)

  # @show length(train_dl)
  # @show size(first(train_dl)[1])
  # @show size(first(train_dl)[1][1])
  #
  # display(Plots.heatmap(ncp.cell.wiring.sens_mask))
  # display(Plots.heatmap(ncp.cell.wiring.syn_mask))
  #
  # display(Plots.heatmap(ncp.cell.wiring.sens_pol))
  # display(Plots.heatmap(ncp.cell.wiring.syn_pol))

  #@show sum([length(p) for p in θ])
  #@show length(lower)
  # @show length(pp)

  # model.(first(train_dl)[1])
  # model.(first(train_dl)[1])
  # model.(first(train_dl)[1])

  opt = Flux.Optimiser(ClipValue(0.5), ADAM(0.05))

  optfun = OptimizationFunction((θ, p, x, y) -> lg(θ,x,y), GalacticOptim.AutoZygote())
  optprob = OptimizationProblem(optfun, pp)
  #using IterTools: ncycle
  res1 = GalacticOptim.solve(optprob, opt, train_dl, cb = cbg, maxiters = n)
  # return res1

  #my_custom_train!(model, (x,y) -> loss(x,y,model), θ, train_dl, opt; cb, lower, upper)
  # Flux.@epochs n my_custom_train!(model, (x,y) -> loss(x,y,model), θ, train_dl, opt; cb, lower, upper)
end

# @time traintest(1, VCABM(), ZygoteAdjoint())

@time traintest(10)

@time traintest(10, VCABM(), ForwardDiffSensitivity())
@time traintest(10, VCABM(), ReverseDiffAdjoint())
@time traintest(10, VCABM(), QuadratureAdjoint(autojacvec=ReverseDiffVJP(true)))
# 192.343950 seconds (206.01 M allocations: 64.016 GiB, 4.65% gc time, 22.36% compilation time)
# 138.803175 seconds (159.75 M allocations: 51.048 GiB, 5.04% gc time)
# 89.533896 seconds (346.13 M allocations: 138.572 GiB, 23.55% gc time)
@time traintest(10, AutoTsit5(Rosenbrock23()), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

@time traintest(1, VCABM(), InterpolatingAdjoint(checkpointing=true, autojacvec=ZygoteVJP()))
@time traintest(1, VCABM(), TrackerAdjoint())
@time traintest(1, VCABM(), InterpolatingAdjoint(autojacvec=ZygoteVJP()))
@time traintest(1, VCABM(), QuadratureAdjoint(compile=true, abstol=1e-3,reltol=1e-3, autojacvec=ZygoteVJP()))
@time traintest(1, VCABM(), QuadratureAdjoint(compile=true, abstol=1e-3,reltol=1e-3, autojacvec=ReverseDiffVJP(true)))
@time traintest(1, Tsit5(), ZygoteAdjoint())
# 142.623700 seconds (1.43 G allocations: 380.410 GiB, 45.09% gc time, 0.02% compilation time)

@time traintest(1, VCABM(), ReverseDiffAdjoint())
@time traintest(1, CVODE_BDF(linear_solver=:GMRES), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
@time traintest(1, VCABM(), InterpolatingAdjoint(checkpointing=true, autojacvec=ReverseDiffVJP(true)))
# 161.160516 seconds (985.58 M allocations: 201.055 GiB, 22.29% gc time, 9.22% compilation time)

#@time traintest(3, VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))

#@time traintest(3, Rosenbrock23(), InterpolatingAdjoint(checkpointing=false))

#@time traintest(3, ImplicitEuler(), InterpolatingAdjoint(autojacvec=ZygoteVJP()))
#@time traintest(3, VCABM(), InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))


text:
"""
returns Core-type relperms for gas and liquid as a
function of the liquid saturation and their derivatives
with respect to the liquid saturation; it also returns
the stars model for the gas relperm in presence of foam.
inputs [unit] \\(default):
   + swc [-] \\(0.1)
   + sgr [-] \\(0.05)
   + krg0 [-] \\(0.9)
   + ng [-] \\(1.8)
   + krw0 [-] \\(0.2)
   + nw [-] \\(4.2)
   + fmmob [-] \\(25000)
   + epdry [-] \\(10000)
   + fmdry [-] \\(0.29)
Outputs:
   + krw
   + krg
   + dkrwdsw
   + dkrgdsw
   + krgf
   + dkrgfdsw
Usage:
   + (krw, krg, dkrwdsw, dkrgdsw, krgf, dkrgfdsw)=corey_rel_perms()
   + (krw, krg, dkrwdsw, dkrgdsw, krgf, dkrgfdsw)=corey_rel_perms(swc,
sgr, krg0, ng, krw0, nw, fmmob, epdry, fmdry)
"""
function corey_rel_perms(;swc=0.1, sgr=0.05, krg0=0.9, ng=1.8, krw0=0.2,
nw=4.2, fmmob=25000.0, epdry=10000.0, fmdry=0.29)
  sws(sw::Real)=((sw>swc)*(sw<1-sgr)*(sw-swc)/(1-sgr-swc)+(sw>=1-sgr)*1.0)
  sws(sw::Array{Float64})=((sw.>swc).*(sw.<1-sgr).*(sw-swc)/(1-sgr-swc)+(sw.>=1-sgr).*ones(size(sw)))
  krg(sw)=((sw.>=swc).*krg0.*(1-sws(sw)).^ng+(sw.<swc).*(1+(krg0-1)/swc*sw))
  krw(sw)=((sw.<=1-sgr).*krw0.*sws(sw).^nw+(sw.>1-sgr).*(-(1-krw0)/sgr.*(1.0-sw)+1.0))
  fm(sw)=(1+fmmob*(0.5+atan(epdry.*(sw-fmdry))/π))
  krgf(sw)=(krg(sw)./fm(sw))
  dkrwdsw(sw)=(nw*krw0*(1/(1-sgr-swc))*sws(sw).^(nw-1))
  dkrgdsw(sw)=((-krg0*ng*(1-sws(sw)).^(ng-1))/(-swc-sgr+1))
  dfmdsw(sw)=(((epdry*fmmob)./(π*(epdry^2*(sw-fmdry).^2+1))))
  dkrgfdsw(sw)=((dkrgdsw(sw).*fm(sw)-dfmdsw(sw).*krg(sw))./fm(sw).^2)

  return krw, krg, dkrwdsw, dkrgdsw, krgf, dkrgfdsw
end

"""
fractional flow formulation and solution of a WAG process
for Corey-type gas and water relative permeabilities on a
Cartesian coordinate (not radial)
inputs [unit] \\(default):
   + muw [Pa.s] \\(1e-3)
   + mug [Pa.s] \\(2e-5)
   + ut [m/s] \\(1e-5)
   + phi [-] \\(0.2)
   + k [m^2] \\(1e-12)
   + swc [-] \\(0.1)
   + sgr [-] \\(0.05)
   + krg0 [-] \\(0.9)
   + ng [-] \\(1.8)
   + krw0 [-] \\(0.2)
   + nw [-] \\(4.2)
   + sw0 [-] \\(1.0)
   + sw_inj [-] \\(0.0)
   + L [m] \\(domain length = 1.0)
   + pv_inj [-] \\(injected pore volumes = 5)
Outputs:
   + xt_shock (x/t shock position)
   + sw_shock
   + xt_prf
   + sw_prf
   + t [s] \\time
   + p_inj [Pa] \\(injection pressure history)
Usage:
   using roots, Dierckx
   + (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj)=frac_flow_wag()
   + (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj)=frac_flow_wag(muw, mug, ut, phi,
k, swc, sgr, krg0, ng, krw0, nw, sw0, sw_inj, L, pv_inj)
"""
function frac_flow_wag(;muw=1e-3, mug=2e-5, ut=1e-5, phi=0.2,
  k=1e-12, swc=0.1, sgr=0.05, krg0=0.9, ng=1.8, krw0=0.2,
  nw=4.2, sw0=1.0, sw_inj=0.0, L=1, pv_inj=5)

  sws(sw::Real)=((sw>swc)*(sw<1-sgr)*(sw-swc)/(1-sgr-swc)+(sw>=1-sgr)*1.0)
  sws(sw::Array{Float64})=((sw.>swc).*(sw.<1-sgr).*(sw-swc)/(1-sgr-swc)+(sw.>=1-sgr).*ones(size(sw)))
  kr(sw)=((sw.>=swc).*krg0.*(1-sws(sw)).^ng+(sw.<swc).*(1+(krg0-1)/swc*sw))
  krw(sw)=((sw.<=1-sgr).*krw0.*sws(sw).^nw+(sw.>1-sgr).*(-(1-krw0)/sgr.*(1.0-sw)+1.0))
  krg(sw)=kr(sw)
  dkrwdsw(sw)=(nw*krw0*(1/(1-sgr-swc))*sws(sw).^(nw-1))
  dkrdsw(sw)=((-krg0*ng*(1-sws(sw)).^(ng-1))/(-swc-sgr+1))
  dkrgdsw(sw)=dkrdsw(sw)
  fw(sw)=((krw(sw)/muw)./(krw(sw)/muw+krg(sw)/mug))
  dfw(sw)=((dkrwdsw(sw)/muw.*(krw(sw)/muw+krg(sw)/mug)-
      (dkrwdsw(sw)/muw+dkrgdsw(sw)/mug).*krw(sw)/muw)./
      (krg(sw)/mug+krw(sw)/muw).^2)
    eps1=1e-3
    dfw_num(sw) = (fw(sw+eps1)-fw(sw-eps1))/(2*eps1)
    ftot(sw)=krg(sw)/mug+krw(sw)/muw

# solve the nl equation to find the shock front saturation
    f_shock(sw)=(dfw(sw)-(fw(sw0)-fw(sw))./(sw0-sw))
    sw_shock = fzero(f_shock, [swc+eps(),1-sgr-eps()])
    s=collect(linspace(0.0,1.0,100))
    s1 = collect(linspace(sw_inj, sw_shock, 500))
    xt_s1 = ut/phi*dfw(s1)
    xt_s = ut/phi*dfw(s)
    xt_shock = ut/phi*dfw(sw_shock)
    xt_prf=[xt_s1; xt_shock; xt_shock+eps(); 2*xt_shock]
    sw_prf=[s1; sw_shock; sw0; sw0]

# Go through the data first
    i=1
    while(true)
      if (i+1)==length(xt_prf)
        break
      elseif xt_prf[i]==xt_prf[i+1]
        deleteat!(xt_prf, i+1)
        deleteat!(sw_prf, i+1)
      else
        i+=1
      end
    end

# find the injection pressure history
    x = collect(linspace(0,L,1000))
    sw_int = Spline1D([xt_prf; L/eps()], [sw_prf; sw0], k=1);
    t_inj=pv_inj*phi*L/ut
    t = collect(linspace(eps(),t_inj, 100)) # [s] time
    p_inj = zeros(length(t))
    p_inj[1]=trapz(x, ut./(k*(krg(sw0*ones(size(x)))/mug+krw(sw0*ones(size(x)))/muw)))
    for i in 2:length(t)
        xt = x/t[i]
        p_inj[i] = trapz(x, ut./(k*(krg(sw_int(xt))/mug+krw(sw_int(xt))/muw)))
    end

# Return the results
    return xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj
end

"""
Foam flow in porous media: a BL formulation of SAG process.
CMG Stars foam model is used for the gas mobility in
presence of foam, and Corey-type relperms are used for
the gas and water. The equations are solved on a
Cartesian coordinate in one dimension.
inputs [unit] \\(default):
   + fmmob [-] \\(25000)
   + epdry [-] \\(10000)
   + fmdry [-] \\(0.29)
   + muw [Pa.s] \\(=1e-3)
   + mug [Pa.s] \\(2e-5)
   + ut [m/s] \\(1e-5)
   + phi [-] \\(0.2)
   + k [m^2] \\(1e-12)
   + swc [-] \\(0.1)
   + sgr [-] \\(0.05)
   + krg0 [-] \\(0.9)
   + ng [-] \\(1.8)
   + krw0 [-] \\(0.2)
   + nw [-] \\(4.2)
   + sw0 [-] \\(1.0)
   + sw_inj [-] \\(0.0)
Outputs:
   + xt_shock (x/t shock position)
   + sw_shock
   + xt_prf
   + sw_prf
Usage:
   using roots
   + (xt_shock, sw_shock, xt_prf, sw_prf)=frac_flow_sag()
   + (xt_shock, sw_shock, xt_prf, sw_prf)=frac_flow_sag(fmmob, epdry, fmdry, muw, mug, ut, phi,
k, swc, sgr, krg0, ng, krw0, nw, sw0, sw_inj)
"""
function frac_flow_sag(;fmmob=25000.0, epdry=10000.0, fmdry=0.29,
  muw=1e-3, mug=2e-5, ut=1e-5, phi=0.2,
  k=1e-12, swc=0.1, sgr=0.05, krg0=0.9, ng=1.8, krw0=0.2,
  nw=4.2, sw0=1.0, sw_inj=0.0, L=1, pv_inj=5)

    sws(sw::Real)=((sw>swc)*(sw<1-sgr)*(sw-swc)/(1-sgr-swc)+(sw>=1-sgr)*1.0)
    sws(sw::Array{Float64})=((sw.>swc).*(sw.<1-sgr).*(sw-swc)/(1-sgr-swc)+(sw.>=1-sgr).*ones(size(sw)))
    kr(sw)=((sw.>=swc).*krg0.*(1-sws(sw)).^ng+(sw.<swc).*(1+(krg0-1)/swc*sw))
    krw(sw)=((sw.<=1-sgr).*krw0.*sws(sw).^nw+(sw.>1-sgr).*(-(1-krw0)/sgr.*(1.0-sw)+1.0))
    fm(sw)=(1+fmmob*(0.5+atan(epdry.*(sw-fmdry))/π))
    krg(sw)=(kr(sw)./fm(sw))
    dkrwdsw(sw)=(nw*krw0*(1/(1-sgr-swc))*sws(sw).^(nw-1))
    dkrdsw(sw)=((-krg0*ng*(1-sws(sw)).^(ng-1))/(-swc-sgr+1))
    dfmdsw(sw)=(((epdry*fmmob)./(π*(epdry^2*(sw-fmdry).^2+1))))
    dkrgdsw(sw)=((dkrdsw(sw).*fm(sw)-dfmdsw(sw).*kr(sw))./fm(sw).^2)
    fw(sw)=((krw(sw)/muw)./(krw(sw)/muw+krg(sw)/mug))
    dfw(sw)=((dkrwdsw(sw)/muw.*(krw(sw)/muw+krg(sw)/mug)-
      (dkrwdsw(sw)/muw+dkrgdsw(sw)/mug).*krw(sw)/muw)./
      (krg(sw)/mug+krw(sw)/muw).^2)
    eps1=1e-3
    dfw_num(sw) = (fw(sw+eps1)-fw(sw-eps1))/(2*eps1)
    ftot(sw)=krg(sw)/mug+krw(sw)/muw

# solve the nl equation to find the shock front saturation
    f_shock(sw)=(dfw(sw)-(fw(sw0)-fw(sw))./(sw0-sw))
    sw_shock = fzero(f_shock, [swc+eps(),1-sgr-eps()])
    s=collect(linspace(0.0,1.0,100))
    s1 = collect([linspace(sw_inj, swc, 100); linspace(swc+eps(), sw_shock, 1000)])
    xt_s1 = ut/phi*dfw(s1)
    xt_s = ut/phi*dfw(s)
    xt_shock = ut/phi*dfw(sw_shock)
    xt_prf=[xt_s1; xt_shock; xt_shock+eps(); 2*xt_shock]
    sw_prf=[s1; sw_shock; sw0; sw0]

# Go through the data first
    i=1
    while(true)
      if (i+1)==length(xt_prf)
        break
      elseif xt_prf[i]==xt_prf[i+1]
        deleteat!(xt_prf, i+1)
        deleteat!(sw_prf, i+1)
      else
        i+=1
      end
    end

# find the injection pressure history
    x = collect(linspace(0,L,1000))
    sw_int = Spline1D([xt_prf; L/eps()], [sw_prf; sw0], k=1);
    t_inj=pv_inj*phi*L/ut
    t = collect(linspace(eps(),t_inj, 100)) # [s] time
    p_inj = zeros(length(t))
    p_inj[1]=trapz(x, ut./(k*(krg(sw0*ones(size(x)))/mug+krw(sw0*ones(size(x)))/muw)))
    for i in 2:length(t)
        xt = x/t[i]
        p_inj[i] = trapz(x, ut./(k*(krg(sw_int(xt))/mug+krw(sw_int(xt))/muw)))
    end

    return xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj
end

function trapz(x,y)
# copied from https://github.com/hwborchers/NumericalMath.jl/blob/master/src/integrate.jl
  n=length(x)
  if (length(y) != n)
    error("Vectors 'x', 'y' must be of same length")
  end
  r=0
  for i in 2:n
    r+=(x[i] - x[i-1])*(y[i] + y[i-1])
  end
  return r/2.0
end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
using Random

randstring(6)

max_stars_repo_name: loocus/StringGenerator

max_issues_repo_name: loocus/StringGenerator

max_forks_repo_name: loocus/StringGenerator


text:
using Pkg
Pkg.API.develop(Pkg.PackageSpec(name="SyringePumpJl", path="/home/caseykneale/Desktop/SyringePumpJl"))
using SyringePumpJl
using Blink, Interact, Measurements, Unitful, LibSerialPort

# mcu = nothing
# try
#     mcu = open( "/dev/ttyACM0", 9600 )
# catch
#
# end
mcu = open( "/dev/ttyACM0", 9600 )

driver          = InputField( "Move N Steps", 100.0,  nothing);
driver_btn      = Widgets.button("Drive Pump");
c               = Calibration("mm");
calibration_plt = Observable{Any}(nothing);
drive_status    = Observable{Any}(nothing);
map!( x ->  drive_pump( mcu, convert( Int64, value(driver) ),
            convert( Int64, 10000) ), drive_status, driver_btn)

map!( x -> makeplot( x ), calibration_plt, c.steps_and_disps)

calibrate_txt = "Drive the stepper motor N steps, and measure the displacement of the plunger with a caliper." *
                "Add the experimental values to the calibration model." *
                "After numerous runs, you will be given a calibration + uncertainty for dispensing fluids based on a linear model."

calibrationbox = vbox( "Calibrate:", calibrate_txt, hbox(driver.field, driver_btn),
                        calibration_plt, c.widget
                );

units       = UncertainInputField( "Dispense Rate", 50.0, 0.02, u"mm / mL");
volume      = InputField( "Volume", 5.0, u"mL");
duration    = InputField( "Duration", 3.0, u"s");
driver_btn2 = Widgets.button("Drive Pump");
stop_btn    = Widgets.button("Stop Pump");

map!( x -> drive_physical( mcu, c.steps_and_disps[], units, volume, duration ), drive_status, driver_btn2)

driverbox   = vbox( units.field, volume.field, duration.field,
                    hbox(driver_btn2, stop_btn) );
tabbed      = Widgets.tabulator(Observable( Dict( "Calibrate" => calibrationbox, "Driver" => driverbox ) ));

w = Window(); body!(w, tabbed)

println("fin")

max_stars_repo_name: JuliaMakers/MakerProjects

max_issues_repo_name: JuliaMakers/MakerProjects

max_forks_repo_name: JuliaMakers/JuliaMakerProjects