text:
using ByteParsers
using ByteParsers: FixEmpty
using Base.Test
using InferenceUtilities

const sc = UInt8(';')           # semicolon, used as separator in tests

"""
   x ≅ y (type with \\cong)

A less strict comparison of MaybeParsed types: positions have to match (`==`),
values only when parsed (`isequal`). Type parameters don't have to match
otherwise. Useful for testing.
"""
≅(x, y) = isequal(x, y)

"""
    x ≂ y (type with \\eqsim)

A very strict comparison of MaybeParsed types: type parameters and positions
have to match (`==`), values only when parsed (`isequal`). Useful for testing.
"""
≂(x, y) = false

≂(x::MaybeParsed{T}, y::MaybeParsed{T}) where {T} = isequal(x, y)

@testset "MaybeParsed methods" begin
    @test MaybeParsed{Int}(-1) ≂ MaybeParsed{Int}(-1)
    @test MaybeParsed(-1, 1) ≂ MaybeParsed(-1, 2)
    @test MaybeParsed(4, 10) ≂ MaybeParsed(4, 10)
    @test !(MaybeParsed{Int}(-3) ≅ MaybeParsed{Int}(-4))
    @test !(MaybeParsed(2, 3) ≅ MaybeParsed(2, 4))
    @test !(MaybeParsed(3, 3) ≅ MaybeParsed(2, 3))
    @test unsafe_get(MaybeParsed(1, 10)) ≡ 10
    @test_throws UndefRefError unsafe_get(MaybeParsed{String}(1))
    @test getpos(MaybeParsed(5, 10)) == 5
    let mp = MaybeParsed(2,5)
        @test getpos(mp) == 2
        @test isparsed(mp)
    end
    let mp = MaybeParsed(ByteParsers.pos_to_error(7), 19)
        @test getpos(mp) == 7
        @test !isparsed(mp)
    end
end

@testset "integer parsing" begin
    @test parsenext(PosInteger(), "119;", 1, ';') ≂ MaybeParsed(5, 119)
    @test parsenext(PosInteger(), "222;xx;", 5, ';') ≂ MaybeParsed{Int}(-5)
    @test parsenext(PosFixedInteger(2), "22;77", 4, ';') ≂ MaybeParsed(6, 77)
    @test parsenext(PosInteger(), b"111", 1, sc) ≂ MaybeParsed{Int}(-4)
    @test parsenext(PosFixedInteger(3), b"11x", 1, sc) ≂ MaybeParsed{Int}(-3)
    @test parsenext(PosFixedInteger(3), b"11", 1, sc) ≂ MaybeParsed{Int}(-3)
    @test parsenext(PosInteger(), b";", 1, sc) ≂ MaybeParsed{Int}(-1)
    @test @isinferred parsenext(PosInteger(), b"11", 1, sc)
    @test @isinferred parsenext(PosFixedInteger(3), b"11", 1, sc)
    # typed
    for T ∈ [Int8, Int16, Int32, Int64, BigInt]
        @test parsenext(PosInteger(T), b"19;", 1, sc) ≂ MaybeParsed(4, T(19))
    end
    @test_throws InexactError parsenext(PosInteger(Int8), b"300", 1, sc)
    @test repr(PosInteger(Int16, Int32)) ==
        "Parse positive integer into Int16 (accumulating in Int32)"
    @test repr(PosFixedInteger(4, Int16)) ==
        "Parse positive integer of 4 digits into Int16"
end

@testset "skip or verbatim strings" begin
    @test parsenext(Skip(), b"xxx;yyyy;", 5, sc) ≂ MaybeParsed(10, nothing)
    @test parsenext(ViewBytes(), b"xxx;yyyy;", 5, sc) ≅ MaybeParsed(10, b"yyyy"[:])
    @test parsenext(Skip(), b"nosep", 1, sc) ≂ MaybeParsed{Void}(-6)
    @test parsenext(ViewBytes(), b"nosep", 1, sc) ≅ MaybeParsed{String}(-6)
    @test @isinferred parsenext(Skip(), b"something", 1, sc)
    @test @isinferred parsenext(ViewBytes(), b"something", 1, sc)
    @test repr(Skip()) == "skip field (parsed as `nothing`)"
    @test repr(ViewBytes()) == "return field as bytes (SubArray)"
end

@testset "dates" begin
    @test parsenext(DateYYYYMMDD(), b"xxx;19800101;", 5, sc) ≂
        MaybeParsed(14, Date(1980, 1, 1))
    @test parsenext(DateYYYYMMDD(), b"19800100;", 1, sc) ≂ MaybeParsed{Date}(-7)
    @test parsenext(DateYYYYMMDD(), b"19809900;", 1, sc) ≂ MaybeParsed{Date}(-5)
    @test parsenext(DateYYYYMMDD(), b"19800101", 1, sc) ≂ MaybeParsed{Date}(-10)
    @test @isinferred parsenext(DateYYYYMMDD(), b"19800101;", 1, sc)
    @test repr(DateYYYYMMDD()) == "Parse date in YYYYMMDD format"
end

@testset "fix empty values" begin
    p = FixEmpty(-1, PosInteger())
    @test parsenext(p, b";", 1, sc) ≂ MaybeParsed{Int}(2, -1)
    @test parsenext(p, b"2;", 1, sc) ≂ MaybeParsed{Int}(3, 2)
    @test parsedtype(p) ≡ Int
    @test repr(p) == "parse empty fields as -1, otherwise " * repr(PosInteger())
end

@testset "parsed types" begin
    @test parsedtype(PosInteger()) ≡ Int
    @test parsedtype(PosFixedInteger(2)) ≡ Int
    @test parsedtype(Skip()) ≡ Void
    @test parsedtype(ViewBytes()) ≡ typeof(@view b"xx"[1:1])
end

@testset "parsing line" begin
    line1 = b";1212;skipped;kept;"
    mp1 = MaybeParsed(length(line1)+1, (-1, 1212, b"kept")) # value skipped
    parser1 = Line(FixEmpty(-1, PosInteger()), PosInteger(), Skip(), ViewBytes())
    @test length(parser1) == 3
    @test parser1[1] ≅ FixEmpty(-1, PosInteger())
    @test parser1[2] == PosInteger()
    @test parser1[3] == ViewBytes()
    @test parsenext(parser1, line1, 1, sc) ≅ mp1
    @test parsenext(Line(PosInteger()), b"1bad;", 1, sc) ≅
        MaybeParsed{Tuple{Int}}(-2)
    @test @isinferred parsenext(parser1, line1, 1, sc)
    @test repr(parser1) ==
        ("parse line as the following:\n    " *
         repr(FixEmpty(-1, PosInteger())) * "\n    " *
         repr(PosInteger()) * "\n    " *
         repr(Skip()) * "\n    " *
         repr(ViewBytes()))
end


text:
# ---
# title: Statistical visualizations
# cover: assets/statistical_visualizations.png
# description: Violin plot, boxplot, qqplot
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie, PalmerPenguins, DataFrames
set_aog_theme!() #src

penguins = dropmissing(DataFrame(PalmerPenguins.load()))

data(penguins) * visual(Violin) *
    mapping(:species, :bill_depth_mm, color=:sex, dodge=:sex) |> draw

#

plt = data(penguins) * visual(Violin, datalimits=extrema)
plt *= mapping(:species, :bill_depth_mm, color=:sex, side=:sex, dodge=:island)
fg = draw(plt, axis=(limits=((0.5, 3.5), nothing),))

#

data(penguins) * visual(BoxPlot, show_notch=true) *
    mapping(:species, :bill_depth_mm, color=:sex, dodge=:sex) |> draw

#

data(penguins) *
    mapping(:bill_length_mm, :bill_depth_mm, col=:sex) *
    visual(QQPlot, qqline=:fit) |> draw

# save cover image #src
mkpath("assets") #src
save("assets/statistical_visualizations.png", fg) #src

ext: jl


text:
# Playing Cards: Playing Hand
import Base: +, -
import Random: shuffle!

struct Hand
    cards::Array{Card, 1}
    function Hand(; cards::Array{Card, 1} = Card[])
        return new(cards);
    end
end

@inline Base.length(hand::Hand) = length(hand.cards);
@inline Base.firstindex(hand::Hand) = 1;
@inline Base.lastindex(hand::Hand) = length(hand);
@inline Base.getindex(hand::Hand, i::Int64) = hand.cards[i];
@inline Base.getindex(hand::Hand, I::UnitRange{Int64}) = [hand[i] for i in I];
@inline Base.getindex(hand::Hand, I::Array{Int64, 1}) = [hand[i] for i in I];
@inline Base.iterate(hand::Hand, i::Int64 = 1) = (i > length(hand)) ? nothing : (hand[i], i + 1);

function Base.show(io::IO, hand::Hand)
    for card in hand
        print(io, card, " ")
    end
    return nothing;
end

@inline +(hand::Hand, card::Card) = push!(hand.cards, card);
@inline +(hand::Hand, cards::Array{Card, 1}) = append!(hand.cards, cards);
@inline -(hand::Hand, card::Card) = try deleteat!(hand.cards, findfirst(x -> x === card, hand.cards)) catch end;
function -(hand::Hand, cards::Array{Card, 1})
    for card in cards
        hand - card
    end
    return hand;
end

function shuffle(hand::Hand)
    shuffle!(hand.cards)
    return hand;
end

max_stars_repo_name: Joshua-Robison/PlayingCards

max_issues_repo_name: Joshua-Robison/PlayingCards

max_forks_repo_name: Joshua-Robison/PlayingCards


text:
# first example 2D

using Test
using LinearAlgebra
using ApproxManifoldProducts

using TransformUtils
const TU = TransformUtils

# using Gadfly

#
# logmap_SO2(Rl::Matrix{Float64}) = sign(Rl[2,1])*acos(Rl[1,1])
# difftheta(wth1::Float64, wth2::Float64)::Float64 = logmap_SO2(TU.R(wth1)'*TU.R(wth2))

#assume an unwrapped circle

@testset "Test SO2 wrap around minimization Gaussian means..." begin

@testset "opposing binary..." begin

# expect \mu = 5pi/6
μ1 = 0.0
μ2 = -pi+0.0

Λ1 = 1.0
Λ2 = 1.0

Lambdas = [Λ1; Λ2]
mus = [μ1; μ2]

μ = get2DMuMin(mus, Lambdas, diffop=difftheta, initrange=(-pi+0.0,pi+0.0), periodicmanifold=wrapRad)[1]
# TU.wrapRad(μ)

N = 1000

BUPS = Float64[get2DMuMin(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad)[1] for i in 1:N]
# BUPS[:] .= TU.wrapRad.(BUPS)

@test 0.3*N < sum( 0 .< BUPS .< 3pi/2.0)
@test 0.3*N < sum( -3pi/2 .< BUPS .< 0.0)
@test sum( -0.1 .< BUPS .< 0.1) < 0.01*N

# plot(x=BUPS, Geom.histogram)

end

@testset "binary near wrap..." begin

# expect \mu = 5pi/6
μ1 = -2.6
μ2 = 2.6

Λ1 = 10.0
Λ2 = 10.0

Lambdas = [Λ1; Λ2]
mus = [μ1; μ2]

μ = get2DMuMin(mus, Lambdas, diffop=difftheta, initrange=(-pi+0.0,pi+0.0), periodicmanifold=wrapRad)[1]
TU.wrapRad(μ)

# doesn't seem to work
# μ = get2DMu(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad)

N = 1000

BUPS = Float64[get2DMuMin(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad, initrange=(-pi+0.0,pi+0.0))[1] for i in 1:N]
# BUPS[:] .= TU.wrapRad.(BUPS)

@test 0.3*N < sum( 0 .< BUPS .< 3pi/2.0)
@test 0.3*N < sum( -3pi/2 .< BUPS .< 0.0)
@test 0.01*N < sum( -0.1 .< BUPS .< 0.1) < 0.3*N

# plot(x=BUPS, Geom.histogram)

end

@testset "fanning ternary means..." begin

μ1 = 0.0
μ2 = -4pi/6.0
μ3 = 4pi/6.0

Λ1 = 1.0
Λ2 = 1.0
Λ3 = 1.0

Lambdas = [Λ1; Λ2; Λ3]
mus = [μ1; μ2; μ3]

μ = get2DMuMin(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad, initrange=(-pi+0.0,pi+0.0))[1]

# TU.wrapRad(μ)

N = 100
BUPS = Float64[get2DMuMin(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad, initrange=(-pi+0.0,pi+0.0))[1] for i in 1:N]
# BUPS[:] .= TU.wrapRad.(BUPS)

# plot(x=BUPS, Geom.histogram)

@test 0.20*N < sum( pi/2 .< BUPS .< 3pi/2.0)
@test 0.20*N < sum( -3pi/2 .< BUPS .< -pi/2)
@test 0.20*N < sum( -0.1 .< BUPS .< 0.1)

end

end

@testset "basic linear 2D Rotation find mu and Lambdas Gaussian product test..." begin

# expect \mu = 5pi/6
μ1 = 4pi/6.0
μ2 = -pi+0.0
μ3 = -1.05
# μ3 = -0.5
#
# μ1 = 0.0
# μ2 = 2.0*pi/3.0
# μ3 = -2.0*pi/3.0
#

