text:
## Vector/Collections API

function append!(x::RLEVector, y::RLEVector)
    last_x_run = nrun(x)
    length_x = length(x)
    if last(x) == first(y)
        last_x_run = last_x_run - 1
    end
    new_nrun = last_x_run + nrun(y)
    resize!(x.runvalues, new_nrun)
    resize!(x.runends, new_nrun)
    x.runvalues[(last_x_run.+1:end)] = y.runvalues
    x.runends[(last_x_run.+1:end)] = y.runends .+ length_x
    x
end

function vcat(x::RLEVector, y::RLEVector...)
    out = copy(x)
    for yi in y
        append!(out, yi)
    end
    out
end

function pop!(x::RLEVector)
    runcount = nrun(x)
    isempty(x) && throw(ArgumentError("array must be non-empty"))
    item = x.runvalues[end]
    x.runends[end] -= 1
    if x.runends[end] == 0
        deleteat!(x.runvalues, runcount)
        deleteat!(x.runends, runcount)
    end
    return (item)
end

function push!(x::RLEVector{T1,T2}, item) where {T1,T2<:Integer}
    item = convert(T1, item) # Copying how base does it for arrays
    if !isempty(x) && (item == x.runvalues[end])
        x.runends[end] += 1
    else
        push!(x.runvalues, item)
        push!(x.runends, (isempty(x) ? 0 : x.runends[end]) + 1)
    end
    return (x)
end

function popfirst!(x::RLEVector)
    isempty(x) && throw(ArgumentError("array must be non-empty"))
    item = x.runvalues[1]
    x.runends[:] = x.runends .- 1
    if x.runends[1] == 0
        deleteat!(x.runvalues, 1)
        deleteat!(x.runends, 1)
    end
    return (item)
end

function pushfirst!(x::RLEVector{T1,T2}, item) where {T1,T2<:Integer}
    item = convert(T1, item) # Copying how base does it for arrays
    x.runends[:] = x.runends .+ 1
    if item != x.runvalues[1]
        pushfirst!(x.runvalues, item)
        pushfirst!(x.runends, 1)
    end
    return (x)
end

function deleterun!(x::RLEVector, i::Integer)
    x.runends[i:end] = x.runends[i:end] .- widths(x, i)
    if (i > 1 && i < nrun(x) && x.runvalues[i-1] == x.runvalues[i+1])
        deleteat!(x.runvalues, (i-1):i)
        deleteat!(x.runends, (i-1):i)
    else
        deleteat!(x.runvalues, i)
        deleteat!(x.runends, i)
    end
    x
end

function decrement_run!(x::RLEVector, run::Integer)
    if widths(x, run) == 1
        deleterun!(x, run)
    else
        x.runends[run:end] .-= 1
    end
    return (x)
end

function deleteat!(x::RLEVector, i::Integer)
    run = ind2run(x, i)
    decrement_run!(x, run)
end

function insert!(x::RLEVector{T1,T2}, i::Integer, item) where {T1,T2<:Integer}
    if i == 1
        pushfirst!(x, item)
    elseif i == length(x) + 1
        push!(x, item)
    else
        _item = convert(T1, item)
        (run, index_in_run, run_remainder) = ind2runcontext(x, i)
        if x.runvalues[run] == _item
            x.runends[run:end] = x.runends[run:end] .+ 1
        else
            new_nrun = nrun(x) + 2
            resize!(x.runvalues, new_nrun)
            resize!(x.runends, new_nrun)
            x.runvalues[(run+2):end] = x.runvalues[run:(end-2)]
            x.runvalues[run+1] = _item
            x.runends[(run+2):end] = x.runends[run:(end-2)] .+ 1
            x.runends[run+1] = i
            x.runends[run] = i - 1
        end
    end
    x
end

_default_splice = RLEVector(Union{}[], Int64[])
function splice!(x::RLEVector, i::Integer, ins::RLEVector = _default_splice)
    nrun_x = nrun(x)
    (1 <= i <= length(x)) || throw(BoundsError())
    if length(ins) == 0
        run = ind2run(x, i)
        current = x.runvalues[run]
        decrement_run!(x, run)
    else
        (run, index_in_run, run_remainder) = ind2runcontext(x, i)
        current = x.runvalues[run]
        # Splice ins into adjusted x
        widths!(x.runends)
        x.runends[run] = run_remainder
        nrun_out = run + nrun(ins) + ((nrun_x - run) + 1)
        resize!(x.runvalues, nrun_out)
        resize!(x.runends, nrun_out)
        x.runvalues[(run+(nrun_out-nrun_x)):nrun_out] = x.runvalues[run:nrun_x]
        x.runends[(run+(nrun_out-nrun_x)):nrun_out] = x.runends[run:nrun_x]
        x.runvalues[(run+1):(run+nrun(ins))] = ins.runvalues
        x.runends[(run+1):(run+nrun(ins))] = widths(ins)
        x.runends[run] = index_in_run - 1
        cumsum!(x.runends, x.runends)
        ree!(x.runvalues, x.runends)
    end
    return (current)
end

function splice!(x::RLEVector, index::UnitRange, ins::RLEVector = _default_splice)
    nrun_x = nrun(x)
    i_left = first(index)
    i_right = last(index)
    if i_left == i_right
        return (splice!(x, i_left, ins))
    end
    (run_left, run_right, index_in_run_left, run_remainder_right) = ind2runcontext(x, index)
    run_range = run_left:run_right
    if i_left > i_right # Insert without removing
        current = similar(x, 0)
    else
        current = RLEVector(x.runvalues[run_range], x.runends[run_range] .- (i_left - 1))
        current.runends[end] = current.runends[end] - run_remainder_right
    end
    # Splice ins into adjusted x
    # FIXME: factor out code shared with setindex!(rle, rle, unit_range)
    widths!(x.runends)
    x.runends[run_right] = run_remainder_right
    nrun_out = run_left + nrun(ins) + ((nrun_x - run_right) + 1)
    resize!(x.runvalues, nrun_out)
    resize!(x.runends, nrun_out)
    x.runvalues[(run_right+(nrun_out-nrun_x)):nrun_out] = x.runvalues[run_right:nrun_x]
    x.runends[(run_right+(nrun_out-nrun_x)):nrun_out] = x.runends[run_right:nrun_x]
    x.runvalues[(run_left+1):(run_left+nrun(ins))] = ins.runvalues
    x.runends[(run_left+1):(run_left+nrun(ins))] = widths(ins)
    x.runends[run_left] = index_in_run_left - 1
    cumsum!(x.runends, x.runends)
    ree!(x.runvalues, x.runends)
    return (current)
end

function splice!(x::RLEVector, i::Union{Integer,UnitRange}, ins::AbstractArray)
    splice!(x, i, RLEVector(ins))
end

# Appended space initialized with zero unlike base array
function resize!(x::RLEVector, nl::Integer) # Based on base version for array
    l = length(x)
    if nl > l
        push!(x.runends, nl)
        push!(x.runvalues, 0)
    else
        nl < 0 && throw(ArgumentError("new length must be ≥ 0"))
        (run, index_in_run, run_remainder) = ind2runcontext(x, nl)
        resize!(x.runends, run)
        resize!(x.runvalues, run)
        x.runends[end] = x.runends[end] - run_remainder
    end
    return (x)
end

function empty!(x::RLEVector)
    empty!(x.runvalues)
    empty!(x.runends)
    return (x)
end

### Set API
# Mostly handled by convert and promote rules in RLEVector-type.jl

function intersect(x::RLEVector, sets...)
    ok = trues(length(x.runvalues))
    for (i, v) in enumerate(x.runvalues)
        for s in sets
            if !in(v, s)
                ok[i] = false
                break
            end
        end
    end
    RLEVector(x.runvalues[ok], cumsum(widths(x)[ok]))
end


text:
module TestTargetedLearning

using FactCheck
# FactCheck.setstyle(:compact)

is_a(T) = x -> isa(x, T)

using TargetedLearning

include("LReg.jl")
include("automaticics.jl")
include("Qmodels.jl")
include("Parameters.jl")
include("TMLEs.jl")
include("CTMLEs.jl")

FactCheck.exitstatus()
end

ext: jl


text:
struct PointerString <: AbstractString
    ptr::Ptr{UInt8}
    len::Int
end

function Base.hash(s::PointerString, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

import Base: ==
function ==(x::String, y::PointerString)
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), pointer(x), y.ptr, y.len) == 0
end
==(y::PointerString, x::String) = x == y

Base.codeunit(s::PointerString) = UInt8
Base.ncodeunits(s::PointerString) = s.len
@inline function Base.codeunit(s::PointerString, i::Integer)
    @boundscheck checkbounds(s, i)
    unsafe_load(s.ptr + i - 1)
end
Base.String(x::PointerString) = unsafe_string(x.ptr, x.len)

function reverseescapechar(b)
    b == UInt8('"')  && return UInt8('"')
    b == UInt8('\\') && return UInt8('\\')
    b == UInt8('/')  && return UInt8('/')
    b == UInt8('b')  && return UInt8('\b')
    b == UInt8('f')  && return UInt8('\f')
    b == UInt8('n')  && return UInt8('\n')
    b == UInt8('r')  && return UInt8('\r')
    b == UInt8('t')  && return UInt8('\t')
    return 0x00
end

utf16_is_surrogate(c::UInt16) = (c & 0xf800) == 0xd800
utf16_get_supplementary(lead::UInt16, trail::UInt16) = Char(UInt32(lead-0xd7f7)<<10 + trail)

charvalue(b) = (UInt8('0') <= b <= UInt8('9')) ? b - UInt8('0') :
               (UInt8('a') <= b <= UInt8('f')) ? b - (UInt8('a') - 0x0a) :
               (UInt8('A') <= b <= UInt8('F')) ? b - (UInt8('A') - 0x0a) :
               throw(ArgumentError("JSON invalid unicode hex value"))

@noinline invalid_escape(str) = throw(ArgumentError("encountered invalid escape character in json string: \"$(String(str))\""))
@noinline unescaped_control(b) = throw(ArgumentError("encountered unescaped control character in json: '$(escape_string(Base.string(Char(b))))'"))

function unescape(s)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == UInt8('\\')
                i += 1
                i > n && invalid_escape(s)
                b = codeunit(s, i)
                if b == UInt8('u')
                    c = 0x0000
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    i += 1
                    i > n && invalid_escape(s)
                    b = codeunit(s, i)
                    c = (c << 4) + charvalue(b)
                    if utf16_is_surrogate(c)
                        i += 3
                        i > n && invalid_escape(s)
                        c2 = 0x0000
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        i += 1
                        i > n && invalid_escape(s)
                        b = codeunit(s, i)
                        c2 = (c2 << 4) + charvalue(b)
                        ch = utf16_get_supplementary(c, c2)
                    else
                        ch = Char(c)
                    end
                    st = codeunits(Base.string(ch))
                    for j = 1:length(st)-1
                        @inbounds buf[len] = st[j]
                        len += 1
                    end
                    b = st[end]
                else
                    b = reverseescapechar(b)
                    b == 0x00 && invalid_escape(s)
                end
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end


text:
@testset "▶ types/ax                    " begin
    # CONSTRUCTORS
    a = GPlot.Axis("x")
    @test a.prefix == "x"
    @test isnothing(a.textstyle.hei)
    @test isnothing(a.textstyle.font)
    @test isnothing(a.textstyle.color)
    @test isnothing(a.title)
    @test a.off == false
    @test isnothing(a.base)
    @test isnothing(a.lwidth)
    @test a.log == false
    @test isnothing(a.min)
    @test isnothing(a.max)

    f = Figure()
    a = GPlot.Axes2D{GPlot.GLE}(parent=f.id)
    @test a.x2axis.prefix == "x2"
    @test isempty(a.drawings)
    @test isnothing(a.title)
    @test isnothing(a.size)
    @test isnothing(a.legend)
    @test a.math == false
    @test isnothing(a.origin)

    # show
    io = IOBuffer()
    Base.show(io, MIME("text/plain"), a)
    s = String(take!(io))
    isin(s, "GPlot.Axes2D{GLE}")
    isin(s, rpad("Title:", 15) * "none")
    isin(s, rpad("N. drawings:", 15) * "0")
    isin(s, rpad("Math mode:", 15) * "false")
    isin(s, rpad("Layout origin:", 15) * "auto")
end

