text:
# this file has not been completed yet

export BDD2

"""
type alias
"""

const UniqueTable2{Ts} = Dict{Tuple{HeaderID,NodeID,NodeID,Bool},AbstractBDDNode{Ts}}
const BinCache2{Ts} = Dict{Tuple{NodeID,Bool,NodeID,Bool},Tuple{AbstractBDDNode{Ts},Bool}}

"""
Node2
"""

struct BDD2Node{Ts} <: AbstractBDDNode{Ts}
    id::NodeID
    header::BDDNodeHeader{Ts}
    low::AbstractBDDNode{Ts}
    high::AbstractBDDNode{Ts}
    neg::Bool
end

"""
BDD2{Ts}

The structure for BDD. This includes all the nodes and caches.
"""

struct BDD2{Ts}
    manager::BDDManager
    headers::Dict{Ts,BDDNodeHeader{Ts}}
    utable::UniqueTable2{Ts}
    zero::BDDTerminal{Ts}
    andcache::BinCache2{Ts}
    orcache::BinCache2{Ts}
    xorcache::BinCache2{Ts}
end

function BDD2(::Type{Ts} = Symbol) where Ts
    manager = BDDManager()
    zero = BDDTerminal{Ts}(nodeid!(manager))
    BDD2{Ts}(
        manager,
        Dict{Ts,BDDNodeHeader{Ts}}(),
        UniqueTable2{Ts}(),
        zero,
        BinCache2{Ts}(),
        BinCache2{Ts}(),
        BinCache2{Ts}()
    )
end

function var(b::BDD2{Ts}, label::Ts)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    h = get(b.headers, label) do
        h = BDDNodeHeader(headerid!(b.manager), level!(b.manager), label)
        b.headers[label] = h
    end
    _node(b, h, b.zero, false, b.zero, true)
end