Λ1 = 1.0
Λ2 = 1.0
Λ3 = 1.0

mus = [μ1; μ2; μ3]
Lambdas = [Λ1; Λ2; Λ3]

difftheta(μ1, 0.0)

difftheta(μ1, -pi+0.0)
difftheta(μ1, pi+0.0)

difftheta(μ1, -1.04)
difftheta(μ1, -1.05)

@test abs(difftheta(μ2, 0.0)) - pi < 1e-8
@test abs(difftheta(μ2, 0.0001) ) - pi < 1e-3
@test abs(difftheta(μ2, -0.0001)) - pi < 1e-3

difftheta(μ2, -1.04)
difftheta(μ2, -1.05)

@test abs(difftheta(μ2, -pi+0.0)) < 1e-8
@test abs(difftheta(μ2, pi+0.0) ) < 1e-8

@test 1e-5 < abs( difftheta(μ1, 5pi/6.0) )
@test 1e-5 < abs( difftheta(μ2, 5pi/6.0) )
@test abs(difftheta(μ1, 5pi/6.0) + difftheta(μ2, 5pi/6.0)) <1e-8

resid2DLinear(0.0, mus, Lambdas, diffop=difftheta)

# @test abs(resid2DLinear(5pi/6.0, mus, Lambdas, diffop=difftheta)) < 1e-10

resid2DLinear(μ1, mus, Lambdas, diffop=difftheta)

resid2DLinear(pi+0.0, mus, Lambdas, diffop=difftheta)
resid2DLinear(-pi, mus, Lambdas, diffop=difftheta)

μ = get2DMu(mus, Lambdas, diffop=difftheta, periodicmanifold=wrapRad)
Λ = get2DLambda(Lambdas)

# @test abs(μ[1] - 8/5) < 1e-10
#
# @test abs(Λ - 5/4) < 1e-10

end

## plot residual and solution for validation
# using Colors
# using Gadfly
#
# μ1 = 0.0
# μ2 = 2.0*pi/3.0
# μ3 = -2.0*pi/3.0
#
# Λ1 = 1.0
# Λ2 = 1.0
# Λ3 = 1.0
#
# Λ1 = 1.0
# Λ2 = 2.0
# Λ3 = 0.25
#
# mus = [μ1; μ2; μ3]
# Lambdas = [Λ1; Λ2; Λ3]
#
# μ1 = 0.0
# μ2 = -pi+0.0
#
# mus = [μ1; μ2]
# Lambdas = [Λ1; Λ2]
#
# th = range(-pi, stop=pi-1e-5, length=1000)
# obj = map((x) -> resid2DLinear(x, mus, Lambdas, diffop=difftheta), th)
#
# PL = []
#
# # draw resdidual objective for delta angles
# push!( PL, Gadfly.layer(x=th,y=obj, Geom.path, Theme(default_color=colorant"red"))[1])
#
# vals = zeros(1000)
# for i in 1:length(vals)
#   vals[i] = get2DMu(mus, Lambdas, initrange=(-pi+0.0, pi-1e-5), diffop=difftheta, periodicmanifold=wrapRad)[1]
# end
#
# push!(PL, layer(x=vals, Geom.histogram(density=true))[1])
#
# plot(PL...)
#
# 0

#
#
#
#
#
# # Two means 0 and -pi, should have two minima at +- pi/2
# # Three means at 0, +-2pi/3, should have three minima at +-pi/3 and -pi
#
#
#
# μ1 = 0.0
# μ2 = -pi+0.0
# μ2 = 1.0
# mus = [μ1; μ2]
#
# Λ1 = 1.0
# Λ2 = 1.0
# lambdas = [Λ1; Λ2]
#
#
#

# obj .= norm.(obj)
#
#
# X = zeros(length(obj))
# Y = zeros(length(obj))
# global i = 0
# for t in th
#   global i
#   i += 1
#   X[i] = cos(t)*obj[i]
#   Y[i] = sin(t)*obj[i]
# end
#
#
#
# plot(
# Gadfly.layer(x=X, y=Y, Geom.path),
# Gadfly.layer(x=X, y=Y, Geom.path, Theme(default_color=colorant"red")),
# Coord.Cartesian(xmin=-6.0,xmax=6.0,ymin=-6.0,ymax=6.0)
# )
#
# 0
#

#
# # get2DMu(mus, Lambdas, initrange=(2.5+0.0, 2.8))[1, periodicmanifold=wrapRad]
# (x) -> resid2DLinear(x, mus, Lambdas, diffop=difftheta)
# gg2 = (x) -> resid2DLinear(x, mus, Lambdas, diffop=difftheta)
#
#
# using NLsolve
#
# res = zeros(1)
# gg = (res, x) -> solveresid2DLinear(res, x, mus, Lambdas)
# initrange = (-pi+0.0, pi+0.0)
# initr = initrange[2]-initrange[1]
#
# r = NLsolve.nlsolve(gg, [initr*rand()+initrange[1]])
# r = NLsolve.nlsolve(gg, [+1.5])
# r.zero[1]
#
# vals0 = zeros(100)
# vals = zeros(100)
# vals2 = zeros(100)
#
# for i in 1:length(th)
#   res[1] = 0.0
#   # gg(res, th[i])
#   res[1] = resid2DLinear(th[i], mus, Lambdas, diffop=difftheta)
#   vals0[i] = res[1]
#   res[1] = 0.0
#   solveresid2DLinear(res, th[i], mus, Lambdas, diffop=difftheta)
#   vals[i] = res[1]
#   vals2[i] = gg2(th[i])
# end
#
# plot(x=th, y=vals0, Geom.path)
# plot(x=th, y=vals, Geom.path)
# plot(x=th, y=vals2, Geom.path)
#
#
# th = range(-pi, stop=pi, length=100)
#
#
# μ1 = 0.0
# mus = [μ1]
#
# Λ1 = 1.0
# Lambdas = [Λ1]
#
#
# obj = map((x) -> resid2DLinear(x, mus, Lambdas, diffop=difftheta), th)
#
# obj .= norm.(obj)
# # obj2= deepcopy(obj)
# # obj .+= obj2
#
#
# X = zeros(length(obj))
# Y = zeros(length(obj))
# global i = 0
# for t in th
#   global i
#   i += 1
#   X[i] = cos(t)*obj2[i]
#   Y[i] = sin(t)*obj2[i]
# end
#
# plot(x=X, y=Y, Geom.path)
#
# 0

# resid2DLinear(-1.05, mus, Lambdas, diffop=difftheta)

#
# using Gadfly
#
# x = range(-1.1pi, stop=1.1pi, length=1000)
#
# y = difftheta.(0.0,x)
#
# plot(x=x, y=y, Geom.line)
#
# for i in 1:length(x)
#   y[i] = resid2DLinear(x[i], mus, Lambdas, diffop=difftheta)
# end

#


text:
module GeophysicalModelGenerator

using Base: String, show_index, Tuple, FieldDescStorage
using Requires

# Load & export some useful commands/functions from GeoParams:
import GeoParams
using .GeoParams
export
        @u_str, uconvert, upreffered, unit, ustrip, NoUnits,  #  Units
        GeoUnit, GEO_units, SI_units, NO_units, AbstratGeoUnits,
        Nondimensionalize, Nondimensionalize!, Dimensionalize, Dimensionalize!,
        superscript, upreferred, GEO, SI, NONE, isDimensional,
        km, m, cm, mm, Myrs, yr, s, MPa, Pa, Pas, K, C, kg, mol

export ReadCSV_LatLon, meshgrid, voxGrav

# julia standard library packages
using DelimitedFiles, Statistics

# other packages
using   WriteVTK, Colors, MeshIO, FileIO, Interpolations, Geodesy

export vtk_multiblock, vtk_save         # Simplifies writing multiblock files
export LLA
export load

# add files for specific tasks
include("data_types.jl")
include("data_import.jl")
include("coord_conversion.jl")
include("utils.jl")
include("Paraview_output.jl")
include("transformation.jl")
include("voxel_gravity.jl")
include("LaMEM_io.jl")
include("LaMEM_geometry.jl")
include("stl.jl")

# Add plotting routines
function __init__()
        @require GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" begin
                @eval include("./GMT_utils.jl")
        end
end

end


text:
@testset "105.Construct-Binary-Tree-from-Preorder-and-Inorder-Traversal.jl" begin
    @test build_tree_105([3, 9, 20, 15, 7], [9, 3, 15, 20, 7]) ==
        TreeNode{Int}([3, 9, 20, nothing, nothing, 15, 7])
    @test build_tree_105([-1], [-1]) == TreeNode(-1)
    @test build_tree_105([3, 9, 8, 5, 4, 10, 20, 15, 7], [4, 5, 8, 10, 9, 3, 15, 20, 7]) ==
        TreeNode{Int}([
        3, 9, 20, 8, nothing, 15, 7, 5, 10, nothing, nothing, nothing, nothing, 4
    ])
end


text:
"""
    @check cond msg

If `cond` evaluates false throw `ArgumentError` and print evaluation of `cond`.
"""
macro check(cond, msg)
    str = _chkmsg(cond, msg)
    quote
        if !$(esc(cond))
            throw(ArgumentError($(str)))
        end
    end
end

"""
    @checkwarn [showwarn] cond msg

If `cond` evaluates false warn and print evaluation of `cond`.
If optional argument `showwarn=false` hide warning.
"""
macro checkwarn(cond, msg)
    str = _chkmsg(cond, msg)
    quote
        if !$(esc(cond))
            @warn $(str)
        end
    end
end

macro checkwarn(showarn, cond, msg)
    str = _chkmsg(cond, msg)
    quote
        if $(esc(showarn)) && !$(esc(cond))
            @warn $(str)
        end
    end
end

"""
    _chkmsg(cond, msg)

Create an expression which evaluates as a error string which contrains the original msg and
some debug information based on the condition.
"""
function _chkmsg(cond, msg)
    head = lstrip(repr(cond), ':')
    if head[begin] == '(' && head[end] == ')'
        head = head[begin+1:end-1]
    end
    head = head * " evaluated false"

    args = Expr[]

    if cond.args[1] == :(!)
        _expr_repr_list!(args, cond.args[2].args[2:end])
    elseif cond.args[1] == :(⊆)
        push!(args, :("\n   ├ set diff = " * _shortrepr(setdiff($(esc(cond.args[2])), $(esc(cond.args[3]))))))
        _expr_repr_list!(args, cond.args[2:end])
    elseif cond.args[1] == :(⊇)
        push!(args, :("\n   ├ set diff = " * _shortrepr(setdiff($(esc(cond.args[3])), $(esc(cond.args[2]))))))
        _expr_repr_list!(args, cond.args[2:end])
    elseif cond.args[1] == :(uniquenames) || cond.args[1] == :(allunique)
        push!(args, :("\n   ├ duplicates = " * _shortrepr(duplicates($(esc(cond.args[2]))))))
        _expr_repr_list!(args, cond.args[2:end])
    else
        _expr_repr_list!(args, cond.args[2:end])
    end
    return :($(esc(msg)) * "\n  " * $head * $(args...))
end

function _expr_repr_list!(list, expressions)
    for (i,a) in enumerate(expressions)
        lhs = lstrip(repr(a), ':') # remove leading : on symbols
        if lhs[begin] == '(' && lhs[end] == ')'
            lhs = lhs[begin+1:end-1]
        end
        symbol = (i == length(expressions)) ? "└ " : "├ "
        push!(list, :("\n   " * $symbol * $lhs * " = " * _shortrepr($(esc(a)))))
    end
end

_shortrepr(x) = repr(x)
function _shortrepr(@nospecialize x::Set)
    isempty(x) && return "(empty)"
    m = match(r"\[(.*)\]\)$", repr(x)::String)
    return m[1]
end
function _shortrepr(@nospecialize x::AbstractArray)
    isempty(x) && return "(empty)"
    m = match(r"\[(.*)\]$", repr(x)::String)
    return m[1]
end

"""
    remove_namespace(namespace, x)

If first namespace of `x` is `namespace` then remove it.
```
remove_namespace(A, A₊B₊x) -> B₊x
remove_namespace(B, A₊B₊x) -> A₊B₊x
````
"""
remove_namespace(namespace, name::T) where T = T(replace(String(name), Regex("^$(namespace)₊") => ""))
remove_namespace(namespace, x::Sym) = rename(x, remove_namespace(namespace, x.name))
remove_namespace(namespace, x::Term) = rename(x, remove_namespace(namespace, operation(x).name))

"""
    remove_namespace(x)

Strips `x` of its first namespace. `A₊B₊x -> B₊x`
"""
remove_namespace(name::T) where T = T(replace(String(name), Regex("^(.+?)₊") => ""))
remove_namespace(x::Sym) = rename(x, remove_namespace(x.name))
remove_namespace(x::Term) = rename(x, remove_namespace(operation(x).name))

eqsubstitute(eq::Equation, rules) = substitute(eq.lhs, rules) ~ substitute(eq.rhs, rules)

uniquenames(syms) = allunique(getname.(syms))

