text:
using Plots

using NetCDF: ncread, ncwrite, nccreate
using Statistics: mean

# data folders
data_folder = DATA_FOLDER
combined_file = COMBINE_FILE
analysis_folder = ANALYSIS_FOLDER

# check if combined file extist, otherwise create it
if !isfile(combined_file)
  ## get all the filenames in the data directory
  files = sort(readdir(data_folder, join=true))
  files = [file for file in files if occursin(".nc", file)]
  files = [file for file in files if !occursin("AtmosCore", file)]

  ## extract grid
  z    = ncread(files[1], "z")
  time = [0.0]
  for file in files
    push!(time, ncread(file, "simtime")[1])
  end
  grid = Dict("time" => time[2:end], "z" => z)

  ## define variables of interest to aggregate
  vars = ["w", "qt", "ql", "cld_frac", "thl", "ht", "var_w", "w3", "var_qt", "var_thl", "cov_w_qt", "cov_w_thl", "w_qt_sgs", "w_ht_sgs", "cld_cover"]

  ## read and aggregate variables of interest from files
  data = Dict(var => [] for var in vars)
  for file in files
    for var in vars
      temp = ncread(file, var)
      temp = reshape(temp, (1, size(temp)...))
      push!(data[var], temp)
    end
  end
  data = Dict(var => vcat(data[var]...) for var in vars)

  ## write data to disk as combined file
  for var in vars
    nccreate(
      combined_file,
      var,
      "time", grid["time"],
      "z", grid["z"],
    )
    ncwrite(data[var], combined_file, var)
  end
end

# load grid
time = ncread(combined_file, "time")
z = ncread(combined_file, "z")

# load data
vars = ["w", "qt", "ql", "cld_frac", "thl", "ht", "var_w", "w3", "var_qt", "var_thl", "cov_w_qt", "cov_w_thl", "w_qt_sgs", "w_ht_sgs", "cld_cover"]
data = Dict(var => ncread(combined_file, var) for var in vars)
println("raw data dimensions (time, height): ", size(data["w"]))

# create time mean data
avg_start_ind = 60 # for removing the spinup time
data_time_mean = Dict(k => mean(v[avg_start_ind:end, :], dims=1)[end, :] for (k, v) in data)
println("time-averaged data dimensions (height): ", size(data_time_mean["w"]))

## Plot time means
p1 = plot(z, data_time_mean["w"], xlabel="w_time_mean", ylabel="z (km)")
p2 = plot(z, data_time_mean["qt"], xlabel="qt_time_mean", ylabel="z (km)")
p3 = plot(z, data_time_mean["ql"], xlabel="ql_time_mean", ylabel="z (km)")
p4 = plot(z, data_time_mean["cld_frac"], xlabel="cld_frac_time_mean", ylabel="z (km)")
p5 = plot(z, data_time_mean["thl"], xlabel="thl_time_mean", ylabel="z (km)")
p6 = plot(z, data_time_mean["ht"], xlabel="ht_time_mean", ylabel="z (km)")
fig = plot(p1, p2, p3, p4, p5, p6, layout=(3, 2), size=(1500, 800))
savefig(fig, string(analysis_folder, "time_means.png"))

# Plot variance and covariances
p1 = plot(z, data_time_mean["var_w"], xlabel="w_time_mean", ylabel="z (km)")
p2 = plot(z, data_time_mean["w3"], xlabel="w3_time_mean", ylabel="z (km)")
p3 = plot(z, data_time_mean["var_qt"], xlabel="var_qt_time_mean", ylabel="z (km)")
p4 = plot(z, data_time_mean["var_thl"], xlabel="var_thl_frac_time_mean", ylabel="z (km)")
p5 = plot(z, data_time_mean["cov_w_qt"], xlabel="cov_w_qt_time_mean", ylabel="z (km)")
p6 = plot(z, data_time_mean["cov_w_thl"], xlabel="cov_w_thl_time_mean", ylabel="z (km)")
p7 = plot(z, data_time_mean["w_qt_sgs"], xlabel="w_qt_sgs_time_mean", ylabel="z (km)")
p8 = plot(z, data_time_mean["w_ht_sgs"], xlabel="w_ht_sgs_time_mean", ylabel="z (km)")
fig = plot(p1, p2, p3, p4, p5, p6, p7, p8, layout=(4, 2), size=(2000, 800))
savefig(fig, string(analysis_folder, "covariances.png"))

# Plot time series
p1 = plot(time/3600., data["cld_cover"], xlabel="Time (hour)", ylabel="cld_cover")
fig = plot(p1, size=(1000, 400))
savefig(fig, string(analysis_folder, "time_series.png"))


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractResidualSampler <: Sampleable{Multivariate, Continuous} end

Generate zero-mean samples from Gaussian posterior, assuming priors are standard gaussians.

# Example
```julia
rng = MersenneTwister(145)
rs = FullResidualSampler(λ_information, jac_dλ_dθ)
samples = rand(rng, rs, 10)  # generate 10 samples
```
"""
abstract type AbstractResidualSampler <: Sampleable{Multivariate, Continuous} end

"""
    FullResidualSampler(λ_information, jac_dλ_dθ)

Gaussian posterior's covariance approximated by the Fisher Information

**Caution: This sampler constructs the covariance matrix explicitly, so memory
use grows quadratically with the number of model parameters**

The Fisher information in canonical coordinates and Jacobian of the coordinate transformation
are provided as arguments. Since `` J \\cdot I_{canon} \\cdot J^T `` does not depend on the definition of canonical coordinates,
we omit description of what canonical coordinates are.

# Arguments

* `λ_information::LinearMap`: Fisher Information in canonical coordinates. Coordinates are chosen in the way
  that this matrix is very simple, and diagonal for the univariate distributions
* `jac_dλ_fθ::LinearMap`: Coordinate transformation jacobian between canonical coordinates and
  coordinates of the model
"""
struct FullResidualSampler <: AbstractResidualSampler
    λ_information::LinearMap
    jac_dλ_dθ::LinearMap
end

Base.length(rs::FullResidualSampler) = size(rs.jac_dλ_dθ, 2)

function Distributions._rand!(rng::AbstractRNG, s::FullResidualSampler, x::AbstractVector{T}) where T<:Real
    θ_information = s.jac_dλ_dθ' * s.λ_information * s.jac_dλ_dθ
    root_covariance = cholesky(PositiveFactorizations.Positive, inv(Matrix(θ_information) + I)).L
    x[:] = root_covariance * randn(eltype(root_covariance), size(root_covariance, 1))
end

"""
    ImplicitResidualSampler(λ_information, jac_dλ_dθ; cg_params=NamedTuple())

Memory efficient Gaussian posterior's covariance approximated by the Fisher Information

This sampler uses Conjugate Gradients to iteratively invert  the Fisher information,
never instantiating the covariance in memory explicitly.

# Arguments

* `λ_information::LinearMap`: Fisher Information in canonical coordinates.
* `jac_dλ_fθ::LinearMap`: Coordinate transformation jacobian between canonical coordinates and
* `cg_params::NamedTuple`: Keyword arguments passed to `cg`. Useful for enabling verbose mode
  or to set accuracy/number of iterations