function _node(b::BDD2{Ts}, h::BDDNodeHeader{Ts}, low::AbstractBDDNode{Ts}, blow::Bool, high::AbstractBDDNode{Ts}, bhigh::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if blow == bhigh && low.id == high.id
        (low, blow)
    else
        flag = blow != bhigh
        key = (h.id, low.id, high.id, flag)
        node = get(b.utable, key) do
            b.utable[key] = BDD2Node(nodeid!(b.manager), h, low, high, flag)
        end
        (node, blow)
    end
end

### binoperator

function not(b::BDD2{Ts}, f::Tuple{AbstractBDDNode{Ts},Bool}) where Ts
    (f[1], !f[2])
end

function and(b::BDD2{Ts}, f::Vararg{Tuple{AbstractBDDNode{Ts},Bool}}) where Ts
    ans, bans = (b.zero, true)
    for x = f
        ans, bans = _binapply(BDDAnd(), b.andcache, b, ans, bans, x[1], x[2])
    end
    (ans, bans)
end

function or(b::BDD2{Ts}, f::Vararg{Tuple{AbstractBDDNode{Ts},Bool}}) where Ts
    ans, bans = (b.zero, false)
    for x = f
        ans, bans = _binapply(BDDOr(), b.orcache, b, ans, bans, x[1], x[2])
    end
    (ans, bans)
end

function xor(b::BDD2{Ts}, f::Tuple{AbstractBDDNode{Ts},Bool}, g::Tuple{AbstractBDDNode{Ts},Bool}) where Ts
    _binapply(BDDXor(), b.xorcache, b, f, g)
end

function imp(b::BDD2{Ts}, f::Tuple{AbstractBDDNode{Ts},Bool}, g::Tuple{AbstractBDDNode{Ts},Bool}) where Ts
    or(b, bddnot(b, f), g)
end

function ite(b::BDD2{Ts}, f::Tuple{AbstractBDDNode{Ts},Bool}, g::Tuple{AbstractBDDNode{Ts},Bool}, h::Tuple{AbstractBDDNode{Ts},Bool}) where Ts
    or(b, and(b, f, g), and(b, not(b, f), h))
end

### primitive

function getzero(f::BDD2Node{Ts}, bf::Bool) where Ts
    (f.low, bf)
end

function getone(f::BDD2Node{Ts}, bf::Bool) where Ts
    (f.high, bf != f.neg)
end

function _binapply(op::AbstractBDDOperator, cache::BinCache2{Ts}, b::BDD2{Ts}, f::BDD2Node{Ts}, bf::Bool, g::BDD2Node{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    key = (f.id, bf, g.id, bg)
    get(cache, key) do
        if f.header.level > g.header.level
            n0, bn0 = _binapply(op, cache, b, getzero(f, bf)..., g, bg)
            n1, bn1 = _binapply(op, cache, b, getone(f, bf)..., g, bg)
            ans, bans = _node(b, f.header, n0, bn0, n1, bn1)
        elseif f.header.level < g.header.level
            n0, bn0 = _binapply(op, cache, b, f, bf, getzero(g, bg)...)
            n1, bn1 = _binapply(op, cache, b, f, bf, getone(g, bg)...)
            ans, bans = _node(b, g.header, n0, bn0, n1, bn1)
        else
            n0, bn0 = _binapply(op, cache, b, getzero(f, bf)..., getzero(g, bg)...)
            n1, bn1 = _binapply(op, cache, b, getone(f, bf)..., getone(g, bg)...)
            ans, bans = _node(b, f.header, n0, bn0, n1, bn1)
        end
        cache[key] = (ans, bans)
    end
end

## and

function _binapply(::BDDAnd, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDD2Node{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == true
        (g, bg)
    else
        (b.zero, false)
    end
end

function _binapply(::BDDAnd, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDD2Node{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bg == true
        (f, bf)
    else
        (b.zero, false)
    end
end

function _binapply(::BDDAnd, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == true && bg == true
        (b.zero, true)
    else
        (b.zero, false)
    end
end

## or

function _binapply(::BDDOr, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDD2Node{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == true
        (b.zero, true)
    else
        (g, bg)
    end
end

function _binapply(::BDDOr, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDD2Node{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bg == true
        (b.zero, true)
    else
        (f, bf)
    end
end

function _binapply(::BDDOr, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == false && bg == false
        (b.zero, false)
    else
        (b.zero, true)
    end
end

## xor

function _binapply!(::BDDXor, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDD2Node{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == true
        (g, !bg)
    else
        (g, bg)
    end
end

function _binapply!(::BDDXor, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDD2Node{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bg == true
        (f, !bf)
    else
        (f, bf)
    end
end

function _binapply!(::BDDXor, ::BinCache2{Ts}, b::BDD2{Ts}, f::BDDTerminal{Ts}, bf::Bool, g::BDDTerminal{Ts}, bg::Bool)::Tuple{AbstractBDDNode{Ts},Bool} where Ts
    if bf == bg
        (b.zero, false)
    else
        (b.zero, true)
    end
end

##

# function fteval(b::BDD{Ts}, f::Tuple{AbstractNode{Ts},Bool}, env::Dict{Ts,Tx})::Tx where {Ts,Tx}
#     cache = Dict{Tuple{UInt,Bool},Tx}()
#     _fteval(b, f, env, cache)
# end

# function _fteval(b::BDD{Ts}, f::Tuple{Node{Ts},Bool}, env::Dict{Ts,Tx}, cache::Dict{Tuple{UInt,Bool},Tx})::Tx where {Ts,Tx}
#     get(cache, (f[1].id, f[2])) do
#         p = env[f[1].header.label]
#         fprob = (1-p) * _fteval(b, getzero(f), env, cache) + p * _fteval(b, getone(f), env, cache)
#         cache[(f[1].id,!f[2])] = 1 - fprob
#         cache[(f[1].id,f[2])] = fprob
#     end
# end

# function _fteval(b::BDD{Ts}, f::Tuple{Terminal{Ts},Bool}, env::Dict{Ts,Tx}, cache::Dict{Tuple{UInt,Bool},Tx})::Tx where {Ts,Tx}
#     (f[2] == false) ? Tx(0) : Tx(1)
# end

# ##

# mutable struct MinPath
#     len::Int
#     set::Vector{Vector{Bool}}
# end

# function ftmcs(b::BDD{Ts}, f::Tuple{AbstractNode{Ts},Bool}) where Ts
#     result = Vector{Ts}[]
#     r = f
#     while r != (b.zero, false)
#         f, s = _ftmcs(b, r)
#         append!(result, s)
#         r = bddnot(b, bddimp(b, r, f))
#     end
#     result
# end

# function _ftmcs(b::BDD{Ts}, f::Tuple{AbstractNode{Ts},Bool}) where Ts
#     vars = Dict([x.level => var!(b, x.label) for (k,x) = b.headers])
#     path = [false for i = 1:Int(b.totalvarid)]
#     s = MinPath(Int(b.totalvarid), Vector{Bool}[])
#     _ftmcs(b, f, path, s)
#     result = (b.zero, false)
#     result2 = Vector{Ts}[]
#     for x = s.set
#         tmp = (b.zero, true)
#         tmp2 = Ts[]
#         for i = 1:length(x)
#             if x[i] == true
#                 tmp = bddand(b, tmp, vars[i])
#                 push!(tmp2, vars[i][1].header.label)
#             end
#         end
#         result = bddor(b, result, tmp)
#         push!(result2, tmp2)
#     end
#     return result, result2
# end

# function _ftmcs(b::BDD{Ts}, f::Tuple{Node{Ts},Bool}, path::Vector{Bool}, s::MinPath) where Ts
#     if s.len < sum(path)
#         return
#     end
#     path[f[1].header.level] = false
#     _ftmcs(b, getzero(f), path, s)
#     path[f[1].header.level] = true
#     _ftmcs(b, getone(f), path, s)
#     path[f[1].header.level] = false
#     nothing
# end

# function _ftmcs(b::BDD{Ts}, f::Tuple{Terminal{Ts},Bool}, path::Vector{Bool}, s::MinPath) where Ts
#     if f[2] == true
#         if s.len > sum(path)
#             s.len = sum(path)
#             s.set = [copy(path)]
#         elseif s.len == sum(path)
#             push!(s.set, copy(path))
#         end
#     end
#     nothing
# end

##

function todot(b::BDD2{Ts}, f::Tuple{AbstractBDDNode{Ts},Bool}) where Ts
    io = IOBuffer()
    visited = Set{AbstractBDDNode{Ts}}()
    println(io, "digraph { layout=dot; overlap=false; splines=true; node [fontsize=10];")
    _todot!(b, f[1], visited, io)
    println(io, "}")
    return String(take!(io))
end

function _todot!(b::BDD2{Ts}, f::BDDTerminal{Ts}, visited::Set{AbstractBDDNode{Ts}}, io::IO)::Nothing where Ts
    if in(f, visited)
        return
    end
    println(io, "\"obj$(f.id)\" [shape = square, label = \"0\"];")
    push!(visited, f)
    nothing
end

function _todot!(b::BDD2{Ts}, f::BDD2Node{Ts}, visited::Set{AbstractBDDNode{Ts}}, io::IO)::Nothing where Ts
    if in(f, visited)
        return
    end

    println(io, "\"obj$(f.id)\" [shape = circle, label = \"$(f.header.label)\"];")
    _todot!(b, f.low, visited, io)
    println(io, "\"obj$(f.id)\" -> \"obj$(f.low.id)\" [label = \"0\"];")
    _todot!(b, f.high, visited, io)
    if f.neg == false
        println(io, "\"obj$(f.id)\" -> \"obj$(f.high.id)\" [label = \"1\"];")
    else
        println(io, "\"obj$(f.id)\" -> \"obj$(f.high.id)\" [label = \"1\", arrowhead = odot];")
    end
    push!(visited, f)
    nothing
end


text:
@testset "Subgradient Plan" begin
    M = Euclidean(2)
    x = RnPoint([4., 2.])
    x0 = RnPoint([5.,2.])
    o = SubGradientMethodOptions(x0,stopAfterIteration(200), DecreasingStepsize(0.1))
    o.∂ = RnTVector([1., 0.])
    f = y -> distance(M,y,x)
    ∂f = y -> distance(M,x,y) == 0 ? zeroTVector(M,y) : -2*log(M,y,x)/distance(M,x,y)
    p = SubGradientProblem(M,f, ∂f)
    oR = solve(p,o)
    xHat = getSolverResult(p,oR)
    @test getInitialStepsize(p,o) == 0.1
    @test getStepsize!(p,o,1) == 0.1
    @test getLastStepsize(p,o,1) == 0.1
    # Check Fallbacks of Problen
    @test getCost(p,x) == 0.
    @test norm(M,x,getSubGradient(p,x)) == 0
    @test_throws ErrorException getGradient(p,o.x)
    @test_throws ErrorException getProximalMap(p,1.,o.x,1)
end

ext: jl


text:
# This was run in Julia 1.5.0

using BGEN, OrdinalMultinomialModels, Random, DataFrames, CSV

bgenfile = "bgen_test.bgen"
bgenfile = "/Users/christophergerman/.julia/dev/OrdinalGWAS/data/bgen_test.bgen"
b = Bgen(BGEN.datadir("example.8bits.bgen"))

people = n_samples(b)
nsnps = n_variants(b)
sampleids = samples(b)
variants = parse_variants(b; from_bgen_start=false)
#sample 3 and 15 and 25 to be causal
v3 = minor_allele_dosage!(b, variants[3]; T=Float64, mean_impute=true)
v25 = minor_allele_dosage!(b, variants[25]; T=Float64, mean_impute=true)

Random.seed!(310)
sex = rand(0:1, people)
X = [sex v3 v25]
θ = [-1.5; 0; 1.5]
β = [0.5; 0.5; -0.5]
y = rpolr(X, β, θ, LogitLink())

df = DataFrame(y = y, sex = sex)
CSV.write("bgen_ex.csv", df)

using OrdinalGWAS

b_rsids = rsids(b)
reps1, rem1 = divrem(nsnps, 10)
repeats = [fill(10, reps1); fill(rem1, 1)]
genenames = ["gene" * string(i) for i in 1:length(repeats .!= 0)]
annotations = vcat(map((gname, reps) ->
    fill(gname, reps),
    genenames, repeats)...)
dfannote = DataFrame(gene = annotations, rsids = b_rsids)
CSV.write("bgen_snpsetfile.txt", dfannote, delim = " ",
    header = false)

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
using Test

include("basictests.jl")

include("testAprilGridCalibration.jl")


text:
push!(LOAD_PATH,"../../../../src/");
# # Bilambertian Canopy Scattering

# Using packages:
using Plots, Distributions
using CanopyOptics
theme(:ggplot2)

# Compute quadrature points:
μ,w = CanopyOptics.gauleg(20,0.0,1.0)

# Use standard planophile leaf distribution:
LD = CanopyOptics.planophile_leaves2()

# Use a Bilambertian model and specify Reflectance R and Transmission T
BiLambMod = CanopyOptics.BiLambertianCanopyScattering(R=0.4,T=0.2)

# Compute Scattering Matrices 𝐙⁺⁺, 𝐙⁻⁺ for 0th Fourier moment:
𝐙⁺⁺, 𝐙⁻⁺ = CanopyOptics.compute_Z_matrices(BiLambMod, μ, LD, 0)

# Plots matrices (shows single scattering contributions for incoming and outgoing directions)

l = @layout [a  b]
p1 = contourf(μ, μ, 𝐙⁻⁺, title="Z⁻⁺ (Reflection)", xlabel="μꜜ", ylabel="μꜛ")
p2 = contourf(μ, μ, 𝐙⁺⁺, title="Z⁺⁺ (Transmission)", xlabel="μꜛ", ylabel="μꜛ")
plot(p1, p2,  layout = l, margin=5Plots.mm)
plot!(size=(900,350))

# ### Animation over different Beta leaf distributions
steps = 0.1:0.1:5
α = [collect(steps); 5*ones(length(steps))]
β = [ 5*ones(length(steps)); reverse(collect(steps));]
x = 0:0.01:1

anim = @animate for i ∈ eachindex(α)
    LD = CanopyOptics.LeafDistribution(Beta(α[i],β[i]), 2/π)
    Z⁺⁺, Z⁻⁺ = CanopyOptics.compute_Z_matrices(BiLambMod, μ, LD, 0)
    l = @layout [a  b  c]
    p0 = plot(rad2deg.(π * x/2), pdf.(LD.LD,x), legend=false, ylim=(0,3), title="Leaf angle distribution", xlabel="Θ (degrees)")
    p1 = contourf(μ, μ, Z⁻⁺, title="Z⁻⁺ (Reflection)", xlabel="μꜜ", ylabel="μꜛ",clims=(0.5,2))
    p2 = contourf(μ, μ, Z⁺⁺, title="Z⁺⁺ (Transmission)", xlabel="μꜛ", ylabel="μꜛ",clims=(0.5,2))
    plot(p0, p1, p2,  layout = l, margin=5Plots.mm)
    plot!(size=(1100,300))
end

gif(anim, "anim_fps10.gif", fps = 10)


text:
function FixedStepSelector(arg0::jdouble, arg1::TimeScale)
    return FixedStepSelector((jdouble, TimeScale), arg0, arg1)
end

function select_dates(obj::FixedStepSelector, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "selectDates", List, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
# ---
# title: 430. Flatten a Multilevel Doubly Linked List
# id: problem430
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Depth-first Search
# link: <https://leetcode.com/problems/flatten-a-multilevel-doubly-linked-list/description/>
# hidden: true
# ---
#
# You are given a doubly linked list which in addition to the next and previous
# pointers, it could have a child pointer, which may or may not point to a
# separate doubly linked list. These child lists may have one or more children
# of their own, and so on, to produce a multilevel data structure, as shown in
# the example below.
#
# Flatten the list so that all the nodes appear in a single-level, doubly linked
# list. You are given the head of the first level of the list.
#
#
#
# **Example 1:**
#
#
#
#     Input: head = [1,2,3,4,5,6,null,null,null,7,8,9,10,null,null,11,12]
#     Output: [1,2,3,7,8,11,12,9,10,4,5,6]
#     Explanation:
#     The multilevel linked list in the input is as follows:
#
#     ![](https://assets.leetcode.com/uploads/2018/10/12/multilevellinkedlist.png)
#
#     After flattening the multilevel linked list it becomes:
#
#     ![](https://assets.leetcode.com/uploads/2018/10/12/multilevellinkedlistflattened.png)
#
#
# **Example 2:**
#
#
#
#     Input: head = [1,2,null,3]
#     Output: [1,3,2]
#     Explanation: The input multilevel linked list is as follows:
#
#       1---2---NULL
#       |
#       3---NULL
#
#
# **Example 3:**
#
#
#
#     Input: head = []
#     Output: []
#
#
#
#
# **How  multilevel linked list is represented in test case:**
#
# We use the multilevel linked list from **Example 1** above:
#
#
#
#      1---2---3---4---5---6--NULL
#              |
#              7---8---9---10--NULL
#                  |
#                  11--12--NULL
#
# The serialization of each level is as follows:
#
#
#
#     [1,2,3,4,5,6,null]
#     [7,8,9,10,null]
#     [11,12,null]
#
#
# To serialize all levels together we will add nulls in each level to signify no
# node connects to the upper node of the previous level. The serialization
# becomes:
#
#
#
#     [1,2,3,4,5,6,null]
#     [null,null,7,8,9,10,null]
#     [null,11,12,null]
#
#
# Merging the serialization of each level and removing trailing nulls we obtain:
#
#
#
#     [1,2,3,4,5,6,null,null,null,7,8,9,10,null,null,11,12]
#
#
#
# **Constraints:**
#
#   * Number of Nodes will not exceed 1000.
#   * `1 <= Node.val <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CryptoTools

using Base: split, join

export encrypt, decrypt
export OneTimePad, AES, ChaCha20, ShamirSecretSharing

"""
The default number of byes of security (256-bit).
"""
const DEFAULT_SECURITY = 32

include("Utils.jl")
include("KeyGen.jl")
include("BlockCipher.jl")
include("StreamCipher.jl")
include("SecretSharing.jl")

end


text:
function rotate(rot, phrase)
    alf = string('a':'z'...)
    things = " -!@#%^&*()_+=,.'"*string(0:9...)
    rotated = alf[rot+1:end]*alf[1:rot]
    alf = alf*uppercase(alf)*things
    rotated = rotated*uppercase(rotated)*things
    decoder = Dict(zip(alf,rotated))
    if typeof(phrase) == Char
        return Char(decoder[phrase])
    end
    join([decoder[c] for c in phrase])
end

#=
It works, but it so ugly.
Also it would crash if phrease contains any special sing that i didn't include (like $).
So it's not realy working. It's just not well tetested.
Please help.
    =#

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
function has_any_children(node_element)
    if length(collect(child_elements(node_element))) > 0
        res = true
    else
        res = false
    end
    return res
end

function isnumeric(x)
    res = false
    try
        res = typeof(parse(x)) <: Number
    end
    return res
end

function xml2dict(element_node, replace_newline = nothing)
    if has_any_children(element_node)
        ces = collect(child_elements(element_node))
        dict_res = MultiDict{Any, Any}()

        for c in ces
            childname = name(c)
            value = xml2dict(c, replace_newline)

            if isa(value, ASCIIString) && isnumeric(value)
                value = parse(value)
            elseif !isa(replace_newline, Void) && isa(value, ASCIIString)
               value = replace(value, "\n", replace_newline)
            end
            insert!(dict_res, childname, value)
        end
        return dict_res
    else
        return content(element_node)
    end
end


text:
using Test
using KnetMetrics
using Random

@testset "classification-metrics" begin
    Random.seed!(42)
    y_true = rand(1:3, 100)
    Random.seed!(123)
    y_pred = rand(1:3,100)
    c = confusion_matrix(y_true, y_pred, labels=[1,2,3])

    @testset "class-confusion" begin
        @test class_confusion(c,ith_class =2) == class_confusion(c, class_name = 2)
    end

    @testset "condition-positive" begin
        @test condition_positive(c) == [31, 34, 35]
        @test condition_positive(c, ith_class = 2) == 34
        @test condition_positive(c, class_name = 3) == condition_positive(c, ith_class =3)
    end

    @testset "condition-negative" begin
        @test condition_negative(c) == [69, 66, 65]
        @test condition_negative(c,ith_class = 3) == 65
    end

    @testset "predicted-positive" begin
        @test predicted_positive(c) == [28, 39, 33]
    end

    @testset "predicted-negative" begin
        @test predicted_negative(c) == [72, 61, 67]
    end

    @testset "correctly-classified" begin
        @test correctly_classified(c) == [57, 55, 56]
    end

    @testset "incorrectly-classified" begin
        @test incorrectly_classified(c) == [43, 45, 44]
    end

    @testset "sensitivity-score" begin
        @test isapprox(sensitivity_score(c), 0.33756212162284266)
    end

    @testset "recall-score" begin
        @test isapprox(recall_score(c), 0.33756212162284266)
    end

    @testset "specificity-score" begin
        @test isapprox(specificity_score(c), 0.6694267085571433)
    end

    @testset "precision-score" begin
        @test isapprox(precision_score(c), 0.3361083361083361)
    end

    @testset "accuracy-score" begin
        @test isapprox(accuracy_score(c), 0.3361083361083361)
    end

    @testset "balanced-accuracy-score" begin
        @test isapprox(balanced_accuracy_score(c), 0.33756212162284266)
    end

    @testset "negative-predictive-value" begin
        @test isapprox(negative_predictive_value(c), 0.6698010403356623)
    end

    @testset "false-negative-rate" begin
        @test isapprox(false_negative_rate(c), 0.6624378783771573)
    end

    @testset "false-positive-rate" begin
        @test isapprox(false_positive_rate(c), 0.3305732914428567)
    end

    @testset "fbeta-score" begin
        @test isapprox(fbeta_score(c), 0.33589642032805694)
    end

    @testset "prevalence-threshold" begin
        @test isapprox(prevalence_threshold(c), 0.4988856718477456)
    end

    @testset "threat-threshold" begin
        @test isapprox(threat_score(c), 0.20281219832565797)
    end

    @testset "matthews-correlation-coefficient" begin
        @test isapprox(matthews_correlation_coeff(c), 0.0705714539179035)
    end

    @testset "fowlkes-mallows-index" begin
        @test isapprox(fowlkes_mallows_index(c), 0.8186209097563695)
    end

    @testset "informedness" begin
        @test isapprox(informedness(c), 0.006988830179986009)
    end

    @testset "markedness" begin
        @test isapprox(markedness(c), 0.005909376443998397)
    end

    @testset "cohen-kappa-score" begin
        @test isapprox(cohen_kappa_score(c), 0.007668019846639673)
    end

    @testset "hamming-loss" begin
        @test isapprox(hamming_loss(c), 0.66)
    end

    @testset "jaccard-score" begin
        @test isapprox(jaccard_score(c), 0.20281219832565797)
    end

end


text:
using RunMetrics, Test

cd(@__DIR__)

tol = 1e-10

sample_data = joinpath("..", "sample_data", "sample.tcx")

rs = parse_run(sample_data)
urs = RunMetrics.unit_run_sum(rs)

rt = sum(rs.time)
urt = sum(urs.time)

@test sum(rs.time) ≈ sum(urs.time) atol = tol
@test sum(rs.dist) ≈ sum(urs.dist) atol = tol
@test sum(rs.alt) ≈ sum(urs.alt) atol = tol
@test rs.hr[end] ≈ urs.hr[end] atol = tol


text:
# function exp_settings!(as::ArgParseSettings)
#     @add_arg_table as begin
#         "--exp_loc"
#         help="Location of experiment"
#         arg_type=String
#         default="tmp"
#         "--seed"
#         help="Seed of rng"
#         arg_type=Int64
#         default=0
#         "--steps"
#         help="number of steps"
#         arg_type=Int64
#         default=100
#         "--prev_action_or_not"
#         action=:store_true
#         "--verbose"
#         action=:store_true
#         "--working"
#         action=:store_true
#         "--progress"
#         action=:store_true
#     end
# end

include("util/flux.jl")


text:
module kara_ex01
using Test
include("../src/JuliaKara_noGUI.jl"); using .JuliaKara_noGUI

@testset "JuliaKara Example 01" begin
    wtest = World(10,10)
    kara = place_kara(wtest,4,4)

    function nextTest(world,kara)
        turnRight(world,kara)
        move(world,kara)
        turnLeft(world,kara)
    end

    place_tree(wtest,4,5)
    place_tree(wtest,5,6)
    place_mushroom(wtest,5,5)
    place_mushroom(wtest,6,5)
    place_mushroom(wtest,6,6)
    place_mushroom(wtest,7,5)
    place_leaf(wtest,8,5)

    @test_throws JuliaKara_noGUI.ActorInvalidMovementError move(wtest,kara)
    nextTest(wtest,kara)
    @test_throws JuliaKara_noGUI.ActorInvalidMovementError move(wtest,kara)
    nextTest(wtest,kara)
    @test_throws JuliaKara_noGUI.ActorInvalidMultipleMovementError move(wtest,kara)
    nextTest(wtest,kara)
    move(wtest,kara)
    @test (mushroomFront(wtest,kara)) == true
    turnLeft(wtest,kara)
    turnLeft(wtest,kara)
    move(wtest,kara)
    turnLeft(wtest,kara)
    turnLeft(wtest,kara)
    @test (mushroomFront(wtest,kara)) == false
    nextTest(wtest,kara)
    move(wtest,kara)
    @test (onLeaf(wtest,kara)) == true
    removeLeaf(wtest,kara)
    @test (onLeaf(wtest,kara)) == false
    @test_throws JuliaKara_noGUI.ActorGrabNotFoundError removeLeaf(wtest,kara)
    putLeaf(wtest,kara)
    @test (onLeaf(wtest,kara)) == true
    @test (treeLeft(wtest,kara)) == false
    @test (treeRight(wtest,kara)) == false
    @test (treeFront(wtest,kara)) == false
    place_tree(wtest,7,5)
    place_tree(wtest,9,5)
    place_tree(wtest,8,6)
    @test (treeLeft(wtest,kara)) == true
    @test (treeRight(wtest,kara)) == true
    @test (treeFront(wtest,kara)) == true
    # Test get kara
    @test get_kara(wtest) == kara
    place_kara(wtest,10,10)
    @test length(get_kara(wtest)) == 2
end

@testset "JuliaKara Bounds" begin
    w = World(3,3)
    k = place_kara(w,1,1)
    move(w,k)
    move(w,k)
    move(w,k)
    @test k.location == JuliaKara_noGUI.ActorsWorld.Location(1,1)
    turnRight(w,k)
    turnRight(w,k)
    move(w,k)
    @test k.location == JuliaKara_noGUI.ActorsWorld.Location(1,3)
    turnRight(w,k)
    move(w,k)
    @test k.location == JuliaKara_noGUI.ActorsWorld.Location(3,3)
    turnRight(w,k)
    turnRight(w,k)
    move(w,k)
    @test k.location == JuliaKara_noGUI.ActorsWorld.Location(1,3)
end

end


text:
export HContainer, add_object!

import Base: collect, delete!

"""
`HContainer` is a device for holding a collection of hyperbolic objects.
It is like a set, but we have to do a lot of work before adding a new
element because equal hyperbolic objects might differ a tiny amount and
that would mess up hashing.

+ `C = HContainer()` creates a new container.
+ `C = HContainer(items...)` creates a new container with the items.
"""
struct HContainer
    objs::Set{HObject}
    function HContainer()
        A = Set{HObject}()
        new(A)
    end
end

function HContainer(args...)
    C = HContainer()
    add_object!(C,args...)
    return C
end

HContainer(HC::HContainer) = HContainer(collect(HC.objs)...)  # copy constructor

function in(X::HObject, C::HContainer)::Bool
    for Z in C.objs
        if Z==X
            return true
        end
    end
    return false
end

length(C::HContainer) = length(C.objs)

collect(C::HContainer) = collect(C.objs)

"""
`add_object!(C::HContainer, X::HObject)` adds `X` to the container `C`.
"""
function add_object!(C::HContainer, X::HObject)::Bool
    # see if we already have it
    if in(X,C)
        return false
    end

    # not here, so OK to add
    push!(C.objs, X)
    return true
end

function add_object!(C::HContainer, args...)
    for X in args
        add_object!(C,X)
    end
end

"""
`delete!(C::HContainer, X::HObject)` deletes `X` from the
container `C` returning `true` if successful (or `false` if
`X` was not in the container).
"""
function delete!(C::HContainer, X::HObject)::Bool
    # see if we already have it; if so, delete
    for Z in C.objs
        if Z==X
            delete!(C.objs, X)
            return true
        end
    end
    return false # never found
end

function show(io::IO, C::HContainer)
    print(io,"HContainer of size $(length(C))")
end


text:
import InteractiveUtils
InteractiveUtils.versioninfo(; verbose = true)

pkgall = ["IJulia", "Pluto", "PlutoUI"]
push!(pkgall, "Plots", "JuMP", "GLPK", "SCIP")
push!(pkgall, "CSV", "JSON", "XLSX", "SQLite", "DataFrames")
push!(pkgall, "LightGraphs", "GraphPlot")
push!(pkgall, "BenchmarkTools", "JuliaFormatter", "PackageCompiler")

haskey(ENV, "CPLEX_STUDIO_BINARIES") == false && (ENV["CPLEX_STUDIO_BINARIES"] = "/usr/local/lib")
Sys.ARCH == :x86_64 && push!(pkgall, "CPLEX")

haskey(ENV, "GUROBI_HOME") == false && (ENV["GUROBI_HOME"] = "/usr/local")
Sys.ARCH == :x86_64 && push!(pkgall, "Gurobi")

import Pkg
Pkg.add(pkgall)

import IJulia
IJulia.installkernel("Julia", "--threads=auto", "--sysimage=/usr/local/lib/SuperJulia.sys")

import PackageCompiler
PackageCompiler.create_sysimage(Symbol.(pkgall), sysimage_path = "/usr/local/lib/SuperJulia.sys")


text:
# %% Example files
# Add new examples here (only compiles these examples)
NEW_EXAMPLES = [
]
# Once approved, move examples here:
EXISTING_EXAMPLES = [
    # The following are sorted ALPHABETICALLY!!!
    "Bacteria Growth" => "growing_bacteria.jl",
    "Battle Royale" => "battle.jl",
    "Continuous space social distancing" => "social_distancing.jl",
    "Conway's game of life" => "game_of_life_2D_CA.jl",
    "Daisyworld" => "daisyworld.jl",
    "Forest fire" => "forest_fire.jl",
    "Fractal Growth" => "fractal_growth.jl",
    "Hegselmann-Krause opinion dynamics" => "hk.jl",
    "Maze Solver" => "maze.jl",
    "Mountain Runners" => "runners.jl",
    "Opinion spread" => "opinion_spread.jl",
    "Predator-Prey" => "predator_prey_fast.jl",
    "Sugarscape" => "sugarscape.jl",
    "Wealth distribution" => "wealth_distribution.jl",
    "Wright-Fisher model of evolution" => "wright-fisher.jl",
]

# Documentation building code (not to be touched)
# Load packages
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.instantiate()
println("Loading Packages")
println("Documenter...")
using Documenter
println("Agents...")
using Agents
println("Literate...")
import Literate
println("InteractiveDynamics...")
using InteractiveDynamics

ENV["GKS_ENCODING"] = "utf-8"
println("Converting Examples...")

# download the themes
println("Theme-ing")
using DocumenterTools:Themes
import Downloads
for file in (
    "juliadynamics-lightdefs.scss",
    "juliadynamics-darkdefs.scss",
    "juliadynamics-style.scss",
)
    Downloads.download(
        "https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file",
        joinpath(@__DIR__, file),
    )
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header * "\n" * theme)
end
# compile the themes
Themes.compile(
    joinpath(@__DIR__, "juliadynamics-light.scss"),
    joinpath(@__DIR__, "src/assets/themes/documenter-light.css"),
)
Themes.compile(
    joinpath(@__DIR__, "juliadynamics-dark.scss"),
    joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"),
)

# %% Build examples with Literate.jl
indir = joinpath(@__DIR__, "examples")
outdir = joinpath(@__DIR__, "src", "examples")
rm(outdir; force = true, recursive = true) # cleans up previous examples
mkpath(outdir)

examples = isempty(NEW_EXAMPLES) ? EXISTING_EXAMPLES : NEW_EXAMPLES
built_examples = String[]
for (title, file) in examples
    Literate.markdown(joinpath(indir, file), outdir; credit = false)
    push!(built_examples, "examples/"*file[1:end-3]*".md")
end

# %%
println("Documentation Build")
ENV["JULIA_DEBUG"] = "Documenter"
makedocs(
    modules = [Agents, InteractiveDynamics],
    sitename = "Agents.jl Example Zoo",
    authors = "George Datseris and contributors.",
    doctest = false,
    format = Documenter.HTML(
        prettyurls = CI,
        assets = [
            asset(
                "https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap",
                class = :css,
            ),
        ],
        collapselevel = 2,
    ),
    pages = [
        "Introduction" => "index.md",
        "Examples" => built_examples,
    ],
)

@info "Deploying Documentation"
if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/AgentsExampleZoo.jl.git",
        target = "build",
        push_preview = true,
        devbranch = "main",
    )
end

println("Finished boulding and deploying docs.")


text:
for metric in (
               :Δxᶜᶜᵃ,
               :Δxᶜᶠᵃ,
               :Δxᶠᶠᵃ,
               :Δxᶠᶜᵃ,

               :Δyᶜᶜᵃ,
               :Δyᶜᶠᵃ,
               :Δyᶠᶠᵃ,
               :Δyᶠᶜᵃ,

               :Δzᵃᵃᶠ,
               :Δzᵃᵃᶜ,

               :Azᶜᶜᵃ,
               :Azᶜᶠᵃ,
               :Azᶠᶠᵃ,
               :Azᶠᶜᵃ,

               :Axᶜᶜᶜ,
               :Axᶠᶜᶜ,
               :Axᶠᶠᶜ,
               :Axᶜᶠᶜ,
               :Axᶠᶜᶠ,
               :Axᶜᶜᶠ,

               :Ayᶜᶜᶜ,
               :Ayᶜᶠᶜ,
               :Ayᶠᶜᶜ,
               :Ayᶠᶠᶜ,
               :Ayᶜᶠᶠ,
               :Ayᶜᶜᶠ,

               :Vᶜᶜᶜ,
               :Vᶠᶜᶜ,
               :Vᶜᶠᶜ,
               :Vᶜᶜᶠ,
              )

    @eval begin
        import Oceananigans.Operators: $metric
        @inline $metric(i, j, k, ibg::ImmersedBoundaryGrid) = $metric(i, j, k, ibg.grid)
    end
end


text: using Flux, Plots


text:
import CSV
using DataFramesMeta
import StatsBase: countmap
import Images: label_components
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
mdates = pyimport("matplotlib.dates")

include("./plotting.jl")
# Formatting functions:
_remove_ticks(ax) = ax.tick_params(
    axis = "both",
    which = "both",
    bottom = false,
    top = false,
    # labelbottom=false, labelleft=false,
    left = false,
    right = false,
)
_set_figsize(fig, h = 3, w = 3.5) = (fig.set_figheight(h), fig.set_figwidth(w))

function plot_rrc_oil_water()
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["font.family"] = "Helvetica"
    rcParams["font.size"] = 16
    rcParams["font.weight"] = "bold"

    data = readdf("summary_inj_prod.csv")
    oilyears = data[!, :production]
    wateryears = data[!, :injection]
    xs = data[!, :Year]

    fig, axes = plt.subplots(2, 1)
    # fig2, axes2 = plt.subplots(2, 1)
    bg_color = "b"
    active_color = "r"
    colors = fill(bg_color, length(xs))
    colors[end - 1] = active_color

    xticks = [2007, 2011, 2015, 2018]
    # yticks=[0, 0.5, 1, 1.5, 2, 2.5]

    titles = ["Oil Production", "Total Injections"]
    for (title, ax, values) in zip(titles, axes, [oilyears, wateryears])
        vals = values ./ 365 ./ 1e6
        # ax.bar(xs, vals, color=colors)
        ax.plot(xs, vals, "-o", lw = 5, ms = 10)
        ax.set_xticks(xticks)
        # ax.set_yticks(yticks)
        ax.set_ylabel("Million bbl/day")
        ax.grid()
        ax.set_title(title)
    end
    fig.set_figheight(9)
    fig.set_figwidth(4.5)
    fig.savefig("cisr-lines.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

end

function pnas_outlier_figure(
    station_name = "TXMC",
    outlier_color = "r",
    h = 3,
    w = 3.5,
    bins = 40,
)
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    # https://matplotlib.org/tutorials/introductory/customizing.html#a-sample-matplotlibrc-file
    # so obscure
    rcParams["pdf.fonttype"] = 42
    rcParams["ps.fonttype"] = 42
    rcParams["font.family"] = "Helvetica"
    rcParams["font.size"] = 16
    rcParams["font.weight"] = "bold"

    # a. histogram before/after outlier
    fig, ax = plt.subplots()

    unw_vals = h5read("gps_pixels_78.h5", station_name)
    geolist = Sario.load_geolist_from_h5("gps_pixels_78.h5")
    intlist = Sario.load_intlist_from_h5("gps_pixels_78.h5")
    end_date = Date(2018, 1, 1)
    geolist18 = geolist[geolist .> end_date]
    idxs = InsarTimeseries._good_idxs(geolist18, intlist)
    intlist17 = intlist[idxs]
    geolist17 = geolist[geolist .< end_date]
    unw_vals17 = unw_vals[idxs]
    unw_vals17_cm = p2c .* unw_vals17
    vm = maximum(abs.(unw_vals17_cm))
    n1, _, _ = ax.hist(
        unw_vals17_cm,
        range = (-vm, vm),
        bins = bins,
        label = "All data",
        color = "red",
    )
    MAT.matwrite("outliers_hist_vals.mat", Dict("unw_vals_$station_name" => unw_vals17_cm))
    # n1, _, _ = ax.hist(unw_vals17.* p2c, range=(-vm, vm), bins=bins, label="All data", color="red", density=true, alpha=0.5)
    ax.grid()

    means = InsarTimeseries.mean_abs_val(geolist17, intlist17, unw_vals17)
    nsigma = 4
    low, high = InsarTimeseries.two_way_cutoff(means, nsigma)
    bigidx = means .> high
    idxs = InsarTimeseries._good_idxs(geolist17[bigidx], intlist17)
    vals_after_outliers_cm = unw_vals17_cm[idxs]
    n1, _, _ = ax.hist(
        vals_after_outliers_cm,
        range = (-vm, vm),
        bins = bins,
        label = "Outliers removed",
        color = "C0",
    )
    # n1, _, _ = ax.hist(vals_after_outliers .* p2c, range=(-vm, vm), bins=bins, density=true, label="Outliers removed", color="C0", alpha=.5)
    ax.set_xlim((-15, 15))
    ax.set_xticks([-15, -7.5, 0, 7.5, 15])
    ax.set_yticks([0, 50, 100])
    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("outliers_a.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    # b) std dev of insar at GPS stations vs distance to reference
    dists = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in all_stations
    ]

    txmc_cm_78 = Float64.(p2c .* h5read("gps_pixels_78.h5", "TXMC"))
    @show median(txmc_cm_78)

    stds78_cm =
        Float64.(std.(p2c .* [h5read("gps_pixels_78.h5", name) for name in station_name_list78],))
    dists78 = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in station_name_list78
    ]

    stds85_cm =
        Float64.(std.(p2c .* [h5read("gps_pixels_85.h5", name) for name in station_name_list85],))
    dists85 = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in station_name_list85
    ]
    stations_all = vcat(station_name_list78, station_name_list85)
    dists_all = vcat(dists78, dists85)
    stds_all = vcat(stds78_cm, stds85_cm)
    dists_all = dists_all[stations_all .!= "TXKM"]
    stds_all = stds_all[stations_all .!= "TXKM"]
    stations_all = stations_all[stations_all .!= "TXKM"]

    df_dists = DataFrame(dists = dists_all, stds = stds_all, names = stations_all)

    sm = lm(@formula(stds ~ dists), df_dists)      # Linear model
    sm2 = lm(@formula(stds ~ sqrt(dists)), df_dists) # sqrt model
    # a, b = coef(sm)
    a, b = coef(sm2)
    @show a, b * 50

    fig, ax = plt.subplots()
    dists, stds = (df_dists[!, :dists], df_dists[!, :stds])
    ax.scatter(dists, stds, s = 40)
    xs = 1:maximum(dists)
    # ax.plot(xs, a .+ b .* xs, lw=3)  # Linear model
    ax.plot(xs, a .+ b .* sqrt.(xs), lw = 3)    # sqrt(dist) model
    ax.set_ylim((0, 5))
    ax.set_yticks([0, 1, 2, 3, 4])
    ax.set_xticks([0, 100, 200])
    ax.set_ylabel("Noise Level [cm]")
    ax.set_xlabel("Distance [km]")
    ax.grid()
    # _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("outliers_b.pdf", bbox_inches = "tight", transparent = true, dpi = 100)
    MAT.matwrite(
        "outliers_b.mat",
        Dict("a" => a, "b" => b, "xs" => collect(xs), "dists" => dists, "stds" => stds),
    )

    return sm, sm2
    # c) errors in 2, 3, 4 year solution at GPS vs distance, stacking only
    # d) errors in 2, 3, 4 year solution at GPS vs distance, after outlier removal
    df, errors = get_gps_errors(to_cm = true)
    df = @where(df, :station .!= "TXKM")

    years = collect(2016:2018)
    cols1 = [Symbol("p$(n)_stack_$(y)") for y in years, n in [78, 85]]
    cols2 = [Symbol("p$(n)_outliers_$(y)") for y in years, n in [78, 85]]
    colors = ["b", "g", "k"]
    fnames = ["outliers_c.pdf", "outliers_d.pdf"]

    maxy = 3
    for (idx, cols) in enumerate([cols1, cols2])
        fig, ax = plt.subplots()
        # errvals = abs.(vcat([df[!, c] for c in cols]...))
        for yidx = 1:3
            year = years[yidx]
            c1, c2 = cols[yidx, :]
            errvals = abs.(vcat(df[!, c1], df[!, c2]))
            dists = repeat(df[!, :dist_to_txkm], 2)
            dists = dists[.!isnan.(errvals)]
            errvals = errvals[.!isnan.(errvals)]
            @show maximum(errvals)

            ax.scatter(dists, errvals, s = 40, c = colors[yidx])
            ax.set_ylim((0, maxy))

            df1 = DataFrame(dist = dists, errval = errvals)
            sm = lm(@formula(errval ~ dist), df1)
            a, b = coef(sm)
            @show a, b * 50
            ax.plot(xs, a .+ b .* xs, colors[yidx], lw = 3)

            _remove_ticks(ax)
            _set_figsize(fig)
            fname = fnames[idx]
            fig.savefig(fname, bbox_inches = "tight", transparent = true, dpi = 100)
        end
    end
end

_max_label(labels) = sort(
    [tup for tup in countmap(vec(labels)) if tup[1] != 0],
    by = x->x[2],
    rev = true,
)[1][1]

function save_pnas_images(;
    years = [2016, 2017, 2018],
    fnames = ["velocities_$(year)_current.h5" for year in years],
    # fnames=["velocities_$(year)_linear_max800_sigma4.h5" for year in years],
    outdir = ".",
    # zoomlats=(31.6, 30.9), zoomlons=(-103.8, -103.),  # Zoomed just to stripes
    zoomlats = (31.9, 31.0),
    zoomlons = (-103.8, -102.85),  # Covers txkm/txmh
    latbounds = (32.8, 30.7041),
    dset = "velos_shifted/1",
    cmap1 = "seismic_wide_y",
    cmap2 = "seismic_wide_y", # cmap2=rdylbl,
    vm1 = 10,
    vm2 = 10,
    dem_mask_height = nothing,
)
    maxdays = maximum([_num_days(Sario.load_geolist_from_h5(f, dset)) for f in fnames])
    @show maxdays
    vm1 *= (maxdays / 3650)
    vm2 *= (maxdays / 3650)
    @show vm1, vm2

    m1, m2 = nothing, nothing
    for f in fnames
        @show f
        days = _num_days(Sario.load_geolist_from_h5(f, dset))
        @show days

        m1 = MapImage(f, dset)
        m1[m1 .== 0] .= NaN

        if !isnothing(dem_mask_height)
            dem = Sario.load("elevation_looked.dem")
            labels = label_components(dem .< 1200)
            maxl = _max_label(labels)
            m1[labels .!= maxl] .= NaN
        end

        m1 = m1[latbounds, (-110.0, -90.0)] # no lon bounds really
        @show MapImages.grid_extent(m1)

        @show extremanan(m1)
        m1 .*= (days / 3650)  # Cumulative, in cm
        @show extremanan(m1)
        m2 = m1[zoomlats, zoomlons]
        out1 = joinpath(outdir, replace(f, ".h5" => ".tif"))
        # out2 = replace(f, ".h5"  => "_inset.tif")
        out2 = joinpath(outdir, replace(f, ".h5" => "_inset_gps.tif"))

        # Save as both png and tif
        save_img_geotiff(out1, m1, cmap = cmap1, vm = vm1)
        # save_img_geotiff(out2, m2, cmap=cmap2, vm=vm2)
        plt.imsave(
            replace(out1, ".tif" => ".png"),
            m1,
            cmap = cmap2,
            vmin = -vm2,
            vmax = vm2,
            dpi = 150,
        )
        plt.imsave(
            replace(out2, ".tif" => ".png"),
            m2,
            cmap = cmap2,
            vmin = -vm2,
            vmax = vm2,
            dpi = 150,
        )
    end
    # Now save one with a colorbar
    fig, axes = plt.subplots(1, 2)
    axim = axes[1].imshow(m1, cmap = cmap1, vmax = vm1, vmin = -vm1)
    fig.colorbar(axim, ax = axes[1])
    axim = axes[2].imshow(m2, cmap = cmap2, vmax = vm2, vmin = -vm2)
    fig.colorbar(axim, ax = axes[2])

    cbar_name = "colorbars.png"
    save_paper_figure(fig, cbar_name, false, 400)
end

function _errs_to_cm(df)
    # mm/year, to cm for 2,3,4 years
    scales = vcat(4.1 .* ones(4), 3 .* ones(4), 2 .* ones(4)) ./ 10
    for (s, c) in zip(scales, names(df)[3:end])
        df[!, c] = df[!, c] .* s
    end
    return df
end

function get_gps_errors(; outfile = "gps_errors.csv", to_cm = false)
    if isfile(outfile)
        df = CSV.read(outfile)
        if to_cm
            df = _errs_to_cm(df)
        end
        return df, [c for c in eachcol(df)][3:end]
    end

    p78 = "/data1/scott/pecos/path78-bbox2/igrams_looked/"
    p85 = "/data4/scott/path85/stitched/igrams_looked/"

    dists = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in all_stations
    ]
    errors = Any[]
    push!(errors, all_stations)
    push!(errors, dists)
    # errors85_4 = get_gps_error(p85*"velocities_2018_current.h5", all_stations, dset="velos_shifted/1", verbose=true, window=7, shift=-0.)
    # errors78_4 = get_gps_error(p78*"velocities_2018_current.h5", all_stations, dset="velos_shifted/1", verbose=true, window=7, shift=-0.)

    for f in [
        "velocities_2018_stackavg_max800.h5",
        "velocities_2018_current.h5",
        "velocities_2017_stackavg_max800.h5",
        "velocities_2017_current.h5",
        "velocities_2016_stackavg_max800.h5",
        "velocities_2016_current.h5",
    ]
        fname = p85 * f
        @show fname
        push!(
            errors,
            round.(get_gps_error(
                    fname,
                    all_stations,
                    dset = "velos_shifted/1",
                    verbose = false,
                    window = 7,
                    shift = -0.0,
                ),
                digits = 2,),
        )
        fname = p78 * f
        @show fname
        push!(
            errors,
            round.(get_gps_error(
                    fname,
                    all_stations,
                    dset = "velos_shifted/1",
                    verbose = false,
                    window = 7,
                    shift = -0.0,
                ),
                digits = 2,),
        )
    end
    columns = [
        :station,
        :dist_to_txkm,
        :p85_stack_2018,
        :p78_stack_2018,
        :p85_outliers_2018,
        :p78_outliers_2018,
        :p85_stack_2017,
        :p78_stack_2017,
        :p85_outliers_2017,
        :p78_outliers_2017,
        :p85_stack_2016,
        :p78_stack_2016,
        :p85_outliers_2016,
        :p78_outliers_2016,
    ]

    df = DataFrame(errors, columns)
    CSV.write("gps_errors.csv", df)
    return df, errors
end

function pnas_gps_error_table()
    df, errors = get_gps_errors()
    show(stdout, MIME("text/latex"), df)
    @show extrema.(filter.(!isnan, errors[3:end]))
    @show rms.(errors[3:end])
    @show maximum.(filter.(!isnan, [abs.(a) for a in errors[3:end]]))
    return df
end

function write_gps_pixels(station_names, max_temp = 800)
    geolist, intlist, valid_idxs =
        InsarTimeseries.load_geolist_intlist("unw_stack.h5", "geolist_ignore.txt", max_temp)
    gpss = [
        get_stack_vals("unw_stack.h5", s, 3, "stack_flat_shifted", valid_idxs)
        for s in station_names
    ]
    for (idx, s) in enumerate(station_names)
        h5write("gps_pixels.h5", s, gpss[idx])
    end
    Sario.save_intlist_to_h5("gps_pixels.h5", intlist)
    Sario.save_geolist_to_h5("gps_pixels.h5", geolist)
    Sario.save_dem_to_h5("gps_pixels.h5", Sario.load("dem.rsc"))
end

function pnas_outlier_figure_old(
    station_name = "TXMC",
    outlier_color = "r",
    h = 3,
    w = 3.5,
    bins = 40,
)

    # unw_vals = [h5read("gps_pixels.h5", name) for name in station_name_list78];
    unw_vals = h5read("gps_pixels_78.h5", station_name)
    geolist = Sario.load_geolist_from_h5("gps_pixels_78.h5")
    intlist = Sario.load_intlist_from_h5("gps_pixels_78.h5")

    end_date = Date(2018, 1, 1)
    geolist18 = geolist[geolist .> end_date]
    idxs = InsarTimeseries._good_idxs(geolist18, intlist)
    intlist17 = intlist[idxs]

    geolist17 = geolist[geolist .< end_date]
    unw_vals17 = unw_vals[idxs]
    geolist17, intlist17, unw_vals17 = (geolist, intlist, unw_vals)

    # 1. plot grouped by date
    fig, ax = plot_big_days(
        geolist17,
        intlist17,
        unw_vals17,
        legend = false,
        color = outlier_color,
    )
    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("outliers2.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    # 2. cm vs baseline igram outlier plot
    # function plot_grouped_by_day(geo, int, vals, nsigma=0, min_spread=2)
    geo = geolist17
    means = InsarTimeseries.mean_abs_val(geo, intlist17, p2c * unw_vals17)
    println("median val: $(median(means))")
    fig, ax = plt.subplots()
    ax.scatter(geo, means, label = "means by date")
    nsigma = 4
    low, high = InsarTimeseries.two_way_cutoff(means, nsigma)
    ax.plot(geo, ones(length(geo)) * high, label = "$nsigma sigma MAD cutoff")
    bigidx = means .> high
    ax.scatter(geo[bigidx], means[bigidx], c = outlier_color)
    _remove_ticks(ax)
    _set_figsize(fig)  # Double the width here
    fig.savefig("outliers1.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    # 3. line  plot of insar after outliers vs GPS
    Blin = sum(InsarTimeseries.prepB(geolist17, intlist17), dims = 2)
    before = (Blin \ unw_vals17) * p2c * 365 * 10  # Change to mm/year

    start_date = Date(2015, 1, 1)
    ref = nothing
    lw = 5
    fig, ax = plot_gps_los(
        station_name,
        before;
        ref = ref,
        lw = lw,
        end_date = end_date,
        start_date = start_date,
        ylim = (-3, 3),
        title = "",
        bigfont = false,
        offset = false,
    )

    # Now solve without outliers
    dts, los = get_gps_los(
        station_name,
        reference_station = ref,
        end_date = end_date,
        start_date = start_date,
    )
    day_nums = _get_day_nums(dts)
    idxs = InsarTimeseries._good_idxs(geolist17[bigidx], intlist17)
    vals_after_outliers = unw_vals17[idxs]
    after = p2c * (Blin[idxs] \ vals_after_outliers)  # cm/day
    full_defo = after * (dts[end] - dts[1]).value
    # ax.plot(dts, -full_defo/2 .+ day_nums .* after, "c", lw=3, label="Insar")
    ax.plot(dts, day_nums .* after, "C0", lw = lw)
    _remove_ticks(ax)
    _set_figsize(fig, h, 2.3w)
    fig.savefig("outliers3.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    println("Before: $before mm/yr, $(before * (dts[end] - dts[1]).value / 3650) total defo")
    println("After: $(after * 3650) mm/yr, $full_defo total defo")

    # 4. histogram before/after outlier
    fig, ax = plt.subplots()

    # ax = axes[ii, jj]
    # ax = axes[idx]
    vm = maximum(abs.(unw_vals17 * p2c))
    n1, _, _ = ax.hist(
        unw_vals17 .* p2c,
        range = (-vm, vm),
        bins = bins,
        label = "All data",
        color = "red",
    )
    n1, _, _ = ax.hist(
        vals_after_outliers .* p2c,
        range = (-vm, vm),
        bins = bins,
        label = "Outliers removed",
        color = "C0",
    )
    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("outliers4.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    stds78 =
        Float64.(std.([h5read("gps_pixels_78.h5", name) for name in station_name_list78]))
    dists78 = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in station_name_list78
    ]
    df78 = DataFrame(dists = dists78, stds = stds78, names = station_name_list78)

    stds85 =
        Float64.(std.([h5read("gps_pixels_85.h5", name) for name in station_name_list85]))
    dists85 = [
        MapImages.latlon_to_dist(
            MapImages.station_latlon("TXKM"),
            MapImages.station_latlon(n),
        ) for n in station_name_list85
    ]
    df85 = DataFrame(dists = dists85, stds = stds85, names = station_name_list85)

    sm78 = lm(@formula(stds ~ dists), df78)
    a1, b1 = coef(sm78)
    sm85 = lm(@formula(stds ~ dists), df85)
    a2, b2 = coef(sm85)

    fig, ax = plt.subplots()
    ax.scatter(dists78, stds78, s = 40)
    xs = 1:maximum(dists78)
    ax.plot(xs, a1 .+ b1 .* xs, lw = 3)

    ax.scatter(dists85, stds85, s = 40)
    xs = 1:maximum(dists85)
    ax.plot(xs, a2 .+ b2 .* xs, lw = 3)

    ax.set_xlabel("km from TXKM")
    ax.set_ylabel("Std. Dev. of unw values")
    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("outliers5.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    return (geolist17, intlist17, unw_vals17)
end

# Hunjoo figures
#
function dict_to_vars(d::AbstractDict{String,Any})
    for key in keys(d)
        @eval $(Symbol(key)) = d[$(QuoteNode(key))]
    end
end
function get_wells()
    wolfidx = occursin.("WOLFCAMP", ff[!, prodcol])
    delidx = occursin.("DELAWARE", ff[!, prodcol])
end

function hunjoo_figures()
    rcParams["font.family"] = "Helvetica"
    rcParams["font.size"] = 16
    rcParams["font.weight"] = "medium"

    # variables = MAT.matread("/Users/scott/Documents/Learning/PNAS-git-overleaf/data/PNAS_Geomechanics3.mat")
    # variables = MAT.matread("/Users/scott/Desktop/figures/hunjoo/PNAS_Geomechanics7.mat")
    variables = MAT.matread("/Users/scott/Documents/Learning/GRL-git-overleaf-resubmit/hunjoo-revision-0526/GRL_rev1_Final_b0.01_Scott.mat")

    # dict_to_vars(variables)  # Fails cuz of some namespace thing
    # for key in keys(variables)
    #     @eval $(Symbol(key)) = variables[$(QuoteNode(key))]
    # end
    #  Area of Interest
    R_earth = 3959 * 5280 # (miles to feet) Earth Radius

    XD_min = -103.59
    XD_max = -103.44
    YD_min = 31.32
    YD_max = 31.47

    # #  Reservoir properties
    # RD = variables["Reservoir"][:, 1:2]   # RD = Reservoir(:,1),Reservoir(:,2)];
    # RR_R_D1 = (variables["Reservoir"][:, 3] / (5 * 3280.84) * 55).^2      # Reservoir radius to scatter plot (55^2 = 5 km = 5*3280.84 ft)
    # RR_R_D2 = (variables["Reservoir"][:, 3] / (1 * 3280.84) * 35).^2      # Reservoir radius to scatter plot (35^2 = 1 km = 3280.84 ft)

    # Cities
    Pecos = [-103.4932, 31.4229]             # (degrees)

    # Figure Scale
    ZD_min = -15.0
    ZD_max = 15.0
    c_min = -10.0
    c_max = 10.0
    FZ_min = -4.0
    FZ_max = 2.0

    demrsc = MapImages.from_grid(variables["XD"], variables["YD"])

    # Vertical Deformation
    cmap = "seismic_narrow_y"
    vm = 10

    # Pecos Area InSAR 2,3,1
    # TODO: DRAW TRANSECT FROM B B''
    fig, ax = plt.subplots()
    ax.imshow(variables["Z_Disp_VV18"], cmap=cmap, vmax=vm, vmin=-vm)
    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("hunjoo1.pdf", bbox_inches="tight", transparent=true, dpi=100)

    # Fault model: 2,3,2
    fig, ax = plt.subplots()
    # hold on
    # surf(XD_VV18,YD_VV18,Z_Disp_VV18,"EdgeColor","none")
    axim = ax.imshow(variables[ "Z_Disp_F3" ], cmap=cmap, vmax=vm, vmin=-vm, extent=(XD_min, XD_max,YD_min,YD_max))  #,"EdgeColor","none")   # InSAR data
    fig.colorbar(axim)
    ax.plot(variables[ "BB_XD" ][[1, end]], variables[ "BB_YD" ][[1, end]], lw=1.5)
    ax.plot(Pecos[1],Pecos[2], "ko");
    save_img_geotiff("hunjoo2.tif", MapImage(variables[ "Z_Disp_F3" ], demrsc), cmap = cmap, vm = vm)

    # reservoir model: 2,3,3
    fig, ax = plt.subplots()
    # hold on
    # surf(XD_VV18,YD_VV18,Z_Disp_VV18,"EdgeColor","none")
    cur_im = variables[ "Z_Disp_R3" ][:, :, 3]
    axim = ax.imshow(cur_im, cmap=cmap, vmax=vm, vmin=-vm, extent=(XD_min, XD_max,YD_min,YD_max))  #,"EdgeColor","none")   # InSAR data
    fig.colorbar(axim)
    ax.plot(variables[ "BB_XD" ][[1, end]], variables[ "BB_YD" ][[1, end]], lw=1.5)
    ax.plot(Pecos[1],Pecos[2], "ko")
    save_img_geotiff("hunjoo3.tif", MapImage(cur_im, demrsc), cmap = cmap, vm = vm)

    # Fault + reservoir model: 2,3,4
    fig, ax = plt.subplots()
    # surf(XD,YD,Z_Disp_T(:,:,3),"EdgeColor","none")
    im_mod = variables["Z_Disp_T3"][:, :, 3]
    axim = ax.imshow(
        im_mod,
        cmap = cmap,
        vmax = vm,
        vmin = -vm,
        extent = (XD_min, XD_max, YD_min, YD_max),
    )  # ,"EdgeColor","none")   # InSAR data
    # fig.colorbar(axim)
    ax.plot(variables["BB_XD"][[1, end]], variables["BB_YD"][[1, end]], lw = 1.5)
    ax.plot(Pecos[1], Pecos[2], "ko")

    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("hunjoo4.pdf", bbox_inches = "tight", transparent = true, dpi = 100)
    save_img_geotiff("hunjoo4.tif", MapImage(im_mod, demrsc), cmap = cmap, vm = vm)

    # Residual map: 2,3,5
    fig, ax = plt.subplots()
    # surf(XD,YD,Z_Disp_T(:,:,3),"EdgeColor","none")
    im_mod = variables["Res_Z_VV18"]
    axim = ax.imshow(
        im_mod,
        cmap = cmap,
        vmax = vm,
        vmin = -vm,
        extent = (XD_min, XD_max, YD_min, YD_max),
    )  # ,"EdgeColor","none")   # InSAR data
    # fig.colorbar(axim)
    ax.plot(variables["BB_XD"][[1, end]], variables["BB_YD"][[1, end]], lw = 1.5)
    ax.plot(Pecos[1], Pecos[2], "ko")

    _remove_ticks(ax)
    _set_figsize(fig)
    fig.savefig("hunjoo5.pdf", bbox_inches = "tight", transparent = true, dpi = 100)
    save_img_geotiff("hunjoo5.tif", MapImage(im_mod, demrsc), cmap = cmap, vm = vm)

    # Supplement: figure reservoir figure with circles
    fig, ax = plt.subplots()
    # surf(XD,YD,Z_Disp_T(:,:,3),"EdgeColor","none")
    im_mod = variables["UZ"]
    axim = ax.imshow(
        im_mod,
        cmap = cmap,
        vmax = vm,
        vmin = -vm,
        extent = (XD_min, XD_max, YD_min, YD_max),
    )
    _remove_ticks(ax)
    _set_figsize(fig)
    # fig.savefig("hunjoo-supp1.pdf", bbox_inches = "tight", transparent = true, dpi = 100)
    save_img_geotiff("hunjoo-supp1.tif", MapImage(im_mod, demrsc), cmap = cmap, vm = vm)

    # # B-B' cross-section
    # fig, ax1 = plt.subplots()
    # ax1.plot(BB_L', BB_Z_Disp_I18, "bs", label = "2018 InSAR")
    # ax1.plot(BB_L', BB_Z_Disp_R18, "r-.", lw = 1.5, label = "Reservoir")
    # ax1.plot(BB_L', BB_Z_Disp_F18, "r:", lw = 1.5, label = "Fault")
    # ax1.plot(BB_L', BB_Z_Disp_T18, "r-", lw = 1.5, label = "Res. + Fault")
    # ax1.set_xlim((0, maximum(BB_L)))
    # ax1.set_ylim((2 * FZ_min, 2 * FZ_max))
    # ax1.set_xlabel("B-B' (km)")
    # ax1.set_ylabel("Surface Deformation (cm)")

    # # yyaxis right
    # ax2 = ax1.twinx()
    # ax2.plot(BB_FC[:, 1, 1], BB_FC[:, 2, 1], "g-", lw = 2.0)
    # ax2.plot(BB_FC[:, 1, 2], BB_FC[:, 2, 2], "g-", lw = 2.0)
    # ax2.plot(BB_FC[:, 1, 3], BB_FC[:, 2, 3], "g-", lw = 2.0)
    # ax2.plot(BB_FC[:, 1, 4], BB_FC[:, 2, 4], "g-", lw = 2.0)
    # # text(BB_FC(2,1,1),BB_FC(2,2,1),"F1",lw=2.0)  #,"VerticalAlignment","top")
    # # text(BB_FC(2,1,2),BB_FC(2,2,2),"F2",lw=2.0)  #,"VerticalAlignment","top")
    # # text(BB_FC(2,1,3),BB_FC(2,2,3),"F3",lw=2.0)  #,"VerticalAlignment","top")
    # # text(BB_FC(2,1,4),BB_FC(2,2,4),"F4",lw=2.0)  #,"VerticalAlignment","top")
    # ax2.set_ylabel("Fault Depth (km)")
    # ax2.set_xlim((0, maximum(BB_L)))
    # ax2.set_ylim((0.5 * FZ_min, 0.5 * FZ_max))
    # # ax2.legend(["2018 InSAR","Reservoir Model","Fault Model","Res.+Fault Model","Location","northeast"])
    # # set(findobj(gcf,"type","axes"),"FontSize",12,"TickDir","out")
    # #
    # xticks = [0, 7, 14]
    # y1ticks = [-8, -4, 0, 4]
    # y2ticks = [-2, -1, 0, 1]

    # ax1.set_xticks(xticks)
    # ax1.set_yticks(y1ticks)
    # ax2.set_yticks(y2ticks)
    # ax1.grid()
    # fig.legend(frameon = false, fontsize = "small", ncol = 2)

    # _set_figsize(fig, 1.3 .* [3, 3.5]...)
    # # _set_figsize(fig)
    # fig.savefig("hunjoo4.pdf", bbox_inches = "tight", transparent = true, dpi = 100)
end

# Functions for weighting by covariance

import ImageFiltering: imfilter, Kernel
function get_aps(stack_in)
    return imfilter(stack_in .- imfilter(stack_in, (Kernel.gaussian((0, 0, 1)), Kernel.Laplacian((false, false, true)))),  Kernel.gaussian((2, 2, 0)))
end

function try_weight()
    test_sub = get_aps(defo_stack_unreg);
    var_by_day = dropdims(var((defo_stack_unreg .- test_sub); dims = (1, 2)), dims = (1, 2));
    var_by_day2 = abs.(test_sub[400, 400, :])
    Qddfull2 = build_Q(abs.(var_by_day2), geolist, intlist);
    Qinv2 = inv(Qddfull2);
    W = inv(A' * Qinv * A) * A' * Qinv2;
end

function build_Q(var_dates, geolist, intlist, diag_val = 0.5)
    N = length(intlist);
    Q = zeros((N, N));
    for (jdx, ig2) in enumerate(intlist)
        for (idx, ig1) in enumerate(intlist)
            if jdx > idx
            # if Q[jdx, idx] != 0
                Q[idx, jdx] = Q[jdx, idx]
                continue  # dont need to redo, symmetric
            end
            # @show idx, jdx
            if idx == jdx
                v1 = var_dates[findfirst(isequal(ig1[1]), geolist)]
                v2 = var_dates[findfirst(isequal(ig1[2]), geolist)]
                Q[idx, jdx] = v1 + v2
            else
                cm = countmap([ig1..., ig2...])
                if length(Set([ig1..., ig2...])) == 4
                    Q[idx, jdx] = 0
                    continue
                end
                repeated_date = sort(collect(countmap([ig1..., ig2...])), by = x->x[2], rev = true)[1][1]
                # @show repeated_date
                # @show findfirst(isequal(repeated_date), geolist)
                sigma_repeat = var_dates[findfirst(isequal(repeated_date), geolist)]
                # if jdx in (1, 2) && idx == 28
                    # @show idx, jdx, repeated_date, sigma_repeat
                if (repeated_date == ig1[1] && repeated_date == ig2[1]) || (repeated_date == ig1[2] && repeated_date == ig2[2])
                    Q[idx, jdx] = sigma_repeat
                else
                    Q[idx, jdx] = -sigma_repeat
                end
            end
        end
    end
    return Q + diag_val * I
end

_intlsubset(geolist, intlist, N) = intlist[[(intl[1] in geolist[1:N] && intl[2] in geolist[1:N]) for intl in intlist]]
function _intlsubset(geolist, intlist, vallist::AbstractArray, N)
    idxs = [(intl[1] in geolist[1:N] && intl[2] in geolist[1:N]) for intl in intlist]
    intlist[idxs], vallist[idxs]
end

function Qsubset(var_by_day, geolist, intlist, N)
    g = geolist[1:N]
    v = var_by_day[1:N]
    igl = _intlsubset(geolist, intlist, N)
    return build_Q(v, g, igl)
end

function build_Q_diag(var_dates, geolist, intlist)
    N = length(intlist);
    Q1 = zeros((N,));
    for (idx, ig) in enumerate(intlist)
        v1 = var_dates[findfirst(isequal(ig[1]), geolist)]
        v2 = var_dates[findfirst(isequal(ig[2]), geolist)]
        Q1[idx] = v1 + v2
    end
    return diagm(Q1)
end

import Combinatorics
function demo_weighted(bad_var = 100)
    geolist1 = Date(2020, 1, 1) .+ Dates.Day(12) .* (0:3)
    # 4-element Array{Date,1}:
    #  2020-01-01
    #  2020-01-13
    #  2020-01-25
    #  2020-02-06
    intlist1 = [tuple(p...) for p in Combinatorics.combinations(geolist1, 2)]
    # 6-element Array{Tuple{Date,Date},1}:
    #  (2020-01-01, 2020-01-13)
    #  (2020-01-01, 2020-01-25)
    #  (2020-01-01, 2020-02-06)
    #  (2020-01-13, 2020-01-25)
    #  (2020-01-13, 2020-02-06)
    #  (2020-01-25, 2020-02-06)

    # Real phase is 0, 1, 2, 3
    phi_true = [0, 1, 2, 3]
    A = InsarTimeseries.build_A_matrix(geolist1, intlist1)
    @show meas_nonoise = A * phi_true[2:end]
    @show A \ meas_nonoise

    # But with noisy second day, looks like it uplifted 10
    phi_noise = [0, 10, 2, 3]
    @show meas_noise = A * phi_noise[2:end]
    println("No weighting:")
    @show A \ meas_noise

    variances = [1, bad_var, 1, 1]
    a, b, c, d = variances
    # Q = [a + b a -b; a a + c c; -b c b + c] + 0.1 * I
    Q = build_Q(variances, geolist1, intlist1, .1)
    print("Weight LS:")
    @show inv(A' * inv(Q) * A) * A' * inv(Q) * meas_noise

    println("or just diagonal?")
    Q2 = build_Q_diag(variances, geolist1, intlist1)
    println(Q2)
    print("Weight LS:")
    @show inv(A' * inv(Q2) * A) * A' * inv(Q2) * meas_noise

end

function plot_l1_vs_stack(;offset=true, alpha=300, h=3, w=4.5, yy=2018,
                          unwfile="unw_stack_shiftedonly.h5", station="TXOE", s=station,
                          days_smooth=1)
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    # https://matplotlib.org/tutorials/introductory/customizing.html#a-sample-matplotlibrc-file
    # so obscure
    rcParams["pdf.fonttype"] = 42
    rcParams["ps.fonttype"] = 42
    rcParams["font.family"] = "Helvetica"
    rcParams["font.size"] = 16
    rcParams["font.weight"] = "bold"
    years = mdates.YearLocator()   # every year
    months = mdates.MonthLocator()  # every month
    years_fmt = mdates.DateFormatter("%Y")
    # yrs = (Date(2015), Date(2016), Date(2017), Date(2018))

    demrsc = Sario.load("dem.rsc")
    input_dset = "stack_flat_shifted"
    geolist, intlist, igram_idxs = load_geolist_intlist(unwfile, "geolist_ignore.txt", 800, max_date = Date(yy,1,1))
    Blin = sum(InsarTimeseries.prepB(geolist, intlist), dims = 2)
    ms = 4
    timediffs = InsarTimeseries.day_diffs(geolist)
    println("s = $(s)")
    unw_vals = get_stack_vals(unwfile, gps.station_rowcol(s, Dict(demrsc))..., 1, input_dset, igram_idxs)

    # Timeseries: unregularized, with all outliers (noisiest)
    B = InsarTimeseries.prepB(geolist, intlist, false, 0)
    println("Loaded $(length(unw_vals)) unw_vals, B: $(size(B))")
    unregged = InsarTimeseries.PHASE_TO_CM .* InsarTimeseries.integrate_velocities(B \ unw_vals, timediffs)

    # Timeseries: regularized, but with all outliers
    Ba = InsarTimeseries.prepB(geolist, intlist, false, alpha)
    unw_vals_a = InsarTimeseries.augment_zeros(Ba, unw_vals)
    regged = InsarTimeseries.PHASE_TO_CM .* InsarTimeseries.integrate_velocities(Ba \ unw_vals_a, timediffs)

    #
    # Plot timeseries with-outlier cases:
    fig, ax = plot_gps_los(s; end_date = Date(yy, 1, 1), offset=offset,
                           days_smooth = days_smooth,
                           gps_color=matlab_colors()[4], ms=ms)

    ax.plot(geolist, unregged, "-x", lw=3, label="Regularized")
    ax.plot(geolist, regged, "-x", lw=3, label="Unregularized")
    ax.format_xdata = years_fmt
    ax.xaxis.set_major_locator(years)
    ax.xaxis.set_major_formatter(years_fmt)
    # ax.tick_params(axis="x", labelrotation=45)
    y0 = ceil(maximum(abs.(unregged)))
    ax.set_ylim((-y0, y0))
    ax.set_yticks(-y0:2:y0)
    _set_figsize(fig, h, w)
    fig.legend()
    fig.savefig("compare_$(yy)_timeseries.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    # No outlier removal linear cases
    stack, l2, l1 = prunesolve(geolist, intlist, unw_vals, Blin, 1000, shrink = false)
    println("No outlier, linear: $(stack), $(l1)")
    println("Difference: $(abs(stack - l1))")

    # Plot linear with-outlier cases
    fig, ax = plot_gps_los(s; end_date = Date(yy, 1, 1),
                             insar_mm_list=[l1, l2], offset=offset,
                             labels=["L1 linear", "L2 linear"],
                             gps_color=matlab_colors()[3],
                             insar_colors=matlab_colors()[1:2],
                             days_smooth = days_smooth,
                             ms=ms)

    # ax.plot(geolist, regged, "-x", lw=3, label="reg")
    ax.plot(geolist, unregged, "-x", lw=3, c=matlab_colors()[4], label="Unregularized")
    ax.plot(geolist, regged, "-x", lw=3, c=matlab_colors()[5], label="Regularized")
    ax.format_xdata = years_fmt
    ax.xaxis.set_major_locator(years)
    ax.xaxis.set_major_formatter(years_fmt)
    # ax.tick_params(axis="x", labelrotation=45)
    # y0 = ceil(maximum(abs.(regged)))
    y0 = 6
    ax.set_ylim((-y0, y0))
    ax.set_yticks(-y0:2:y0)
    _set_figsize(fig, h, w)
    fig.legend()
    fig.savefig("compare_$(yy)_linear.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

    ###### Outlier remove cases ################
    # timeseries: regularized with outliers removed
    geo_clean, intlist_clean, unw_clean = remove_outliers(geolist, intlist, unw_vals, mean_sigma_cutoff = 4)
    B2 = InsarTimeseries.prepB(geo_clean, intlist_clean, false, 0)
    td_clean = InsarTimeseries.day_diffs(geo_clean)
    unregged2 = InsarTimeseries.PHASE_TO_CM .* InsarTimeseries.integrate_velocities(B2 \ unw_clean, td_clean)
    Ba2 = InsarTimeseries.prepB(geo_clean, intlist_clean, false, alpha)
    unw_vals_a2 = InsarTimeseries.augment_zeros(Ba2, unw_clean)
    regged2 = InsarTimeseries.PHASE_TO_CM .* InsarTimeseries.integrate_velocities(Ba2 \ unw_vals_a2, td_clean)

    # linear solves with outlier removal
    stack2, l22, l12 = prunesolve(geolist, intlist, unw_vals, Blin, 4, shrink = false)

    # PLOT:
    fig, ax = plot_gps_los(s; end_date = Date(yy, 1, 1),
                             insar_mm_list=[l12, l22], offset=offset,
                             labels=["L1 linear", "L2 linear"],
                             gps_color=matlab_colors()[3],
                             insar_colors=matlab_colors()[1:2],
                             days_smooth = days_smooth,
                             ms=ms)
    # ax.plot(geo_clean, regged2, "-x", lw=3, label="reg")
    ax.plot(geo_clean, unregged2, "-x", c=matlab_colors()[4], lw=3, label="Unregularized")
    ax.plot(geo_clean, regged2, "-x", lw=3, c=matlab_colors()[5], label="Regularized")
    # y0 = ceil(maximum(abs.(regged2)))
    # y0 = 4
    ax.set_ylim((-y0, y0))
    ax.set_yticks(-y0:2:y0)
    ax.format_xdata = years_fmt
    ax.xaxis.set_major_locator(years)
    ax.xaxis.set_major_formatter(years_fmt)
    _set_figsize(fig, h, w)
    fig.legend()
    fig.savefig("compare_$(yy)_removed.pdf", bbox_inches = "tight", transparent = true, dpi = 100)

end


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots
using DelimitedFiles
using Plots.Measures

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if titleArg == "3 Isoforms" && path == "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    elseif titleArg == "Real"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,left_margin = 50mm,
        right_margin = 50mm)
    end
    #vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end

#for countsMatrix in scnorm_counts
# "Weibull", "Random", "UniformObserved", "CellVariable"
Weibull_1 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Weibull_real = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_real_isoforms_array.txt", "Real")

Random_1 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Random_real = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_real_isoforms_array.txt", "Real")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
UniformObserved_real = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_real_isoforms_array.txt", "Real")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
CellVariable_real = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_real_isoforms_array.txt", "Real")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       Weibull_real, Random_real,UniformObserved_real, CellVariable_real,
       layout = (5,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure3.png")

###############################################################
using JLD2
using DelimitedFiles
using HypothesisTests
using DataFrames

Weibull_1 = readdlm("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
Weibull_2 = readdlm("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_2_isoforms_array.txt")
Weibull_3 = readdlm("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_3_isoforms_array.txt")
Weibull_4 = readdlm("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_4_isoforms_array.txt")
Weibull_real = readdlm("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_real_isoforms_array.txt")

Random_1 = readdlm("figures/figure5_data/Random_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
Random_2 = readdlm("figures/figure5_data/Random_H9_24_scnorm_counts.txt_2_isoforms_array.txt")
Random_3 = readdlm("figures/figure5_data/Random_H9_24_scnorm_counts.txt_3_isoforms_array.txt")
Random_4 = readdlm("figures/figure5_data/Random_H9_24_scnorm_counts.txt_4_isoforms_array.txt")
Random_real = readdlm("figures/figure5_data/Random_H9_24_scnorm_counts.txt_real_isoforms_array.txt")

UniformObserved_1 = readdlm("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
UniformObserved_2 = readdlm("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_2_isoforms_array.txt")
UniformObserved_3 = readdlm("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_3_isoforms_array.txt")
UniformObserved_4 = readdlm("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_4_isoforms_array.txt")
UniformObserved_real = readdlm("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_real_isoforms_array.txt")

CellVariable_1 = readdlm("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
CellVariable_2 = readdlm("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_2_isoforms_array.txt")
CellVariable_3 = readdlm("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_3_isoforms_array.txt")
CellVariable_4 = readdlm("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_4_isoforms_array.txt")
CellVariable_real = readdlm("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_real_isoforms_array.txt")

df = DataFrame(IsoformsSimulated=Int64[], Comparison = String[], pValue = Float64[])

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["All"], pValue=[pvalue(KSampleADTest(vec(Weibull_1), vec(Random_1), vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_2), vec(Random_2), vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_3), vec(Random_3), vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_4), vec(Random_4), vec(UniformObserved_4), vec(CellVariable_4)))]))

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest( vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_4), vec(CellVariable_4)))]))

jldopen("figures/figure4_data/Supfig3_AD_test_results.jld2", "w") do file
    file["df"] = df
end

###########################################
# overlap
#############################################

function make_hist(path, titleArg)

     Arr = readdlm(path)
     p1 = histogram(Arr, title = titleArg, normalize=:pdf, bins = 100,
     ylims = (0,4), color = palette(:default)[2],
     xlims = (0,1),
     axis = (font(150)), titlefontsize=150, yticks = 0:2:4,
     xticks = 0:0.5:1,
     margin = 35mm,left_margin = 50mm,
     right_margin = 50mm)
     return p1

end

Weibull_1 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H9_24_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

Random_1 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H9_24_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H9_24_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H9_24_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       layout = (4,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure3_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
# include the package
using LatticePhysics

# obtain the unitcell of the honeycomb lattice
unitcell = getUnitcellHoneycomb()

# build a lattice for visual confirmation that this is the correct unitcell
lattice = getLatticeInBox(unitcell, [10.0, 8.0])
# plot the lattice for visual confirmation
plotLattice2D(lattice, openfile=true)


text:
include(joinpath(@__DIR__, "objective.jl"))
include(joinpath(@__DIR__, "constraints.jl"))
include(joinpath(@__DIR__, "data.jl"))
include(joinpath(@__DIR__, "rollout.jl"))
include(joinpath(@__DIR__, "augmented_lagrangian.jl"))
include(joinpath(@__DIR__, "derivatives.jl"))
include(joinpath(@__DIR__, "backward_pass.jl"))
include(joinpath(@__DIR__, "forward_pass.jl"))
include(joinpath(@__DIR__, "solve.jl"))


text:
const GAEVF = JuMP.GenericAffExpr{Float64, VariableRef}
const GQEVF = JuMP.GenericQuadExpr{Float64, VariableRef}

function moi_tests(op_model::OperationModel,
                   params::Bool,
                   vars::Int64,
                   interval::Int64,
                   lessthan::Int64,
                   greaterthan::Int64,
                   equalto::Int64,
                   binary::Bool)

    JuMPmodel = op_model.canonical.JuMPmodel
    @test (:params in keys(JuMPmodel.ext)) == params
    @test JuMP.num_variables(JuMPmodel) == vars
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.Interval{Float64}) == interval
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.LessThan{Float64}) == lessthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.GreaterThan{Float64}) == greaterthan
    @test JuMP.num_constraints(JuMPmodel, GAEVF, MOI.EqualTo{Float64}) == equalto
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(JuMPmodel)) == binary

    return

end

function psi_constraint_test(op_model::OperationModel, constraint_names::Vector{Symbol})

    for con in constraint_names
        @test !isnothing(get(op_model.canonical.constraints, con, nothing))
    end

    return

end

function psi_checkbinvar_test(op_model::OperationModel, bin_variable_names::Vector{Symbol})

    for variable in bin_variable_names
        for v in op_model.canonical.variables[variable]
            @test JuMP.is_binary(v)
        end
    end

    return

end

function psi_checkobjfun_test(op_model::OperationModel, exp_type)

    @test JuMP.objective_function_type(op_model.canonical.JuMPmodel) == exp_type

    return

end

function moi_lbvalue_test(op_model::OperationModel, con_name::Symbol, value::Number)

    for con in op_model.canonical.constraints[con_name]
        @test JuMP.constraint_object(con).set.lower == value
    end

    return

end

function psi_checksolve_test(op_model::OperationModel, status)
    JuMP.optimize!(op_model.canonical.JuMPmodel)
    @test termination_status(op_model.canonical.JuMPmodel) in status
end

function psi_checksolve_test(op_model::OperationModel, status, expected_result, tol = 0.0)
    res = solve_op_model!(op_model)
    @test termination_status(op_model.canonical.JuMPmodel) in status
    @test isapprox(res.total_cost[:OBJECTIVE_FUNCTION], expected_result, atol = tol)
end


text:
using OffsetArrays
using Test
using DelimitedFiles
using OffsetArrays: IdentityUnitRange

@test isempty(detect_ambiguities(OffsetArrays, Base, Core))

# Basics
for n = 0:5
    for z in (OffsetArray(ones(Int,ntuple(d->1,n)), ntuple(x->x-1,n)),
              fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)), 1),
              fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)...), 1),
              fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)), 1),
              fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)...), 1))
        @test length(LinearIndices(z)) == 1
        @test axes(z) == ntuple(x->x:x, n)
        @test z[1] == 1
    end
end
a0 = reshape([3])
a = OffsetArray(a0)
@test axes(a) == ()
@test ndims(a) == 0
@test a[] == 3

# missing and nothing constructors
for (T, t) in ((Missing, missing), (Nothing, nothing))
    @test !isassigned(OffsetArray{Union{T,Vector{Int}}}(undef, -1:1, -1:1), -1, -1)
    @test OffsetArray{Union{T,Vector{Int}}}(t, -1:1, -1:1)[-1, -1] === t
    @test !isassigned(OffsetVector{Union{T,Vector{Int}}}(undef, -1:1), -1)
    @test OffsetVector{Union{T,Vector{Int}}}(t, -1:1)[-1] === t
end

y = OffsetArray{Float64}(undef, -1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
@test axes(y) == (-1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
y[-1,-7,-128,-5,-1,-3,-2,-1] = 14
y[-1,-7,-128,-5,-1,-3,-2,-1] += 5
@test y[-1,-7,-128,-5,-1,-3,-2,-1] == 19

r = -2:5
y = OffsetArray(r, r)
@test axes(y) == (r,)
y = OffsetArray(r, (r,))
@test axes(y) == (r,)

y = OffsetArray{Float32}(undef, (IdentityUnitRange(-1:1),))
@test axes(y) === (IdentityUnitRange(-1:1),)

A0 = [1 3; 2 4]
A = OffsetArray(A0, (-1,2))                   # IndexLinear
S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))   # IndexCartesian
@test axes(A) == axes(S) == (0:1, 3:4)
@test size(A) == size(A0)
@test size(A, 1) == size(A0, 1)
@test length(A) == length(A0)
@test A == OffsetArray(A0, 0:1, 3:4)
@test_throws DimensionMismatch OffsetArray(A0, 0:2, 3:4)
@test_throws DimensionMismatch OffsetArray(A0, 0:1, 2:4)

# Scalar indexing
@test @inferred(A[0,3]) == @inferred(A[0,3,1]) == @inferred(A[1]) == @inferred(S[0,3]) == @inferred(S[0,3,1]) == @inferred(S[1]) == 1
@test A[1,3] == A[1,3,1] == A[2] == S[1,3] == S[1,3,1] == S[2] == 2
@test A[0,4] == A[0,4,1] == A[3] == S[0,4] == S[0,4,1] == S[3] == 3
@test A[1,4] == A[1,4,1] == A[4] == S[1,4] == S[1,4,1] == S[4] == 4
@test @inbounds(A[0,3]) == @inbounds(A[0,3,1]) == @inbounds(A[1]) == @inbounds(S[0,3]) == @inbounds(S[0,3,1]) == @inbounds(S[1]) == 1
@test @inbounds(A[1,3]) == @inbounds(A[1,3,1]) == @inbounds(A[2]) == @inbounds(S[1,3]) == @inbounds(S[1,3,1]) == @inbounds(S[2]) == 2
@test @inbounds(A[0,4]) == @inbounds(A[0,4,1]) == @inbounds(A[3]) == @inbounds(S[0,4]) == @inbounds(S[0,4,1]) == @inbounds(S[3]) == 3
@test @inbounds(A[1,4]) == @inbounds(A[1,4,1]) == @inbounds(A[4]) == @inbounds(S[1,4]) == @inbounds(S[1,4,1]) == @inbounds(S[4]) == 4
@test_throws BoundsError A[1,1]
@test_throws BoundsError S[1,1]
@test_throws BoundsError A[0,3,2]
@test_throws BoundsError A[0,3,0]
Ac = copy(A)
Ac[0,3] = 10
@test Ac[0,3] == 10
Ac[0,3,1] = 11
@test Ac[0,3] == 11
@inbounds Ac[0,3,1] = 12
@test Ac[0,3] == 12

# Vector indexing
@test A[:, 3] == S[:, 3] == OffsetArray([1,2], (A.offsets[1],))
@test A[:, 4] == S[:, 4] == OffsetArray([3,4], (A.offsets[1],))
@test_throws BoundsError A[:, 1]
@test_throws BoundsError S[:, 1]
@test A[0, :] == S[0, :] == OffsetArray([1,3], (A.offsets[2],))
@test A[1, :] == S[1, :] == OffsetArray([2,4], (A.offsets[2],))
@test_throws BoundsError A[2, :]
@test_throws BoundsError S[2, :]
@test A[0:1, 3] == S[0:1, 3] == [1,2]
@test A[[1,0], 3] == S[[1,0], 3] == [2,1]
@test A[0, 3:4] == S[0, 3:4] == [1,3]
@test A[1, [4,3]] == S[1, [4,3]] == [4,2]
@test A[:, :] == S[:, :] == A

# CartesianIndexing
@test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
@test A[CartesianIndex((0,3)),1] == S[CartesianIndex((0,3)),1] == 1
@test @inbounds(A[CartesianIndex((0,3))]) == @inbounds(S[CartesianIndex((0,3))]) == 1
@test @inbounds(A[CartesianIndex((0,3)),1]) == @inbounds(S[CartesianIndex((0,3)),1]) == 1
@test_throws BoundsError A[CartesianIndex(1,1)]
@test_throws BoundsError A[CartesianIndex(1,1),0]
@test_throws BoundsError A[CartesianIndex(1,1),2]
@test_throws BoundsError S[CartesianIndex(1,1)]
@test_throws BoundsError S[CartesianIndex(1,1),0]
@test_throws BoundsError S[CartesianIndex(1,1),2]
@test eachindex(A) == 1:4
@test eachindex(S) == CartesianIndices(IdentityUnitRange.((0:1,3:4)))

# view
S = view(A, :, 3)
@test S == OffsetArray([1,2], (A.offsets[1],))
@test S[0] == 1
@test S[1] == 2
@test_throws BoundsError S[2]
@test axes(S) === (IdentityUnitRange(0:1),)
S = view(A, 0, :)
@test S == OffsetArray([1,3], (A.offsets[2],))
@test S[3] == 1
@test S[4] == 3
@test_throws BoundsError S[1]
@test axes(S) === (IdentityUnitRange(3:4),)
S = view(A, 0:0, 4)
@test S == [3]
@test S[1] == 3
@test_throws BoundsError S[0]
@test axes(S) === (Base.OneTo(1),)
S = view(A, 1, 3:4)
@test S == [2,4]
@test S[1] == 2
@test S[2] == 4
@test_throws BoundsError S[3]
@test axes(S) === (Base.OneTo(2),)
S = view(A, :, :)
@test S == A
@test S[0,3] == S[1] == 1
@test S[1,3] == S[2] == 2
@test S[0,4] == S[3] == 3
@test S[1,4] == S[4] == 4
@test_throws BoundsError S[1,1]
@test axes(S) === IdentityUnitRange.((0:1, 3:4))

# iteration
let a
    for (a,d) in zip(A, A0)
        @test a == d
    end
end

# show
@test sprint(show, A) == "[1 3; 2 4]"
@test sprint(show, S) == "[1 3; 2 4]"
strs = split(strip(sprint(show, MIME("text/plain"), A)), '\n')
@test strs[2] == " 1  3"
@test strs[3] == " 2  4"
v = OffsetArray(rand(3), (-2,))
@test sprint(show, v) == sprint(show, parent(v))
io = IOBuffer()
function cmp_showf(printfunc, io, A)
    ioc = IOContext(io, :limit=>true, :compact=>true)
    printfunc(ioc, A)
    str1 = String(take!(io))
    printfunc(ioc, parent(A))
    str2 = String(take!(io))
    @test str1 == str2
end
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,5), (10,-9)))       # rows&cols fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,5), (10,-9)))    # columns fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,10^3), (10,-9)))    # rows fit
cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neither fits