function duplicates(X)
    X = collect(X)
    dups = Set{eltype(X)}()
    for (i, x) in enumerate(X)
        for j in i+1:lastindex(X)
            if isequal(x, X[j])
                push!(dups, x)
                break
            end
        end
    end
    return dups
end

"""
    eq_type(eq::Equation)

Checks the type of the equation. Returns:
- `(:diffeq, lhs_variable)` for differential equations
- `(:explicit_algebraic, lhs_variable)` for explicit algebraic equations
- `(:implicit_algebraic, lhs_variable)` for implicit algebraic equations
"""
function eq_type(eq::Equation)
    if eq.lhs isa Term && operation(eq.lhs) isa Differential
        vars = get_variables(eq.lhs)
        @check length(vars) == 1 "Diff. eq $eq has more than on variable in lhs!"
        return (:diffeq, vars[1])
    elseif eq.lhs isa Sym || eq.lhs isa Term
        vars = get_variables(eq.lhs)
        @check length(vars) == 1 "Algebraic eq $eq has more than on variable in lhs!"
        if vars[1] ∈ Set(get_variables(eq.rhs))
            return (:implicit_algebraic, vars[1])
        else
            return (:explicit_algebraic, vars[1])
        end
    elseif isequal(eq.lhs, 0)
        return (:implicit_algebraic, nothing)
    else
        throw(ArgumentError("Uknnown equation type $eq"))
    end
end

"""
    lhs_var(eq::Equation)

Returns the variable on the lhs of the equation for equations.
"""
lhs_var(eq::Equation) = eq_type(eq)[2]

"""
    recusive_substitute(term, rules::Dict)

Apply substitutions until there is no more change in `term`.
"""
function recursive_substitute(term, rules::Dict)
    new_term = substitute(term, rules)
    while !isequal(new_term, term)
        term = new_term
        new_term = substitute(term, rules)
    end
    return new_term
end

"""
    rhs_differentials(iob::IOBlock)

Return Set of all differentials which are present in the rhs of the system.
"""
function rhs_differentials(iob)
    diffs = Set{SymbolicUtils.Symbolic}()
    for eq in equations(iob)
        _collect_differentials!(diffs, eq.rhs)
    end
    return diffs
end

_collect_differentials(ex) = _collect_differentials!(Set{SymbolicUtils.Symbolic}(), ex)

function _collect_differentials!(found, ex)
    if SymbolicUtils.istree(ex)
        if operation(ex) isa Differential
            push!(found, ex)
        else
            for arg in arguments(ex)
                _collect_differentials!(found, arg)
            end
        end
    end
    return found
end


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

"""
    DiscretizationMethod

A method for discretizing spatial regions.
"""
abstract type DiscretizationMethod end

"""
    discretize(region, method)

Discretize spatial `region` with discretization `method`.
"""
function discretize end

# ----------------
# IMPLEMENTATIONS
# ----------------
include("discretizing/block.jl")