"""
struct ImplicitResidualSampler <: AbstractResidualSampler
    λ_information::LinearMap
    jac_dλ_dθ::LinearMap
    cg_params::NamedTuple
end

const _default_cg_params=NamedTuple()

ImplicitResidualSampler(λ_information::LinearMap,
                        jac_dλ_dθ::LinearMap;
                        cg_params::NamedTuple=_default_cg_params) = ImplicitResidualSampler(λ_information,
                                                                                jac_dλ_dθ,
                                                                                cg_params)

Base.length(rs::ImplicitResidualSampler) = size(rs.jac_dλ_dθ, 2)

function _implicit_rand_impl_args(s::ImplicitResidualSampler)
    ℐ_λ = s.λ_information
    dλ_dθ = s.jac_dλ_dθ
    num_λs = size(dλ_dθ, 1)
    num_θs = size(dλ_dθ, 2)
    sqrt_Id = cholesky_L(ℐ_λ)
    Σ⁻¹_θ_est = dλ_dθ' * ℐ_λ * dλ_dθ + I
    (num_λs, num_θs, sqrt_Id, Σ⁻¹_θ_est)
end

function _implicit_rand_impl!(rng::AbstractRNG, s::ImplicitResidualSampler, x::AbstractVector{T},
                              num_λs, num_θs, root_Id, Σ⁻¹_θ_est) where T<:Real
    dλ_dθ = s.jac_dλ_dθ
    sample_n = randn(rng, num_λs)
    sample_eta = randn(rng, num_θs)
    Δφ = dλ_dθ' * (root_Id * sample_n) + sample_eta
    x[:] = cg(Σ⁻¹_θ_est, Δφ; s.cg_params...)  # Δξ
end

function Distributions._rand!(rng::AbstractRNG, s::ImplicitResidualSampler, x::AbstractVector{T}) where T<:Real
    args = _implicit_rand_impl_args(s)
    _implicit_rand_impl!(rng, s, x, args...)
end

function Distributions._rand!(rng::AbstractRNG, s::ImplicitResidualSampler, A::DenseMatrix{T}) where T<:Real
    args = _implicit_rand_impl_args(s)
    for i = 1:size(A,2)
        _implicit_rand_impl!(rng, s, view(A,:,i), args...)
    end
    return A
end

export FullResidualSampler,
       ImplicitResidualSampler

ext: jl


text:
using DReal
using Base.Test

"""
(set-logic QF_NRA)
(declare-fun x () Int)
(declare-fun y () Real)
(assert (<= 10 x))
(assert (<= x 20))
(assert (<= 30 y))
(assert (<= y 66))
(assert (<= (- (sin x) y) 0.3))
(check-sat)
(exit)
"""

# x = Var(Int, "x", 10,20)
# y = Var(Float64, "y", 30.0, 66.)
# add!(sin(x) - y <= 0.3)
# @test is_satisfiable()
# @show model(x,y)
begin
reset_global_ctx!()

x = Var(Float64, "x", -10.0, 10.0)
y = Var(Int, "y", -10, 10)
z = Var(Float64, "z", 0.0, 10.0)
add!(x-y>z)
@show is_satisfiable()
@show model(x,y,z)
end


text:
using Base.Test
using PyCall

roundtrip(T, x) = convert(T, PyObject(x))
roundtrip(x) = roundtrip(PyAny, x)
roundtripeq(T, x) = roundtrip(T, x) == x
roundtripeq(x) = roundtrip(x) == x

@test roundtripeq(17)
@test roundtripeq(0x39)
@test roundtripeq(true) && roundtripeq(false)
@test roundtripeq(3.14159)
@test roundtripeq(1.3+4.5im)
@test roundtripeq(nothing)
@test roundtripeq("Hello world")
@test roundtripeq("Hëllö")
@test roundtripeq(Symbol, :Hello)
@test roundtripeq(C_NULL) && roundtripeq(convert(Ptr{Void}, 12345))
@test roundtripeq([1,3,4,5]) && roundtripeq([1,3.2,"hello",true])
@test roundtripeq([1 2 3;4 5 6]) && roundtripeq([1. 2 3;4 5 6])
@test roundtripeq((1,(3.2,"hello"),true)) && roundtripeq(())
@test roundtripeq(Int32)
@test roundtripeq([1 => "hello", 2 => "goodbye"]) && roundtripeq(Dict())

@test pycall(PyObject(x -> x + 1), PyAny, 314158) == 314159
@test roundtrip(x -> x + 1)(314158) == 314159

if PyCall.npy_initialized
    @test PyArray(PyObject([1. 2 3;4 5 6])) == [1. 2 3;4 5 6]
end
@test PyVector(PyObject([1,3.2,"hello",true])) == [1,3.2,"hello",true]
@test PyDict(PyObject([1 => "hello", 2 => "goodbye"])) == [1 => "hello", 2 => "goodbye"]

@test roundtripeq({1 2 3; 4 5 6})
@test roundtripeq([])
@test convert(Array{PyAny,1}, PyObject({1 2 3; 4 5 6})) == {{1,2,3},{4,5,6}}
if PyCall.npy_initialized
    @test roundtripeq(begin A = Array(Int); A[1] = 3; A; end)
end
@test convert(PyAny, PyObject(begin A = Array(Any); A[1] = 3; A; end)) == 3

array2py2arrayeq(x) = PyCall.py2array(Float64,PyCall.array2py(x)) == x
@test array2py2arrayeq(rand(3))
@test array2py2arrayeq(rand(3,4))
@test array2py2arrayeq(rand(3,4,5))

@test roundtrip(1:10) == [1:10]

@pyimport math
@test_approx_eq math.sin(3) sin(3)


text:
# ---
# title: 148. Sort List
# id: problem148
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Sort
# link: <https://leetcode.com/problems/sort-list/description/>
# hidden: true
# ---
#
# Given the `head` of a linked list, return _the list after sorting it in
# **ascending order**_.
#
# **Follow up:** Can you sort the linked list in `O(n logn)` time and `O(1)`
# memory (i.e. constant space)?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/14/sort_list_1.jpg)
#
#
#
#     Input: head = [4,2,1,3]
#     Output: [1,2,3,4]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/14/sort_list_2.jpg)
#
#
#
#     Input: head = [-1,5,3,4,0]
#     Output: [-1,0,3,4,5]
#
#
# **Example 3:**
#
#
#
#     Input: head = []
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is in the range `[0, 5 * 104]`.
#   * `-105 <= Node.val <= 105`
#
#
## @lc code=start
using LeetCode

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


text:
module ClimateTasks

using Distributed, DistributedArrays, SparseArrays
using MeshArrays, MITgcmTools, FortranFiles
using JLD, YAML

include("Main.jl")

export StartWorkers, TaskDriver, task1_loop
export MITgcm

end # module


text:
const timestamp_regex = r"\/(\d\d\d\d-\d\d-\d\d-\d\d-\d\d-\d\d-\d\d\d)[\/\-]"

function delete_old_pull_request_branches(AUTOMERGE_INTEGRATION_TEST_REPO, older_than)
    with_cloned_repo(AUTOMERGE_INTEGRATION_TEST_REPO) do git_repo_dir
        cd(git_repo_dir) do
            all_origin_branches = list_all_origin_branches(git_repo_dir)::Vector{String}
            branches_to_delete = String[]
            for branch_name in all_origin_branches
                if occursin(timestamp_regex, branch_name)
                    commit = strip(read(`git rev-parse origin/$(branch_name)`, String))
                    age = get_age_of_commit(commit)
                    if age >= older_than
                        push!(branches_to_delete, branch_name)
                    end
                end
            end
            unique!(branches_to_delete)
            m = 50
            for _ in 1:(ceil(Int, length(branches_to_delete)/m) + 1)
                k = min(m, length(branches_to_delete))
                current_branches = branches_to_delete[1:k]
                branches_to_delete = branches_to_delete[(k+1):end]
                if !isempty(current_branches)
                    try
                        run(`git push origin --delete $(current_branches)`)
                    catch ex
                        @info "Encountered an error while trying to delete multiple branches" exception=(ex, catch_backtrace()) current_branches
                    end
                end
            end
        end
    end
    with_cloned_repo(AUTOMERGE_INTEGRATION_TEST_REPO) do git_repo_dir
        cd(git_repo_dir) do
            all_origin_branches = list_all_origin_branches(git_repo_dir)::Vector{String}
            for branch_name in all_origin_branches
                if occursin(timestamp_regex, branch_name)
                    commit = strip(read(`git rev-parse origin/$(branch_name)`, String))
                    age = get_age_of_commit(commit)
                    if age >= older_than
                        try
                            run(`git push origin --delete $(branch_name)`)
                        catch ex
                            @info "Encountered an error while trying to delete branch" exception=(ex, catch_backtrace()) branch_name
                        end
                    end
                end
            end
        end
    end
    return nothing
end

function now_localzone()
    return TimeZones.now(TimeZones.localzone())
end

function get_age_of_commit(commit)
    commit_date_string = strip(read(`git show -s --format=%cI $(commit)`, String))
    commit_date = TimeZones.ZonedDateTime(commit_date_string, "yyyy-mm-ddTHH:MM:SSzzzz")
    now = TimeZones.ZonedDateTime(TimeZones.now(), TimeZones.localzone())
    age = max(now - commit_date, Dates.Millisecond(0))
    return age
end

function list_all_origin_branches(git_repo_dir)
    result = Vector{String}(undef, 0)
    original_working_directory = pwd()
    cd(git_repo_dir)
    a = try
        read(`git branch -a`, String)
    catch
        ""
    end
    b = split(strip(a), '\n')
    b_length = length(b)
    c = Vector{String}(undef, b_length)
    for i = 1:b_length
        c[i] = strip(strip(strip(b[i]), '*'))
        c[i] = first(split(c[i], "->"))
        c[i] = strip(c[i])
    end
    my_regex = r"^remotes\/origin\/(.*)$"
    for i = 1:b_length
        if occursin(my_regex, c[i])
            m = match(my_regex, c[i])
            if m[1] != "HEAD"
                push!(result, m[1])
            end
        end
    end
    cd(original_working_directory)
    return result
end

function with_cloned_repo(f, repo_url)
    return mktempdir() do tmp_dir
        return cd(tmp_dir) do
            git_repo_dir = joinpath(tmp_dir, "REPO")
            run(`git clone $(repo_url) REPO`)
            return cd(git_repo_dir) do
                return f(git_repo_dir)
            end
        end
    end
end

function empty_git_repo(git_repo_dir::AbstractString)
    original_working_directory = pwd()
    cd(git_repo_dir)
    for x in readdir(git_repo_dir)
        if x != ".git"
            path = joinpath(git_repo_dir, x)
            rm(path; force = true, recursive = true)
        end
    end
    cd(original_working_directory)
    return nothing
end

function _generate_branch_name(name::AbstractString)
    sleep(5)
    _now = now_localzone()
    _now_utc_string = utc_to_string(_now)
    b = "integration/$(_now_utc_string)/$(rand(UInt32))/$(name)"
    sleep(5)
    return b
end

function generate_branch(name::AbstractString,
                         path_to_content::AbstractString,
                         parent_branch::AbstractString = "master";
                         repo_url)
    original_working_directory = pwd()
    b = _generate_branch_name(name)
    with_cloned_repo(repo_url) do git_repo_dir
        cd(git_repo_dir)
        run(`git checkout $(parent_branch)`)
        run(`git branch $(b)`)
        run(`git checkout $(b)`)
        empty_git_repo(git_repo_dir)
        for x in readdir(path_to_content)
            src = joinpath(path_to_content, x)
            dst = joinpath(git_repo_dir, x)
            rm(dst; force = true, recursive = true)
            cp(src, dst; force = true)
        end
        cd(git_repo_dir)
        CompatHelper.my_retry(() -> run(`git add -A`))
        CompatHelper.my_retry(() -> run(`git commit -m "Automatic commit - CompatHelper integration tests"`))
        CompatHelper.my_retry(() -> run(`git push origin $(b)`))
        cd(original_working_directory)
        rm(git_repo_dir; force = true, recursive = true)
    end
    return b
end

function generate_master_branch(path_to_content::AbstractString,
                                parent_branch::AbstractString = "master";
                                repo_url)
    name = "master"
    b = generate_branch(name, path_to_content, parent_branch; repo_url = repo_url)
    return b
end

function templates(parts...)
    this_filename = @__FILE__
    test_directory = dirname(this_filename)
    templates_directory = joinpath(test_directory, "templates")
    result = joinpath(templates_directory, parts...)
    return result
end

function username(auth::GitHub.Authorization)
    sleep(5)
    user_information = GitHub.gh_get_json(GitHub.DEFAULT_API,
                                          "/user";
                                          auth = auth)
    sleep(5)
    return user_information["login"]::String
end

function utc_to_string(zdt::TimeZones.ZonedDateTime)
    zdt_as_utc = TimeZones.astimezone(zdt, TimeZones.tz"UTC")
    year = TimeZones.Year(zdt_as_utc.utc_datetime).value
    month = TimeZones.Month(zdt_as_utc.utc_datetime).value
    day = TimeZones.Day(zdt_as_utc.utc_datetime).value
    hour = TimeZones.Hour(zdt_as_utc.utc_datetime).value
    minute = TimeZones.Minute(zdt_as_utc.utc_datetime).value
    second = TimeZones.Second(zdt_as_utc.utc_datetime).value
    millisecond = TimeZones.Millisecond(zdt_as_utc.utc_datetime).value
    result = Printf.@sprintf "%04d-%02d-%02d-%02d-%02d-%02d-%03d" year month day hour minute second millisecond
    return result
end

function with_master_branch(f::Function,
                            path_to_content::AbstractString,
                            parent_branch::AbstractString;
                            repo_url)
    b = generate_master_branch(path_to_content,
                               parent_branch;
                               repo_url = repo_url)
    result = f(b)
    return result
end


text:
module TestIdealPoints
    using Base.Test

    import RollCallDataIO
    import IdealPoints
    import FiniteDiff

    io = open(joinpath("data", "senate", "112.ord"), "r")
    ord_file = read(RollCallDataIO.ORDFile, io)
    close(io)

    roll_call = convert(
        RollCallDataIO.SparseRollCall,
        RollCallDataIO.roll_call(ord_file),
    )

    n_legislators, n_bills = roll_call.n_legislators, roll_call.n_bills

    for d in 1:2
        a, b, c, misc = IdealPoints.ideal_points(ord_file, d)

        θ = fill(
            0.0,
            d * n_legislators + d * n_bills + n_bills
        )

        nlp = IdealPoints.make_nlp(
            d,
            1.0,
            n_legislators,
            n_bills,
            roll_call.legislators,
            roll_call.bills,
            roll_call.votes,
        )

        nlp_gr! = IdealPoints.make_nlp_gr(
            d,
            1.0,
            n_legislators,
            n_bills,
            roll_call.legislators,
            roll_call.bills,
            roll_call.votes,
        )

        θstar = vcat(a, b, c)
        @test nlp(θstar) < nlp(θ)

        gr = copy(θ)
        nlp_gr!(θ, gr)
        @test norm(FiniteDiff.gradient(nlp, θ) - gr) < 1e-4

        nlp_gr!(θstar, gr)
        @test norm(FiniteDiff.gradient(nlp, θstar) - gr) < 1e-4
        @test norm(gr) < 1e-4

        # Check in a region around θstar that θstar produces a smaller nlp.
        for _ in 1:10
            for i in 1:length(θ)
                θ[i] = θstar[i] + 0.01 * randn()
            end
            @test nlp(θstar) < nlp(θ)

            nlp_gr!(θ, gr)
            n1 = norm(gr)

            nlp_gr!(θstar, gr)
            n2 = norm(gr)
            @test n2 < n1
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function get_test_2d_model()
    X = Dict{Int64, Vector{Float64}}(
        7 => [0.0, 1.0],
        8 => [5/4, 1.0],
        9 => [2.0, 1.0],
        10 => [0.0, 1.0],
        11 => [3/4, 1.0],
        12 => [2.0, 1.0])
    mel1 = Element(Seg2, [7, 8])
    mel2 = Element(Seg2, [8, 9])
    sel1 = Element(Seg2, [10, 11])
    sel2 = Element(Seg2, [11, 12])
    update!([mel1, mel2, sel1, sel2], "geometry", X)
    update!([sel1, sel2], "master elements", [sel1, sel2])
    calculate_normals!([sel1, sel2], 0.0, Val{1})
    return [sel1, sel2], [mel1, mel2]
end

@testset "calculate flat 2d projection from slave to master" begin
    (sel1, sel2), (mel1, mel2) = get_test_2d_model()

    time = 0.0
    X1 = sel1("geometry", [-1.0], time)
    n1 = sel1("normal", [-1.0], time)
    xi2 = project_from_slave_to_master(mel1, X1, n1, time)
    @test isapprox(xi2, -1.0)

    X1 = sel1("geometry", [1.0], time)
    n1 = sel1("normal", [1.0], time)
    xi2 = project_from_slave_to_master(mel1, X1, n1, time)
    @test isapprox(xi2, 0.2)

    X2 = mel1("geometry", xi2, time)
    @test isapprox(X2, [3/4, 1.0])
end

@testset "calculate flat 2d projection from master to slave" begin
    (sel1, sel2), (mel1, mel2) = get_test_2d_model()
    time = 0.0
    x2 = mel1("geometry", [-1.0], time)
    xi1 = project_from_master_to_slave(sel1, x2, time)
    @test isapprox(xi1, -1.0)
    x2 = mel1("geometry", [1.0], time)
    xi1 = project_from_master_to_slave(sel1, x2, time)
    X1 = sel1("geometry", xi1, time)
    @test isapprox(X1, [5/4, 1.0])
end

@testset "calculate flat 2d projection rotated 90 degrees" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [0.0, 1.0],
        3 => [0.0, 1.0],
        4 => [0.0, 0.0])
    sel1 = Element(Seg2, [1, 2])
    mel1 = Element(Seg2, [3, 4])
    update!([sel1, mel1], "geometry", X)
    time = 0.0
    calculate_normals!([sel1], time, Val{1})

    X2 = mel1("geometry", [-1.0], time)
    xi = project_from_master_to_slave(sel1, X2, time)
    @test isapprox(xi, 1.0)

    X2 = mel1("geometry", [1.0], time)
    xi = project_from_master_to_slave(sel1, X2, time)
    @test isapprox(xi, -1.0)

    X1 = sel1("geometry", [-1.0], time)
    n1 = sel1("normal", [-1.0], time)
    xi = project_from_slave_to_master(mel1, X1, n1, time)
    @test isapprox(xi, 1.0)

    X1 = sel1("geometry", [1.0], time)
    n1 = sel1("normal", [1.0], time)
    xi = project_from_slave_to_master(mel1, X1, n1, time)
    @test isapprox(xi, -1.0)
end


text:
# It appears that PythonCall and PyCall cannot co-exist in a module.
# see https://github.com/cjdoris/PythonCall.jl/issues/140

using SparseArrays

import PythonCall: Py, PyList
import PythonCall: pyimport, pyconvert, pyconvert_return, pyconvert_add_rule

"""
Convert KVXOPT.spmatrix to SparseMatrixCSC using PythonCall.
"""
function convert_spmatrix(S, x::Py)
    if pyconvert(String, x.typecode) == "d"
        S = SparseMatrixCSC(pyconvert(Int64, x.size[0]), pyconvert(Int64, x.size[1]),
                            pyconvert(Vector{Int64}, PyList(x.CCS[0])) .+ 1,
                            pyconvert(Vector{Int64}, PyList(x.CCS[1])) .+ 1,
                            pyconvert(Vector{Float64}, PyList(x.CCS[2])))::SparseMatrixCSC{Float64,Int64}

    elseif pyconvert(String, x.typecode) == "z"
        S = SparseMatrixCSC(pyconvert(Int64, x.size[0]), pyconvert(Int64, x.size[1]),
                            pyconvert(Vector{Int64}, PyList(x.CCS[0])) .+ 1,
                            pyconvert(Vector{Int64}, PyList(x.CCS[1])) .+ 1,
                            pyconvert(Vector{ComplexF64}, PyList(x.CCS[2])))::SparseMatrixCSC{ComplexF64,ComplexF64}
    end
    return pyconvert_return(S)
end

pc = pyimport("andes")
ss = pc.run(pc.get_case("kundur/kundur_full.xlsx"), no_output=true)

convert_spmatrix(Any, ss.dae.gy)

pyconvert_add_rule("kvxopt.base:spmatrix", SparseMatrixCSC, convert_spmatrix)

@time pyconvert(SparseMatrixCSC, ss.dae.gy)


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

function nnCostFunction(nn_params, input_layer_size, hidden_layer_size,
                        num_labels, X, y, lambda)
  # Reshape nn_params back into the parameters Theta1 and Theta2, the weight matrices
  # for our 2 layer neural network
  Theta1 = reshape(nn_params[1:hidden_layer_size * (input_layer_size + 1)],
                   hidden_layer_size, (input_layer_size + 1))

  Theta2 = reshape(nn_params[(1 + (hidden_layer_size * (input_layer_size + 1))):end],
                   num_labels, (hidden_layer_size + 1))

  # Setup some useful variables
  m = size(X, 1)

  # You need to return the following variables correctly
      [J grad] = NNCOSTFUNCTON(nn_params, hidden_layer_size, num_labels,
    X, y, lambda)
  Theta1_grad = zeros(size(Theta1))
  Theta2_grad = zeros(size(Theta2))

  # ====================== YOUR CODE HERE ======================
  # Instructions: You should complete the code by working through the
  #               following parts.
  #
  # Part 1: Feedforward the neural network and return the cost in the
  #         variable J. After implementing Part 1, you can verify that your
  #         cost function computation is correct by verifying the cost
  #         computed in ex4.m
  #
  # Part 2: Implement the backpropagation algorithm to compute the gradients
  #         Theta1_grad and Theta2_grad. You should return the partial derivatives of
  #         the cost function with respect to Theta1 and Theta2 in Theta1_grad and
  #         Theta2_grad, respectively. After implementing Part 2, you can check
  #         that your implementation is correct by running checkNNGradients
  #
  #         Note: The vector y passed into the function is a vector of labels
  #               containing values from 1..K. You need to map this vector into a
  #               binary vector of 1's and 0's to be used with the neural network
  #               cost function.
  #
  #         Hint: We recommend implementing backpropagation using a for-loop
  #               over the training examples if you are implementing it for the
  #               first time.
  #
  # Part 3: Implement regularization with the cost function and gradients.
  #
  #         Hint: You can implement this around the code for
  #               backpropagation. That is, you can compute the gradients for
  #               the regularization separately and then add them to Theta1_grad
  #               and Theta2_grad from Part 2.
  #
  # =========================================================================

  # Unroll gradients
  grad = [Theta1_grad[:] ; Theta2_grad[:]]
  return (J, grad)
end


text:
using STMMRunner
using Documenter

DocMeta.setdocmeta!(STMMRunner, :DocTestSetup, :(using STMMRunner); recursive = true)

makedocs(;
    modules = [STMMRunner],
    authors = "Gabriel Wu <wuzihua@pku.edu.cn> and contributors",
    repo = "https://github.com/JuliaRemoteSensing/STMMRunner.jl/blob/{commit}{path}#{line}",
    sitename = "STMMRunner.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://JuliaRemoteSensing.github.io/STMMRunner.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Configurations" => "config.md",
        "MSTM v4" => "mstmv4.md",
        "MSTM v3" => "mstmv3.md",
    ],
)

deploydocs(; repo = "github.com/JuliaRemoteSensing/STMMRunner.jl", devbranch = "main")


text:
using Oceananigans: instantiated_location

#####
##### Algorithm for adding fluxes associated with non-trivial flux boundary conditions.
##### Value and Gradient boundary conditions are handled by filling halos.
#####

"""
    apply_x_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