# Similar
B = similar(A, Float32)
@test isa(B, OffsetArray{Float32,2})
@test axes(B) === axes(A)
B = similar(A, (3,4))
@test isa(B, Array{Int,2})
@test size(B) == (3,4)
@test axes(B) === (Base.OneTo(3), Base.OneTo(4))
B = similar(A, (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test axes(B) === IdentityUnitRange.((-3:3, 1:4))
B = similar(parent(A), (-3:3,1:4))
@test isa(B, OffsetArray{Int,2})
@test axes(B) === IdentityUnitRange.((-3:3, 1:4))
@test isa([x for x in [1,2,3]], Vector{Int})
@test similar(Array{Int}, (0:0, 0:0)) isa OffsetArray{Int, 2}
@test similar(Array{Int}, (1, 1)) isa Matrix{Int}
@test similar(Array{Int}, (Base.OneTo(1), Base.OneTo(1))) isa Matrix{Int}

# Reshape
B = reshape(A0, -10:-9, 9:10)
@test isa(B, OffsetArray{Int,2})
@test parent(B) === A0
@test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
B = reshape(A, -10:-9, 9:10)
@test isa(B, OffsetArray{Int,2})
@test pointer(parent(B)) === pointer(A0)
@test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
b = reshape(A, -7:-4)
@test axes(b) == (IdentityUnitRange(-7:-4),)
@test isa(parent(b), Vector{Int})
@test pointer(parent(b)) === pointer(parent(A))
@test parent(b) == A0[:]
a = OffsetArray(rand(3,3,3), -1:1, 0:2, 3:5)
# Offset axes are required for reshape(::OffsetArray, ::Val) support
b = reshape(a, Val(2))
@test isa(b, OffsetArray{Float64,2})
@test pointer(parent(b)) === pointer(parent(a))
@test axes(b) == IdentityUnitRange.((-1:1, 1:9))
b = reshape(a, Val(4))
@test isa(b, OffsetArray{Float64,4})
@test pointer(parent(b)) === pointer(parent(a))
@test axes(b) == (axes(a)..., IdentityUnitRange(1:1))

# Indexing with OffsetArray axes
i1 = OffsetArray([2,1], (-5,))
i1 = OffsetArray([2,1], -5)
b = A0[i1, 1]
@test axes(b) === (IdentityUnitRange(-4:-3),)
@test b[-4] == 2
@test b[-3] == 1
b = A0[1,i1]
@test axes(b) === (IdentityUnitRange(-4:-3),)
@test b[-4] == 3
@test b[-3] == 1
v = view(A0, i1, 1)
@test axes(v) === (IdentityUnitRange(-4:-3),)
v = view(A0, 1:1, i1)
@test axes(v) === (Base.OneTo(1), IdentityUnitRange(-4:-3))

# logical indexing
@test A[A .> 2] == [3,4]

# copyto!
a = OffsetArray{Int}(undef, (-3:-1,))
fill!(a, -1)
copyto!(a, (1,2))   # non-array iterables
@test a[-3] == 1
@test a[-2] == 2
@test a[-1] == -1
fill!(a, -1)
copyto!(a, -2, (1,2))
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copyto!(a, 1, (1,2))
fill!(a, -1)
copyto!(a, -2, (1,2,3), 2)
@test a[-3] == -1
@test a[-2] == 2
@test a[-1] == 3
@test_throws BoundsError copyto!(a, -2, (1,2,3), 1)
fill!(a, -1)
copyto!(a, -2, (1,2,3), 1, 2)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2

b = 1:2    # copy between AbstractArrays
bo = OffsetArray(1:2, (-3,))
@test_throws BoundsError copyto!(a, b)
fill!(a, -1)
copyto!(a, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
fill!(a, -1)
copyto!(a, -2, bo)
@test a[-3] == -1
@test a[-2] == 1
@test a[-1] == 2
@test_throws BoundsError copyto!(a, -4, bo)
@test_throws BoundsError copyto!(a, -1, bo)
fill!(a, -1)
copyto!(a, -3, b, 2)
@test a[-3] == 2
@test a[-2] == a[-1] == -1
@test_throws BoundsError copyto!(a, -3, b, 1, 4)
am = OffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
fill!(am, -1)
copyto!(am, b)
@test am[1] == 1
@test am[2] == 2
@test am[3] == -1
@test am[1,7] == 1
@test am[1,8] == 2
@test am[1,9] == -1

dest = similar(am)
map!(+, dest, am, am)
@test dest[1,7] == 2
@test dest[1,8] == 4
@test dest[1,9] == -2

A = OffsetArray(rand(4,4), (-3,5))
@test maximum(A) == maximum(parent(A))
@test minimum(A) == minimum(parent(A))
@test extrema(A) == extrema(parent(A))
C = similar(A)
cumsum!(C, A, dims = 1)
@test parent(C) == cumsum(parent(A), dims = 1)
@test parent(cumsum(A, dims = 1)) == cumsum(parent(A), dims = 1)
cumsum!(C, A, dims = 2)
@test parent(C) == cumsum(parent(A), dims = 2)
R = similar(A, (1:1, 6:9))
maximum!(R, A)
@test parent(R) == maximum(parent(A), dims = 1)
R = similar(A, (-2:1, 1:1))
maximum!(R, A)
@test parent(R) == maximum(parent(A), dims = 2)
amin, iamin = findmin(A)
pmin, ipmin = findmin(parent(A))
@test amin == pmin
@test A[iamin] == amin
@test amin == parent(A)[ipmin]
amax, iamax = findmax(A)
pmax, ipmax = findmax(parent(A))
@test amax == pmax
@test A[iamax] == amax
@test amax == parent(A)[ipmax]

v  = OffsetArray([1,1e100,1,-1e100], (-3,))*1000
v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000
@test isa(v, OffsetArray)
cv  = OffsetArray([1,1e100,1e100,2], (-3,))*1000
cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000
# @test isequal(cumsum_kbn(v), cv)
# @test isequal(cumsum_kbn(v2), cv2)
# @test isequal(sum_kbn(v), sum_kbn(parent(v)))

io = IOBuffer()
writedlm(io, A)
seek(io, 0)
@test readdlm(io, eltype(A)) == parent(A)

amin, amax = extrema(parent(A))
@test clamp.(A, (amax+amin)/2, amax) == OffsetArray(clamp.(parent(A), (amax+amin)/2, amax), axes(A))

@test unique(A, dims=1) == OffsetArray(parent(A), 0, first(axes(A, 2)) - 1)
@test unique(A, dims=2) == OffsetArray(parent(A), first(axes(A, 1)) - 1, 0)
v = OffsetArray(rand(8), (-2,))
@test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
@test sortslices(A; dims=1) == OffsetArray(sortslices(parent(A); dims=1), A.offsets)
@test sortslices(A; dims=2) == OffsetArray(sortslices(parent(A); dims=2), A.offsets)
@test sort(A, dims = 1) == OffsetArray(sort(parent(A), dims = 1), A.offsets)
@test sort(A, dims = 2) == OffsetArray(sort(parent(A), dims = 2), A.offsets)

@test mapslices(v->sort(v), A, dims = 1) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 1), A.offsets)
@test mapslices(v->sort(v), A, dims = 2) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 2), A.offsets)