text:
@testset "more $N" for (N, FR, Ti) = (("FastQ32", FastQ32, Int32), ("FastQ64", FastQ64, Int64))

    @testset "constructors $N" begin
        @test_throws DivideError FR(1, 0)
        @test_throws DomainError FR(Ti(1), typemin(Ti))
        @test_throws DomainError FR(1, typemin(Ti))
        mx = FR(typemin(Ti))
        @test_throws DivideError inv(mx)
        @test FR(typemin(Ti) + 10000) - FR(10000) == mx
        my = FR(typemin(Ti) + 10000) - FR(10000)
        @test_throws DivideError Ti(1) / my
    end

    @testset "avoid overflow during add/subtract $N" begin
        x =  div(typemax(Ti),6)
        @test FR(17, 2x) + FR(17, 3x) == FR(85, 6x)
        @test FR(17, 2x) - FR(17, 3x) == FR(17, 6x)
    end

    @testset "constructors1 $N" begin
        @test FastRational{Ti}(Ti(2),Ti(4),Val(true)).num == 2
        @test FastRational{Ti}(Ti(2),Ti(4),Val(true)).den == 4
        @test FastRational{Ti}(2,4).num == 1
        @test FastRational{Ti}(2,4).den == 2
        @test FR((2,4)).num == 1
        @test FR((2,4)).den == 2
    end

    @testset "convert from Rational $N" begin
        @test FR(1//2) == FR(1, 2)
        @test FR((2, 4)) == FR(1, 2)
    end

    @testset "inverse matrix $N" begin
        B = FR.([(1,2) (1,3); (1,3) (1,4)])
        @test inv(B) == FR.([(18,1) (-24,1); (-24,1) (36,1)])
    end

end


text:
#####
##### Boundary conditions on tendency terms are derived from the boundary
##### conditions on their repsective fields.
#####

TendencyBC(::BC)   = BoundaryCondition(Flux, nothing)
TendencyBC(::PBC)  = PeriodicBC()
TendencyBC(::NPBC) = NoPenetrationBC()

TendencyCoordinateBCs(bcs) =
    CoordinateBoundaryConditions(TendencyBC(bcs.left), TendencyBC(bcs.right))

TendencyFieldBCs(field_bcs) =
    FieldBoundaryConditions(Tuple(TendencyCoordinateBCs(bcs) for bcs in field_bcs))

TendenciesBoundaryConditions(solution_bcs) =
    NamedTuple{propertynames(solution_bcs)}(Tuple(TendencyFieldBCs(bcs) for bcs in solution_bcs))


text:
using DeutschlandDaten
using Test

@testset "DeutschlandDaten.jl" begin

    @test [deaths(deathdat; alter=a, jahr=2020, monat=1, geschlecht="Ins")
           for a in DeutschlandDaten.Altersgruppen_Monate] ==
               [ 272, 329, 248, 374, 486, 856, 1943, 3159, 4414, 5820, 6653, 11018,  16556, 15577+12282+ 4993]

    @test [deaths(deathdat; alter=a, jahr=2018, monat=6, geschlecht="Ins")
           for a in DeutschlandDaten.Altersgruppen_Monate] ==
        [ 322, 345, 210, 317, 418, 873, 1845, 2748, 3863, 4842, 5555, 10125, 12784, 12196+ 9168+ 3717 ]

    @test [deaths(deathdat; alter=a, jahr=2019, monat=3, geschlecht="Männlich")
           for a in DeutschlandDaten.Altersgruppen_Monate] ==
               [ 181, 228, 149, 217, 338, 585, 1295, 2045, 2877, 3640, 3991, 6926, 8337, 6668 + 3916 + 1037 ]

    @test [ deaths(deathdat; alter=a, jahr=2020, kw=3, geschlecht="Ins")
           for a in DeutschlandDaten.Altersgruppen_KW ] ==
               [ 130, 62, 79, 102, 190, 449, 723, 1003, 1259, 1459, 2424, 3701, 3526 + 2746 + 1100 ]

    @test [ deaths(deathdat; alter=a, jahr=2018, kw=6, geschlecht="Männlich")
           for a in DeutschlandDaten.Altersgruppen_KW ] ==
               [ 98, 44, 43, 77, 122, 273, 463, 617, 793, 992, 1738, 1874, 1511 + 862 + 259 ]

    @test [ deaths(deathdat; alter=a, jahr=2017, bundesland="Schleswig-Holstein", kw=1, geschlecht="Männlich")
            for a in DeutschlandDaten.Altersgruppen_Bundesländer ] ==
                [  61, 75, 101, 80  ]

    @test [ deaths(deathdat; alter=a, jahr=2018, bundesland="Bayern", kw=6, geschlecht="Männlich")
           for a in DeutschlandDaten.Altersgruppen_Bundesländer ] ==
               [ 250, 256, 487, 406  ]

    @test [ deaths(deathdat; alter=a, jahr=2021, bundesland="Niedersachsen", kw=6, geschlecht="Weiblich")
           for a in DeutschlandDaten.Altersgruppen_Bundesländer ] ==
               [ 97, 124, 324, 503 ]

    @test [ deaths(deathdat; alter=a, jahr=2016, bundesland="Baden-Württemberg", kw=52, geschlecht="Weiblich")
           for a in DeutschlandDaten.Altersgruppen_Bundesländer ] ==
               [ 94, 124, 350, 595 ]

    # Write your tests here.
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
"""
Hidden Markov Models for Julia.

[Documentation](https://maxmouchet.github.io/HMMBase.jl/stable/).
[Issues](https://github.com/maxmouchet/HMMBase.jl/issues).
"""
module HMMBase

using ArgCheck
using Clustering
using Distributions
using Hungarian
using LinearAlgebra

using Base: OneTo
using Random: AbstractRNG, GLOBAL_RNG

# Extended functions
import Base: ==, copy, rand, size
import Distributions: MixtureModel, fit_mle, loglikelihood

export
    # hmm.jl
    AbstractHMM,
    HMM,
    TimeVaryingHMM,
    fit_mle!,
    copy,
    rand,
    size,
    nparams,
    permute,
    statdists,
    istransmat,
    likelihoods,
    # messages.jl
    forward,
    backward,
    posteriors,
    loglikelihood,
    # likelihoods.jl
    loglikelihoods,
    # mle.jl
    fit_mle,
    # viterbi.jl
    viterbi,
    # utilities.jl,
    gettransmat,
    randtransmat,
    remapseq

include("hmm.jl")
include("mle.jl")
include("mle_init.jl")
include("messages.jl")
include("viterbi.jl")
include("likelihoods.jl")
include("utilities.jl")
include("experimental.jl")

include("timevaryinghmm.jl")

# To be removed in a future version
# ---------------------------------
#!format: off

# < v1.1

export likelihoods

function likelihoods(args...; logl = false, kwargs...)
    @warn "`likelihoods(...)` is deprecated, use `loglikelihoods(...)` or `exp.(loglikelihoods(...))` instead."
    logl ? loglikelihoods(args...; kwargs...) : exp.(loglikelihoods(args...; kwargs...))
end

function deprecate_kwargs(name)
    @warn "`$(name)` keyword argument is deprecated."
end

# < v1.0

@deprecate n_parameters(hmm) nparams(hmm)
@deprecate log_likelihoods(hmm, observations) loglikelihoods(hmm, observations)

@deprecate forward_backward(init_distn, trans_matrix, log_likelihoods) posteriors(init_distn, trans_matrix, log_likelihoods, logl = true)
@deprecate messages_forwards(init_distn, trans_matrix, log_likelihoods) forward(init_distn, trans_matrix, log_likelihoods, logl = true)
@deprecate messages_backwards(init_distn, trans_matrix, log_likelihoods) backward(init_distn, trans_matrix, log_likelihoods, logl = true)

@deprecate forward_backward(hmm, observations) posteriors(hmm, observations, logl = true)
@deprecate messages_forwards(hmm, observations) forward(hmm, observations, logl = true)
@deprecate messages_backwards(hmm, observations) backward(hmm, observations, logl = true)

@deprecate messages_forwards_log(init_distn, trans_matrix, log_likelihoods) log.(forward(init_distn, trans_matrix, log_likelihoods, logl = true)[1])
@deprecate messages_backwards_log(trans_matrix, log_likelihoods) log.(backward(init_distn, trans_matrix, log_likelihoods, logl = true)[1])

@deprecate compute_transition_matrix(seq) gettransmat(seq, relabel = true)
@deprecate rand_transition_matrix(K, α = 1.0) randtransmat(K, α)

end


text:
module IR

using ..MLIR

import Base: ==, insert!, push!, getindex

include("utils.jl")

include("ir/context.jl")
include("ir/type.jl")
include("ir/value.jl")
include("ir/attribute.jl")
include("ir/location.jl")
include("ir/operation.jl")
include("ir/identifier.jl")
include("ir/block.jl")
include("ir/region.jl")
include("ir/dialect.jl")
include("ir/module.jl")
include("ir/builders.jl")

include("pass_manager.jl")

end # module


text:
using BinaryBuilder

name = "GLFW"
version = v"3.3.6"

# Collection of sources required to build glfw
sources = [
    ArchiveSource("https://github.com/glfw/glfw/releases/download/$(version)/glfw-$(version).zip",
                  "45537305d44c0a9f3612d4ec4a48414547cf854bff3ed613078f7ec648a12781")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glfw-*/
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DBUILD_SHARED_LIBS=ON \
    -DCMAKE_BUILD_TYPE=Release \
    -DGLFW_BUILD_EXAMPLES=false \
    -DGLFW_BUILD_TESTS=false \
    -DGLFW_BUILD_DOCS=OFF
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libglfw", "glfw3"], :libglfw)
]

x11_platforms = filter(p ->Sys.islinux(p) || Sys.isfreebsd(p), platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"; platforms=x11_platforms),
    Dependency("Libglvnd_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXcursor_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXi_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXinerama_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXrandr_jll"; platforms=x11_platforms),
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
for i = 1 : (37 + 5)
	j = 2 * idx
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
push!(LOAD_PATH, joinpath(expanduser("~/repos/InsarTimeseries.jl/src/")))
import InsarTimeseries
using InsarTimeseries:
    PHASE_TO_CM,
    STACK_FLAT_SHIFTED_DSET,
    STACK_FLAT_DSET,
    STACK_DSET,
    CC_FILENAME,
    UNW_FILENAME
using Dates
using HDF5
using LinearAlgebra
using Statistics: mean
# using Plots
# pyplot()
using PyCall: pyimport
import Polynomials
import NaNMath
using Printf: @printf
import PyPlot
plt = PyPlot

nm = NaNMath
# import Convex
# import ECOS
# import SCS

p2c = InsarTimeseries.PHASE_TO_CM
p2mm = InsarTimeseries.PHASE_TO_CM * 365 * 10

# TODO: find the shift to match insar to gps best
# insar_arr_txmc = [365 * 10 * PHASE_TO_CM * solve_insar_ts(station_name, 5, "TXMC")[3][1] for station_name in  station_name_list]
# [ rms(gps_arr78 - (insar_arr_txmc .+ cc)) for cc in range(1, stop=3, length=10)]
#
YLIMS = (-12, 12)
# TODO: maybe get these names automatically?
# PATH 78
station_name_list = [
    "NMHB",
    "TXAD",
    "TXBG",
    "TXBL",
    "TXCE",
    "TXFS",
    "TXKM",
    "TXL2",
    "TXMC",
    "TXMH",
    "TXOE",
    "TXOZ",
    "TXS3",
    "TXSO",
]
#
# PATH 85
# station_name_list = ["TXKM", "TXMH", "TXFS", "TXAL"]
# BAD: MDO1 (nothing 2014-2017), TXVH (started 2018), TXPC (ended 2017)

# REFERENCE_STATION = "TXKM"
# REFERENCE_STATION = "TXAD"
REFERENCE_STATION = nothing

latlon = pyimport("apertools.latlon")
sario = pyimport("apertools.sario")
gps = pyimport("apertools.gps")

# const DEFO_FILENAME = "deformation_unreg_maxtemp400.h5"
# const DEFO_FILENAME_LINEAR = "deformation_linear_maxtemp400.h5"
# const DEFO_FILENAME_L1 = "deformation_linear_maxtemp400_l1.h5"
# println("Defo filename $DEFO_FILENAME, defo linear: $DEFO_FILENAME_LINEAR")

const ignore_geo_file = "geolist_ignore.txt"
# max_temporal_baseline = 500
max_temporal_baseline = 800
# max_temporal_baseline = 1200
unw_filename = "unw_stack.h5"

# NOTE: not using the ignore file for now- using all .geos
GEOLIST, INTLIST, VALID_IGRAM_INDICES =
    InsarTimeseries.load_geolist_intlist(unw_filename, "", max_temporal_baseline);
timediffs = InsarTimeseries.day_diffs(GEOLIST)

# For checking / plotting specific interestion points:
# Sinkhole
row, col = 1838, 1367
# row, col = [361, 257] .+ [1, 1]
# well uplift
row, col = 1223, 943
# row, col = [244, 188] .+ [1, 1]
# Pecos subs.
# row, col = [368, 131]

function get_stack_vals(
    unw_stack_file::String,
    station_name::String,
    window = 5,
    dset = STACK_FLAT_DSET,
    valid_indices = VALID_IGRAM_INDICES;
    reference_station = nothing,
)
    dem_rsc = sario.load("dem.rsc")
    lon, lat = gps.station_lonlat(station_name)
    row, col =
        map(x -> convert(Int, x), latlon.nearest_pixel(dem_rsc, lon = lon, lat = lat))

    # println("Getting data at $station_name: row $row, col $col, lon $lon, lat $lat")
    unw_vals = _read_vals(unw_stack_file, row, col, window, dset, valid_indices)
    return _subtract_reference(
        unw_vals,
        reference_station,
        unw_stack_file,
        window,
        dset,
        valid_indices,
    )
end

function get_stack_vals(
    unw_stack_file::String,
    row::Int,
    col::Int,
    window = 5,
    dset = STACK_FLAT_DSET,
    valid_indices = VALID_IGRAM_INDICES;
    reference_station = nothing,
)
    unw_vals = _read_vals(unw_stack_file, row, col, window, dset, valid_indices)
    return _subtract_reference(
        unw_vals,
        reference_station,
        unw_stack_file,
        window,
        dset,
        valid_indices,
    )
end

function _read_vals(
    unw_stack_file,
    row,
    col,
    window = 5,
    dset = STACK_FLAT_DSET,
    valid_indices = VALID_IGRAM_INDICES,
)
    # println("Loading $row, $col from $dset, avging window $window")
    halfwin = div(window, 2)
    # Note: swapping the col and row in h5read since julia is col-major
    unw_depth =
        h5read(unw_stack_file, dset, (col-halfwin:col+halfwin, row-halfwin:row+halfwin, :))
    unw_vals_all = vec(mean(unw_depth, dims = (1, 2)))
    unw_vals = unw_vals_all[valid_indices]
end

function _subtract_reference(
    unw_vals,
    reference_station,
    unw_stack_file,
    window,
    dset,
    valid_indices,
)
    isnothing(reference_station) && return unw_vals
    return unw_vals -
           get_stack_vals(unw_stack_file, reference_station, window, dset, valid_indices)
end

function solve_insar_ts(
    station_name::String,
    window::Int = 5,
    reference_station = nothing;
    cutoff = false,
)
    unw_vals = get_stack_vals(
        UNW_FILENAME,
        station_name,
        window,
        reference_station = reference_station,
    )
    return solve_insar_ts(unw_vals, window, cutoff = cutoff)
end

function solve_insar_ts(
    row::Int,
    col::Int,
    window::Int = 5,
    reference_station = nothing;
    cutoff = false,
)
    unw_vals = get_stack_vals(
        UNW_FILENAME,
        row,
        col,
        window,
        reference_station = reference_station,
    )
    return solve_insar_ts(unw_vals, window, cutoff = cutoff)
end

function solve_insar_ts(unw_vals::Array{<:AbstractFloat,1}, window::Int = 5; cutoff = false)

    B = InsarTimeseries.build_B_matrix(GEOLIST, INTLIST)
    Blin = sum(B, dims = 2)

    if cutoff
        v_linear_l1, v_linear_lstsq = InsarTimeseries.solve_after_cutoff(
            GEOLIST,
            INTLIST,
            unw_vals,
            Blin,
            in_mm_yr = false,
        )
        return v_linear_lstsq, 0, v_linear_l1, 0
    end

    v_linear_lstsq = Blin \ unw_vals
    v_unreg_lstsq = B \ unw_vals

    # # solver = SCS.SCSSolver()
    # solver = ECOS.ECOSSolver(verbose=0)
    # # Using functions in module:
    # var_linear = Convex.Variable(1)
    # v_linear_l1 = InsarTimeseries.invert_pixel(unw_vals, Blin, var_linear)
    # var_unreg = Convex.Variable(size(B, 2))
    # v_unreg_l1 = InsarTimeseries.invert_pixel(unw_vals, B, var_unreg)

    # Using Huber Loss in ADMM
    v_linear_l1 = InsarTimeseries.invert_pixel(unw_vals, Blin)
    v_unreg_l1 = InsarTimeseries.invert_pixel(unw_vals, Blin)

    return v_linear_lstsq, v_unreg_lstsq, v_linear_l1, v_unreg_l1
end

function integrate_velos(v_linear_lstsq, v_unreg_lstsq, v_linear_l1, v_unreg_l1)
    linear_lstsq =
        PHASE_TO_CM .* InsarTimeseries.integrate_velocities(v_linear_lstsq, timediffs)
    unreg_lstsq =
        PHASE_TO_CM .* InsarTimeseries.integrate_velocities(v_unreg_lstsq, timediffs)
    linear_l1 = PHASE_TO_CM .* InsarTimeseries.integrate_velocities(v_linear_l1, timediffs)
    unreg_l1 = PHASE_TO_CM .* InsarTimeseries.integrate_velocities(v_unreg_l1, timediffs)
    return linear_lstsq, unreg_lstsq, linear_l1, unreg_l1
end

function plot_insar(geolist, insar_linear_ts, insar_unreg_ts; title = "", ylims = YLIMS)
    # TODO: fix plotting to pyplot
    plt.figure()
    p = plt.plot(
        geolist,
        insar_linear_ts,
        label = "linear insar",
        linewidth = 3,
        ylims = ylims,
    )
    plt.ylabel("cm")
    # plt.plot(geolist, insar_unreg_ts, label="unreg insar", marker=:o, linealpha=0.0)
    return p
end

function plot_gps(dts, gps_los_data, gps_poly)
    day_nums = _get_day_nums(dts)
    # TODO: fix plotting to pyplot
    plt.plot(dts, gps_los_data, "gx", markersize = 1, label = "gps")
    # plt.plot(dts, gps_poly(day_nums),  color=:green, linewidth=3, label="gps line fit")
end

_get_day_nums(dts) = [(d - dts[1]).value for d in dts]

function process_pixel(;
    station_name = nothing,
    plotting = false,
    reference_station = REFERENCE_STATION,
    window = 5,
    verbose = true,
    cutoff = false,
)
    total_interval_days = (GEOLIST[end] - GEOLIST[1]).value

    # First solve for the velocities in L1 vs L2 to compare to GPS
    reference_station = "TXKM"
    v_linear_lstsq, v_unreg_lstsq, v_linear_l1, v_unreg_l1 =
        solve_insar_ts(station_name, window, reference_station, cutoff = cutoff)

    # NOTE: CURRENLT IGNORING THE REFERENCE STATION AND FORCING IT TO BE NOTHING
    slope_gps_mm_yr = InsarTimeseries.solve_gps_ts(station_name, nothing)

    slope_insar_l2_mm_yr = (365*10*PHASE_TO_CM*v_linear_lstsq)[1] # solution currently in phase
    slope_insar_l1_mm_yr = (365*10*PHASE_TO_CM*v_linear_l1)[1]

    l2_error = slope_gps_mm_yr - slope_insar_l2_mm_yr
    l1_error = slope_gps_mm_yr - slope_insar_l1_mm_yr

    if verbose
        @show station_name reference_station
        println("GPS velocity (mm / year): $slope_gps_mm_yr")
        println("InSAR L2 linear velocity (mm / year): $slope_insar_l2_mm_yr")
        println("InSAR L1 linear velocity (mm / year): $slope_insar_l1_mm_yr")
        println("==="^10)
        println("Difference (mm / yr): GPS - Insar L2 solution: $l2_error")
        println("Difference: (mm / yr) GPS - Insar L1 solution: $l1_error")
        println("==="^10)
    end

    if plotting
        # Now get the full time series of insar by integrating
        linear_lstsq, unreg_lstsq, linear_l1, unreg_l1 =
            integrate_velos(v_linear_lstsq, v_unreg_lstsq, v_linear_l1, v_unreg_l1)
        # Plot the solutions vs gps
        plot_insar(GEOLIST, linear_lstsq, unreg_lstsq, title = "L2 least squares solution")
        plot_gps(dts, gps_los_data, gps_poly)

        plot_insar(GEOLIST, linear_l1, unreg_l1, title = "L1 norm minimization")
        plot_gps(dts, gps_los_data, gps_poly)

        plot(p1, p2)
        println("Saving $station_name.png")
        png(station_name)
    end
    return l1_error, l2_error
end

function calc_error_matrix(station_name_list, plotting = false)
    println("Calculating errors at all station/reference combos")
    l1_error_matrix = zeros(length(station_name_list), length(station_name_list))
    l2_error_matrix = zeros(length(station_name_list), length(station_name_list))

    for (i, ref_name) in enumerate(station_name_list)
        println("Calculating errors for $ref_name")
        for (j, station_name) in enumerate(station_name_list)
            # println("Processing $station_name")
            l1_error, l2_error = process_pixel(
                station_name = station_name,
                plotting = plotting,
                reference_station = ref_name,
                verbose = false,
            )
            l1_error_matrix[i, j] = l1_error
            l2_error_matrix[i, j] = l2_error
        end
    end
    return l1_error_matrix, l2_error_matrix
end

# Reminder: error matrix rows correspond to 1 single reference station
function print_matrix_stats(l1_error_matrix, l2_error_matrix, station_name_list)
    println("Errors compared to GPS when using <station name> as the reference point")
    for (mat, lname) in zip((l2_error_matrix, l1_error_matrix), ("L2", "L1"))
        println("$lname RESULTS:")
        println("Name | RMS error | Total abs | Max abs")
        for i = 1:length(station_name_list)
            row = mat[i, :]
            a, b, c, d =
                station_name_list[i], rms(row), total_abs_error(row), maximum(abs.(row))
            @printf("%s :    %.2f     |   %.2f   |    %.2f \n", a, b, c, d)
        end
    end
end

# Functions for error evaluation
rms(arr::AbstractArray{Any,1}) = sqrt(mean(arr .^ 2))
rms(arr::AbstractArray{Any,2}) = [rms(arr[i, :]) for i = 1:size(arr, 1)]
total_abs_error(arr::AbstractArray{Any,1}) = sum(abs.(arr))
total_abs_error(arr::AbstractArray{Any,2}) = sum(abs.(arr), dims = 2)

plotting = false

############
# Try and remove a day, check all GPS errors
############

num_geos = length(GEOLIST)
num_stations = length(station_name_list)

l1_diff_matrix = zeros(num_stations, num_geos)
base_errors = zeros(num_stations)

for (idx, station_name) in enumerate(station_name_list)
    window = 5
    # I'm assuming TXKM is the best based on other analysis here
    ref_stat = "TXKM"
    # unw_vals = get_stack_vals(UNW_FILENAME, station_name, window, reference_station=ref_stat)

    # The "diff" being positive means an improvement (new error off GPS is lower), while
    # negative means the new err from GPS is bigger than before
    #
    # l1_diffs, l2_diffs, base_l1_error = InsarTimeseries.compare_solutions_with_gps(GEOLIST, INTLIST, unw_vals, station_name)

    # l1_diff_matrix[idx, :] = l1_diffs
    # base_errors[idx] = base_l1_error
end

# # top row is RMS, bottom row is max
# new_errs = zeros(2, num_geos)
# for idx in 1:num_geos
#     # if l1_diff_matrix entry is positive, it is an improvement,
#     # so we subtract from base.
#     new_err = abs.(abs.(base_errors) .- l1_diff_matrix[:, idx])
#     # reminder: l1_diffs[idx] = abs(base_l1_error) - abs(l1d)
#     # which means we are getting back to `l1d = l1 - slope_gps_mm_yr`
#
#     # println("$(rms(new_err)) , $(maximum(new_err))")
#     new_errs[:, idx] = [rms(new_err) , maximum(new_err)]
# end

# bases = [rms(abs.(base_errors)), maximum(abs.(base_errors))]
# err_diffs = new_errs .- bases;

function plot_errs(err_diffs, geolist)
    scatter(
        geolist,
        err_diffs[1, :],
        label = "RMS",
        title = "difference in errors (negative=improvement)",
    )
    scatter!(geolist, err_diffs[2, :], label = "maximum")
end

######

l1_errors, l2_errors = [], []

# Compute all combos:
# @time l1_error_matrix, l2_error_matrix = calc_error_matrix(station_name_list, plotting)
# print_matrix_stats(l1_error_matrix, l2_error_matrix, station_name_list)

# station_name = station_name_list[1]
# station_name = "TXOZ"
# station_name_list = ["TXMH"]

# l1_error, l2_error = process_pixel(station_name, UNW_FILENAME, GEOLIST, INTLIST, VALID_IGRAM_INDICES, plotting=plotting)
# append!(l2_errors, l2_error)
# append!(l1_errors, l1_error)

if false
    for station_name in station_name_list
        l1_error, l2_error =
            process_pixel(station_name = station_name, plotting = plotting, cutoff = true)
        append!(l2_errors, l2_error)
        append!(l1_errors, l1_error)
    end

    println("TOTAL ERRORS (all in mm / year of velocity):")
    println("L2 RMS error for all stations: $(rms(l2_errors))")
    println("L1 RMS error for all stations: $(rms(l1_errors))")
    println("L2 sum of abs errors for all stations: $(total_abs_error(l2_errors))")
    println("L1 sum of abs errors for all stations: $(total_abs_error(l1_errors))")
    println("L2 maximum errors : $(maximum(abs.(l2_errors)))")
    println("L1 maximum errors : $(maximum(abs.(l1_errors)))")

end

# function check_convex_l2(Blin, unw_vals)
#     # To check that l2 norm min matches the backslash result
#     v2linear_test = Convex.Variable(1)
#     prob_linear_test2 = minimize(norm(Blin*v2linear_test - unw_vals, 2))
#     Convex.solve!(prob_linear_test2, ECOSSolver())
# end

# # The values already solved for:
# ts_linear = h5read(DEFO_FILENAME_LINEAR, "stack", (col, row, :))[1, 1, :];
# ts_unreg = h5read(DEFO_FILENAME, "stack", (col, row, :))[1, 1, :];
# p1 = plot(GEOLIST, ts_linear)
# plot!(p1, GEOLIST, ts_unreg)

function plotunws(B1, B2, unw1, unw2)
    yh = maximum(vcat(unw1, unw2)) + 1
    yl = minimum(vcat(unw1, unw2)) - 1
    p1 = plot(B1, unw1, line = :nothing, marker = :o, color = :blue, ylims = (yl, yh))
    p2 = plot(B2, unw2, line = :nothing, marker = :o, color = :green, ylims = (yl, yh))
    plot(p1, p2)
end

function load_multi_temp(baselines...; station_name = nothing, rowcol = nothing, window = 1)
    intlists, idxs, vals, Bs, ccvals = [], [], [], [], []
    for baseline in baselines
        geolist, cur_intlist, cur_valid_idxs =
            InsarTimeseries.load_geolist_intlist(UNW_FILENAME, ignore_geo_file, baseline)
        if isnothing(station_name)
            cur_unw_vals = get_stack_vals(
                UNW_FILENAME,
                rowcol...,
                window,
                STACK_FLAT_DSET,
                cur_valid_idxs,
                reference_station = "TXKM",
            )
            # TODO: add back in when the stack sizes match up again
            # ccval = get_stack_vals(UNW_FILENAME, rowcol..., window, STACK_DSET, cur_valid_idxs)
        else
            cur_unw_vals = get_stack_vals(
                UNW_FILENAME,
                station_name,
                window,
                STACK_FLAT_DSET,
                cur_valid_idxs,
                reference_station = "TXKM",
            )
            # TODO: add back in when the stack sizes match up again
            # ccval = get_stack_vals(UNW_FILENAME, station_name, window, STACK_DSET, cur_valid_idxs)
        end
        B = InsarTimeseries.build_B_matrix(geolist, cur_intlist)

        push!(intlists, cur_intlist)
        push!(idxs, cur_valid_idxs)
        push!(vals, cur_unw_vals)
        push!(Bs, sum(B, dims = 2))
        # push!(ccvals, ccval)
    end
    return intlists, idxs, vals, Bs # , ccvals
end

function plot_multi_temp(Bs, vals, temps, title = ".unw vals at different baselines")
    n = length(Bs)

    sort!(temps, rev = true)
    sort!(Bs, by = x -> length(x), rev = true)
    sort!(vals, by = x -> length(x), rev = true)

    ym = maximum(maximum(abs.(v)) for v in vals) + 3
    # fig, axes = plt.subplots(1, n)
    fig, axes = plt.subplots()
    for ii = 1:n
        # axes[ii].plot(Bs[ii], vals[ii], ".")
        axes.plot(Bs[ii], vals[ii], ".", label = temps[ii])
    end
    axes.set_ylim((-ym, ym))
    axes.legend()
    axes.set_title(title)
    return fig, axes
end

function calc_exclude_dates(fname, dset)
    geolist = Sario.load_geolist_from_h5(fname, dset)
    excls = Sario.load(fname, "excluded/1")
    excl_map = InsarTimeseries.decode_excluded(excls, geolist)
    return geolist, get_excl_hist(excl_map, geolist)
end

# @time excl_map = InsarTimeseries.decode_excluded(excls, geolist)
function get_excl_hist(excl_map, geolist; as_pct = true)
    counts = zeros(length(geolist))
    for idx in eachindex(excl_map)
        was_excluded = .!isnothing.(indexin(geolist, excl_map[idx]))
        counts .+= was_excluded
    end
    return as_pct ? counts ./ length(excl_map) : counts
end


text:
using NeutronTransport
using Documenter
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

makedocs(bib;
    modules=[NeutronTransport],
    authors="Ramiro Vignolo <ramirovignolo@gmail.com>",
    repo="https://github.com/rvignolo/NeutronTransport.jl/blob/{commit}{path}#L{line}",
    sitename="NeutronTransport.jl",
    format=Documenter.HTML(;
        prettyurls=false,
        canonical="https://rvignolo.github.io/NeutronTransport.jl/stable",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Method of Characteristics" => "moc.md",
        "References" => "references.md"
    ],
)

deploydocs(;
    repo="github.com/rvignolo/NeutronTransport.jl",
    devbranch = "main"
)


text:
include("try.jl")

# solvers_list = [:Mosek, :Ipopt, :Knitro_Ip, :Knitro_IpCG, :Knitro_AS, :Knitro_SQP, :Cvxopt]

solvers_list = [:Mosek]

tmpFloatDatatype = Float64

list= String[]
for i in 1:5
    push!(list, "../PDFs/Data/gauss-$i")
end
size =[1,0.75,0.5,0.4]
list_instances=String[]
for l in 1:length(list)
    for s in size
        push!(list_instances, string(list[l],"-$s.dens"))
    end#^n
end#^s

for solver in solvers_list
    for instance in list_instances
        println("///////////////////////////")
        println("Starting computations of the distribution: ", instance," with ", solver)
        println("///////////////////////////")
        test(solver, instance, instance,  tmpFloatDatatype)
        println("///////////////////////////")
        println("Finishing computations of the distribution: ", instance," with ", solver)
        println("///////////////////////////")
    end
end

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
@testset "pml" begin

@testset "PMLParam" begin
    pml = PMLParam()
    @test pml.m===4.0 && pml.R===exp(-16) && pml.κmax===1.0 && pml.amax===0.0 && pml.ma===4.0
end  # @testset "PMLParam"

@testset "1D" begin
    Npmln = 5
    Npmlp = 7
    Npml = ([Npmln], [Npmlp])
    N = 10
    M = N + Npmln + Npmlp
    ∆ldual = rand(M)
    L = sum(∆ldual)
    l₀ = L / 2
    lprim = cumsum([-l₀; ∆ldual])
    ldual = StaggeredGridCalculus.movingavg(lprim)

    lpml = ([lprim[1+Npmln]], [lprim[end-Npmlp]])
    Lpml = ([lpml[nN][1]-lprim[1]], [lprim[end]-lpml[nP][1]])

    bounds = ([lprim[1]], [lprim[end]])
    lprim = lprim[1:end-1]

    @test StaggeredGridCalculus.pml_loc(tuple(lprim), bounds, Npml) == (lpml, Lpml)

    ω = rand()
    pml = PMLParam()
    σmax = (-(pml.m+1) * log(pml.R) / 2 ./ Lpml[nN], -(pml.m+1) * log(pml.R) / 2 ./ Lpml[nP])  # see calc_stretch_factor
    f(l) = l≤lpml[nN][1] ? (1 + σmax[nN][1] * ((lpml[nN][1]-l)/Lpml[nN][1])^pml.m / (im*ω)) : (l≥lpml[nP][1] ? (1 + σmax[nP][1] * ((l-lpml[nP][1])/Lpml[nP][1])^pml.m / (im*ω)) : 1)
    sprim = f.(lprim)
    sdual = f.(ldual)

    l = (tuple(lprim), tuple(ldual))
    sfactor = StaggeredGridCalculus.create_sfactor(ω, l, lpml, Lpml)
    @test sfactor == (tuple(sprim), tuple(sdual))

    ∆lprim = diff(ldual)
    pushfirst!(∆lprim, ldual[1]+L-ldual[end])  # periodic boundary

    ∆l = (tuple(∆lprim), tuple(∆ldual))
    s∆lprim = sprim .* ∆lprim
    s∆ldual = sdual .* ∆ldual

    s∆l = create_stretched_∆l(ω, ∆l, l, bounds, Npml)
    @test s∆l == (tuple(s∆lprim), tuple(s∆ldual))

    isbloch = [true]
    g1 = Grid(tuple([lprim..., bounds[nP][1]]), isbloch)
    @test s∆l ≈ create_stretched_∆l(ω, g1, Npml)
end  # @testset "1D"

@testset "1D, Npml = 0" begin
    Npml = ([0],[0])
    lprim = [0.0, 1.0]  # with ghost point
    ldual = [0.5]  # without ghost point

    lpml = ([0.0], [1.0])
    Lpml = ([0.0], [0.0])

    bounds = ([lprim[1]], [lprim[end]])
    lprim = lprim[1:end-1]

    @test StaggeredGridCalculus.pml_loc(tuple(lprim), bounds, Npml) == (lpml, Lpml)

    ω = rand()
    l = (tuple(lprim), tuple(ldual))
    sfactor = StaggeredGridCalculus.create_sfactor(ω, l, lpml, Lpml)
    @test sfactor == (tuple([1.0]), tuple([1.0]))

    ∆lprim  = ∆ldual = [1.0]
    ∆l = (tuple(∆lprim), tuple(∆ldual))

    s∆l = create_stretched_∆l(ω, ∆l, l, bounds, Npml)
    @test s∆l == (tuple([1.0]), tuple([1.0]))

    isbloch = [true]
    g1 = Grid(tuple([0.0,1.0]), isbloch)
    @test s∆l ≈ create_stretched_∆l(ω, g1, Npml)
end  # @testset "1D, Npml = 0"

@testset "3D" begin
    Npml = ([10,5,1], [9,6,4])
    N = [10, 12, 1]
    M = sum(Npml) + N
    ∆ldual = ntuple(d->rand(M[d]), length(N))

    L = SVec(sum.(∆ldual))  # SFloat{3}
    l₀ = L ./ 2  # SFloat{3}
    lprim = map((v,s)->v.-s, map(x->[0; cumsum(x)], ∆ldual), (l₀...,))  # tuple of vectors
    ldual = StaggeredGridCalculus.movingavg.(lprim)

    lpml = (SVec{3}((w->lprim[w][1+Npml[nN][w]]).(1:3)), SVec{3}((w->lprim[w][end-Npml[nP][w]]).(1:3)))
    Lpml = (lpml[nN].-SVec{3}((w->lprim[w][1]).(1:3)), SVec{3}((w->lprim[w][end]).(1:3)).-lpml[nP])

    bounds = ([v[1] for v = lprim], [v[end] for v = lprim])
    lprim = map(v->v[1:end-1], lprim)

    @test StaggeredGridCalculus.pml_loc(lprim, bounds, Npml) == (lpml, Lpml)

    ω = rand()
    pml = PMLParam()
    σmax = (-(pml.m+1) * log(pml.R) / 2 ./ Lpml[nN], -(pml.m+1) * log(pml.R) / 2 ./ Lpml[nP])  # see calc_stretch_factor
    f(l, nw) = l≤lpml[nN][nw] ? (1 + σmax[nN][nw] * ((lpml[nN][nw]-l)/Lpml[nN][nw])^pml.m / (im*ω)) : (l≥lpml[nP][nw] ? (1 + σmax[nP][nw] * ((l-lpml[nP][nw])/Lpml[nP][nw])^pml.m / (im*ω)) : 1)
    sprim = (f.(lprim[1],1), f.(lprim[2],2), f.(lprim[3],3))
    sdual = (f.(ldual[1],1), f.(ldual[2],2), f.(ldual[3],3))

    l = (lprim, ldual)
    sfactor = StaggeredGridCalculus.create_sfactor(ω, l, lpml, Lpml)
    @test sfactor == (sprim, sdual)

    ∆lprim = diff.(ldual)
    for k = 1:3
        pushfirst!(∆lprim[k], ldual[k][1]+L[k]-ldual[k][end])  # periodic boundary
    end
    ∆l = (∆lprim, ∆ldual)
    s∆lprim = map((x,y)->x.*y, sfactor[nPR], ∆l[nPR])
    s∆ldual = map((x,y)->x.*y, sfactor[nDL], ∆l[nDL])

    s∆l = create_stretched_∆l(ω, ∆l, l, bounds, Npml)
    @test s∆l == (s∆lprim, s∆ldual)

    isbloch = [true,true,true]
    g3 = Grid(ntuple(k->[lprim[k]..., bounds[nP][k]], 3), isbloch)
    @test s∆l ≈ create_stretched_∆l(ω, g3, Npml)
end  # @testset "3D"

@testset "3D, Npml = 0" begin
    Npml = ([0,0,0], [0,0,0])
    lprim = ([0.0,1.0], [0.0,1.0], [0.0,1.0])  # with ghost points
    ldual = ([0.5], [0.5], [0.5])  # without ghost points

    lpml = ([0.0,0.0,0.0], [1.0,1.0,1.0])
    Lpml = ([0.0,0.0,0.0], [0.0,0.0,0.0])

    bounds = ([v[1] for v = lprim], [v[end] for v = lprim])
    lprim = map(v->v[1:end-1], lprim)

    @test StaggeredGridCalculus.pml_loc(lprim, bounds, Npml) == (lpml, Lpml)

    ω = rand()
    l = (lprim, ldual)
    sfactor = StaggeredGridCalculus.create_sfactor(ω, l, lpml, Lpml)
    @test sfactor == (([1.0],[1.0],[1.0]), ([1.0],[1.0],[1.0]))

    ∆lprim  = ∆ldual = ([1.0], [1.0], [1.0])
    ∆l = (∆lprim, ∆ldual)

    s∆l = create_stretched_∆l(ω, ∆l, l, bounds, Npml)
    @test s∆l == (([1.0],[1.0],[1.0]), ([1.0],[1.0],[1.0]))

    isbloch = [true,true,true]
    g3 = Grid(([0.0,1.0], [0.0,1.0], [0.0,1.0]), isbloch)
    @test s∆l ≈ create_stretched_∆l(ω, g3, Npml)
end  # @testset "3D, Npml = 0"

end  # @testset "pml"


text:
using TimerOutputs
using Base.Test

import TimerOutputs: sections, ncalls

to = TimerOutput()

@timeit to "sleep" sleep(0.1)

@test "sleep" in keys(sections(to))

@timeit to "multi statement" begin
1+1
sleep(0.1)
end

@test "multi statement" in keys(sections(to))

@timeit to "sleep" sleep(0.1)
@timeit to "sleep" sleep(0.1)
@timeit to "sleep" sleep(0.1)

@test ncalls(sections(to)["sleep"]) == 4

# Check reset works
reset_timer!(to)

@test length(keys(sections(to))) == 0

# Check return values get propagated
function foo(a)
    a+a
end

to2 = TimerOutput()

a = @timeit to2 "foo" foo(5)

@test a === 10
@test "foo" in keys(sections(to2))

# Test nested
c = @timeit to2 "nest 1" begin
    sleep(0.1)
    @timeit to2 "nest 2" sleep(0.2)
    @timeit to2 "nest 1" sleep(0.2)
    5
end

@test ncalls(sections(to2)["nest 1"]) == 2
@test ncalls(sections(to2)["nest 2"]) == 1
@test c === 5

# test throws

function foo2(v)
    @timeit to "throwing" begin
        sleep(1)
        print(v[6]) # OOB
    end
end

try
    foo2(rand(5))
catch e
    isa(e, BoundsError) || rethrow(e)
end

@test "throwing" in keys(sections(to))

@test_throws ArgumentError @timeit to 5 sleep(1)

function foo3()
    to = TimerOutput()

    to_tester = TimerOutput()
    @timeit to "runtime" for i in 1:10^6
        @timeit to_tester "test" 1+1
    end

    return to, to_tester
end

to3, to4 = foo3()

@test "runtime" in keys(sections(to3))
@test "test" in keys(sections(to4))

function foo4()
    to = TimerOutput()

    enter_section(to, "sec1")
    sleep(0.25)
    enter_section(to, "sec2")
    sleep(0.25)
    exit_section(to)

    exit_section(to, "sec1")

    return to
end

to5 = foo4()

@test "sec1" in keys(sections(to5))
@test "sec2" in keys(sections(to5))

print(to5)


text:
using JuliaSets

Q(z) = z.^4 + 0.6 + 1im * 0.55
n = 1000

xs = [linspace(-2, 2, n);]
ys = copy(xs)

A = escape_time(Q, xs, ys, 20, 100.0)

# map to [0,2 \pi]
A = atan(A)

using Gadfly

draw(PNG("example.png", 20cm, 20cm), spy(A, Guide.xticks(ticks=nothing), Guide.yticks(ticks=nothing), Guide.xlabel(nothing), Guide.ylabel(nothing), Theme(key_position = :none)))


text:
function forward!(
    α::AbstractArray,
    c::AbstractArray,
    a::AbstractArray,
    A::AbstractArray,
    L::AbstractArray,
    n::Int,
)
    _, K, N = size(L)
    # (T == 0) && return
    α[1, :, n] = a .* L[1, :, n]
    c[1, n] = sum(α[1, :, n])
    α[1, :, n] /= c[1, n]
    T = length(filter(!isnothing, L[:, 1, n]))
    @inbounds for t = 2:T
        for j2 = 1:K
            @simd for j1 = 1:K
                α[t, j2, n] += α[t-1, j1, n] * A[j1, j2]
            end
        end
        for j = 1:K
            α[t, j, n] *= L[t, j, n]
            c[t, n] += α[t, j, n]
        end
        for j = 1:K
            α[t, j, n] /= c[t, n]
        end
    end
    @inbounds for t = T+1:length(L[:, 1, n])
        α[t, :, n] .= nothing
        c[t, n] = nothing
    end
end

function backward!(
    β::AbstractArray,
    c::AbstractArray,
    A::AbstractArray,
    L::AbstractArray,
    n::Int,
)
    _, K, N = size(L)
    T = length(filter(!isnothing, L[:, 1, n]))
    @inbounds for j = 1:K
        β[T, j, n] = 1.0
    end
    @inbounds for t in reverse(1:T-1)
        for j1 = 1:K
            @simd for j2 = 1:K
                β[t, j1, n] += β[t+1, j2, n] * A[j1, j2] * L[t+1, j2, n]
            end
        end
        for j = 1:K
            β[t, j, n] /= c[t+1, n]
        end
    end
    @inbounds for t = T+1:length(L[:, 1, n])
        β[t, :, n] .= nothing
    end
end

function posterior!(
    γ::AbstractArray,
    α::AbstractArray,
    β::AbstractArray,
    L::AbstractArray,
    n::Int,
)
    @argcheck size(α) == size(β)
    _, K, N = size(L)
    T = length(filter(!isnothing, L[:, 1, n]))
    @inbounds for t = 1:T
        @simd for j = 1:K
            γ[t, j, n] = α[t, j, n] * β[t, j, n]
        end
    end
end

function update_ξ!(
    ξ::AbstractArray,
    α::AbstractArray,
    β::AbstractArray,
    c::AbstractArray,
    A::AbstractMatrix,
    L::AbstractArray,
    n::Int,
)
    @argcheck size(α) == size(β)
    _, K, N = size(L)
    T = length(filter(!isnothing, L[:, 1, n]))
    @inbounds for t = 1:T-1
        for j1 = 1:K
            @simd for j2 = 1:K
                ξ[t, j1, j2, n] += (α[t, j1, n] * A[j1, j2] * L[t+1, j2, n]) / c[t+1, n]
            end
        end
    end
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
# we are missing Rayleigh Scattering

using Interpolations: LinearInterpolation, Throw
using StaticArrays: SA

using ..ContinuumAbsorption: hydrogenic_bf_absorption, hydrogenic_ff_absorption, ionization_energies

const _H_I_ion_energy = ionization_energies[1][1] # not sure if this is a good idea
const _H⁻_ion_energy = 0.7552 # eV

_H_I_bf(ν, T, nH_I_div_partition, ion_energy = _H_I_ion_energy, nmax_explicit_sum = 8,
        integrate_high_n = true) =
           hydrogenic_bf_absorption(ν, T, 1, nH_I_div_partition, ion_energy, nmax_explicit_sum,
                                    integrate_high_n)
"""
    H_I_bf(ν, T, nH_I_div_partition, [ion_energy], [nmax_explicit_sum], [integrate_high_n];
           kwargs...)

Compute the bound-free linear absorption coefficient contributed by all energy states of a
neutral Hydrogen atom.

# Required Arguments
- `ν::AbstractVector{<:Real}`: sorted frequency vector in Hz
- `T`: temperature in K
- `nH_I_div_partition` is the total number density of neutral Hydrogen divided by the its
   partition function.

# Optional Arguments
- `ion_energy`: The ionization energy of Hydrogen. By default, this is set to the values loaded
  into the global `ionization_energies` list.
- `nmax_explicit_sum`: The highest energy level whose absorption contribution is included
   in the explicit sum. The contributions from higher levels are included in an integral.
- `integrate_high_n::bool`: When this is `false`, bf absorption from higher energy states are not
   estimated at all. Default is `true`.

For a description of the kwargs, see [Continuum Absorption Kwargs](@ref).

# Notes
This function wraps [`hydrogenic_ff_absorption`](@ref). See that function for implementation
details.
"""
H_I_bf = bounds_checked_absorption(_H_I_bf; ν_bound = Interval(0, Inf),
                                   temp_bound = Interval(0, Inf))

_H_I_ff(ν, T, nH_II, ne) = hydrogenic_ff_absorption(ν, T, 1, nH_II, ne)

"""
    H_I_ff(ν, T, nH_II, ne; kwargs...)

Compute the H I free-free linear absorption coefficient α.

The naming scheme for free-free absorption is counter-inutitive. This actually refers to the
reaction:  `photon + e⁻ + H II -> e⁻ + H II`.

# Arguments
- `ν::AbstractVector{<:Real}`: sorted frequency vector in Hz
- `T`: temperature in K
- `nH_II`: the number density of ionized Hydrogen in cm⁻³.
- `ne`: the number density of free electrons.

For a description of the kwargs, see [Continuum Absorption Kwargs](@ref).

# Notes
This function wraps [`hydrogenic_ff_absorption`](@ref). See that function for implementation
details.
"""
H_I_ff = bounds_checked_absorption(_H_I_ff; ν_bound = λ_to_ν_bound(_gauntff_λ_bounds),
                                   temp_bound = _gauntff_T_bounds)

"""
    _ndens_Hminus(nH_I_div_partition, ne, T, ion_energy = _H⁻_ion_energy)

Compute the number density of H⁻ (implements eqn 5.10 of Kurucz 1970). This is an application of
the saha equation where the "ground state" is H⁻ and the "first ionization state" is H I. The
partition function of H⁻ is 1 at all temperatures.
"""
function _ndens_Hminus(nH_I_div_partition, ne, T, ion_energy = _H⁻_ion_energy)
    # The value of _H_I_ion_energy energy doesn't matter at all becasue the exponential evaluates
    # to 1 for n = 1
    nHI_groundstate = ndens_state_hydrogenic(1, nH_I_div_partition, T, _H_I_ion_energy)

    # should should probably factor this out
    if T < 1000
        # this isn't a hard requirement at all. This is just a sanity check in case an argument got
        # forgotten (which is allowed since ion_energy has a default value)
        # In reality, we should use the lower temperature limit in which the exponential evaluates
        # to Inf. That should probably be our domain limit
        throw(DomainError(T, "Unexpected T value."))
    end

    # coef = (h^2/(2*π*m))^1.5
    coef = 3.31283018e-22 # cm³*eV^1.5
    β = 1.0/(kboltz_eV*T)
    0.25 * nHI_groundstate * ne * coef * β^1.5 * exp(ion_energy * β)
end

const _Hminus_bf_table = SA[1250.0 5.431;   1500.0 6.512;   1750.0 7.918;   2000.0 9.453;
                            2250.0 11.08;   2500.0 12.75;   2750.0 14.46;   3000.0 16.19;
                            3250.0 17.92;   3500.0 19.65;   3750.0 21.35;   4000.0 23.02;
                            4250.0 24.65;   4500.0 26.24;   4750.0 27.77;   5000.0 29.23;
                            5250.0 30.62;   5500.0 31.94;   5750.0 33.17;   6000.0 34.32;
                            6250.0 35.37;   6500.0 36.32;   6750.0 37.17;   7000.0 37.91;
                            7250.0 38.54;   7500.0 39.07;   7750.0 39.48;   8000.0 39.77;
                            8250.0 39.95;   8500.0 40.01;   8750.0 39.95;   9000.0 39.77;
                            9250.0 39.48;   9500.0 39.06;   9750.0 38.53;   10000.0 37.89;
                            10250.0 37.13;  10500.0 36.25;  10750.0 35.28;  11000.0 34.19;
                            11250.0 33.01;  11500.0 31.72;  11750.0 30.34;  12000.0 28.87;
                            12250.0 27.33;  12500.0 25.71;  12750.0 24.02;  13000.0 22.26;
                            13250.0 20.46;  13500.0 18.62;  13750.0 16.74;  14000.0 14.85;
                            14250.0 12.95;  14500.0 11.07;  14750.0 9.211;  15000.0 7.407;
                            15250.0 5.677;  15500.0 4.052;  15750.0 2.575;  16000.0 1.302;
                            16100.0 0.8697; 16200.0 0.4974; 16300.0 0.1989]

const _Hminus_bf_cross_section_interp = LinearInterpolation(view(_Hminus_bf_table, :, 1),
                                                            view(_Hminus_bf_table, :, 2),
                                                            extrapolation_bc=Throw())

"""
    _Hminus_bf_cross_section(λ, ion_energy_H⁻)

Compute the H⁻ bound-free cross-section at a given wavelength (specified in Å). The cross-section
has units of megabarns per H⁻ particle and does NOT include a correction for stimulated emission.

This function linearly interpolates the table provided in Wishart (1979).
"""
function _Hminus_bf_cross_section(λ, ion_energy_H⁻)
    Å_per_eV = 1e8 * (hplanck_eV * c_cgs)
    max_λ_ionize = Å_per_eV/ion_energy_H⁻

    last_table_λ = _Hminus_bf_table[end, 1]
    last_table_cross_section = _Hminus_bf_table[end, 2]

    @assert max_λ_ionize > last_table_λ

    if λ < _Hminus_bf_table[1,1]
        throw(DomainError(λ, "λ must be ≥ $(_Hminus_bf_table[1,1]) Å"))
    elseif λ <= last_table_λ
        _Hminus_bf_cross_section_interp(λ)
    elseif λ <= max_λ_ionize
        # linearly interpolate b/t (last_table_λ, last_table_cross_section) and (max_λ_ionize, 0.0)
        m = last_table_cross_section/(last_table_λ - max_λ_ionize)
        m*(λ - max_λ_ionize)
    else
        0.0
    end

end

function _Hminus_bf(ν::Real, T::Real, nH_I_div_partition::Real, ne::Real,
                    ion_energy_H⁻::Real = _H⁻_ion_energy)
    λ = c_cgs*1e8/ν # in ångstroms
    cross_section = _Hminus_bf_cross_section(λ, ion_energy_H⁻) # in units of megabarn
    # convert from megabarn to cm² and include contributions from stimulated emission  1e-18
    cross_section *= (1 - exp(-hplanck_cgs*ν/(kboltz_cgs*T))) * 1e-18
    _ndens_Hminus(nH_I_div_partition, ne, T, _H⁻_ion_energy) * cross_section
end

"""
    Hminus_bf(ν, T, nH_I_div_partition, ne, [ion_energy_H⁻]; kwargs...)

Compute the H⁻ bound-free linear absorption coefficient α

# Arguments
- `ν::AbstractVector{<:Real}`: sorted frequency vector in Hz
- `T`: temperature in K
- `nH_I_div_partition`: the total number density of H I divided by its partition function.
- `ne`: the electron number density
- `ion_energy_H⁻`: Specifies the ionization energy of the single state of H⁻ in eV. This is
   roughly 0.7552 eV.

For a description of the kwargs, see [Continuum Absorption Kwargs](@ref).

# Notes
This function assumes that n(H⁻) ≪ n(H I) + n(H II). The number density of n(H⁻) should not be
pre-computed (instead it's computed internally by this function).

This function is adapted from equation 8.12 from Grey (2005). That equation specifies the
absorption coefficient per H I atom (uncorrected by stimulated emission) as:

``8.316 \\times 10^{-10} \\sigma_{bf}(H^-) P_e \\theta^{2.5} 10^{(\\theta \\chi_{\\mathrm{H}^-})} U(H⁻,T)/ U(H I,T)``

where:
- ``\\sigma_{bf}(\\mathrm{H}^-)`` is the photo dissociation cross-section (Grey uses the variable
  ``\\alpha_{bf}`` in place of ``\\sigma_{bf}``). We estimate this by linearly interpolating data
  from [Wishart (1979)].
- ``\\chi_{\\mathrm{H}^-}`` is the ionization enrgy of H⁻.
- θ = log10(e)/(k*T) or θ = 5040/T in units of eV⁻¹
- U(H⁻,T) is the partition function of H⁻ at temperature T. This is always 1
- U(H I,T) is the partition function of H I at temperature T. This is 2 at low to intermediate T.
Eqn 8.12 of Grey (2005) implicitly assumes that (U(H⁻,T)/ U(H I,T)) is always equal to 0.5.

This expression is simply a rewritten form of ``\\sigma_{bf}(H^-) n(H^-)/n(H I),`` where
``n(H^-)/n(H I)`` has been replaced with the expanded form of the saha equation, in which
n₀ = n(H⁻) and n₁ = n(H I).

Equations 8.18 and 8.19 of Gray (2005), indicate that the H⁻ bound-free opacity (with the
stimulated emission correction) is given by:

``\\kappa_\\nu = \\sigma_{bf}(H^-) \\frac{n(H^-)}{n(H I)} (1 - \\exp \\left( \\frac{-h\\nu}{k T}\\right)) \\frac{n(H I)}{n(H I) + n(H II)} \\frac{n(H I) + n(H II)}{\\rho}``

In other words, the linear absorption coefficient is: ``\\alpha_\\nu = \\sigma_{bf}(H^-) n(H⁻) (1 - \\exp \\left( \\frac{-h\\nu}{k T}\\right))``

Wishart (1979) expects the tabulated data to have better than 1% percent accuracy. Mathisen (1984)
suggests that this data has better than 3% accuracy.
"""
Hminus_bf = bounds_checked_absorption(
    _Hminus_bf;
    ν_bound = λ_to_ν_bound( Interval(1.25e-5, Inf; exclusive_lower=false, exclusive_upper=true) ),
    temp_bound = Interval(0, Inf)
)

function _Hminus_ff(ν::Real, T::Real, nH_I_div_partition::Real, ne::Real)
    λ = c_cgs*1e8/ν # in Angstroms

    logλ = log10(λ)
    log2λ = logλ * logλ
    log3λ = log2λ * logλ
    log4λ = log3λ * logλ

    f0 =  -2.2763 -   1.6850 * logλ +  0.76661 * log2λ -  0.053346 * log3λ
    f1 = +15.2827 -   9.2846 * logλ +  1.99381 * log2λ -  0.142631 * log3λ
    f2 = -197.789 + 190.266  * logλ - 67.9775  * log2λ + 10.6913   * log3λ - 0.625151 * log4λ

    logθ = log10(5040.0/T)
    αff_H⁻ = 1e-26 * 10.0^(f0 + f1 * logθ + f2 * (logθ * logθ))
    # Pₑ * α_ff(H⁻) gives the absorption coefficient in units of cm² per ground state H I atom
    Pₑ = ne * kboltz_cgs * T

    # Account for the fact that n(H I, n=1) might be slightly smaller than the entire number
    # density of H I. There is only really a difference at the highest temperatures. For the
    # temperature range where this approximation is valid, less than 0.23% of all H I atoms are not
    # in the ground state.

    # this calculation could reduce to nHI_gs = 2.0*nH_I_div_partition
    nHI_gs = ndens_state_hydrogenic(1, nH_I_div_partition, T, _H_I_ion_energy)

    return αff_H⁻ * Pₑ * nHI_gs
end

"""
    Hminus_ff(ν, T, nH_I_div_partition, ne; kwargs...)

Compute the H⁻ free-free linear absorption coefficient α

The naming scheme for free-free absorption is counter-inutitive. This actually refers to the
reaction:  photon + e⁻ + H I -> e⁻ + H I.

# Arguments
- `ν::AbstractVector{<:Real}`: sorted frequency vector in Hz
- `T`: temperature in K
- `nH_I_div_partition::Flt`: the total number density of H I divided by its partition function.
- `ne`: the number density of free electrons.

For a description of the kwargs, see [Continuum Absorption Kwargs](@ref).

# Notes
This is taken from equation 8.13 of Gray (2005). The equation uses a polynomial fig against Table 1
of [Bell & Berrington's (1987)](https://doi.org/10.1088/0022-3700/20/4/019), which tabulates values
for "the H⁻ absorption coefficient" (the values include the correction for stimulated emission). We
denote this quantity as `f(H⁻)` (Grey instead denotes it with the variable, ``\\alpha_{ff}(H^-)``).
The polynomial fits `f(H⁻)` in the range 2520 K ≤ T ≤ 10080 K and 2604 Å ≤ λ ≤ 113918 Å.
According to Grey, the polynomial fit to the tabulated data typically has 1% precision. We find
that at worst, the discrepancy never exceeds 2.25%.

Based on equations 8.13, 8.18, and 8.19 of Gray (2005), the free-free linear absorption coefficient
of H⁻ is given by: `α = f(H⁻) * Pₑ * n(H I)`.

Based on Section 5.3 from Kurucz (1970), I'm fairly confident that the "more correct" version of
this equation should actually read `α = f(H⁻) * Pₑ * n(H I, n = 1)`, and that the form inferred
from Gray (2005) implicitly assumes that `n(H I, n = 1) ≈ n(H I)`. In detail, the Boltzmann
equation relates `n(H I, n = 1)` and `n(H I)` as:
```
 n(H I, n = 1) = n(H I)*gₙ₌₁/U(T)*exp(-Eₙ₌₁/(k*T)) = n(H I) * 2/U(T)*exp(0) = n(H I) * 2/U(T).
 ```
Based on the values of the partition function in Table D.2 of Grey, this approximation seems fairly
good. It only introduces a ≤0.5% overestimate in n(H I, n = 1) over the temperature range where the
polynomial is valid.

We also considered the polynomial fit in Section 5.3 from Kurucz (1970). Unfortunately, it seems
to be wrong (it gives lots of negative numbers).
"""
Hminus_ff = bounds_checked_absorption(_Hminus_ff,
                                      ν_bound = λ_to_ν_bound(closed_interval(2.604e-5,1.13918e-3)),
                                      temp_bound = closed_interval(2520, 10080))

function _H2plus_bf_and_ff(ν::Real, T::Real, nH_I::Real, nH_II::Real)
    λ = c_cgs*1e8/ν # in ångstroms
    β_eV = 1.0/(kboltz_eV * T)
    stimulated_emission_correction = (1 - exp(-hplanck_eV*ν*β_eV))

    K = Stancil1994.K_H2plus(T) # n(H I) n(H II) / n(H₂⁺)
    σbf = Stancil1994.σ_H2plus_bf(λ, T)
    σff = Stancil1994.σ_H2plus_ff(λ, T)

    (σbf/K + σff) * nH_I * nH_II * stimulated_emission_correction
end

"""
    H2plus_bf_and_ff(ν, T, nH_I_div_partition, n_HII; kwargs...)

Compute the combined H₂⁺ bound-free and free-free linear absorption coefficient α using the tables
from [Stancil 1994](https://ui.adsabs.harvard.edu/abs/1994ApJ...430..360S/abstract).

# Arguments
- `ν::AbstractVector{<:Real}`: sorted frequency vector in Hz
- `T`: temperature in K
- `nH_I`: the total number density of H I divided by its partition
   function.
- `nH_II`: the number density of H II (not of H₂⁺).

For a description of the kwargs, see [Continuum Absorption Kwargs](@ref).

# Notes
This computes the H₂⁺ number density from those of H I and H II, since ionized molecules
are not included in the molecular equlibrium calculation.  Stancil provides approximate equilibrium
constants, K, for the molecule, but the Barklem and Collet values used elsewhere by Korg may be more
reliable.  Once ionized molecules are fully supported, those values should be used instead. While
cross sections are tabulated down to only 3150 K, the cross sections could be linearly interpolated
1000 K or so lower, if reliable K values are availble.

Because n(H₂⁺) is computed on the fly, the "cross-sections" used (internally) by this function have units of
cm^-5, since they must be multiplied by n(H I) and n(H II) to obtain absorption coefficients.

Stancil, [Bates (1952)](https://ui.adsabs.harvard.edu/abs/1952MNRAS.112...40B/abstract), and
Gray (2005) all use n(H I), but Kurucz (1970) notes in section 5.2 that the photodisociation of H₂⁺
produces a ground state H atom and that we should therefore use n(H I, n=1) instead. This difference
causes a descrepancy of a fraction of a percent (at most) for T ≤ 1.2e4 K. Here we use n(H I).
"""
H2plus_bf_and_ff = bounds_checked_absorption(
    _H2plus_bf_and_ff;
    ν_bound = λ_to_ν_bound(closed_interval(7e-6, 2e-3)),
    temp_bound = closed_interval(3150, 25200)
)


text:
function showresult(value::AbstractString, opts; under = false, html = false)
  raise(opts[:editor], "julia.result",
        @d("value" => value,
          "bounds" => opts[:bounds],
          "under"  => under,
          "html"   => html,
          "info"   => opts[:info]))
end

function showexception(req, value::AbstractString, bounds)
  raise(req, "julia.error",
        @d("value" => value,
           "start" => bounds[1],
           "end"   => bounds[2]))
end

function showexception(editor, e, bt, bounds)
  notify_error(sprint(showerror, e))
  showexception(
    editor,
    sprint(showerror_html, e, bt, :include_string),
    bounds)
end


text:
# See https://en.wikipedia.org/wiki/Tower_of_Hanoi
# Original code from https://rosettacode.org/wiki/Towers_of_Hanoi#Julia

function hanoi(n, from, to, via)
    if n != 1
        hanoi(n - 1, from, via, to)
        hanoi(1, from, to, via)
        hanoi(n - 1, via, to, from)
    end
    true
end

expressions = [
    :(hanoi(4, 1, 2, 3))
    :(hanoi(24, 1, 2, 3))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
#! /usr/bin/julia

# Rosetta Code, Abundant, deficient and perfect number classifications

# A natural number, n, can be written as a product of powers of it's
# prime factors, n = Prod(p_i^(a_i)).  The sum of all of n's divisors,
# including n itself, is dsum = Prod((p_i^(a_i+1) - 1)/(p_i - 1)).

const L = 2*10^4

function pcontrib(p::Int64, a::Int64)
    # (p^(a+1) - 1)/(p_i - 1) = p^a + p^(a-1) + ... + p + 1, Perhaps
    # this could be made more efficient by caching results to avoid
    # repeated calculations.
    n = one(p)
    pcon = one(p)
    for i in 1:a
        n *= p
        pcon += n
    end
    return pcon
end

function divisorsum(n::Int64)
    #the "proper divisor sum", or aliquot sum of n's divisors,
    #excluding n itself
    dsum = one(n)
    for (p, a) in factor(n)
        dsum *= pcontrib(p, a)
    end
    dsum -= n
end

# Work with a three element array, rather than three separate
# variables, and take advantage of the fact that the sign of dsum - n
# depends upon its class to increment this array of classifications.

iclasslabel = ["Deficient", "Perfect", "Abundant"]
iclass = zeros(Int64, 3)
iclass[1] = one(Int64) #by convention 1 is deficient

for n in 2:L
    if isprime(n)
        iclass[1] += 1
    else
        iclass[sign(divisorsum(n)-n)+2] += 1
    end
end

println("Classification of integers from 1 to ", L)
for i in 1:3
    println("   ", iclasslabel[i], ", ", iclass[i])
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
## TODO: simplify forward function through forward! (etc...)

export PlanarFlow

import Base: eltype, size, length

@doc raw"""
The PlanarFlow function is defined as

```math
f({\bf{x}}) = {\bf{x}} + {\bf{u}} \tanh({\bf{w}}^\top {\bf{x}} + b)
```

with input and output dimension ``D``. Here ``{\bf{x}}\in \mathbb{R}^D`` represents the input of the function. Furthermore ``{\bf{u}}\in \mathbb{R}^D``, ``{\bf{w}}\in \mathbb{R}^D`` and ``b\in\mathbb{R}`` represent the parameters of the function. The function contracts and expands the input space.

This function has been introduced in:

Rezende, Danilo, and Shakir Mohamed. "Variational inference with normalizing flows." _International conference on machine learning._ PMLR, 2015.
"""
mutable struct PlanarFlow{T1, T2 <: Real} <: AbstractCouplingFlow
    u       :: T1
    w       :: T1
    b       :: T2
    function PlanarFlow(u::T1, w::T1, b::T2) where { T1, T2 <: Real}
        @assert length(u) == length(w) "The parameters u and w in the PlanarFlow structure should have equal length."
        return new{T1,T2}(u, w, float(b))
    end
end
struct PlanarFlowEmpty{N} <: AbstractCouplingFlowEmpty end
PlanarFlowEmpty(dim::Int) = PlanarFlowEmpty{dim}()
struct PlanarFlowPlaceholder <: AbstractCouplingFlowPlaceholder end
@doc raw"""
The `PlanarFlow()` function creates a planar flow object. Its dimensionality is automatically set when wrapping this object inside a model. Its parameters are initialized during compilation.
"""
PlanarFlow() = PlanarFlowPlaceholder()
function prepare(dim::Int, flow::PlanarFlowPlaceholder)
    return PlanarFlowEmpty(dim)
end

# compile placeholder
compile(f::PlanarFlowEmpty{1})          = PlanarFlow(randn(), randn(), randn())
compile(f::PlanarFlowEmpty)             = PlanarFlow(randn(getdim(f)), randn(getdim(f)), randn())
compile(f::PlanarFlowEmpty{1}, params)  = PlanarFlow(params[1], params[2], params[3])
compile(f::PlanarFlowEmpty, params)     = PlanarFlow(params[1:getdim(f)], params[1+getdim(f):2*getdim(f)], params[2*getdim(f)+1])

@doc raw"""
The `PlanarFlow(dim::Int64)` function creates a mutable `PlanarFlow` structure with parameters corresponding to input of dimensions `dim`. The parameters are each random sampled from a standard (multivariate) normal distribution.
"""
function PlanarFlow(dim::Int64)
    if dim == 1
        return PlanarFlow(randn(), randn(), randn())
    else
        return PlanarFlow(randn(dim), randn(dim), randn())
    end
end

# number of parameters
nr_params(flow::PlanarFlow)      = length(flow.u) + length(flow.w) + length(flow.b)
nr_params(flow::PlanarFlowEmpty{N}) where {N} = 2*N + 1

# get-functions for the PlanarFlow structure.
getu(f::PlanarFlow)                         = return f.u
getw(f::PlanarFlow)                         = return f.w
getb(f::PlanarFlow)                         = return f.b
getall(f::PlanarFlow)                       = return f.u, f.w, f.b
getdim(f::PlanarFlowEmpty{N}) where { N }   = return N

# set-functions for the PlanarFlow structure
function setu!(f::PlanarFlow{T1,T2}, u::T1) where { T1, T2 <: Real}
    @assert length(f.u) == length(u) "The dimensionality of the current value of u and its new value do not match."
    f.u = u
end

function setw!(f::PlanarFlow{T1,T2}, w::T1) where { T1, T2 <: Real}
    @assert length(f.w) == length(w) "The dimensionality of the current value of w and its new value do not match."
    f.w = w
end

function setb!(f::PlanarFlow{T1,T2}, b::T2) where { T1, T2 <: Real }
    f.b = b
end

# custom Base function for the PlanarFlow structure
eltype(f::PlanarFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                 = promote_type(T1, T2)
eltype(f::PlanarFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}       = promote_type(eltype(T1), T2)
eltype(::Type{PlanarFlow{T1,T2}}) where { T1 <: Real, T2 <: Real}            = promote_type(T1, T2)
eltype(::Type{PlanarFlow{T1,T2}}) where { T1 <: AbstractVector, T2 <: Real}  = promote_type(eltype(T1), T2)

size(f::PlanarFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                   = 1
size(f::PlanarFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}         = length(f.u)
size(f::PlanarFlowEmpty{N}) where { N }                                      = return N

length(f::PlanarFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                 = 1
length(f::PlanarFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}       = length(f.u)
length(f::PlanarFlowEmpty{N}) where { N }                                    = return N

# forward pass through the PlanarFlow function (multivariate input)
function _forward(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real }

    # fetch values
    u, w, b = getall(f)

    # calculate result
    result = u.*ones(size(u))
    result .*= tanh(dot(w, input) + b)
    result .+= input

    # return result
    return result

end
forward(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::PlanarFlow{T1,T2}, input::AbstractVector{T1}) where { T1, T2 <: Real } = broadcast(_forward, Ref(f), input)

# forward pass through the PlanarFlow function (univariate input)
function _forward(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real }

    # fetch values
    u, w, b = getall(f)

    # calculate result
    result = copy(u)
    result *= tanh(dot(w, input) + b)
    result += input

    # return result
    return result

end
forward(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::PlanarFlow{T1,T2}, input::AbstractVector{ <: Real }) where { T1 <: Real, T2 <: Real } = broadcast(_forward, Ref(f), input)

function _forward(f::PlanarFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real }
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the PlanarFlow flow."
    return forward(f, input[1])
end
forward(f::PlanarFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::PlanarFlow{T1,T2}, input::AbstractVector) where { T1 <: Real, T2 <: Real} = broadcast(_forward, Ref(f), input)

# inplace forward pass through the PlanarFlow function (multivariate input)
function forward!(output::T1, f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real }

    # check dimensionality
    @assert length(output) == length(input) "The length of the preallocated vector does not seem to match the length of the input vector."

    # fetch values
    u, w, b = getall(f)

    # calculate result
    output .= u
    output .*= tanh(dot(w, input) + b)
    output .+= input

end

# jacobian of the PlanarFlow function (multivariate input)
function _jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real}

    # fetch values
    u, w, b = getall(f)

    # calculate result
    result = u*w'
    result .*= dtanh(dot(w, input) + b)
    @inbounds for k = 1:length(input)
        result[k,k] += 1.0
    end

    # return result
    return result

end
jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::PlanarFlow{T1,T2}, input::AbstractVector{T1}) where { T1, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

# jacobian of the PlanarFlow function (univariate input)
function _jacobian(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real }

    # fetch values
    u, w, b = getall(f)

    # calculate result (optimized)
    result = u * w * dtanh(w * input + b) + 1

    # return result
    return result

end
jacobian(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::PlanarFlow{T1,T2}, input::AbstractVector{ <: Real}) where { T1 <: Real, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

function _jacobian(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 }
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the PlanarFlow flow."
    return jacobian(f, input[1])
end
jacobian(f::PlanarFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::PlanarFlow{T1,T2}, input::AbstractVector) where { T1 <: Real, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

# inplace jacobian of the PlanarFlow function (multivariate input)
function jacobian!(output::AbstractMatrix{T2}, f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real}

    # check whether the dimensionality is correct
    @assert size(output) == (length(input), length(f.u)) "The dimensionality of the preallocated jacobian matrix seems incorrect."

    # fetch values
    u, w, b = getall(f)

    # calculate result
    for ku = 1:length(u)
        for kw = 1:length(w)
            output[ku,kw] = u[ku]*w[kw]
        end
    end
    output .*= dtanh(dot(w, input) + b)
    @inbounds for k = 1:length(input)
        output[k,k] += 1.0
    end

end

# determinant of the jacobian of the PlanarFlow function (multivariate input)
det_jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real} = det(jacobian(f, input))

# determinant of the jacobian of the PlanarFlow function (univariate input)
function det_jacobian(f::PlanarFlow{T, T}, input::T) where { T <: Real }

    # fetch values
    u, w, b = getall(f)

    # calculate result
    result = 1 + dot(u, w)*dtanh(dot(w, input) + b)

    # return result
    return result

end

# extra utility function (multivariate)
inv_jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real }       = inv(jacobian(f, input))
absdet_jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real }    = abs(det_jacobian(f, input))
logabsdet_jacobian(f::PlanarFlow{T1,T2}, input::T1) where { T1, T2 <: Real}  = logabsdet(jacobian(f, input))

# extra utility functions (univariate)
inv_jacobian(f::PlanarFlow{T,T}, input::T) where { T <: Real }               = 1.0 / jacobian(f, input)
absdet_jacobian(f::PlanarFlow{T,T}, input::T) where { T <: Real }            = abs(det_jacobian(f, input))
logabsdet_jacobian(f::PlanarFlow{T,T}, input::T) where { T <: Real}          = log(absdet_jacobian(f, input))

function inv_jacobian(f::PlanarFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real }
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the PlanarFlow flow."
    return inv_jacobian(f, input[1])
end