@testset "▶ /ax                         " begin
    f = Figure()
    # axis limits
    # -- without kwargs
    xlim(0, 3); ylim(0, 3)
    x2lim(0, 3); y2lim(0, 3)
    @test f.axes[1].xaxis.min  == 0.0
    @test f.axes[1].xaxis.max  == 3.0
    @test f.axes[1].x2axis.min == 0.0
    @test f.axes[1].x2axis.max == 3.0
    @test f.axes[1].yaxis.min  == 0.0
    @test f.axes[1].yaxis.max  == 3.0
    @test f.axes[1].y2axis.min == 0.0
    @test f.axes[1].y2axis.max == 3.0
    # -- with kwargs
    xlim(min=1, max=2); ylim(min=1, max=2)
    x2lim(min=1, max=2); y2lim(min=1, max=2)
    @test f.axes[1].xaxis.min  == 1.0
    @test f.axes[1].xaxis.max  == 2.0
    @test f.axes[1].x2axis.min == 1.0
    @test f.axes[1].x2axis.max == 2.0
    @test f.axes[1].yaxis.min  == 1.0
    @test f.axes[1].yaxis.max  == 2.0
    @test f.axes[1].y2axis.min == 1.0
    @test f.axes[1].y2axis.max == 2.0
    # scale
    erase!(f)
    xscale("log");
    @test f.axes[1].xaxis.log == true
    @test f.axes[1].x2axis.log == false
    x2scale("log"); yscale("log"); y2scale("log")
    @test f.axes[1].x2axis.log == true
    erase!(f)
    yscale("log")
    @test f.axes[1].xaxis.log == false
    @test f.axes[1].yaxis.log == true
    erase!(f); x2scale("log"); erase!(f); y2scale("log")
    @test f.axes[1].y2axis.log == true

    # diverse extra ones prompted by codecov
    f = Figure()
    xlim(min=0, max=2)
    @test f.axes[1].xaxis.min == 0.0
    @test f.axes[1].xaxis.max == 2.0
    xlim(0, 2)
    @test f.axes[1].xaxis.min == 0.0
    @test f.axes[1].xaxis.max == 2.0

    erase!(gcf())
    xscale("log")
    @test f.axes[1].xaxis.log
    @test_throws G.OptionValueError yscale("las")
    xscale("lin")
    @test f.axes[1].xaxis.log == false
    xscale("log")
    @test f.axes[1].xaxis.log

    # xaxis, yaxis, ... with shorthands
    erase!(gcf())
    xaxis("log")
    @test f.axes[1].xaxis.log
    xaxis("lin")
    @test !f.axes[1].xaxis.log
    y2axis("off")
    @test f.axes[1].y2axis.off
    y2axis("on")
    @test !f.axes[1].y2axis.off
    cla()
    yaxis(title="blah", log=true)
    @test f.axes[1].yaxis.log
    @test f.axes[1].yaxis.title.text == "blah"
    @test_throws G.OptionValueError x2axis("somethingweird")
    @test_throws G.OptionValueError grid(axis=["akdjf"])

    clf()
    grid(lw=0.5)
    @test f.axes[1].xaxis.ticks.linestyle.lwidth == 0.5
    @test f.axes[1].yaxis.ticks.linestyle.lwidth == 0.5

    # math
    clf()
    math()
    @test f.axes[1].math

    # more misc
    @test_throws ArgumentError xlim(5, 2)

    # test of clo
    cla()
    text("blah", (0.5,0.5))
    clo()
    @test isempty(gca().objects)

    # grid
    cla()
    grid()
    @test gca().xaxis.ticks.grid && gca().yaxis.ticks.grid
    grid("off")
    @test !gca().xaxis.ticks.grid && !gca().yaxis.ticks.grid
    grid("on")
    @test gca().xaxis.ticks.grid && gca().yaxis.ticks.grid

    # shorthands for axes
    cla()
    axis("math")
    @test gca().math
end

@testset "▶ set_prop/ax                 " begin
    f = Figure()
    yaxis(title="blah", base=0.3, min=0, max=4, log=true)
    @test f.axes[1].yaxis.title.text  == "blah"
    @test f.axes[1].yaxis.base == 0.3
    @test f.axes[1].yaxis.min == 0.0
    @test f.axes[1].yaxis.max == 4.0
    @test f.axes[1].yaxis.log
end

@testset "▶ apply_gle/ax                " begin
    # AXES2D (see also apply_title, apply_drawings, apply_legend)
    g = G.GLE()
    f = G.Figure();
    G.add_axes2d!()
    G.apply_axes!(g, f.axes[1], f.id, 1)
    s = String(take!(g))
    isin(s, "begin graph")
    isin(s, "scale auto")
    isin(s, "end graph")
    f.axes[1].math = true
    f.axes[1].size = (10.,8.)
    G.apply_axes!(g, f.axes[1], f.id, 1)
    s = String(take!(g))
    isin(s, "math")
    isin(s, "size 10.0 8.0")

    # AXIS (see also apply_ticks, apply_textstyle)
    erase!(f)
    G.add_axes2d!()
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    # XXX by default subticks are off (may change)
    isin(s, "xsubticks off")
    isin(s, "ysubticks off")
    isin(s, "x2subticks off")
    isin(s, "y2subticks off")
    f.axes[1].xaxis.off = true
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    isin(s, "xaxis off")
    f.axes[1].xaxis.off = false
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    notisin(s, "xaxis off")
    f.axes[1].xaxis.base = 0.1
    f.axes[1].yaxis.lwidth = 2.0
    f.axes[1].x2axis.ticks.grid = true
    f.axes[1].y2axis.log = true
    f.axes[1].xaxis.min = 0.0
    f.axes[1].yaxis.max = 2.0
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    isin(s, "xaxis base 0.1 min 0.0")
    isin(s, "x2axis grid")
    isin(s, "yaxis lwidth 2.0 max 2.0")
    isin(s, "y2axis log")
end


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : false
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 512
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 512
###
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    # CUDA.device!(6) # select GPU
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, LinearAlgebra, MAT

macro innH3(ix,iy)       esc(:( H[$ix+1,$iy+1] * H[$ix+1,$iy+1] * H[$ix+1,$iy+1] )) end
macro av_xi_H3(ix,iy)    esc(:( 0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]) * 0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]) * 0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]) )) end
macro av_yi_H3(ix,iy)    esc(:( 0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]) * 0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]) * 0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]) )) end
macro av_xi_Re(ix,iy)    esc(:( π + sqrt(π*π + max_lxy2 / @av_xi_H3($ix,$iy) * _dt) )) end
macro av_yi_Re(ix,iy)    esc(:( π + sqrt(π*π + max_lxy2 / @av_yi_H3($ix,$iy) * _dt) )) end
macro Re(ix,iy)          esc(:( π + sqrt(π*π + max_lxy2 / @innH3($ix,$iy)    * _dt) )) end
macro av_xi_θr_dτ(ix,iy) esc(:( max_lxy / Vpdτ / @av_xi_Re($ix,$iy) * Resc )) end
macro av_yi_θr_dτ(ix,iy) esc(:( max_lxy / Vpdτ / @av_yi_Re($ix,$iy) * Resc )) end
macro dτ_ρ(ix,iy)        esc(:( Vpdτ * max_lxy / @innH3($ix,$iy) / @Re($ix,$iy) * Resc )) end

@parallel_indices (ix,iy) function compute_flux!(qHx, qHy, H, Vpdτ, Resc, _dt, max_lxy, max_lxy2, _dx, _dy)
    if (ix<=size(qHx,1) && iy<=size(qHx,2))  qHx[ix,iy]  = (qHx[ix,iy] * @av_xi_θr_dτ(ix,iy) - @av_xi_H3(ix,iy) * _dx * (H[ix+1,iy+1] - H[ix,iy+1]) ) / (1.0 + @av_xi_θr_dτ(ix,iy))  end
    if (ix<=size(qHy,1) && iy<=size(qHy,2))  qHy[ix,iy]  = (qHy[ix,iy] * @av_yi_θr_dτ(ix,iy) - @av_yi_H3(ix,iy) * _dy * (H[ix+1,iy+1] - H[ix+1,iy]) ) / (1.0 + @av_yi_θr_dτ(ix,iy))  end
    return
end

@parallel_indices (ix,iy) function compute_update!(H, Hold, qHx, qHy, Vpdτ, Resc, _dt, max_lxy, max_lxy2, _dx, _dy, size_innH_1, size_innH_2)
    if (ix<=size_innH_1 && iy<=size_innH_2)  H[ix+1,iy+1] = (H[ix+1,iy+1] + @dτ_ρ(ix,iy) * (_dt * Hold[ix+1,iy+1] - (_dx * (qHx[ix+1,iy] - qHx[ix,iy]) + _dy * (qHy[ix,iy+1] - qHy[ix,iy])) )) / (1.0 + _dt * @dτ_ρ(ix,iy))  end
    return
end

@parallel_indices (ix,iy) function compute_flux_res!(qHx2, qHy2, H, _dx, _dy)
    if (ix<=size(qHx2,1) && iy<=size(qHx2,2))  qHx2[ix,iy] = -@av_xi_H3(ix,iy) * _dx * (H[ix+1,iy+1] - H[ix,iy+1])  end
    if (ix<=size(qHy2,1) && iy<=size(qHy2,2))  qHy2[ix,iy] = -@av_yi_H3(ix,iy) * _dy * (H[ix+1,iy+1] - H[ix+1,iy])  end
    return
end

@parallel_indices (ix,iy) function check_res!(ResH, H, Hold, qHx2, qHy2, _dt, _dx, _dy)
    if (ix<=size(ResH,1) && iy<=size(ResH,2))  ResH[ix,iy] = -_dt * (H[ix+1,iy+1] - Hold[ix+1,iy+1]) - (_dx * (qHx2[ix+1,iy] - qHx2[ix,iy]) + _dy * (qHy2[ix,iy+1] - qHy2[ix,iy]))  end
    return
end

@parallel_indices (ix,iy) function assign!(Hold, H)
    if (ix<=size(H,1) && iy<=size(H,2))  Hold[ix,iy] = H[ix,iy]  end
    return
end