@test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]])
@test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]])
@test reverse(A, dims = 1) == OffsetArray(reverse(parent(A), dims = 1), A.offsets)
@test reverse(A, dims = 2) == OffsetArray(reverse(parent(A), dims = 2), A.offsets)

@test A.+1 == OffsetArray(parent(A).+1, A.offsets)
@test 2*A == OffsetArray(2*parent(A), A.offsets)
@test A+A == OffsetArray(parent(A)+parent(A), A.offsets)
@test A.*A == OffsetArray(parent(A).*parent(A), A.offsets)

B = fill(5, 1:3, -1:1)
@test axes(B) == (1:3,-1:1)
@test all(B.==5)

# @inbounds
a = OffsetArray(zeros(7), -3:3)
unsafe_fill!(x) = @inbounds(for i in axes(x,1); x[i] = i; end)
function unsafe_sum(x)
    s = zero(eltype(x))
    @inbounds for i in axes(x,1)
        s += x[i]
    end
    s
end
unsafe_fill!(a)
for i = -3:3
    @test a[i] == i
end
@test unsafe_sum(a) == 0

a = OffsetArray([1 2; 3 4], -1:0, 5:6)
@test summary(a) == "OffsetArray(::Array{$(Int),2}, -1:0, 5:6) with eltype $(Int) with indices -1:0×5:6"
@test summary(view(a, :, 5)) == "view(OffsetArray(::Array{$(Int),2}, -1:0, 5:6), :, 5) with eltype $(Int) with indices -1:0"
a = OffsetArray(reshape([1]))
@test summary(a) == "0-dimensional OffsetArray(::Array{$(Int),0}) with eltype $(Int)"