apply_x_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :yz, _apply_x_bcs!, Gc, instantiated_location(Gc), c.grid,
                                                  c.boundary_conditions.x.left, c.boundary_conditions.x.right, args...,
                                                  dependencies=dep)

"""
    apply_y_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
apply_y_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :xz, _apply_y_bcs!, Gc, instantiated_location(Gc), c.grid,
                                                  c.boundary_conditions.y.left, c.boundary_conditions.y.right, args...,
                                                  dependencies=dep)

"""
    apply_z_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the top and bottom.
"""
apply_z_bcs!(Gc, c, arch, dep, args...) = launch!(arch, c.grid, :xy, _apply_z_bcs!, Gc, instantiated_location(Gc), c.grid,
                                                  c.boundary_conditions.z.left, c.boundary_conditions.z.right, args...,
                                                  dependencies=dep)

apply_x_bcs!(::Nothing, args...) = nothing
apply_y_bcs!(::Nothing, args...) = nothing
apply_z_bcs!(::Nothing, args...) = nothing

"""
    _apply_x_bcs!(Gc, grid, west_bc, east_bc, args...)

Apply a west and/or east boundary condition to variable `c`.
"""
@kernel function _apply_x_bcs!(Gc, loc, grid, west_bc, east_bc, args...)
    j, k = @index(Global, NTuple)
    apply_x_west_bc!(Gc, loc, west_bc, j, k, grid, args...)
    apply_x_east_bc!(Gc, loc, east_bc, j, k, grid, args...)
end

"""
    _apply_y_bcs!(Gc, grid, south_bc, north_bc, args...)