@views function diffusion_2D_()
    # Physics
    lx, ly  = 10.0, 10.0    # domain size
    ttot    = 0.4           # total simulation time
    dt      = 0.2           # physical time step
    # Numerics
    # nx, ny  = 2*256, 2*256  # numerical grid resolution
    tol     = 1e-8          # tolerance
    itMax   = 1e2#1e5           # max number of iterations
    nout    = 2000            # tol check
    CFL     = 1 / sqrt(2)   # CFL number
    Resc    = 1 / 1.2       # iteration parameter scaling
    # Derived numerics
    dx, dy  = lx / nx, ly / ny  # grid size
    Vpdτ    = CFL * min(dx, dy)
    max_lxy = max(lx, ly)
    max_lxy2= max_lxy^2
    xc, yc  = LinRange(dx/2, lx - dx/2, nx), LinRange(dy/2, ly - dy/2, ny)
    _dx, _dy, _dt = 1.0/dx, 1.0/dy, 1.0/dt
    # Array allocation
    qHx     = @zeros(nx-1, ny-2)
    qHy     = @zeros(nx-2, ny-1)
    qHx2    = @zeros(nx-1, ny-2)
    qHy2    = @zeros(nx-2, ny-1)
    ResH    = @zeros(nx-2, ny-2)
    # Initial condition
    H0      = Data.Array(exp.(-(xc .- lx/2).^2 .- ((yc .- ly/2)').^2))
    Hold    = @ones(nx,ny) .* H0
    H       = @ones(nx,ny) .* H0
    size_innH_1, size_innH_2 = size(H,1)-2, size(H,2)-2
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot/dt)); t_tic = 0.0; niter = 0
    # Physical time loop
    while it < nt
        iter = 0; err = 2 * tol
        # Pseudo-transient iteration
        while err > tol && iter < itMax
            if (it==1 && iter==0) t_tic = Base.time(); niter = 0 end
            @parallel compute_flux!(qHx, qHy, H, Vpdτ, Resc, _dt, max_lxy, max_lxy2, _dx, _dy)
            @parallel compute_update!(H, Hold, qHx, qHy, Vpdτ, Resc, _dt, max_lxy, max_lxy2, _dx, _dy, size_innH_1, size_innH_2)
            iter += 1;  niter += 1
            if iter % nout == 0
                @parallel compute_flux_res!(qHx2, qHy2, H, _dx, _dy)
                @parallel check_res!(ResH, H, Hold, qHx2, qHy2, _dt, _dx, _dy)
                err = norm(ResH) / sqrt(length(ResH))
                if isnan(err) error("NaN") end
            end
        end
        ittot += iter; it += 1; t += dt
        @parallel assign!(Hold, H)
    end
    t_toc = Base.time() - t_tic
    A_eff = (2*3+2)/1e9*nx*ny*sizeof(Data.Number) # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                           # Execution time per iteration [s]
    T_eff = A_eff/t_it                            # Effective memory throughput [GB/s]
    @printf("PERF: Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=3), niter)
    @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, nx, ittot)
    # Visualise
    if do_viz display(heatmap(xc, yc, Array(H'), aspect_ratio=1, framestyle=:box, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), xlabel="lx", ylabel="ly", c=:viridis, clims=(0, 1), title="linear diffusion (nt=$it, iters=$ittot)")) end
    if do_save
        !ispath("../output") && mkdir("../output")
        open("../output/out_diff_2D_nonlin_perf.txt","a") do io
            println(io, "$(nx) $(ny) $(ittot) $(t_toc) $(A_eff) $(t_it) $(T_eff)")
        end
    end
    if do_save_viz
        !ispath("../out_visu") && mkdir("../out_visu")
        matwrite("../out_visu/diff_2D_nonlin_perf.mat", Dict("H_2D"=> Array(H), "xc_2D"=> Array(xc), "yc_2D"=> Array(yc)); compress = true)
    end
    return xc, yc, H
end

if use_return
    xc, yc, H = diffusion_2D_();
else
    diffusion_2D = begin diffusion_2D_(); return; end
end


text:
using CategoricalArrays
using DataFrames
using Codex
using Test

@testset "Utils" begin
    A = ["1", 2, missing]
    @test dirparent("/a/b/c") == "/a/b"
    @test dirparent("/a/b/c/") == "/a/b"
    @test dirparent("/a/b/c", 2) == "/a"
    @test has_duplicates([1]) == false
    @test has_duplicates([1,1]) == true

    @test endswith(Codex.PROJECT_ROOT, "Codex.jl")

    fns = [
        df -> transform(df, :A => :A2),
        df -> select(df, :A2)
    ]
    @test apply(fns, DataFrame(:A => 1)) == DataFrame(:A2 => 1)
    @test map(apply([x -> x + 1]), [1]) == [2]

    df = DataFrame(from = [1], to = ["one"])
    @test Codex.map_by_df([1, 2], df, :from, :to) == ["one", 2]

    @test Codex.rescale(5, 1, 5, 0, 6) == 6

    df = DataFrame(group = [4, 9, 9])
    @test Codex.nrow_per_group(df, :group; col1="group", col2="n") == DataFrame(group = [4, 9], n = [1, 2])

    cmd = `echo lorem`
    f(out, err) = pipeline(cmd, stdout=out, stderr=err)
    expected = Output(0, "lorem\n", "")
    @test Codex.output(f) == expected
    @test Codex.output(cmd) == expected

    df = DataFrame(; group=[:A, :A, :B, :B], data=[1, 2, 3, 4])
    @test Codex.split_data(df, :group, :data, [:A, :B]) == ([1, 2], [3, 4])
end


text:
using LinearAlgebra, Distributions
using Test
using DistributionVectors
using StatsBase # coef

@testset "cormatrix_for_acf" begin
    acf0 = [1,0.4,0.1]
    nsum = 5
    Sigma = cormatrix_for_acf(nsum, acf0);
    @test diag(Sigma) == fill(1.0, nsum)
    @test diag(Sigma,1) == diag(Sigma,-1) == fill(0.4, nsum - 1)
    @test diag(Sigma,2) == diag(Sigma,-2) == fill(0.1, nsum - 2)
    @test diag(Sigma,3) == diag(Sigma,-3) == fill(0.0, nsum - 3)
    #dmn = MvNormal(ones(nsum), Symmetric(Sigma));
end;

@testset "AutoCorrelationFunction" begin
    acfvec = [1,0.4,0.1]
    acf = @inferred AutoCorrelationFunction(acfvec)
    @test @inferred coef(acf) == acfvec
    @test @inferred coef(acf,1) == acfvec[2] # vector starts with lag 0
end;


text:
#=- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 Optical Pump Probe (OPP) Spectroscopy Data Structures
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -=#

export TimeDelay, DiffReflect

struct TimeDelay{T} <: AbstractSpectra{T}
    td::Vector{T}

    TimeDelay(td::Vector{T}) where T<:Real = new{T}(td)

    function TimeDelay(pth::String)
        raw = readdlm(pth); nrow, ncol = size(raw); T = eltype(raw)

        if !isone(ncol)
            error("Strange TimeDelay data file of ncol = $ncol")
        end

        td = Vector{T}(undef, nrow)

        @inbounds @simd for i in eachindex(td)
            td[i] = raw[i]
        end

        return TimeDelay(td)
    end
end

struct DiffReflect{N,T,S} <: AbstractSpectra{T}
    ΔR::Matrix{T}; spec::NTuple{N, S}

    DiffReflect(ΔR::Matrix{T}, spec::NTuple{N,S}) where {N,T<:Real,S} = new{N,T,S}(ΔR, spec)

    function DiffReflect(pth::String, spec::NTuple{N,S}; ifnorm::Bool=true) where {N,S}
        ΔR = readdlm(pth)

        if ifnorm
            nrow, ncol = size(ΔR)
            @inbounds begin
                for j in 1:ncol
                    abs_minΔR = abs(minimum(view(ΔR, :, j)))
                    @simd for i in 1:nrow
                        ΔR[i, j] /= abs_minΔR
                    end
                end
            end
        end

        return DiffReflect(ΔR, spec)
    end
end


text:
# #absolutely not important! Temporary file

import Pkg
Pkg.activate(normpath(@__DIR__))

import ColorTypes: RGB
import ProjectMoana: HdrImage
import ProjectMoana: _read_line
import ProjectMoana: _read_float
import ProjectMoana: read_pfm_image
import ProjectMoana: Vec
function main(ARGS)
    a=Vec(1,2,3.4)
    b=Vec(1,2,3.4)
    print(a+b)
end

main(ARGS)


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>1.0, :alg=>GLGM06(δ=1e-4)),
    "X0"   => [(1, Singleton([3.1, 4.0, 0.0, 0.0, 0.0]))],
    "info" => Dict(
        "name" => "linear_switching",
        "dim" => 1,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using BlackBoxOptim: set_candidate!, candidate

# a 2d optimization problem with an optimum at (3.14, 7.2)
fixed_optimum_prob(x) = (x[1] - 3.14)^2 + (x[2] - 7.2)^4
const FixedOptimum = [3.14, 7.2]
const FitnessOptimum = fixed_optimum_prob(FixedOptimum)

@testset "set_candidate! for DE population optimizers" begin
    for m in [:de_rand_1_bin, :de_rand_2_bin, :de_rand_1_bin_radiuslimited, :de_rand_2_bin_radiuslimited,
                :adaptive_de_rand_1_bin, :adaptive_de_rand_1_bin_radiuslimited]
        PopSize = 10
        b = bbsetup(fixed_optimum_prob; Method = m, MaxFuncEvals = 10*PopSize, # Give it a chance to be sampled so best fitness is in archive
            NumDimensions = 2, SearchRange = (-10.0, 10.0), PopulationSize = PopSize)
        set_candidate!(b.optimizer, FixedOptimum)
        inds = population(b.optimizer).individuals
        # There is at least one position in the population that is set to the FixedOptimum:
        @test any(i -> inds[:, i] == FixedOptimum, 1:PopSize)

        res = bboptimize(b; TraceMode = :silent)
        @test isapprox(best_fitness(res), FitnessOptimum)
    end
end

@testset "set_candidate! for stepping optimizers" begin
    for m in [:generating_set_search, :resampling_memetic_search, :resampling_inheritance_memetic_search,
                :simultaneous_perturbation_stochastic_approximation]
        b = bbsetup(fixed_optimum_prob; Method = m, MaxFuncEvals = 2,
            NumDimensions = 2, SearchRange = (-10.0, 10.0), PopulationSize = 5)
        set_candidate!(b.optimizer, FixedOptimum)
        @test candidate(b.optimizer) == FixedOptimum

        # Note that we exclude simultaneous_perturbation_stochastic_approximation since it updates even
        # if new candidates give worse fitness => not guaranteed to find the optimum. We should debug this at some point.
        if m != :simultaneous_perturbation_stochastic_approximation
            res = bboptimize(b; TraceMode = :silent)
            @test isapprox(best_fitness(res), FitnessOptimum)
        end
    end
end

@testset "set_candidate! for NES optimizers" begin
    for m in [:dxnes, :xnes, :separable_nes]
        b = bbsetup(fixed_optimum_prob; Method = m, MaxFuncEvals = 100,
            NumDimensions = 2, SearchRange = (-10.0, 10.0), PopulationSize = 5)
        set_candidate!(b.optimizer, FixedOptimum)
        @test candidate(b.optimizer) == FixedOptimum

        res = bboptimize(b; TraceMode = :silent)
        @test isapprox(best_fitness(res), FitnessOptimum, atol = 1e-1)
    end
end

@testset "set_candidate! for BORG" begin
    # Best aggregated fitness should be for [0.5, 0.5, 0.5]:
    fitness_2obj(x) = (sum(abs2, x), sum(abs2, x .- 1.0))
    x0 = 0.5 * ones(3)

    PopSize = 10
    b = bbsetup(fitness_2obj; Method=:borg_moea,
            FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true),
            MaxFuncEvals = 10*PopSize, PopulationSize = PopSize,
            SearchRange=(-10.0, 10.0), NumDimensions=3, ϵ=0.05);

    set_candidate!(b.optimizer, x0)

    inds = b.optimizer.population.individuals
    # There is at least one position in the population that is set to x0
    @test any(i -> inds[:, i] == x0, 1:PopSize)

    # Now ensure we actually get back (sum(abs2, x0), sum(abs2, x0 .- 1.0)) as the best fitness (since it is the best aggregated one).
    res = bboptimize(b; TraceMode = :silent)
    @test best_fitness(res) == fitness_2obj(x0)
end


text:
# Parameter helper functions
const gas_molpL = 22.4

watts_to_light_mol(watts) = watts * 4.57e-6
light_mol_to_watts(light_mol) = light_mol / 4.57e-6
water_content_to_mols_per_litre(wc) = wc * 55.5 # L/L of water to mol/L
fraction_per_litre_gas_to_mols(frac) = frac / 22.4

" Variables for carbon assimilation "
@description @limits @units @udefault_kw @flattenable mutable struct CarbonVars{MoMS,MoL,KPA}
    J_L_F::MoMS              | false | watts_to_light_mol(800.0) | mol*m^-2*s^-1 | watts_to_light_mol.([0.0, 2000.0]) | "flux of useful photons"
    X_C::MoL                 | false | 400.0*1e-6                | mol*L^-1      | [200.0, 700.0] .* 1e-6             | "carbon dioxide concentration in air"
    X_O::MoL                 | false | 0.21 * gas_molpL          | mol*L^-1      | [0.1, 0.3] .* gas_molpL            | "oxygen concentration in air"
    soilwaterpotential::KPA  | false | -100.0                    | kPa           | [0.0, -10000.0]                    | "soil water potential"
end

" Variables for nitgroen assimilation "
@description @limits @units @udefault_kw @flattenable mutable struct NitrogenVars{F,KPA,MoL}
    # TODO work out the naming conventions here
    soilwaterpotential::KPA  | false | -100.0 | kPa      | [0.0, -10000.0] | "soil water potential"
    soilwaterconent::F       | false | -100.0 | _        | [0.0, -10000.0] | _
    X_NH::MoL                | false | 0.005  | mol*L^-1 | [0.0, 0.1]      | "concentration of ammonia"
    X_NO::MoL                | false | 0.01   | mol*L^-1 | [0.0, 0.1]      | "concentration of nitrate see e.g. [_@crawford1998molecular]"
    X_H::MoL                 | false | 10.0   | mol*L^-1 | [0.0, 20.0]     | _
end

" Assimilation "
abstract type AbstractAssim end

" Parent of all Carbon assimilation types"
abstract type AbstractCAssim <: AbstractAssim end

" Parent of all Nitrogen assimilation types"
abstract type AbstractNAssim <: AbstractAssim end

" Parent of Ammonia/Nitrate assimilation types"
abstract type AbstractNH4_NO3Assim <: AbstractNAssim end

@columns struct ConstantCAssim{μMoMS} <: AbstractCAssim
    # Field       | Def | Unit              | Prior            | Limits      | Log | Description
    uptake::μMoMS | 0.1 | μmol*mol^-1*s^-1  | Gamma(2.0, 2.0)  | [0.0, 10.0] | _   | "constant rate of C uptake"
end

@columns struct ConstantNAssim{μMoMS} <: AbstractNAssim
    uptake::μMoMS | 0.1  | μmol*mol^-1*s^-1 | Gamma(2.0, 2.0)  | [0.0, 0.5]  | _   | "constant rate of N uptake"
end

@mix @columns struct SLA{MG}
    SLA::MG       | 24.0 | m^2*kg^-1        | Gamma(10.0, 1.0) | [5.0, 30.0] | _   | "Specific leaf Area. Ferns: 17.4, Forbs: 26.2, Graminoids: 24.0, Shrubs: 9.10, Trees: 8.30"
end

" Uses FvCB photosynthesis model from Photosynthesis.jl "
@SLA struct FvCBPhotosynthesis{P,V} <: AbstractCAssim
    vars::V        | Photosynthesis.PhotoVars           | _ | _ | _ | _ | _
    photoparams::P | Photosynthesis.FvCBEnergyBalance() | _ | _ | _ | _ | _
end

@mix @SLA @columns struct KooijmanPhoto{μMoMoS,MoL,μMoMS,MoMS,V}
    #Field              | Default           | Unit             | Prior           | Limits           | Log  | Description
    vars::V             | CarbonVars()      | _                | _               | _                | _    | _
    k_C_binding::μMoMoS | 10000.0           | μmol*mol^-1*s^-1 | Gamma(2.0, 2.0) | [1e-5, 2000.0]   | true | "Scaling rate for carbon dioxide"
    k_O_binding::μMoMoS | 10000.0           | μmol*mol^-1*s^-1 | Gamma(2.0, 2.0) | [1e-5, 2000.0]   | true | "Scaling rate for oxygen"
    K_C::MoL            | 50*1e-6/gas_molpL | mol*L^-1         | Gamma(2.0, 2.0) | [1e-7, 100.0]    | true | "Half-saturation concentration of carbon dioxide"
    K_O::MoL            | 0.0021/gas_molpL  | mol*L^-1         | Gamma(2.0, 2.0) | [1e-7, 100.0]    | true | "Half-saturation concentration of oxygen"
    J_L_K::MoMS         | 2000.0            | mol*m^-2*s^-1    | Gamma(2.0, 2.0) | [1e-3, 100000.0] | true | "Half-saturation flux of useful photons"
    j_L_Amax::μMoMS     | 100.01            | μmol*m^-2*s^-1   | Gamma(2.0, 2.0) | [1e-4, 10000.0]  | true | "Max spec uptake of useful photons"
    j_C_Amax::μMoMS     | 20.0              | μmol*m^-2*s^-1   | Gamma(2.0, 2.0) | [5.0,  100.0]    | true | "Max spec uptake of carbon dioxide"
    j_O_Amax::μMoMS     | 0.1               | μmol*m^-2*s^-1   | Gamma(2.0, 2.0) | [0.01,  50.0]    | true | "Max spec uptake of oxygen"
end

abstract type AbstractKooijmanPhoto <: AbstractCAssim end

" Parameters for simple photosynthesis module. With specific leaf area to convert area to mass "
@KooijmanPhoto struct KooijmanSLAPhotosynthesis{} <: AbstractKooijmanPhoto end

@KooijmanPhoto struct KooijmanWaterPotentialPhotosynthesis{PL} <: AbstractKooijmanPhoto
    potential_modifier::PL | Photosynthesis.ZhouPotentialDependence() | _ | _ | _ | _ | "modify photosynthesis with a water potential model"
end

# @columns struct WaterPotentialCutoff{KPA}
    # cutoff::KPA    | -250.0  | kPa  | Gamma(2.0, 2.0) | [0.0, -500.0] | _ | "Max spec uptake of oxygen"
# end

" Parameters for Ammonia/Nitrate assimilation "
@columns struct KooijmanNH4_NO3Assim{μMoMS,F,MoMo,MoL,V} <: AbstractNH4_NO3Assim
    vars::V          | NitrogenVars() | _ | _ | _ | _ | _
    j_NH_Amax::μMoMS | 50.0    | μmol*mol^-1*s^-1 | Gamma(50.0, 1.0) | [0.1, 1000.0] | _ | "Max spec uptake of ammonia"
    j_NO_Amax::μMoMS | 50.0    | μmol*mol^-1*s^-1 | Gamma(50.0, 1.0) | [0.1, 1000.0] | _ | "Max spec uptake of nitrate"
    ρNO::F           | 0.7     | _                | Beta(0.7, 1.0)   | [1e-4, 1.0]   | _ | "Weights preference for nitrate relative to ammonia." # 1 or less but why?
    y_E_CH_NH::MoMo  | 1.25    | mol*mol^-1       | Gamma(1.25, 1.0) | [1e-4, 4.0]   | _ | "From roots C-reserve to reserve using ammonia"
    K_NH::MoL        | 0.01    | mol*L^-1         | Gamma(0.01, 1.0) | [1e-4, 10.0]  | _ | "Half-saturation concentration of ammonia"
    K_NO::MoL        | 0.01    | mol*L^-1         | Gamma(0.01, 1.0) | [1e-4, 10.0]  | _ | "Half-saturation concentration of nitrate"
    K_H::MoL         | 10.0    | mol*L^-1         | Gamma(10.0, 1.0) | [5.0, 20.0]   | _ | "Half-saturation concentration of water"
end

" Parameters for lumped Nitrogen assimilation "
@columns struct NAssim{μMoS,MoL,V} <: AbstractNAssim
    vars::V          | NitrogenVars() | _ | _ | _ | _ | _
    j_N_Amax::μMoS   | 50.0    | μmol*mol^-1*s^-1 | Gamma(50.0, 1.0) | [0.1, 1000.0] | _ | "Max spec uptake of ammonia"
    K_N::MoL         | 0.01    | mol*L^-1         | Gamma(0.01, 1.0) | [1e-4, 1.0]   | _ | "Half-saturation concentration of nitrate"
    K_H::MoL         | 10.0    | mol*L^-1         | Gamma(10.0, 1.0) | [1e-2, 100.0] | _ | "Half-saturation concentration of water"
end

"""
    assimilation!(o, u)
Runs assimilation methods, depending on formulation and state.
"""
assimilation!(organs::Tuple, u) = apply(assimilation!, organs, u)
assimilation!(o, u) = begin
    is_germinated(o, u) && assimilation!(has_reserves(o), assimilation_pars(o), o, u)
    nothing
end
assimilation!(::Nothing, x, o, u, env) = nothing

"""
    assimilation!(f::AbstractAassim, o, u)
Runs nitrogen uptake, and combines N with translocated C.
"""
assimilation!(p::HasCNE, f::AbstractCAssim, o, u) = begin
    J = flux(o)
    c_uptake = photosynthesis(f, o, u) * u.V * shape(o)
    n_tra = J[:N,:tra]
    # Merge rejected N from root and photosynthesized C into reserves
    J[:C,:ass], J[:N,:tra], J[:E,:ass] = stoich_merge(c_uptake, n_tra, y_E_CH_NO(o), y_E_EN(o))

    # lc, ln = stoich_merge_losses(c_uptake, n_tra, J[:C,:ass], J[:N,:tra], J[:E,:ass],
                                 # 1, 1, n_N_E(o))
    # J1[:C,:los] += lc
    # J1[:N,:los] += ln
end
assimilation!(::HasCN, f::AbstractCAssim, o, u) = begin
    flux(o)[:C,:ass] = photosynthesis(f, o, u) * u.V * shape(o)
end

"""
    assimilation!(f::AbstractNH4_NO3Assim, o, u)
Runs nitrogen uptake for nitrate and ammonia, and combines N with translocated C.
Unused ammonia is discarded.
"""
assimilation!(p::HasCN, f::AbstractNH4_NO3Assim, o, u) = begin
    J = flux(o)
    J_N_ass, J_NO_ass, J_NH_ass = nitrogen_uptake(f, o, u) .* u.V * shape(o)

    θNH = J_NH_ass/J_N_ass                          # Fraction of ammonia in arriving N-flux
    θNO = 1 - θNH                                   # Fraction of nitrate in arriving N-flux
    y_E_EC = θNH * f.y_E_CH_NH + θNO * shared(o).y_E_EC  # Yield coefficient from C-reserve to reserve

    C_tra = J[:C,:tra]

    # Merge rejected C from shoot and uptaken N into reserves
    J[:C,:tra], J[:N,:ass], J[:E,:ass] = stoich_merge(C_tra, J_N_ass, y_E_EC, 1/n_N_E(o))
    # stoich_merge_losses(c_tra, J_N_ass, J[:C,:tra], J[:N,:ass], J[:E,:ass], 1, 1, n_N_E(o))

    # Unused NH₄ remainder is lost so we recalculate N assimilation for NO₃ only
    J[:N,:ass] = (J_NO_ass - θNO * n_N_E(o) * J[:E,:ass]) * 1/n_N_EN(o)
end

"""
    assimilation!(f::AbstractNAssim, o, u)
Runs nitrogen uptake, and combines N with translocated C.
"""
assimilation!(::HasCNE, f::AbstractNAssim, o, u) = begin
    J = flux(o)
    J_N_assim = nitrogen_uptake(f, o, u) * u.V * shape(o)
    c_tra = J[:C,:tra]

    # Merge rejected C from shoot and uptaken N into reserves
    # treating N as N reserve now carbon has been incorporated.
    J[:C,:tra], J[:N,:ass], J[:E,:ass] = stoich_merge(c_tra, J_N_assim, y_E_CH_NO(o), y_E_EN(o))
    # lc, ln = stoich_merge_losses(c_tra, J_N_assim, J[:C,:tra], J[:N,:ass], J[:E,:ass], 1, 1, n_N_E(o))
    # J1[:C,:los] += lc
    # J1[:N,:los] += ln
end
assimilation!(::HasCN, f::AbstractNAssim, o, u) = begin
    J = flux(o)
    max_N = nitrogen_uptake(f, o, u) * u.V * shape(o)
    unused_C, unused_N, J[:N,:ass] = stoich_merge(su_pars(o), J[:C,:tra], max_N, 1.0, 1.0)
    J[:C,:ass] = unused_C - J[:C,:tra]
end

"""
    photosynthesis(f::ConstantCAssim, o, u)
Returns a constant rate of carbon assimilation.
"""
photosynthesis(f::ConstantCAssim, o, u) = f.uptake

"""
    photosynthesis(f::FvCBPhotosynthesis, o, u)
Returns carbon assimilated in mols per time.
"""
photosynthesis(f::FvCBPhotosynthesis, o, u) = f.vars.aleaf * f.SLA * w_V(o)

"""
    photosynthesis(f::KooijmanSLAPhotosynthesis, o, u)
Returns carbon assimilated in mols per time.
"""
photosynthesis(f::KooijmanSLAPhotosynthesis, o, u) = begin
    v = vars(o); va = f.vars
    mass_area_coef = w_V(o) * f.SLA

    # Photon flux is not temperature dependent, but O and C flux is.
    # See "Stylized facts in microalgal growth" Lorena, Marques, Kooijman, & Sousa.
    j1_l = half_saturation(f.j_L_Amax, f.J_L_K, va.J_L_F) * mass_area_coef / 2
    j1_c = half_saturation(f.j_C_Amax, f.K_C, va.X_C) * mass_area_coef * tempcorrection(o)
    j1_o = half_saturation(f.j_O_Amax, f.K_O, va.X_O) * mass_area_coef * tempcorrection(o)

    # photorespiration.
    bound_o = j1_o/f.k_O_binding # mol/mol
    bound_c = j1_c/f.k_C_binding # mol/mol

    # C flux
    j_c_intake = (j1_c - j1_o)
    j1_co = j1_c + j1_o
    co_l = j1_co/j1_l - j1_co/(j1_l + j1_co)

    j_c_intake / (1 + bound_c + bound_o + co_l)
end

photosynthesis(f::KooijmanWaterPotentialPhotosynthesis, o, u) = begin
    va = f.vars
    mass_area_coef = w_V(o) * f.SLA

    j1_l = half_saturation(f.j_L_Amax, f.J_L_K, va.J_L_F) * mass_area_coef / 2

    # Modify CO2 and O2 intake by water availability to simulate stomatal closure.
    potentialcorrection = potential_dependence(f.potential_modifier, va.soilwaterpotential)
    j1_c = half_saturation(f.j_C_Amax, f.K_C, va.X_C) * mass_area_coef * potentialcorrection * tempcorrection(o)
    j1_o = half_saturation(f.j_O_Amax, f.K_O, va.X_O) * mass_area_coef * potentialcorrection * tempcorrection(o)

    # photorespiration.
    bound_o = j1_o/f.k_O_binding # mol/mol
    bound_c = j1_c/f.k_C_binding # mol/mol

    # C flux
    j_c_intake = (j1_c - j1_o)

    j1_co = j1_c + j1_o
    co_l = j1_co/j1_l - j1_co/(j1_l + j1_co)

    j_c_intake / (1 + bound_c + bound_o + co_l)
end

# photosynthesis(f::KooijmanWaterPotentialCutoffPhotosynthesis, o, u) =
    # if f.vars.soilwaterpotential > assimilation_pars(o).potential_cutoff
        # kooijman_photosynthesis(f, o, u)
    # else
        # zero(kooijman_photosynthesis(f, o, u))
    # end

"""
    nitrogen_uptake(f::ConstantNAssim, o, u)
Returns constant nitrogen assimilation.
"""
nitrogen_uptake(f::ConstantNAssim, o, u) = f.uptake * tempcorrection(o)

"""
    nitrogen_uptake(f::KooijmanNH4_NO3Assim, o, u)
Returns total nitrogen, nitrate and ammonia assimilated in mols per time.
"""
function nitrogen_uptake(f::KooijmanNH4_NO3Assim, o, u)
    va = f.vars

    K1_NH = half_saturation(f.K_NH, f.K_H, va.X_H) # Ammonia saturation. va.X_H was multiplied by ox.scaling. But that makes no sense.
    K1_NO = half_saturation(f.K_NO, f.K_H, va.X_H) # Nitrate saturation
    J1_NH_ass = half_saturation(f.j_NH_Amax, K1_NH, va.X_NH) * tempcorrection(o) # Arriving ammonia mols.mol⁻¹.s⁻¹
    J_NO_ass = half_saturation(f.j_NO_Amax, K1_NO, va.X_NO) * tempcorrection(o) # Arriving nitrate mols.mol⁻¹.s⁻¹

    J_N_ass = J1_NH_ass + f.ρNO * J_NO_ass # Total arriving N flux
    return (J_N_ass, J_NO_ass, J1_NH_ass)
end

"""
    nitrogen_uptake(f::NAssim, o, u)
Returns nitrogen assimilated in mols per time.
"""
function nitrogen_uptake(f::NAssim, o, u)
    va = f.vars
    # Ammonia proportion in soil water
    K1_N = half_saturation(f.K_N, f.K_H, va.X_H)
    # Arriving ammonia in mol mol^-1 s^-1
    half_saturation(f.j_N_Amax, K1_N, va.X_NO) * tempcorrection(o)
end


text:
"""
Methods that read in all draws of forecast output:
```
reverse_transform(m, input_type, cond_type, class, product, vars;
                           forecast_string = "", fourquarter = false, verbose = :low)

reverse_transform(path, class, product, var, rev_transform;
                           pop_growth = [], fourquarter = false, verbose = :low)
```

Methods that take a single draw of a forecast output that is already in memory:
```
reverse_transform(m::AbstractDSGEModel, untransformed::Matrix, start_date,
                           row_labels, class; fourquarter = false,
                           verbose = :low)

reverse_transform(m, untransformed::DataFrame, class; fourquarter = false, verbose = :low)
```

Lowest level:
```
reverse_transform(y, rev_transform; fourquarter = false, y0 = NaN, y0s = [],
    pop_growth = [])
```

The purpose of the `reverse_transform` function is to transform raw
forecast outputs from model units into final output units. We
provide methods at various levels of abstraction.

At the highest level, the first two methods listed above read in raw
forecast output from a file and transform **all draws** corresponding to
the given `vars` or `var`. The model-agnostic function takes the name of the
file containing the forecast output and a reverse transform
function. The model-dependent function determines the filepath and
appropriate reverse transform from the information contained in the
model.

The second two methods take in raw forecast output that is already in
memory corresponding to a **single draw**. In the first method, the
user must specify `row_labels`, which label the rows of the `untransformed`
matrix. Each element of `row_labels` must be a valid observable or
pseudoobservable key that corresponds to the row it labels. These
methods return a `DataFrame` of all transformed series.

The lowest-level method requires the caller to provide additional
information. It is applied to **one or many draws** of a single
series. This method is called by the higher-level functions as well as
`compute_means_bands`.

## Inputs
- `m::AbstractDSGEModel`
- `input_type::Symbol`: see `?forecast_one`
- `cond_type::Symbol`: see `?forecast_one`
- `class::Symbol`: Indicates whether the reverse transform should be
                   applied to observables or pseudo observables
- `product::Symbol`: Indicates the product the reverse transform should be applied to (e.g. `:hist`,
                    `:forecast`)
- `vars::Vector{Symbol}`: Names of specific series for which to apply the reverse transform
- `untransformed`: `nperiods` x `nseries` `Matrix` or `DataFrame` of
                   series in model units.
- `y`:`ndraws` x `nperiods` array of a series in model units.
- `rev_transform::Function`: a function used to to transform a variable

## Keyword arguments
- `fourquarter::Bool`: produce four-quarter results instead of annualized quarter-to-quarter results
- `verbose::Symbol`: level of verbosity: `:high`, `:low`, `:none`
- `y0::AbstractFloat`: value of `y` in the period prior to the
  first period. This is used to compute growth rates for the first
  period.
- `y0s::Vector{AbstractFloat}`: contains the value of `y` in the `n`
  periods prior to the first period that `y` contains.
- `pop_growth::Vector{AbstractFloat}`: a vector of population growth
  data. This is required when the reverse transform associated with
  `untransformed` or `y` adjusts for population.
"""
function reverse_transform(m::AbstractDSGEModel, untransformed::AbstractArray, start_date::Dates.Date,
                           vars::Vector{Symbol}, class::Symbol;
                           fourquarter::Bool = false,
                           verbose::Symbol = :low)
    df = DataFrame()
    nperiods = size(untransformed, 2)
    end_date = iterate_quarters(start_date, nperiods - 1)
    df[!, :date] = quarter_range(start_date, end_date)

    for (ind, var) in enumerate(vars)
        series = untransformed[ind, :]
        df[!, var] = series
    end

    reverse_transform(m, df, class; fourquarter = fourquarter, verbose = verbose)
end

function reverse_transform(m::AbstractDSGEModel, untransformed::DataFrame, class::Symbol;
                           fourquarter::Bool = false, verbose::Symbol = :low)
    # Dates
    @assert (:date in names(untransformed)) "untransformed must have a date column"
    date_list  = convert(Vector{Dates.Date}, untransformed[!,:date])
    start_date = date_list[1]
    end_date   = date_list[end]

    # Get mapping from variable name to Observable or PseudoObservable instance
    if class == :obs
        dict = m.observable_mappings
    elseif class == :pseudo
        dict = m.pseudo_observable_mappings
    else
        error("Class $class does not have reverse transformations")
    end

    # Prepare population series
    population_data, population_forecast = load_population_growth(m, verbose = verbose)
    population_series = get_population_series(:population_growth, population_data,
                                              population_forecast, start_date, end_date)
    population_series = convert(Vector{Float64}, population_series)
    if fourquarter
        population_series = vcat(fill(NaN,3), population_series)
    end

    # Apply reverse transform
    transformed = DataFrame()
    transformed[!,:date] = untransformed[!,:date]

    vars = setdiff(names(untransformed), [:date])
    for var in vars
        y = convert(Vector{Float64}, untransformed[!,var])
        rev_transform = dict[var].rev_transform
        if fourquarter
            rev_transform = get_transform4q(rev_transform)
        end
        transformed[!,var] = reverse_transform(y, rev_transform; fourquarter = fourquarter,
                                             pop_growth = population_series)
    end
    return transformed
end

function reverse_transform(y::AbstractArray, rev_transform::Function;
                           fourquarter::Bool = false,
                           y0 = missing, y0s::Vector = Vector(undef, 0),
                           pop_growth::Vector = Vector(undef, 0)) where T<:AbstractFloat
    if fourquarter
        if rev_transform in [loggrowthtopct_4q_percapita, loggrowthtopct_4q,
                             loggrowthtopct_4q_approx]
            # Sum growth rates y_{t-3}, y_{t-2}, y_{t-1}, and y_t
            y0s = isempty(y0s) ? fill(missing, 3) : y0s
            if rev_transform == loggrowthtopct_4q_percapita
                rev_transform(y, pop_growth, y0s)
            else
                rev_transform(y, y0s)
            end
        elseif rev_transform in [logleveltopct_4q_percapita, logleveltopct_4q,
                                 logleveltopct_4q_approx]
            # Divide log levels y_t by y_{t-4}
            y0s = isempty(y0s) ? fill(missing, 4) : y0s
            if rev_transform == logleveltopct_4q_percapita
                rev_transform(y, pop_growth, y0s)
            else
                rev_transform(y, y0s)
            end
        elseif rev_transform in [quartertoannual, identity]
            rev_transform(y)
        else
            error("Invalid 4-quarter reverse transform: $rev_transform")
        end
    else
        if rev_transform in [logleveltopct_annualized_percapita]
            rev_transform(y, pop_growth, y0)
        elseif rev_transform in [loggrowthtopct_annualized_percapita, loggrowthtopct_percapita]
            rev_transform(y, pop_growth)
        elseif rev_transform in [logleveltopct_annualized, logleveltopct_annualized_approx]
            rev_transform(y, y0)
        elseif rev_transform in [loggrowthtopct_annualized, loggrowthtopct, quartertoannual, identity]
            rev_transform(y)
        else
            error("Invalid reverse transform: $rev_transform")
        end
    end
end

# function reverse_transform(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
#                            class::Symbol, product::Symbol, vars::Vector{Symbol};
#                            forecast_string = "", fourquarter::Bool = false, verbose::Symbol = :low)
#     # Figure out dates
#     start_date, end_date  = if product in [:hist]
#         date_mainsample_start(m), date_mainsample_end(m)
#     elseif product in [:forecast, :bddforecast, :forecast4q, :bddforecast4q]
#         date_forecast_start(m), date_forecast_end(m)
#     elseif product in [:shockdec, :trend, :dettrend]
#         date_shockdec_start(m), date_shockdec_end(m)
#     else
#         error("reverse_transform not supported for product $product")
#     end

#     # Load population series
#     population_data, population_forecast = load_population_growth(m,
#                                                                   verbose = verbose)
#     population_series = get_population_series(:population_growth, population_data,
#                                               population_forecast,
#                                               start_date, end_date)
#     if sum(ismissing.(population_series)) > 0
#         error("Population series being loaded has a missing value.")
#     else
#         population_series = Vector{Float64}(population_series)
#     end

#     if class == :obs
#         dict = m.observable_mappings
#     elseif class == :pseudo
#         dict = m.pseudo_observable_mappings
#     else
#         error("Class $class does not have reverse transformations")
#     end

#     # Determine which file to read in
#     path       = get_forecast_filename(m, input_type, cond_type, Symbol(product, class);
#                                  forecast_string = forecast_string)

#     # Apply reverse transform to each desired variable
#     results    = Dict{Symbol, Array{Float64,2}}()
#     for var in vars
#         series = dict[var]
#         transformed = reverse_transform(path, class, product,
#                                         var, series.rev_transform,
#                                         pop_growth = population_series,
#                                         fourquarter = fourquarter)
#         results[var] = transformed
#     end

#     return results
# end

# function reverse_transform(path::String, class::Symbol, product::Symbol,
#                            var::Symbol, rev_transform::Function;
#                            pop_growth::Vector{Float64} = Vector{Float64}(),
#                            fourquarter::Bool = false)

#     # Read in the draws for this variable
#     draws = JLD2.jldopen(path, "r") do jld
#         read_forecast_output(jld, class, product, var)
#     end

#     # Transform and return
#     transformed = reverse_transform(draws, rev_transform, pop_growth = pop_growth,
#                                     fourquarter = fourquarter)

#     return transformed
# end


text:
# ============================================================
# Mathematical Methods
# ============================================================

export linspace,
       heaviside,
       fortsign,
       mat_split
export central_diff,
       central_diff!,
       central_diff2,
       central_diff2!,
       upwind_diff,
       upwind_diff!,
       unstruct_diff
export maxwell_boltzmann,
       maxwell_boltzmann_prime,
       maxwell_boltzmann_dual,
       maxwell_boltzmann_dual_prime,
       kinetic_entropy
export optimize_closure, realizable_reconstruct
export ylm, rlylm, eval_spherharmonic
export basis_size, eval_sphermonomial

include("math_general.jl")
include("math_entropy.jl")
include("math_closure.jl")
include("sphere_harmonics.jl")
include("sphere_monomials.jl")


text:
using Documenter, Vann

makedocs(
  format = :html,
  sitename = "Vann.jl",
  authors = "Jan Magnusson",
  pages = [
      "Home" => "index.md",
      "Hydrological models" => "hyd_models.md",
      "Snow models" => "snow_models.md",
      "Model calibration" => "calibration.md",
      "Data assimilation" => "data_assim.md"
  ]
)

deploydocs(
  repo  = "github.com/jmgnve/Vann.jl.git",
  target = "build",
  deps = nothing,
  make = nothing,
  julia = "0.5"
)

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
"""Calculates the 2nd order centred gradient in x-direction on any grid (u,v,T or q).
The size of dudx must be m-1,n compared to m,n = size(u)"""
function ∂x!(dudx::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
    m,n = size(dudx)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n, i ∈ 1:m
        dudx[i,j] = u[i+1,j] - u[i,j]
    end
end

"""Calculates the 2nd order centred gradient in y-direction on any grid (u,v,T or q).
The size of dudy must be m,n-1 compared to m,n = size(u)."""
function ∂y!(dudy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
    m,n = size(dudy)
    @boundscheck (m,n+1) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n, i ∈ 1:m
            dudy[i,j] = u[i,j+1] - u[i,j]
    end
end

""" ∇² is the 2nd order centred Laplace-operator ∂/∂x^2 + ∂/∂y^2.
The 1/Δ²-factor is omitted and moved into the viscosity coefficient."""
function ∇²!(du::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
    m, n = size(du)
    @boundscheck (m+2,n+2) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            #        1
            # the 1 -4 1  -stencil in low-precision resilient form
            #        1
            ui1j1 = u[i+1,j+1]
            du[i,j] = ((u[i,j+1] - ui1j1) + (u[i+2,j+1] - ui1j1)) + ((u[i+1,j] - ui1j1) + (u[i+1,j+2] - ui1j1))
        end
    end
end

"""∂x is the 2nd order centred Gradient-operator ∂/∂x with grid spacing Δ (default 1)."""
function ∂x(u::Array{T,2},Δx::Real) where {T<:AbstractFloat}

    m,n = size(u)

    dudx = Array{T,2}(undef,m-1,n)
    one_over_dx = T(1.0/Δx)

    @inbounds for j ∈ 1:n
        for i ∈ 1:m-1
            dudx[i,j] = one_over_dx*(u[i+1,j] - u[i,j])
        end
    end

    return dudx
end

"""∂y is the 2nd order centred Gradient-operator ∂/∂y with grid spacing Δ (default 1)."""
function ∂y(u::Array{T,2},Δy::Real=1) where {T<:AbstractFloat}

    m,n = size(u)

    dudy = Array{T,2}(undef,m,n-1)
    one_over_dy = T(1.0/Δy)

    @inbounds for j ∈ 1:n-1
        for i ∈ 1:m
            dudy[i,j] = one_over_dy*(u[i,j+1] - u[i,j])
        end
    end

    return dudy
end

""" ∇² is the 2nd order centred Laplace-operator ∂/∂x^2 + ∂/∂y^2 with grid spacing Δ (default 1)."""
function ∇²(u::Array{T,2},Δ::Real=1) where {T<:AbstractFloat}

    m, n = size(u)
    du = Array{T,2}(undef,m-2,n-2)

    minus_4 = T(-4.0)
    one_over_dx² = T(1/Δ^2)

    @inbounds for j ∈ 2:n-1
        for i ∈ 2:m-1
            du[i-1,j-1] = one_over_dx²*(minus_4*u[i,j] + u[i,j-1] + u[i,j+1] + u[i-1,j] + u[i+1,j])
        end
    end
    return du
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a2a72c58-9dd5-4e35-8fdd-a23c81b42739
begin
	using AeroMDAO
	using Plots
	using PlutoUI
	plotlyjs()
end

# ╔═╡ 88d02510-af57-11eb-333b-73e3d9b214d9
md"""
# Trimmed drag polars for the whole aircraft
"""

# ╔═╡ 910f1a2d-f977-4a37-b54a-b3e465cfc996
md"### Fuselage drag"

# ╔═╡ af55a6ef-096c-4d7b-bc2a-9365ea1bfb58
CD0_fuse(C_f,f_LD,f_M,S_wetF,S)=C_F*F_LD*f_M*S_wet/S

# ╔═╡ 5e9c31a1-ec37-46b5-add9-7398d3bf68df
begin
	#C_f= #skin friction coef.
	#f_LD= #factor relating to length-to-diameter ratio
	#f_M= #factor relating to Mach no.
	#S_wetF= #wetted area of fuselage
	#S= #surface area
end

# ╔═╡ b5ffb132-c25f-48a5-8f1a-3d0230b3e787
function compute_skin_friction_drag(rho,V,L,miu)
	Re(rho,V,L,miu)=rho*V*L/miu
	#L is length of fuselage for fuselage, for tail it is the chord length

	if Re<200000
		C_friction=1.327/(sqrt(Re))
	end

	if Re>2000000
		C_friction=0.455/((log10(Re))^2.58)
	end

	C_friction
end

# ╔═╡ e0ebfab6-fca0-4b9a-8d0f-5106a35f4242
f_LtoD(L_f,D_f)=1+(60/((L_f/D_f)^3))+0.0025*(L_f/D_f)

# ╔═╡ 3a2567a1-29bd-4a44-8ce8-f1b983c7ec83
f_Mach(M)=1-0.08*M^1.45

# ╔═╡ 0d98227c-401e-4e56-bf54-c20b8d7bd6fa
md"### Wing, horizontal, vertical tail drags"

# ╔═╡ 06594da2-3cea-474d-9797-eafc0c4dfe81
C_D0(C_f, f_M, f_t, Swet, C_dmin, S)=C_f*f_M*f_t*(Swet/S)*((C_dmin/0.004)^0.4)         #obtain CDmin from Cd-Cl curve of airfoil, skin friction is similar to fuselage

# ╔═╡ 51df3b06-5d94-4693-84e5-061e13de1902
f_t(t_cmax)=1+(2.7*t_cmax)+(100*(t_cmax^4)) #obtain

# ╔═╡ d896f8e2-4d42-487d-82b9-392842ed8e00
#find CD of LERX --> can just make it an additional wing

# ╔═╡ 23381ccb-ed20-4d5b-b3e0-f034efaf5a4b
md"### Landing gear drag"

# ╔═╡ 2e758033-337a-49c6-90d2-babccbdda401
#fixed may incure 50% extra drag, fairings can reduce drag

# ╔═╡ 9c3b3453-0d4d-4bc8-9dee-d1132e34a62f
CD0_LG(CD_LG,S_LG,S)=CD_LG*S_LG/S

# ╔═╡ e9d7447f-c6a4-4f7c-811d-1a86f17b94b4
begin
	CD_LG=0.15#with fairings
end

# ╔═╡ 3142c8fc-3f2e-4dc1-b77e-c0aa81d28aee
md"### Strut drag"

# ╔═╡ 694cdd4c-c042-418b-932a-682bec0daf51
CD0_strut(CD0_Si,S_s,S)=CD0_Si*S_s/S

# ╔═╡ b67a3dad-29bc-45c8-b1e1-92a5c8eb8935
#struts with fairings (airfoil section) decreases drag significantly, but at a high cost

# ╔═╡ 4d5b962e-9301-47d0-a27d-6ad775c854b6
#refer to L13 P21 for CD of different shapes

# ╔═╡ a7865faf-383c-467e-83a5-b62615617a1c
md"### Drag polar at clean configuration"

# ╔═╡ abfe110e-dca1-4d5f-9046-0d06912ed401
#C_Dclean=CD0_clean+K*(CL_c)^2

# ╔═╡ 33dadd70-7e82-46a6-8100-91daf073f574
#CL_c is CL at cruise

# ╔═╡ Cell order:
# ╟─88d02510-af57-11eb-333b-73e3d9b214d9
# ╠═a2a72c58-9dd5-4e35-8fdd-a23c81b42739
# ╟─910f1a2d-f977-4a37-b54a-b3e465cfc996
# ╠═af55a6ef-096c-4d7b-bc2a-9365ea1bfb58
# ╠═5e9c31a1-ec37-46b5-add9-7398d3bf68df
# ╠═b5ffb132-c25f-48a5-8f1a-3d0230b3e787
# ╠═e0ebfab6-fca0-4b9a-8d0f-5106a35f4242
# ╠═3a2567a1-29bd-4a44-8ce8-f1b983c7ec83
# ╟─0d98227c-401e-4e56-bf54-c20b8d7bd6fa
# ╠═06594da2-3cea-474d-9797-eafc0c4dfe81
# ╠═51df3b06-5d94-4693-84e5-061e13de1902
# ╠═d896f8e2-4d42-487d-82b9-392842ed8e00
# ╟─23381ccb-ed20-4d5b-b3e0-f034efaf5a4b
# ╠═2e758033-337a-49c6-90d2-babccbdda401
# ╠═9c3b3453-0d4d-4bc8-9dee-d1132e34a62f
# ╠═e9d7447f-c6a4-4f7c-811d-1a86f17b94b4
# ╟─3142c8fc-3f2e-4dc1-b77e-c0aa81d28aee
# ╠═694cdd4c-c042-418b-932a-682bec0daf51
# ╠═b67a3dad-29bc-45c8-b1e1-92a5c8eb8935
# ╠═4d5b962e-9301-47d0-a27d-6ad775c854b6
# ╟─a7865faf-383c-467e-83a5-b62615617a1c
# ╠═abfe110e-dca1-4d5f-9046-0d06912ed401
# ╠═33dadd70-7e82-46a6-8100-91daf073f574

max_stars_repo_path: Drag Polar.jl

max_stars_repo_name: rahmanfahd/aircraft-design

max_issues_repo_path: Drag Polar.jl

max_issues_repo_name: rahmanfahd/aircraft-design

max_forks_repo_path: Drag Polar.jl

max_forks_repo_name: rahmanfahd/aircraft-design


text:
# Define the new measure data type with <: AbstractMeasureData
struct NewMeasureData <: AbstractMeasureData
    attr1::String # REPLACE WITH ACTUAL ATTRIBUTE
    attr2::DiscreteMeasureData # REPLACE WITH ACTUAL ATTRIBUTE
    # ADD MORE ATTRIBUTES AS NEEDED
    # constructor
    function NewMeasureData(attr1::String, attr2::DiscreteMeasureData)
        # INSERT CHECKS AND OTHER CONSTRUCTOR METHODS HERE
        return new(attr1, attr2) # REPLACE WITH ACTUAL ATTRIBUTES
    end
end

# Extend parameter_refs to return the parameter(s) being measured by a measure using NewMeasureData
function InfiniteOpt.parameter_refs(data::NewMeasureData)::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}}
    return data.attr2.parameter_refs # REPLACE WITH ACTUAL PARAMETER LOCATION
end

# Extend expand_measure to return the finite reformulation of a measure using NewMeasureData
function InfiniteOpt.expand_measure(expr::JuMP.AbstractJuMPScalar,
                                    data::NewMeasureData,
                                    write_model::JuMP.AbstractModel
                                    )::JuMP.AbstractJuMPScalar
    # INSERT APPROPRIATE METHODS HERE
    # USING make_point_variable_ref AND make_reduced_variable_ref MAY BE USEFUL
    return expand_measure(expr, data.attr2, write_model) # REPLACE ACTUAL RESULT
end

# Extend add_supports_to_parameters to add supports to parameters according to the NewMeasureData
# This is only optional if the measure is always analytic
function InfiniteOpt.add_supports_to_parameters(data::NewMeasureData)
    return add_supports_to_parameters(data.attr2) # REPLACE WITH ACTUAL SUPPORT ADDITION STEPS
end

# Extend supports to return any infinite parameter supports employed by NewMeasureData
# This is only optional if the new abstraction doesn't use supports at all
function InfiniteOpt.supports(data::NewMeasureData)::Vector
    return data.attr2.supports # REPLACE WITH ACTUAL LOCATION
end

# Extend measure_data_in_hold_bounds to determine if NewMeasureData is in the
# domain of hold variable bounds. (Enables hold variable error checking)
function InfiniteOpt.measure_data_in_hold_bounds(data::NewMeasureData,
                                                 bounds::ParameterBounds)::Bool
    # INSERT ACTUAL CHECK HERE
    in_bounds = measure_data_in_hold_bounds(data.attr2, bounds) # REPLACE WITH ACTUAL RESULT
    return in_bounds
end

# Extend coefficients to return the coefficients stored in NewMeasureData if appropriate
# This is optional (returns empty vector otherwise)
function InfiniteOpt.coefficients(data::NewMeasureData)::Vector{Float64}
    return data.attr2.coefficients # REPLACE WITH ACTUAL LOCATION
end

# Extend weight_function to return the weight function stored in NewMeasureData if appropriate
# This is optional (returns default_weight otherwise)
function InfiniteOpt.weight_function(data::NewMeasureData)::Function
    return data.attr2.weight_function # REPLACE WITH ACTUAL LOCATION
end

# Make a convenient measure constructor function for our new measure type
# This should employ measure(expr, data)
function new_measure(expr::JuMP.AbstractJuMPScalar, param::GeneralVariableRef,
                     lb::Number, ub::Number; name::String = "NewMeas",
                     num_supports::Int = 10)::GeneralVariableRef # REPLACE ARGS WITH ACTUAL DESIRED
    # INSERT RELAVENT CHECKS AND OPERATIONS HERE
    # REPLACE BELOW WITH ACTUAL CONSTRUCTION
    attr2 = DiscreteMeasureData(param, ones(num_supports), [lb + (ub - lb) / (num_supports - 1) * i for i in 1:num_supports]) # just an example
    data = NewMeasureData(name, attr2) # REPLACE WITH ACTUAL
    # built the measure using the built-in constructor
    return measure(expr, data)
end


text:
module HTMexperiments

using HierarchicalTemporalMemory
using Random, Statistics
import Lazy: @>
using UnPack

include("metrics.jl")
export sparseness, entropy, noiserobustness, stabilitySP
include("makedata.jl")
export randomSDR

end

max_stars_repo_name: Oblynx/HTMexperiments

max_issues_repo_name: Oblynx/HTMexperiments

max_forks_repo_name: Oblynx/HTMexperiments


text:
#Test matrix with cyclones

matrixwithcyclones = zeros(Float64,20,20)
matrixwithcyclones[2,1] = -2 #Center
matrixwithcyclones[8:10,5:7] .= -1
matrixwithcyclones[9,6] = -2 #Center
matrixwithcyclones[17:19,17:19] .= -1
matrixwithcyclones[18,18] = -2

## Fake pressure perturbation

@testset "Cyclone detection tools" begin
    @testset "Detect centers" begin
        centerstest =  AvailablePotentialEnergyFramework.findcyclonecenters_aspressureminima(matrixwithcyclones,0)
        @test length(centerstest) == 2
        @test centerstest[1] == CartesianIndex(9, 6)
        @test centerstest[2] == CartesianIndex(18, 18)
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mobile_analytics
using AWS.Compat
using AWS.UUIDs

"""
    put_events(events, x-amz-_client-_context)
    put_events(events, x-amz-_client-_context, params::Dict{String,<:Any})

The PutEvents operation records one or more events. You can have up to 1,500 unique custom
events per app, any combination of up to 40 attributes and metrics per custom event, and
any number of attribute or metric values.

# Arguments
- `events`: An array of Event JSON objects
- `x-amz-_client-_context`: The client context including the client ID, app title, app
  version and package name.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"x-amz-Client-Context-Encoding"`: The encoding used for the client context.
"""
function put_events(
    events, x_amz_Client_Context; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mobile_analytics(
        "POST",
        "/2014-06-05/events",
        Dict{String,Any}(
            "events" => events,
            "headers" => Dict{String,Any}("x-amz-Client-Context" => x_amz_Client_Context),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_events(
    events,
    x_amz_Client_Context,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mobile_analytics(
        "POST",
        "/2014-06-05/events",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "events" => events,
                    "headers" =>
                        Dict{String,Any}("x-amz-Client-Context" => x_amz_Client_Context),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
export StaticGaussianDriver

type StaticGaussianDriver{A} <: DriverModel{A,IntegratedContinuous}
    C::IntegratedContinuous
    P::MvNormal
end

get_name(::StaticGaussianDriver) = "StaticGaussian"
action_context(model::StaticGaussianDriver) = model.C
function Base.rand{A}(model::StaticGaussianDriver{A})
    a = rand(model.P)
    convert(A, a)
end
Distributions.pdf{A}(model::StaticGaussianDriver{A}, a::A) = pdf(model.P, convert(Vector{Float64}, a))
Distributions.logpdf{A}(model::StaticGaussianDriver{A}, a::A) = logpdf(model.P, convert(Vector{Float64}, a))


text:
###############################################################################
"""

# General

This module contains Julia code for performing uncertainty analysis of camera
calbiration and georeferencing, as described in Schweitzer & Cowen, 2021.

# References
Schweitzer, S. A., & Cowen, E. A. (2021). "Instantaneous river-wide water surface
velocity field measurements at centimeter scales using infrared quantitative
image velocimetry". Water Resources Research, 57,
e2020WR029279. https://doi.org/10.1029/2020WR029279

# License
This software is released under the MIT license, however, the authors request to
be credited if you use the software, preferably by citing the above-mentioned
Schweitzer & Cowen (2021) paper.

MIT License

Copyright (c) 2022 Seth A. Schweitzer

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

"""
module GeoreferencingUncertaintySimulations
###############################################################################

###############################################################################
using Statistics
using ProgressMeter
#
import ImageRectification; Rect = ImageRectification
import ExtrinsicCameraCalibration; Calib = ExtrinsicCameraCalibration
import RemoveImageDistortion; RadialDist = RemoveImageDistortion
###############################################################################

###############################################################################
export SimulateDLTFromGCPsWithAddedNoise
###############################################################################

###############################################################################
function GeorefGCPSimulationWrapper(TrueGCPxyz, TrueGCPuv, TrueWaterSurfaceElev, u, v, χinitial, ι;
                                     xyzNoiseFactor = 0.0,
                                     uvNoiseFactor = 0.0,
                                     wseNoiseFactor = 0.0,
                                     AddRadialDistortion = false,
                                     nSimulations = 500,
                                     MaxIterations = 50)
    InfoString = "n Simulations = $nSimulations\nAdding the following types of noise:"
    InfoString = (xyzNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nXYZ $xyzNoiseFactor"])
    InfoString = (uvNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nUV $uvNoiseFactor"])
    InfoString = (wseNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nWSE $wseNoiseFactor"])
    InfoString = join([InfoString, "\nMaxIterations $MaxIterations"])
    @info(InfoString)
    #
    # Calculate the calibration without adding any noise:
    χcalculated = Calib.HollandEtAlExtrinsicCalibration(TrueGCPuv, TrueGCPxyz, χinitial, ι,
                                                            MaxIterations=MaxIterations)
    # Calculate the DLTs without adding any noise:
    DLTcalculated = Calib.CalculateDLTcoefficients(χcalculated, ι)
    # Calculate the x,y coordinates without adding any noise:
    Calculated_x, Calculated_y = Rect.FindXYFromUVZ(u,v, TrueWaterSurfaceElev, DLTcalculated)
    # Generate a set of DLTs that by running the calibraiotn with noise added to
    # the GCP uv and/or xyz coordinates:
    DLTsimulated, GCP_perturb_uv_xyz = SimulateDLTFromGCPsWithAddedNoise(TrueGCPxyz, TrueGCPuv, u, v, χinitial, ι;
                                                                         xyzNoiseFactor = xyzNoiseFactor,
                                                                         uvNoiseFactor = uvNoiseFactor,
                                                                         AddRadialDistortion = AddRadialDistortion,
                                                                         nSimulations = nSimulations,
                                                                         MaxIterations = MaxIterations)
    # Calculate the x,y coordinates of pixels that would be found using the
    # simulated DLTs (i.e., with noise added to the GCP coordinates). Optionally
    # add more noise related to uncertainty in determining the water surface
    # elevation at this step:
    Modified_x, Modified_y, Perturb_wse = CalculateModifiedxyCoordinates(u,v,
                                                                         TrueWaterSurfaceElev, DLTsimulated;
                                                                         wseNoiseFactor=wseNoiseFactor)
    Perturbations = (GCPxyz=GCP_perturb_uv_xyz[1], GCPuv=GCP_perturb_uv_xyz[2], wse=Perturb_wse)
    return (;Calculated_x, Calculated_y, Modified_x, Modified_y, InfoString, Perturbations)
end
###############################################################################

###############################################################################
"""
    function SimulateDLTFromGCPsWithAddedNoise(TrueGCPxyz, TrueGCPuv, u, v, z, χinitial, ι;
                                       xyzNoiseFactor = 1.0, uvNoiseFactor = 1.0,
                                       AddRadialDistortion = false,
                                       nSimulations = 500,
            GeorefSim.GeorefGCPSimulationWrapper2                           MaxIterations = 50)

    Calculates extrinsic calibration for a set of inputs with random noise.

    The function adds normally distributed noise with mean zero and input-dependent
    standard deviation to the pixel and/or physical coordinates of the GCPs.
    """
function SimulateDLTFromGCPsWithAddedNoise(TrueGCPxyz, TrueGCPuv, u, v, χinitial, ι;
                                   xyzNoiseFactor = 1.0, uvNoiseFactor = 1.0,
                                   AddRadialDistortion = false,
                                   nSimulations = 500,
                                   MaxIterations = 50)
    xyzNoiseFactor = isa(xyzNoiseFactor, Number) ? repeat([xyzNoiseFactor], 3) : xyzNoiseFactor
    uvNoiseFactor = isa(uvNoiseFactor, Number) ? repeat([uvNoiseFactor], 2) : uvNoiseFactor
    # DLTs = fill(NaN, 11, nSimulations)
    DLTs = Array{Any}(undef, nSimulations)
    GCP_xyz_perturb = [ hcat(randn(size(TrueGCPxyz,1)) .* xyzNoiseFactor[1],
                             randn(size(TrueGCPxyz,1)) .* xyzNoiseFactor[2],
                             randn(size(TrueGCPxyz,1)) .* xyzNoiseFactor[3]) for i in 1:nSimulations ]
    GCP_uv_perturb = [ hcat(randn(size(TrueGCPuv,1)) .* uvNoiseFactor[1],
                            randn(size(TrueGCPuv,1)) .* uvNoiseFactor[2]) for i in 1:nSimulations ]
    InputGCPuv = ( AddRadialDistortion ?
        hcat(RadialDist.SimulateRadialDistortion(TrueGCPuv[:,1], TrueGCPuv[:,2], ι)...) :
        TrueGCPuv )
    @showprogress 1 "simulating DLT uncertainty..." for i in 1:nSimulations
        χmodified = Calib.HollandEtAlExtrinsicCalibration(
            InputGCPuv .+ GCP_uv_perturb[i],
            TrueGCPxyz .+ GCP_xyz_perturb[i],
            χinitial,
            ι;
            MaxIterations=MaxIterations,
            PrintoutStats=false)
        DLTs[i] = Calib.CalculateDLTcoefficients(χmodified, ι)
    end
    return DLTs, (GCP_xyz_perturb=GCP_xyz_perturb, GCP_uv_perturb=GCP_uv_perturb)
end
###############################################################################

###############################################################################
""" Calculate the uv -> xyz coordinate conversion using each of the sets of DLTs """
function CalculateModifiedxyCoordinates(u,v, WaterSurfaceElev, DLTs;
                                        wseNoiseFactor = 0.0)
    Perturb_wse = randn(size(DLTs,2)) .* wseNoiseFactor
    Modified_x, Modified_y = Rect.FindXYFromUVZ_DLT_z_sensitivity(u,v,
                                                                  WaterSurfaceElev .+ Perturb_wse,
                                                                  DLTs)
    return Modified_x, Modified_y, Perturb_wse
end
###############################################################################

###############################################################################
"""
Calculate the distance b/w the simulated and "true" xy coordinates of each pixel
footprint:
"""
function CalculateDistanceFromOriginal(Modified_x, Calculated_x, Modified_y, Calculated_y)
    return hypot.((Modified_x .- Calculated_x), (Modified_y .- Calculated_y))
end
###############################################################################

###############################################################################
function GeorefExtrinsicParamsSimulationWrapper(χinitial, ι,
                                                TrueWaterSurfaceElev,
                                                u, v;
                                                φNoiseFactor = 0.0,
                                                τNoiseFactor = 0.0,
                                                σNoiseFactor = 0.0,
                                                x_cNoiseFactor = 0.03,
                                                y_cNoiseFactor = 0.03,
                                                z_cNoiseFactor = 0.05,
                                                wseNoiseFactor = 0.03,
                                                AddRadialDistortion = false,
                                                nSimulations = 500,
                                                MaxIterations = 50)
    InfoString = "n Simulations = $nSimulations\nAdding the following types of noise:"
    roundedrad2deg(t) = round(rad2deg(t), digits=3)
    InfoString = (φNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nφ $(roundedrad2deg(φNoiseFactor)) [deg]"])
    InfoString = (τNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nτ $(roundedrad2deg(τNoiseFactor)) [deg]"])
    InfoString = (σNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nσ $(roundedrad2deg(σNoiseFactor)) [deg]"])
    InfoString = (z_cNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nz_c $z_cNoiseFactor"])
    InfoString = (wseNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nWSE $wseNoiseFactor"])
    InfoString = join([InfoString, "\nMaxIterations $MaxIterations"])
    @info(InfoString)
    #
    DLTcalculated = Calib.CalculateDLTcoefficients(χinitial, ι)
    # Calculate the x,y coordinates without adding any noise:
    Calculated_x, Calculated_y = Rect.FindXYFromUVZ(u,v, TrueWaterSurfaceElev, DLTcalculated)
    # Generate a set of DLTs that by running the calibraiotn with noise added to
    # the GCP uv and/or xyz coordinates:
    DLTsimulated, ExtrinsicParamsPerturbs = SimulateDLTFromExtrinsicParamsWithAddedNoise(χinitial, ι, u, v;
                                                                                φNoiseFactor = φNoiseFactor,
                                                                                τNoiseFactor = τNoiseFactor,
                                                                                σNoiseFactor = σNoiseFactor,
                                                                                z_cNoiseFactor = z_cNoiseFactor,
                                                                                nSimulations = nSimulations,
                                                                                MaxIterations = MaxIterations)
    # Simulate radial distortion by changing the pixel coordinates of the input grid:
    u_sim, v_sim = AddRadialDistortion ? RadialDist.SimulateRadialDistortion(u, v, ι) : (u, v)
    # Calculate the x,y coordinates of pixels that would be found using the
    # simulated DLTs (i.e., with noise added to the extrinsic
    # calibration). Optionally add more noise related to uncertainty in
    # determining the water surface elevation at this step:
    Modified_x, Modified_y, Perturb_wse = CalculateModifiedxyCoordinates(u_sim, v_sim,
                                                                         TrueWaterSurfaceElev,
                                                                         DLTsimulated;
                                                                         wseNoiseFactor=wseNoiseFactor)
    # Calculate the horizontal difference in position between the set of xy
    # coordinates calculated with random noise added, and the set of coordinates
    # found from the original input (no additional noise added to GCP or water
    # surface elevation coordinates):
    DeltaDist = CalculateDistanceFromOriginal(Modified_x, Calculated_x, Modified_y, Calculated_y)
    # return DeltaDist, InfoString, (; ExtrinsicParamsPerturbs..., Perturb_wse)
    return (; Calculated_x, Calculated_y, Modified_x, Modified_y, InfoString, ExtrinsicParamsPerturbs..., Perturb_wse)
end
###############################################################################

###############################################################################
function SimulateDLTFromExtrinsicParamsWithAddedNoise(
    χinitial, ι, u, v;
    φNoiseFactor = 0.0,
    τNoiseFactor = 0.0,
    σNoiseFactor = 0.0,
    x_cNoiseFactor = 0.05,
    y_cNoiseFactor = 0.05,
    z_cNoiseFactor = 0.1,
    nSimulations = 500,
    MaxIterations = 50)
    #
    DLTs = Array{Any}(undef, nSimulations)
    φ_perturb = randn(nSimulations) .* φNoiseFactor
    τ_perturb = randn(nSimulations) .* τNoiseFactor
    σ_perturb = randn(nSimulations) .* σNoiseFactor
    x_c_perturb = randn(nSimulations) .* x_cNoiseFactor
    y_c_perturb = randn(nSimulations) .* y_cNoiseFactor
    z_c_perturb = randn(nSimulations) .* z_cNoiseFactor
    @showprogress 1 "simulating DLT uncertainty from extrinsic params..." for i in 1:nSimulations
        χmodified = ExtrinsicCameraCalibration.ExtrinsicCalibrationParameters(
            χinitial.φ + φ_perturb[i],
            χinitial.τ + τ_perturb[i],
            χinitial.σ + σ_perturb[i],
            χinitial.f,
            χinitial.x_c + x_c_perturb[i],
            χinitial.y_c + y_c_perturb[i],
            χinitial.z_c + z_c_perturb[i])
        #
        DLTs[i] = Calib.CalculateDLTcoefficients(χmodified, ι)
    end
    return DLTs, (; φ_perturb, τ_perturb, σ_perturb, x_c_perturb, y_c_perturb, z_c_perturb)
end
###############################################################################

###############################################################################
function SimulateAndReturnXYcoordinates(TrueGCPxyz, TrueGCPuv, TrueWaterSurfaceElev, u, v, χinitial, ι;
                                        xyzNoiseFactor = 0.0,
                                        uvNoiseFactor = 0.0,
                                        wseNoiseFactor = 0.0,
                                        AddRadialDistortion = false,
                                        nSimulations = 500,
                                        MaxIterations = 50)
    InfoString = "n Simulations = $nSimulations\nAdding the following types of noise:"
    InfoString = (xyzNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nXYZ $xyzNoiseFactor"])
    InfoString = (uvNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nUV $uvNoiseFactor"])
    InfoString = (wseNoiseFactor == 0.0) ? InfoString : join([InfoString, "\nWSE $wseNoiseFactor"])
    InfoString = join([InfoString, "\nMaxIterations $MaxIterations"])
    @info(InfoString)
    #
    # Calculate the calibration results without adding any noise:
    χcalculated = Calib.HollandEtAlExtrinsicCalibration(TrueGCPuv, TrueGCPxyz, χinitial, ι,
                                                            MaxIterations=MaxIterations)
    # Calculate the DLTs without adding any noise:
    # DLTcalculated = Calib.CalculateDLTcoefficients(χcalculated, ι)
    DLTcalculated = Calib.CalculateDLTcoefficients(χcalculated, ι)
    # Calculate the x,y coordinates without adding any noise:
    Calculated_x, Calculated_y = Rect.FindXYFromUVZ(u,v, TrueWaterSurfaceElev, DLTcalculated)
    # Generate a set of DLTs that by running the calibraiotn with noise added to
    # the GCP uv and/or xyz coordinates:
    DLTsimulated, GCP_perturb_uv_xyz = SimulateDLTFromGCPsWithAddedNoise(TrueGCPxyz, TrueGCPuv, u, v, χinitial, ι;
                                                                         xyzNoiseFactor = xyzNoiseFactor,
                                                                         uvNoiseFactor = uvNoiseFactor,
                                                                         AddRadialDistortion = AddRadialDistortion,
                                                                         nSimulations = nSimulations,
                                                                         MaxIterations = MaxIterations)
    # Calculate the x,y coordinates of pixels that would be found using the
    # simulated DLTs (i.e., with noise added to the GCP coordinates). Optionally
    # add more noise related to uncertainty in determining the water surface
    # elevation at this step:
    Modified_x, Modified_y, Perturb_wse = CalculateModifiedxyCoordinates(u,v,
                                                                         TrueWaterSurfaceElev, DLTsimulated;
                                                                         wseNoiseFactor=wseNoiseFactor)
    return Calculated_x, Calculated_y, Modified_x, Modified_y# , ρ, θ
end
###############################################################################

###############################################################################
randinrange(minval, maxval) = (minval + rand() * (maxval-minval))
randinrange(minval, maxval, n) = [ randinrange(minval, maxval) for _ in 1:n ]
blurintegers(x::Array) = x .+ rand(size(x)) .- 0.5 # adds a random value b/w (-0.5,0.5), to avoind integer locking
blurintegers(x::Number) = x .+ rand() .- 0.5 # adds a random value b/w (-0.5,0.5), to avoind integer locking
function GenerateGCPuvCoordinatesForSimulation(nGCPs,
                                               u,v,
                                               GCPzrange;
                                               EdgeWidth = 50, # how far in from the edges can GCPs be?
                                               IdealEdge = 10, # distance from edge for fixed, "ideal" GCPs
                                               GCPdistribution = "perimeter") # how to distribute the GCPs throughout the frame
    #
    ### Define GCPs:
    #
    # We'll randomly select u,v,z coordinates for the GCPs, possibly with some
    # constraints, and then calculate the true xyz coordinates from these
    # values. This will be the basis for the uncertainty calculations.
    #
    # GCPdistribution = "perimeter" # "left_and_right" # "corners", "perimeter", "anywhere", "left_and_right", "top_and_bottom"
    if GCPdistribution == "anywhere"
        # # ---------- Randomly distributed across image: ----------
        GCPu = randinrange(extrema(u)..., nGCPs)
        GCPv = randinrange(extrema(v)..., nGCPs)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "corners"
        # # ---------- Randomly in corners of image: ----------
        # First four GCPs, one in each corner:
        GCPu = [rand(min(u...):EdgeWidth),
                rand((max(u...)-EdgeWidth):max(u...)),
                rand(min(u...):EdgeWidth),
                rand((max(u...)-EdgeWidth):max(u...))]
        GCPv = [rand(min(v...):EdgeWidth),
                rand(min(v...):EdgeWidth),
                rand((max(v...)-EdgeWidth):max(v...)),
                rand((max(v...)-EdgeWidth):max(v...))]
        # Any additional GCPs get distributed among the corners randomly:
        if nGCPs > 4
            GCPu = vcat(GCPu, rand(union(min(u...):EdgeWidth, (max(u...)-EdgeWidth):max(u...)), nGCPs-4))
            GCPv = vcat(GCPv, rand(union(min(v...):EdgeWidth, (max(v...)-EdgeWidth):max(v...)), nGCPs-4))
        end
        GCPu =  blurintegers(GCPu)
        GCPv =  blurintegers(GCPv)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        # #
        # trueGCPuvz = hcat(Random.shuffle!(GCPu), Random.shuffle!(GCPv), Random.shuffle!(GCPz))
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "perimeter"
        # ---------- Randomly along edges of image: ----------
        # set up half of the GCPs along the top/bottom edges, and half along the
        # left/right edges. If an odd number of GCPs is requested, generate
        # ceil(nGCPs/2) of each, leadaing to a total of nGCPs+1, then randomly
        # select nGCPs out of the set.
        #
        # left/right
        GCPu_lr = rand(union(min(u...):EdgeWidth, (max(u...)-EdgeWidth):max(u...)), ceil(Int, nGCPs/2)) |> blurintegers
        GCPv_lr = randinrange(extrema(v)..., ceil(Int, nGCPs/2))
        # top/bottom:
        GCPu_tb = randinrange(extrema(u)..., ceil(Int, nGCPs/2))
        GCPv_tb = rand(union(min(v...):EdgeWidth, (max(v...)-EdgeWidth):max(v...)), ceil(Int, nGCPs/2)) |> blurintegers
        #
        GCPu = vcat(GCPu_lr, GCPu_tb)
        GCPv = vcat(GCPv_lr, GCPv_tb)
        GCPz = randinrange(extrema(GCPzrange)..., length(GCPu))
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)[StatsBase.sample(1:size(GCPu,1), nGCPs; replace=false),:]
    elseif GCPdistribution == "perimeter_ideal"
        # Place one GCP in the center of each side of the image (IdealEdge pixels
        # from edge), then add additional GCPs up to nGCPs in the corners:
        GCPu = [round(Int,mean(u)), round(Int,mean(u)),   IdealEdge,          maximum(u)-IdealEdge]
        GCPv = [IdealEdge,          maximum(v)-IdealEdge, mean(v),            mean(v)             ]
        if nGCPs > 4
            GCPu = [GCPu..., IdealEdge]
            GCPv = [GCPv..., IdealEdge]
        end
        if nGCPs > 5
            GCPu = [GCPu..., maximum(u)-IdealEdge]
            GCPv = [GCPv..., maximum(v)-IdealEdge]
        end
        if nGCPs > 6
            GCPu = [GCPu..., IdealEdge]
            GCPv = [GCPv..., maximum(v)-IdealEdge]
        end
        if nGCPs > 7
            GCPu = [GCPu..., maximum(u)-IdealEdge]
            GCPv = [GCPv..., IdealEdge]
        end
        GCPu = round.(Int, GCPu)
        GCPv = round.(Int, GCPv)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "left_and_right"
        # ---------- Randomly along left and right edges of image: ----------
        # left/right
        # GCPu = rand(union(min(u...):EdgeWidth, (max(u...)-EdgeWidth):max(u...)), nGCPs) |> blurintegers
        GCPu = union( rand(min(u...):EdgeWidth, floor(Int,nGCPs/2)),
                      rand((max(u...)-EdgeWidth):max(u...), ceil(Int,nGCPs/2))) |> blurintegers
        GCPv = randinrange(extrema(v)..., nGCPs)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "right"
        # ---------- Randomly along right edge of image: ----------
        # left/right
        GCPu = rand( (max(u...)-EdgeWidth):max(u...), nGCPs) |> blurintegers
        GCPv = randinrange(extrema(v)..., nGCPs)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "left"
        # ---------- Randomly along left edge of image: ----------
        # left/right
        GCPu = rand(min(u...):EdgeWidth, nGCPs) |> blurintegers
        GCPv = randinrange(extrema(v)..., nGCPs)
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "top_and_bottom"
        # ---------- Randomly along top and bottom edges of image: ----------
        GCPu = randinrange(extrema(u)..., nGCPs)
        GCPv = union(rand(min(v...):EdgeWidth, floor(Int,nGCPs/2)),
                     rand( (max(v...)-EdgeWidth):max(v...), ceil(Int, nGCPs/2))) |> blurintegers
        # GCPv = rand(union(min(v...):EdgeWidth, (max(v...)-EdgeWidth):max(v...)), nGCPs) |> blurintegers
        #
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "top"
        # ---------- Randomly along top edge of image: ----------
        GCPu = randinrange(extrema(u)..., nGCPs)
        GCPv = rand((max(v...)-EdgeWidth):max(v...), nGCPs) |> blurintegers
        #
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif GCPdistribution == "bottom"
        # ---------- Randomly along bottom edge of image: ----------
        GCPu = randinrange(extrema(u)..., nGCPs)
        GCPv = rand(min(v...):EdgeWidth, nGCPs) |> blurintegers
        #
        GCPz = randinrange(extrema(GCPzrange)..., nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    elseif lowercase(GCPdistribution) == "none"
        GCPu = fill(NaN, nGCPs)
        GCPv = fill(NaN, nGCPs)
        GCPz = fill(NaN, nGCPs)
        trueGCPuvz = hcat(GCPu, GCPv, GCPz)
    end
    return (trueGCPuv=trueGCPuvz[:,1:2], trueGCPz=trueGCPuvz[:,3])
end
###############################################################################

###############################################################################
end                             # end of module
###############################################################################

max_stars_repo_name: saschweitzer/camera-calibration

max_issues_repo_name: saschweitzer/camera-calibration

max_forks_repo_name: saschweitzer/camera-calibration