@testset "Resizing OffsetVectors" begin
    local a = OffsetVector(rand(5),-3)
    axes(a,1) == -2:2
    length(a) == 5
    resize!(a,3)
    length(a) == 3
    axes(a,1) == -2:0
    @test_throws ArgumentError resize!(a,-3)
end

@testset "OffsetVector constructors" begin
    local v = rand(5)
    @test OffsetVector(v, -2) == OffsetArray(v, -2)
    @test OffsetVector(v, -2:2) == OffsetArray(v, -2:2)
    @test typeof(OffsetVector{Float64}(undef, -2:2)) == typeof(OffsetArray{Float64}(undef, -2:2))
end


text:
export
    matches_keywords,
    get_files_matching

"""
    matches_keywords(p,keywords)

Check if any of the keywords occurs in p.
"""
function matches_keywords(p,keywords)
    if isempty(keywords)
        return true
    end
    for k in keywords
        if findfirst(k,p) != nothing
            return true
        end
    end
    return false
end

"""
    get_files_matching(base_path,ext,keywords=[])

Get all files in a directory that match the extension `ext`, and (optionally)
which contain any of `keywords`
"""
function get_files_matching(base_path,ext,keywords=[])
    paths = String[]
    for (root,dirs,files) in walkdir(base_path)
        for f in files
            if matches_keywords(f,keywords) && splitext(f)[end] == ext
                push!(paths,joinpath(root,f))
            end
        end
    end
    return paths
end

export redirect_to_files

"""
	redirect_to_files(dofunc, outfile, errfile)

redirects output of `stdout` and `stderr` to `outfile` and `errfile`,
respectively.
Usage:
	```
	redirect_to_files(prefix * ".log", prefix * ".err") do
    	compute(...)
	end
	```
"""
function redirect_to_files(dofunc, outfile, errfile)
    open(outfile, "w") do out
        open(errfile, "w") do err
            redirect_stdout(out) do
                redirect_stderr(err) do
                    dofunc()
                end
            end
        end
    end
end


text:
using BasisFunctions
using MoreBasisFunctions
using Test

include("vandermonde_tests.jl")
include("bernstein_tests.jl")
include("lagrange_tests.jl")