Apply a south and/or north boundary condition to variable `c`.
"""
@kernel function _apply_y_bcs!(Gc, loc, grid, south_bc, north_bc, args...)
    i, k = @index(Global, NTuple)
    apply_y_south_bc!(Gc, loc, south_bc, i, k, grid, args...)
    apply_y_north_bc!(Gc, loc, north_bc, i, k, grid, args...)
end

"""
    _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)

Apply a top and/or bottom boundary condition to variable `c`.
"""
@kernel function _apply_z_bcs!(Gc, loc, grid, bottom_bc, top_bc, args...)
    i, j = @index(Global, NTuple)
    apply_z_bottom_bc!(Gc, loc, bottom_bc, i, j, grid, args...)
       apply_z_top_bc!(Gc, loc, top_bc,    i, j, grid, args...)
end

# Avoid some computation / memory accesses for Value, Gradient, Periodic, NormalFlow,
# and zero-flux boundary conditions --- every boundary condition that does *not* prescribe
# a non-trivial flux.
const NotFluxBC = Union{VBC, GBC, PBC, NFBC, ZFBC}

# Fall back functions for boundary conditions that are not of type Flux.
@inline apply_x_east_bc!(  Gc, args...) = nothing
@inline apply_x_west_bc!(  Gc, args...) = nothing
@inline apply_y_north_bc!( Gc, args...) = nothing
@inline apply_y_south_bc!( Gc, args...) = nothing
@inline apply_z_top_bc!(   Gc, args...) = nothing
@inline apply_z_bottom_bc!(Gc, args...) = nothing

# Shortcuts for 'zero' flux boundary conditions.
@inline apply_x_east_bc!(  Gc, loc, ::ZFBC, args...) = nothing
@inline apply_x_west_bc!(  Gc, loc, ::ZFBC, args...) = nothing
@inline apply_y_north_bc!( Gc, loc, ::ZFBC, args...) = nothing
@inline apply_y_south_bc!( Gc, loc, ::ZFBC, args...) = nothing
@inline apply_z_top_bc!(   Gc, loc, ::ZFBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, loc, ::ZFBC, args...) = nothing

@inline flip(::Center) = Face()
@inline flip(::Face) = Center()

"""
    apply_x_west_bc!(Gc, west_flux::BC{<:Flux}, j, k, grid, args...)

Add the flux divergence associated with a west flux boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive west flux is associated with an *increase* in `Gc` near the west boundary.
If `west_bc.condition` is a function, the function must have the signature

    `west_bc.condition(j, k, grid, boundary_condition_args...)`

The same logic holds for south and bottom boundary conditions in `y`, and `z`, respectively.
"""
@inline function apply_x_west_bc!(Gc, loc, west_flux::BC{<:Flux}, j, k, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[1, j, k] += getbc(west_flux, j, k, grid, args...) * Ax(1, j, k, grid, flip(LX), LY, LZ) / volume(1, j, k, grid, LX, LY, LZ)
    return nothing
end

@inline function apply_y_south_bc!(Gc, loc, south_flux::BC{<:Flux}, i, k, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[i, 1, k] += getbc(south_flux, i, k, grid, args...) * Ay(i, 1, k, grid, LX, flip(LY), LZ) / volume(i, 1, k, grid, LX, LY, LZ)
    return nothing
end

@inline function apply_z_bottom_bc!(Gc, loc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) * Az(i, j, 1, grid, LX, LY, flip(LZ)) / volume(i, j, 1, grid, LX, LY, LZ)
    return nothing
end

"""
    apply_x_east_bc!(Gc, top_flux::BC{<:Flux}, j, k, grid, args...)

Add the part of flux divergence associated with a east boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive east flux is associated with a *decrease* in `Gc` near the east boundary.
If `east_bc.condition` is a function, the function must have the signature

    `east_bc.condition(i, j, grid, boundary_condition_args...)`

The same logic holds for north and top boundary conditions in `y`, and `z`, respectively.
"""
@inline function apply_x_east_bc!(Gc, loc, east_flux::BC{<:Flux}, j, k, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[grid.Nx, j, k] -= getbc(east_flux, j, k, grid, args...) * Ax(grid.Nx+1, j, k, grid, flip(LX), LY, LZ) / volume(grid.Nx, j, k, grid, LX, LY, LZ)
    return nothing
end

@inline function apply_y_north_bc!(Gc, loc, north_flux::BC{<:Flux}, i, k, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[i, grid.Ny, k] -= getbc(north_flux, i, k, grid, args...) * Ay(i, grid.Ny+1, k, grid, LX, flip(LY), LZ) / volume(i, grid.Ny, k, grid, LX, LY, LZ)
    return nothing
end

@inline function apply_z_top_bc!(Gc, loc, top_flux::BC{<:Flux}, i, j, grid, args...)
    LX, LY, LZ = loc
    @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) * Az(i, j, grid.Nz+1, grid, LX, LY, flip(LZ)) / volume(i, j, grid.Nz, grid, LX, LY, LZ)
    return nothing
end


text:
# Script executed by Biowulf to run code

@everywhere include("/home/carsonc/StochasticGene/src/StochasticGene.jl")

include("/home/carsonc/StochasticGene/runfiles/scriptfunctions.jl")
@time fit_rna(parse(Int,ARGS[1]),ARGS[2],ARGS[3],ARGS[15],ARGS[16],ARGS[4],parse(Int,ARGS[5]),parse(Float64,ARGS[6]),ARGS[7],ARGS[8],ARGS[9],ARGS[10],ARGS[11],parse(Int,ARGS[12]),parse(Bool,ARGS[13]),parse(Bool,ARGS[14]))

# fit_rna(nchains::Int,gene::String,cell::String,fittedparam::String,fixedeffects::String,datacond,G::Int,maxtime::Float64,infolder::String,resultfolder::String,datafolder,inlabel,label,nsets,runcycle::Bool=false,transient::Bool=false,samplesteps::Int=100000,warmupsteps=20000,annealsteps=0,temp=1.,tempanneal=100.,root = "/home/carsonc/scrna/")

# Arguments
# 1: nchains
# 2: gene
# 3: cell
# 4: cond
# 5: G
# 6: maxtime
# 7: infolder
# 8: resultfolder
# 9: datafolder
# 10: inlabel
# 11: label
# 12: nsets (number of rate sets)
# 13: runcycle (bool)
# 14: transient (bool)
# 15: fittedparams (string), e.g. "1-2-4-7"
# 16: fixedeffects (string), e.g. "off"

max_stars_repo_name: nih-niddk-mbs/Transcription

max_issues_repo_name: nih-niddk-mbs/Transcription

max_forks_repo_name: nih-niddk-mbs/Transcription


text:
# Methods for the Semidefinite interface

@compat abstract type AbstractSDModel <: MPB.AbstractMathProgModel end
export AbstractSDModel

MPB.@define_interface begin
    SDModel
    setconstrB!
    setconstrentry!
    setobjentry!
end


text:
"""
# OptionalData

This package provides the `@OptionalData` macro and the corresponding `OptData`
type which is a thin wrapper around a nullable value (of type `Union{T, Nothing} where T`).
It allows you to load and access globally available data at runtime in a type-stable way.

## Usage

*OptionalData* has the following use cases:

1. Parts of your package depend on data from the internet while other parts do not.
   In the case of a network outage the package should offer a degraded experience but
   the independent parts should still function.
2. Your package requires manual initialisation steps, e.g. loading data from a
   user-supplied file, and you do not want to repeat yourself writing code that
   checks for the availability of the data.

You declare optional global data with the `@OptionalData` macro:

```julia
using OptionalData

# @OptionalData name type [error_msg]
@OptionalData OPT_FLOAT Float64 "Forgot to load it?"

# this expands to
const OPT_FLOAT = OptData{Float64}(string(:OPT_FLOAT), "Forgot to load it?")
```

You access its value with `get` and check whether it is available with `isavailable`:

```julia
# This will throw a `NoDataError` because OPT_FLOAT does not contain a value, yet.
get(OPT_FLOAT)
# ERROR: OPT_FLOAT is not available. Forgot to load it?
isavailable(OPT_FLOAT) == false
```

Use `push!` to load the data:

```julia
push!(OPT_FLOAT, 3.0)
isavailable(OPT_FLOAT) == true
get(OPT_FLOAT) == 3.0
```
"""
module OptionalData

export @OptionalData, OptData, NoDataError, show, push!, isavailable, get

import Base: push!, get, show

mutable struct OptData{T}
    data::Union{T, Nothing}
    name::String
    msg::String
end
OptData{T}(name, msg="") where {T} = OptData{T}(nothing, name, msg)

"""
    NoDataError(opt::OptData)

`opt` does not contain data.
"""
struct NoDataError <: Exception
    opt::OptData
end

function Base.showerror(io::IO, err::NoDataError)
    print(io, "`", err.opt.name, "` is not available.")
    if !isempty(err.opt.msg)
        print(io, " ", err.opt.msg)
    end
end

function show(io::IO, opt::OptData{T}) where T
    val = isavailable(opt) ? get(opt) : ""
    print(io, "OptData{$T}($val)")
end

"""
    push!(opt::OptData{T}, data::T) where T

Push `data` of type `T` to `opt`.
"""
function push!(opt::OptData{T}, data::T) where T
    opt.data = data
    opt
end

"""
    push!(opt::OptData, ::Type{T}, args...) where T

Construct an object of type `T` from `args` and push it to `opt`.
"""
function push!(opt::OptData, ::Type{T}, args...) where {T}
    push!(opt, T(args...))
    opt
end
push!(opt::OptData{T}, args...) where {T} = push!(opt, T, args...)

"""
    get(opt::OptData)

Get data from `opt`. Throw an exception if no data has been pushed to `opt` before.
"""
function get(opt::OptData{T}) where T
    !isavailable(opt) && throw(NoDataError(opt))
    opt.data::T
end

"""
    @OptionalData name type msg=""

Initialise a constant `name` with type `OptData{type}`.
An exception with the custom error message `msg` is thrown when `name` is
accessed before data has been pushed to it.

# Example

```julia
@OptionalData OPT_FLOAT Float64
```
"""
macro OptionalData(name, typ, msg="")
    :(const $(esc(name)) = OptData{$(esc(typ))}($(string(name)), $msg))
end

"""
    isavailable(opt::OptData)

Check whether data has been pushed to `opt`.
"""
isavailable(opt::OptData) = opt.data !== nothing

end # module

max_stars_repo_name: helgee/OptionalData

max_issues_repo_name: helgee/OptionalData

max_forks_repo_name: helgee/OptionalData


text:
## FALL-BACK METHODS FOR COMPOSITE MODELS (EXPORTED LEARNING NETWORKS)

# *Note.* Be sure to read Note 4 in src/operations.jl to see see how
# fallbacks are provided for operations acting on Composite models.

caches_data_by_default(::Type{<:Composite}) = true

# builds on `fitted_params(::CompositeFitresult)` defined in
# composition/learning_networks/machines.jl:
fitted_params(::Union{Composite,Surrogate}, fitresult::CompositeFitresult) =
    fitted_params(glb(fitresult))

function update(model::M,
                verbosity::Integer,
                fitresult::CompositeFitresult,
                cache,
                args...) where M <: Composite

    # This method falls back to `fit` to force rebuilding the
    # underlying learning network if, since the last fit:
    #
    # (i) Any hyper-parameter associated with a model in the learning network
    #     has been replaced with a new model instance (and not merely
    #     mutated), OR

    # (ii) Any OTHER hyper-parameter has changed it's value (in the sense
    # of `==`).

    # Otherwise, a "smart" fit is carried out by calling `fit!` on a
    # greatest lower bound node for nodes in the signature of the
    # underlying learning network machine. For this it is necessary to
    # temporarily "de-anonymize" the source nodes.

    network_model_names = cache.network_model_names
    old_model = cache.old_model

    glb_node = glb(fitresult) # greatest lower bound

    if fallback(model, old_model, network_model_names, glb_node)
        return fit(model, verbosity, args...)
    end

    # return data to source nodes for fitting:
    sources, data = cache.sources, cache.data
    for k in eachindex(sources)
        rebind!(sources[k], data[k])
    end

    fit!(glb_node; verbosity=verbosity)
    update!(fitresult) # updates report_additions

    # anonymize data again:
    for s in sources
        rebind!(s, nothing)
    end

    # record current model state:
    cache = (sources=cache.sources,
             data=cache.data,
             network_model_names=cache.network_model_names,
             old_model = deepcopy(model))

    return (fitresult,
            cache,
            merge(report(glb_node),
                  report_additions(fitresult)))

end

# helper for preceding method (where logic is explained):
function fallback(model::M, old_model, network_model_names, glb_node) where M
    # check the hyper-parameters corresponding to models:
    network_models = MLJBase.models(glb_node)
    for j in eachindex(network_models)
        name = network_model_names[j]
        name === nothing ||
            objectid(network_models[j])===objectid(getproperty(model, name)) ||
            return true
    end
    # check any other hyper-parameter:
    for name in propertynames(model)
        if !(name in network_model_names)
            old_value = getproperty(old_model, name)
            value = getproperty(model, name)
            value == old_value || return true
        end
    end
    return false
end


text:
function _lookup(terms)
    mapreduce(intersect, _tokenize(terms)) do token
        get(_INDEX, token, Set{String}())
    end
end


text:
using GLAbstraction, Colors, GeometryTypes, GLVisualize, Reactive

if !isdefined(:runtests) # We use the examples to generate docs and tests
    window = glscreen()
end

vx = -10:0.1:10;
vy = -8:0.1:8;
x = [i for i in vx, j in vy];
y = [j for i in vx, j in vy];
#x,y = meshgrid(vx, vy)
z = (x .* y .* (x.^2 - y.^2) ./ (x.^2 .+ y.^2 + eps())) / 5.;
xyz  = map(Point3f0, zip(x, y, z));
# ma
lines3d = visualize(
	xyz, :lines,
	color=RGBA(0.0f0,0.74736935f0,1.0f0,0.6f0),
	thickness=4f0
)
_view(lines3d, window, camera=:perspective)

if !isdefined(:runtests)
    @async renderloop(window)
end


text:
include(joinpath("toolkit","all.jl"))
include(joinpath("geometric_center","all.jl"))
include(joinpath("distance","all.jl"))
include(joinpath("principal_axes","all.jl"))
include(joinpath("angle","all.jl"))


text:
export reproject, convertmode, userval, userbounds

# These methods are only available if ArchGDAL is loaded.
# Otherwise Projected selector crs field must be `nothing`,
# and no reprojection can occur.

sel2indices(mode::Projected, dim::Dimension, sel::Contains{<:Number}) = begin
    selval = reproject(usercrs(mode), crs(mode), dim, val(sel))
    DD.contains(dim, rebuild(sel, selval))
end
sel2indices(mode::Projected, dim::Dimension, sel::At{<:Number}) = begin
    selval = reproject(usercrs(mode), crs(mode), dim, val(sel))
    DD.at(dim, rebuild(sel, selval))
end
sel2indices(mode::Projected, dim::Dimension, sel::Between) = begin
    selval = map(v -> reproject(usercrs(mode), crs(mode), dim, v), val(sel))
    DD.between(dim, rebuild(sel, selval))
end

"""
`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
a single dimension at a time.
"""
reproject(source, target, dim::Dimension, val) = val
reproject(source::GeoFormat, target::GeoFormat, dim::Lon, val::Number) =
    ArchGDAL.reproject((Float64(val), 0.0), source, target; order=:trad)[1]
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, val::Number) =
    ArchGDAL.reproject((0.0, Float64(val)), source, target; order=:trad)[2]

reproject(source::GeoFormat, target::GeoFormat, ::Lon, vals::AbstractArray) =
    [r[1] for r in ArchGDAL.reproject([(Float64(v), 0.0) for v in vals], source, target; order=:trad)]
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, vals::AbstractArray) =
    [r[2] for r in ArchGDAL.reproject([(0.0, Float64(v)) for v in vals], source, target; order=:trad)]

reproject(source::GeoFormat, target::GeoFormat, ::Lon, vals::Tuple) =
    Tuple(r[1] for r in ArchGDAL.reproject([(Float64(v), 0.0) for v in vals], source, target; order=:trad))
reproject(source::GeoFormat, target::GeoFormat, dim::Lat, vals::Tuple) =
    Tuple(r[2] for r in ArchGDAL.reproject([(0.0, Float64(v)) for v in vals], source, target; order=:trad))

convertmode(dstmode::Type{Converted}, srcmode::Type{Projected}, dim::Dimension) where M = begin
    m = mode(dim)
    newval = reproject(crs(m), usercrs(m), dim, val(dim))
    newbounds = reproject(crs(m), usercrs(m), dim, bounds(dim))
    newmode = Converted(order(m), Irregular(newbounds), sampling(m), crs(m), usercrs(m))
    rebuild(dim; val=newval, mode=newmode)
end
convertmode(dstmode::Type{Projected}, srcmode::Type{Converted}, dim::Dimension) where M = begin
    m = mode(dim)
    start, stop = reproject(dimcrs(m), crs(m), dim, [first(dim), last(dim)])
    newval = LinRange(start, stop, length(dim))
    newmode = Projected(order(m), Regular(step(newval)), sampling(m), crs(m), dimcrs(m))
    rebuild(dim; val=newval, mode=newmode)
end

# Add Lat/Lon methods to reproject bounds and val
userbounds(dim::Union{Lat,Lon}) = userbounds(mode(dim), dim)
userbounds(::Projected, dim::Union{Lat,Lon}) =
    reproject(crs(dim), usercrs(dim), dim, bounds(dim))
userbounds(::IndexMode, dim::Union{Lat,Lon}) = bounds(dim)

userval(dim::Union{Lat,Lon}) = userval(mode(dim), dim)
userval(::Projected, dim::Union{Lat,Lon}) =
    reproject(crs(dim), usercrs(dim), dim, val(dim))
userval(::IndexMode, dim::Union{Lat,Lon}) = val(dim)


text:
ds_small = GDAL.open("data/utmsmall.tif", GDAL.GA_ReadOnly)
ds_point = GDAL.openex("data/point.geojson", GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)

# GDALInfo
# infooptionsnew checks if the options are valid
@test_throws GDAL.GDALError GDAL.infooptionsnew(["-novalidoption"], C_NULL)
# check not only that a GDALError is thrown, but also its contents
try
    GDAL.infooptionsnew(["-novalidoption"], C_NULL)
catch err
    @test err.class === GDAL.CE_Failure
    @test err.code === Cint(6)
    @test err.msg == "Unknown option name '-novalidoption'"
end

options = GDAL.infooptionsnew(["-checksum"], C_NULL)
infostr = GDAL.info(ds_small, options)
@test contains(infostr, "Checksum=50054")
GDAL.infooptionsfree(options)

# GDALTranslate
# resample to a 5×5 ascii grid
optvec = ["-of","AAIGrid","-r","cubic","-tr","1200","1200"]
options = GDAL.translateoptionsnew(optvec, C_NULL)
ds_tiny_asc = GDAL.translate("data/utmtiny.asc", ds_small, options, C_NULL)
GDAL.translateoptionsfree(options)
GDAL.close(ds_tiny_asc)
@test replace(readstring("data/utmtiny.asc"), "\r", "") == """
    ncols        5
    nrows        5
    xllcorner    440720.000000000000
    yllcorner    3745320.000000000000
    cellsize     1200.000000000000
     128 126 161 185 193
     171 164 175 206 205
     127 148 177 205 209
     93 114 164 172 181
     83 101 140 128 122
    """
# also create a Byte GeoTIFF version for later use
optvec = ["-r","cubic","-tr","1200","1200","-ot","Byte"]
options = GDAL.translateoptionsnew(optvec, C_NULL)
ds_tiny = GDAL.translate("data/utmtiny.tif", ds_small, options, C_NULL)
GDAL.translateoptionsfree(options)

tinydata = UInt8[
    128 126 161 185 193;
    171 164 175 206 205;
    127 148 177 205 209;
     93 114 164 172 181;
     83 101 140 128 122]

band = GDAL.getrasterband(ds_tiny, 1)
read_data = zeros(tinydata)
GDAL.rasterio(band, GDAL.GF_Read, 0, 0, 5, 5,
              read_data, 5, 5, GDAL.GDT_Byte, 0, 0)
@test read_data' == tinydata

# GDALWarp

# cannot reproject file
# ERROR: LoadError: LoadError: GDALError (Failure, code 6):
#        Unable to load PROJ.4 library (libproj.so), creation of OGRCoordinateTransformation failed.
# when using
# options = GDAL.warpappoptionsnew(["-t_srs","EPSG:4326"], C_NULL)
# should we require PROJ.4 or not?
options = GDAL.warpappoptionsnew(["-of","MEM"], C_NULL)  # in memory
ds_warped = GDAL.warp("data/utmsmall.mem", Ptr{GDAL.GDALDatasetH}(C_NULL), 1, [ds_small], options, C_NULL)
GDAL.warpappoptionsfree(options)
GDAL.close(ds_warped)

# GDALDEMProcessing
# calculate hillshade
options = GDAL.demprocessingoptionsnew(["-of","AAIGrid"], C_NULL)
ds_dempr = GDAL.demprocessing("data/utmtiny-hillshade.asc", ds_tiny, "hillshade", C_NULL, options, C_NULL)
GDAL.demprocessingoptionsfree(options)
GDAL.close(ds_dempr)
@test replace(readstring("data/utmtiny-hillshade.asc"), "\r", "") == """
    ncols        5
    nrows        5
    xllcorner    440720.000000000000
    yllcorner    3745320.000000000000
    cellsize     1200.000000000000
    NODATA_value 0
     0 0 0 0 0
     0 183 184 183 0
     0 180 182 181 0
     0 181 181 177 0
     0 0 0 0 0
    """
rm("data/utmtiny-hillshade.asc")
rm("data/utmtiny-hillshade.prj")

# GDALNearblack
# not ascii because it doesn't support `create`
options = GDAL.nearblackoptionsnew(["-of","GTiff","-color","0"], C_NULL)
ds_nearblack = GDAL.nearblack("data/utmtiny-nearblack.tif", Ptr{GDAL.GDALDatasetH}(C_NULL), ds_tiny, options, C_NULL)
GDAL.nearblackoptionsfree(options)
# do rasterio to read result
band = GDAL.getrasterband(ds_nearblack, 1)
read_data = zeros(tinydata)
GDAL.rasterio(band, GDAL.GF_Read, 0, 0, 5, 5,
              read_data, 5, 5, GDAL.GDT_Byte, 0, 0)
nbdata = read_data'
# by default the outer two pixels are zeroed
nbexpected = zeros(tinydata)
nbexpected[3,3] = tinydata[3,3]
@test nbdata == nbexpected

GDAL.close(ds_nearblack)
rm("data/utmtiny-nearblack.tif")

# GDALGrid
options = GDAL.gridoptionsnew(["-of","MEM","-outsize","3","10","-txe","100","100.3","-tye","0","0.1"], C_NULL)
ds_grid = GDAL.grid("data/point-grid.mem", ds_point, options, C_NULL)
GDAL.gridoptionsfree(options)
geotransform = zeros(6)
GDAL.getgeotransform(ds_grid, geotransform)
@test geotransform ≈ [100.0,0.1,0.0,0.0,0.0,0.01]
GDAL.close(ds_grid)

# GDALRasterize
options = GDAL.rasterizeoptionsnew(["-of","MEM","-tr","0.05","0.05"], C_NULL)
ds_rasterize = GDAL.rasterize("data/point-rasterize.mem", Ptr{GDAL.GDALDatasetH}(C_NULL), ds_point, options, C_NULL)
GDAL.rasterizeoptionsfree(options)
geotransform = zeros(6)
GDAL.getgeotransform(ds_rasterize, geotransform)
@test geotransform ≈ [99.975,0.05,0.0,0.1143,0.0,-0.05]
GDAL.close(ds_rasterize)

# GDALBuildVRT
options = GDAL.buildvrtoptionsnew([], C_NULL)
ds_buildvrt = GDAL.buildvrt("data/utmtiny.vrt", 1, [ds_tiny], C_NULL, options, C_NULL)
GDAL.buildvrtoptionsfree(options)
GDAL.close(ds_buildvrt)
@test contains(readstring("data/utmtiny.vrt"), "<SourceFilename relativeToVRT=\"1\">utmtiny.tif</SourceFilename>")
rm("data/utmtiny.vrt")

# GDALVectorTranslate
# convert a GeoJSON to CSV with X and Y columns
options = GDAL.vectortranslateoptionsnew(["-f","CSV","-lco","GEOMETRY=AS_XY"], C_NULL)
ds_csv = GDAL.vectortranslate("data/point.csv", Ptr{GDAL.GDALDatasetH}(C_NULL), 1, [ds_point], options, C_NULL)
GDAL.vectortranslateoptionsfree(options)
GDAL.close(ds_csv)
dstcsv = """
    X,Y,FID,pointname
    100,0,2,point-a
    100.2785,0.0893,3,point-b
    100,0,0,a
    100.2785,0.0893,3,b
    """
@test replace(readstring("data/point.csv"), "\r", "") == dstcsv
rm("data/point.csv")

GDAL.close(ds_small)
GDAL.close(ds_tiny)
GDAL.close(ds_point)

rm("data/utmtiny.asc.aux.xml")
rm("data/utmtiny.prj")
rm("data/utmtiny.asc")
rm("data/utmtiny.tif")


text:
using QuantumESPRESSO
using Test

@testset "QuantumESPRESSO.jl" begin
end


text:
module test_octo_predicates

using Test
using Octo.Adapters.SQL # from to_sql
using Octo: Enclosed

@test to_sql([Enclosed(Any[1, true])]) == "(1, true)"

struct User
end

u = from(User)
enc = Enclosed([])
@test to_sql([u.val == enc]) == "val = ()"

@test to_sql([u.name IS NOT NULL]) == "name IS NOT NULL"
@test to_sql([30 >= u.age AND u.age >= 20]) == "30 >= age AND age >= 20"

struct PGStatioUserTable
end
Schema.model(PGStatioUserTable, table_name="pg_statio_user_tables")
p = from(PGStatioUserTable)
q = [SELECT
     (as(SUM(p.heap_blks_read), :heap_read), as(SUM(p.heap_blks_hit), :heap_hit), as(SUM(p.heap_blks_hit) / (SUM(p.heap_blks_hit) + SUM(p.heap_blks_read)), :ratio))
     FROM
     p]
@test to_sql(q) == """SELECT SUM(heap_blks_read) AS heap_read, SUM(heap_blks_hit) AS heap_hit, SUM(heap_blks_hit) / (SUM(heap_blks_hit) + SUM(heap_blks_read)) AS ratio FROM pg_statio_user_tables"""

q = [SELECT
     (p.relname, as(100 * p.idx_scan / (p.seq_scan + p.idx_scan), :percent_of_times_index_used), as(p.n_live_tup, :rows_in_table))
     FROM
     p
     WHERE
     p.seq_scan + p.idx_scan > 0
     ORDER
     BY
     p.n_live_tup
     DESC]
@test to_sql(q) == "SELECT relname, 100 * idx_scan / (seq_scan + idx_scan) AS percent_of_times_index_used, n_live_tup AS rows_in_table FROM pg_statio_user_tables WHERE (seq_scan + idx_scan) > 0 ORDER BY n_live_tup DESC"

@test to_sql([u.name != NULL]) == "name != NULL"

end # module test_octo_predicates


text:
"""
    Bond

Defines a bond defintion in lattice.

# Fields
$(TYPEDFIELDS)
"""
struct Bond

    "Dimension of system bond is in."
    D::Int

    "Initial and final orbital species respectively."
    orbitals::Vector{Int}

    "Displacement in unit cells."
    displacement::Vector{Int}
end

"""
    Bond(orbitals::AbstractVector{Int}, displacement::AbstractVector{Int})

Constrcut a [`Bond`](@ref)
"""
function Bond(orbitals::AbstractVector{Int},displacement::AbstractVector{Int})

    D = length(displacement)
    return  Bond(D,orbitals,displacement)
end

Bond(; orbitals, displacement) = Bond(orbitals,displacement)

"""
    Base.show(io::IO, bond::Bond)
    Base.show(io::IO, ::MIME"text/plain", bond::Bond)

Show lattice.
"""
Base.show(io::IO, bond::Bond) = print(io,"Bond(", (bond.D),", ",bond.o,", ",bond.Δl)
function Base.show(io::IO, ::MIME"text/plain", bond::Bond)

    (; D, orbitals, displacement) = bond
    print(io, "Bond:\n")
    print(io, " • D = ", D, "\n")
    print(io, " • orbitals = ", orbitals, "\n")
    print(io, " • displacement = ", displacement)
    return nothing
end

"""
    Base.:==(b₁::Bond,b₂::Bond)

Tests if two bonds `b₁` and `b₂` are equivalent.
"""
function Base.:(==)(b₁::Bond,b₂::Bond)

    D₁  = b₁.D
    Δl₁ = b₁.displacement
    o₁  = b₁.orbitals
    D₂  = b₂.D
    Δl₂ = b₂.displacement
    o₂  = b₂.orbitals

    equal = false
    if D₁==D₂
        if Δl₁==Δl₂ && o₁==o₂
            equal = true
        elseif all(Δl₁[i]==-Δl₂[i] for i in 1:D₁) && o₁[1]==o₂[2] && o₁[2]==o₂[1]
            equal = true
        end
    end

    return equal
end


text:
@testset "SVDD_util_test" begin

    @testset "merge pool" begin
        pools = labelmap([:U, :Lin, :Lout, :U, :U, :Lin])

        @test [] == SVDD.merge_pools(pools, [])

        @test [1,4,5] == SVDD.merge_pools(pools, :U)
        @test [1,4,5] == sort(SVDD.merge_pools(pools, :U, :U))
        @test [2,6] == SVDD.merge_pools(pools, :Lin)
        @test [3] == SVDD.merge_pools(pools, :Lout)

        @test [2,3,6] == sort(SVDD.merge_pools(pools, :Lin, :Lout))
        @test collect(1:6) == sort(SVDD.merge_pools(pools, :U, :Lin, :Lout))

        @test_throws ArgumentError SVDD.merge_pools(pools, :A)
    end

    @testset "classify" begin
        scores = [2.0, 0.0, 0.01, -0.1, -5.0]
        expected = [:outlier, :inlier, :outlier, :inlier, :inlier]
        actual = SVDD.classify.(scores)
        @test MLLabelUtils.islabelenc(actual, SVDD.class_label_enc)
        @test expected == actual
    end

    @testset "adjust kernel" begin
        K1 = [2 -1 0; -1 2 -1; 0 -1 2]
        @assert all(eigen(K1).values .> 0)
        @test K1 ≈ SVDD.adjust_kernel_matrix(K1)

        K2 = [1 2; 2 1]
        @assert any(eigen(K2).values .< 0)
        K2_adjusted = SVDD.adjust_kernel_matrix(K2, warn_threshold = 2)
        @test !(K2 ≈ K2_adjusted)
        @test all(eigen(K2_adjusted).values .>= 0.0)

        Random.seed!(42)
        dummy_data, _ = generate_mvn_with_outliers(2, 100, 42, true, true)
        model = SVDD.VanillaSVDD(dummy_data)
        init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(4), 0.1)
        SVDD.initialize!(model, init_strategy)
        K_old = copy(model.K)
        @assert any(eigen(model.K).values .< 0)

        @test_throws ArgumentError SVDD.adjust_kernel_matrix([3 -2; 4 -1])
    end
end


text:
function initiate_list(data, p)
    first_atom = zeros(Int64,data.nc,data.nc)
    next_atom = zeros(Int64, length(p))
    return first_atom, next_atom
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
using Test

#=
To test with only one solver:
julia> using Pkg
julia> Pkg.test("REoptLite"; test_args=["Cbc"])

nlaws 200721: only running Cbc tests here b/c cannot get CPLEX and Xpress licences on to Github
    servers (runtests.jl is automated in Github Actions with ci.yml).

TODO: combine tests into one file and pass in Solver
=#

@testset "REoptLite.jl" begin
    # if isempty(ARGS) || "all" in ARGS
    #     all_tests = true
    # else
    #     all_tests = false
    # end
    # if all_tests || "CPLEX" in ARGS
    #     @testset "test_with_cplex" begin
    #         include("test_with_cplex.jl")
    #     end
    # end
    # if all_tests || "Xpress" in ARGS
    #     @testset "test_with_xpress" begin
    #         include("test_with_xpress.jl")
    #     end
    # end
    # if all_tests || "Cbc" in ARGS
        @testset "test_with_cbc" begin
            include("test_with_cbc.jl")
        end
    # end
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
# prediction from single tree - assign each observation to its final leaf
function predict_gpu!(pred::AbstractMatrix{T}, tree::Tree_gpu{T,S}, X::AbstractMatrix) where {T,S}
    @inbounds @threads for i in 1:size(X,1)
        K = length(tree.nodes[1].pred)
        id = 1
        x = view(X, i, :)
        @inbounds while tree.nodes[id].split
            if x[tree.nodes[id].feat] < tree.nodes[id].cond
                id = tree.nodes[id].left
            else
                id = tree.nodes[id].right
            end
        end
        @inbounds for k in 1:K
            pred[i,k] += tree.nodes[id].pred[k]
        end
        # @views performance much improived in Julia 1.5
        # @views pred[i,:] .+= tree.nodes[id].pred
    end
end

# prediction from single tree - assign each observation to its final leaf
function predict_gpu(tree::Tree_gpu{T,S}, X::AbstractMatrix, K) where {T,S}
    pred = zeros(T, size(X, 1), K)
    predict_gpu!(pred, tree, X)
    return pred
end

# prediction from single tree - assign each observation to its final leaf
function predict_gpu(model::GBTree_gpu{T,S}, X::AbstractMatrix) where {T,S}
    K = length(model.trees[1].nodes[1].pred)
    pred = zeros(T, size(X, 1), K)
    for tree in model.trees
        predict_gpu!(pred, tree, X)
    end
    if typeof(model.params.loss) == Logistic
        @. pred = sigmoid(pred)
    elseif typeof(model.params.loss) == Poisson
        @. pred = exp(pred)
    elseif typeof(model.params.loss) == Gaussian
        pred[:,2] = exp.(pred[:,2])
    elseif typeof(model.params.loss) == Softmax
        pred = transpose(reshape(pred, model.K, :))
        for i in 1:size(pred,1)
            pred[i,:] .= softmax(pred[i,:])
        end
    end
    return pred
end

# prediction in Leaf - GradientRegression
function pred_leaf_gpu(loss::S, node::TrainNode_gpu{T}, params::EvoTypes, δ²) where {S<:GradientRegression,T}
    - params.η .* node.∑δ ./ (node.∑δ² .+ params.λ .* node.∑𝑤)
end

# prediction in Leaf - GaussianRegression
function pred_leaf_gpu(loss::S, node::TrainNode_gpu{T}, params::EvoTypes, δ²) where {S<:GaussianRegression,T}
    - params.η .* node.∑δ ./ (node.∑δ² .+ params.λ .* node.∑𝑤)
end


text:
# things that don't directly have anything to do with this package but are necessary

@generated unsafe_initialize(::Type{T}, args::Tuple) where {T} = Expr(:splatnew, :T, :args)

# Val wraps the number of axes to retain
naxes(A::AbstractArray, v::Val) = naxes(axes(A), v)
naxes(axs::Tuple, v::Val{N}) where {N} = _naxes(axs, N)
@inline function _naxes(axs::Tuple, i::Int)
    if i === 0
        return ()
    else
        return (first(axs), _naxes(tail(axs), i - 1)...)
    end
end

@inline function _naxes(axs::Tuple{}, i::Int)
    if i === 0
        return ()
    else
        return (SimpleAxis(1), _naxes((), i - 1)...)
    end
end

known_offset1(::Type{T}) where {T} = first(known_offsets(T))

offset1(::Type{T}) where {T} = first(offsets(T))

function same_root_offset(::Type{A}, ::Type{I}) where {A<:SimpleAxis,I}
    offset_axis1(A) === offset_axis1(I)
end
function same_root_offset(::Type{A}, ::Type{I}) where {A<:AbstractAxis,I}
    return same_root_offset(parent_type(A), I)
end

is_dynamic(x) = is_dynamic(typeof(x))
function is_dynamic(::Type{T}) where {T}
    if can_change_size(T) || ismutable(T)
        return true
    elseif parent_type(T) <: T
        return false
    else
        return is_dynamic(parent_type(T))
    end
end

function is_dynamic(::Type{T}) where {K,Ks,T<:Axis{K,Ks}}
    return is_dynamic(Ks) || is_dynamic(parent_type(T))
end

function assign_indices(axis, inds)
    if can_change_size(axis) && !((known_length(inds) === nothing) || known_length(inds) === known_length(axis))
        return unsafe_reconstruct(axis, inds)
    else
        return axis
    end
end

function throw_offset_error(@nospecialize(axis))
    throw("Cannot wrap axis $axis due to offset of $(first(axis))")
end


text:
module NDPrototype
using LightGraphs
using OrderedCollections
using Unrolled

include("utils.jl")
include("edge_coloring.jl")
include("component_functions.jl")
include("network_structure.jl")
include("construction.jl")
include("coreloop.jl")

end


text:
"""
    q = query_reference2image(sectiondataset, x, y, z)

Given a `sectiondataset`, find a point in a 2d slice image closest to
the specified `x`, `y`, `z` location in the reference atlas (specified
in microns).
"""
function query_reference2image(sectiondataset::Dict, x, y, z)
    id = sectiondataset["id"]
    refspace = sectiondataset["reference_space_id"]
    xr, yr, zr = inmicrons(x), inmicrons(y), inmicrons(z)
    str = "http://api.brain-map.org/api/v2/reference_to_image/$refspace.json?x=$xr&y=$yr&z=$zr&section_data_set_ids=$id"
    rq = HTTP.request("GET", str)
    data = JSON.parse(String(rq.data))
    data["msg"][1]["image_sync"]
end

"""
    p, i, r = query_image2reference(sectionimageid, x, y)

Given a particular slice `sectionimageid` and pixel location `x`, `y`,
return the posterior, inferior, right coordinates `p`, `i`, `r` in microns.
"""
function query_image2reference(sectionimageid::Integer, x, y)
    str = "http://api.brain-map.org/api/v2/image_to_reference/$sectionimageid.json?x=$x&y=$y"
    rq = HTTP.request("GET", str)
    data = JSON.parse(String(rq.data))
    ret = data["msg"]["image_to_reference"]
    ret["x"], ret["y"], ret["z"]
end

"""
    sectionimage(secdata, w, h, filename; downsample=0)

Download a section image specified by `secdata`, where `secdata` comes
from [`query_reference2image`](@ref). Width `w` and height `h` are expressed
in pixels prior to downsampling (shrinking horizontal and vertical by
a factor of `2^downsample`). Image is saved to `filename`.
"""
function sectionimage(secdata::Dict, w, h, filename; downsample=0)
    id = secdata["section_image_id"]
    x, y = secdata["x"], secdata["y"]
    download("http://api.brain-map.org/api/v2/image_download/$id?downsample=$downsample&left=$x&top=$y&width=$w&height=$h", filename)
end

"""
    splice_sectionimages(dirname, Δplane)

Combine a series of section images into a 3d volume.  Adjacent planes
will be spaced by `Δplane`, expressed in microns. The resolution
within the slice plane is determined by the raw images.

Uses interpolation along the slicing axis where appropriate, but
otherwise this takes a conservative approach and marks pixels with NaN
whenever data are missing. Data can be missing for at least two
reasons:

* Adjacent available slices are separated by more than `2Δplane`, so
  that some intermediate plane is more distant than `Δplane` from any
  available image;

* Some individual images do not cover the entire region expressed by
  the within-plane components of the BoundingBox used to extract the
  slices (stored in the metadata file in `dirname`). This can happen
  when the relevant structure lies near the tissue boundary.

## Note

Currently this assumes that the image is in PIR coordinates and that
these are coronal slices, so that the first axis is the slicing axis.
"""
function splice_sectionimages(dirname, Δplane)
    cd(dirname) do
        md = load("metadata.jld2")
        pos = md["sectionposition"]
        ds = md["downsample"]
        fn = md["filenames"]
        bb = md["boundingbox"]
        res = md["sectiondataset"]["section_images"][1]["resolution"]*μm*2^ds
        ivx, ivy, ivz = bb.intervals               # in "Allen Coordinates"
        xrange = minimum(ivx):Δplane:maximum(ivx)  # the slice axis (assumes coronal)
        sz2d = map(x->round(Int, x/res), IntervalSets.width.((ivy, ivz)))
        img1 = load(first(fn))
        frames = Vector{typeof(img1)}(length(fn))
        frames[1] = img1
        for i = 2:length(fn)
            frames[i] = load(fn[i])
        end
        C = base_colorant_type(eltype(img1)){Float32}
        img = Array{C}(sz2d..., length(xrange))
        fill!(img, nan(RGB{Float32}))
        computeframes!(img, bb, xrange, frames, pos, ds)
        imgaa = AxisArray(img, (:I, :R, :P), (res, res, Δplane))
        save("spliced.nhdr", imgaa; props=Dict("datafile"=>"spliced.raw"))
    end
end

function computeframes!(img, bb, xrange, frames, pos, ds)
    slicepos = [p[1]*μm for p in pos]
    for (j, x) in enumerate(xrange)
        # Find the nearest slices to this particular "x" (slice) position
        im = findfirst(y->y>=x, slicepos) # nearest slice "below"
        ip = im+1                         # nearest slice "above"
        # The available data might be too far away to reasonably
        # interpolate, in which case it's better to leave this plane
        # as NaN to mark it as missing data.
        m_inrange = isinrange(im, slicepos, x, step(xrange))
        p_inrange = isinrange(ip, slicepos, x, step(xrange))
        m_inrange || p_inrange || continue
        if !m_inrange
            computeframe!(view(img, :, :, j), bb, ds, frames[ip], pos[ip])
        elseif !p_inrange
            computeframe!(view(img, :, :, j), bb, ds, frames[im], pos[im])
        else  # both are in range, use linear interpolation along the slice dimension
            f = (x - slicepos[im])/(slicepos[ip] - slicepos[im])
            computeframe!(view(img, :, :, j), bb, ds, frames[im], pos[im], frames[ip], pos[ip], f)
        end
    end
    img
end

"Copy frame data, compensating if needed for a shift in the slice plane"
function computeframe!(dest, bb, downsample, src, pos)
    sa = extrapolant(eltype(dest), bb, downsample, src, pos)
    for i in CartesianRange(indices(dest))
        dest[i] = sa[i]
    end
    dest
end

"Interpolate frame data, compensating if needed for a shift in the slice plane"
function computeframe!(dest, bb, downsample, src1, pos1, src2, pos2, f)
    sa1 = extrapolant(eltype(dest), bb, downsample, src1, pos1)
    sa2 = extrapolant(eltype(dest), bb, downsample, src2, pos2)
    for i in CartesianRange(indices(dest))
        dest[i] = (1-f)*sa1[i] + f*sa2[i]
    end
    dest
end

"""
Correct for an integer-pixel shift in the frame. We could use subpixel
interpolation, but given that the resolution in the slice plane is so
much better than the resolution along the slice axis, there seems
little point in being more sophisticated.
"""
function extrapolant(::Type{T}, bb, downsample, src, pos) where T
    d1, d2 = inmicrons(minimum(bb.intervals[2])), inmicrons(minimum(bb.intervals[3]))
    s1, s2 = pos[2], pos[3]
    mag = 2^downsample
    Δ1, Δ2 = round(Int, (s1-d1)/mag), round(Int, (s2-d2)/mag)
    inds1, inds2 = offsetindex(indices(src, 1), Δ1), offsetindex(indices(src, 2), Δ2)
    extrapolate(interpolate(view(src, inds1, inds2), NoInterp(), OnGrid()), nan(T))
end

offsetindex(inds, Δ) = OffsetArray(inds, inds+Δ)

function isinrange(idx, pos, x, Δx)
    1 <= idx <= length(pos) || return false
    abs(pos[idx] - x) <= Δx
end

function imagesize(filename)
    wand = ImageMagick.MagickWand()
    ImageMagick.readimage(wand, filename)
    size(wand)
end