text:
include("csv_parser.jl")

function methods_return(selected_method::String = "None")
    # Returns an Array of all methods in the database directory if input is "None"
    # If input is any other string, it will return an Array of just that string
    if selected_method == "None"
        return filter(x-> isdir(joinpath(dirname(pathof(OpenSAFT)), "../database", x)), readdir(joinpath(dirname(pathof(OpenSAFT)), "../database")))
    else
        if selected_method in filter(x-> isdir(joinpath(dirname(pathof(OpenSAFT)), "../database",x)), readdir(joinpath(dirname(pathof(OpenSAFT)), "../database")))
            return [selected_method]
        else
            error("Database for selected method $selected_method not found.")
        end
    end
end

function createfilepath(selected_method, database_type; variant="None")
    return joinpath(dirname(pathof(OpenSAFT)), "../database", selected_method, variant=="None" ? "data_$(selected_method)_$(database_type).csv" : "$variant/data_$(selected_method)_$(variant)_$(database_type).csv")
end

function customdatabase_check(selected_method, database_type, customdatabase_filepath; variant="None")
    if !isfile(customdatabase_filepath)
        error("Custom database $customdatabase_filepath not found.")
    elseif selected_method == "None"
        error("Method has to be specified for custom databases.")
    end
    database_filepath = createfilepath(selected_method, database_type; variant=variant)
    database_header = parseline(database_filepath, 3)
    customdatabase_header = parseline(customdatabase_filepath, 3)
    if Set(database_header) != Set(customdatabase_header)
        error("Custom database $customdatabase_filepath header not consistent with format in $selected_method (variant: $variant) for type $database_type.")
    end
end

function searchdatabase_like(components::Array{String, 1}, selected_method="None"; customdatabase_filepath="None", variant="None")
    # Returns a dictionary of components containing another dictionary
    # where the key is a method in which its database contains that parameter,
    # with the found line number as value
    if customdatabase_filepath != "None"
        customdatabase_check(selected_method, "like", customdatabase_filepath, variant=variant)
    end
    methods = methods_return(selected_method)
    found_methods = Dict(Set([component]) => Dict{String, Int64}() for component in components)
    for component in components
        for method in methods
            filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "like"; variant=variant) : customdatabase_filepath
            if isfile(filepath)
                found_lines = findmatches(filepath, component, "species"; header_row=3)
                if length(found_lines) > 1
                    println("The species $component is not unique in database for method $method. Selecting most recent entry.")
                elseif length(found_lines) == 0
                    if length(found_lines) == 0
                        if selected_method != "None"
                            error("The species $component cannot be found in database for method $method.")
                        end
                    end
                end
                if length(found_lines) != 0
                    found_methods[Set([component])][method] = found_lines[end]
                end
            end
        end
    end
    return found_methods
end

function searchdatabase_unlike(components::Array{String, 1}, selected_method="None"; customdatabase_filepath="None", variant="None")
    # Returns a dictionary of pairs (no same pair) containing another dictionary
    # where the key is a method in which its database contains that parameter,
    # with the found line number as value
    if customdatabase_filepath != "None"
        customdatabase_check(selected_method, "unlike", customdatabase_filepath, variant=variant)
    end
    methods = methods_return(selected_method)
    pairs = [Set(i) for i in collect(combinations(components, 2))]
    found_methods = Dict{Set{String}, Dict{String, Int64}}()
    for pair in pairs
        for method in methods
            filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "unlike"; variant=variant) : customdatabase_filepath
            if isfile(filepath)
                components = [i for i in pair]
                found_lines = findmatches_pair(filepath, components[1], components[2], "species1", "species2"; header_row=3)
                if length(found_lines) > 1
                    println("The pair (" * join([i for i in pair], ", ") * ") is not unique in database for method " * method * ". Selecting most recent entry.")
                end
                if length(found_lines) != 0
                    found_methods[pair] = Dict()
                    found_methods[pair][method] = found_lines[end]
                end
            end
        end
    end
    return found_methods
end

function searchdatabase_assoc(components::Array{String, 1}, selected_method="None"; customdatabase_filepath="None", variant="None")
    # Returns a dictionary of pairs (same pair allowed) containing another dictionary
    # where the key is a method in which its database contains that parameter,
    # with the found line number as value
    if customdatabase_filepath != "None"
        customdatabase_check(selected_method, "assoc", customdatabase_filepath, variant=variant)
    end
    methods = methods_return(selected_method)
    pairs = vcat([Set([i]) for i in components], [Set(i) for i in collect(Combinatorics.combinations(components, 2))])
    found_methods = Dict{Set{String}, Dict{String, Array{Int64,1}}}()
    for pair in pairs
        for method in methods
            filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "assoc"; variant=variant) : customdatabase_filepath
            if isfile(filepath)
                if length(pair) == 1
                    component = [i for i in pair]
                    components = [component component]
                else
                    components = [i for i in pair]
                end
                found_lines = findmatches_pair(filepath, components[1], components[2], "species1", "species2"; header_row=3)
                if length(found_lines) > 0
                    found_methods[pair] = Dict()
                    found_methods[pair][method] = found_lines
                end
            end
        end
    end
    return found_methods
end

function retrieveparams_like(components::Array{String, 1}, selected_method; customdatabase_filepath="None", variant="None", redirect="None")
    # Returns a dictionary of all columns in the like database for a selected method
    # for each component
    filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "like"; variant=variant) : customdatabase_filepath
    if redirect != "None"
        selected_method = redirect
    end
    header = parseline(filepath, 3)
    found_method = searchdatabase_like(components, selected_method; customdatabase_filepath=customdatabase_filepath, variant=variant)
    found_params = Dict{Set{String}, Dict{String, Any}}()
    for component in keys(found_method)
        found_params[component] = Dict(zip(header, parseline(filepath, found_method[component][selected_method])))
    end
    return found_params
end

function retrieveparams_unlike(components::Array{String, 1}, selected_method; customdatabase_filepath="None", variant="None", redirect="None")
    # Returns a dictionary of all columns in the unlike database for a selected method
    # for each pair
    filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "unlike"; variant=variant) : customdatabase_filepath
    if redirect != "None"
        selected_method = redirect
    end
    header = parseline(filepath, 3)
    found_method = searchdatabase_unlike(components, selected_method; customdatabase_filepath=customdatabase_filepath, variant=variant)
    found_params = Dict{Set{String}, Dict{String, Any}}()
    for pair in keys(found_method)
        found_params[pair] = Dict(zip(header, parseline(filepath, found_method[pair][selected_method])))
    end
    return found_params
end

function retrieveparams_assoc(components::Array{String, 1}, selected_method; customdatabase_filepath="None", variant="None", redirect="None")
    # Returns a dictionary of all columns in the assoc database for a selected method
    # for each pair (same pair allowed)
    filepath = customdatabase_filepath == "None" ? createfilepath(selected_method, "assoc"; variant=variant) : customdatabase_filepath
    if redirect != "None"
        selected_method = redirect
    end
    header = parseline(filepath, 3)
    found_method = searchdatabase_assoc(components, selected_method; customdatabase_filepath=customdatabase_filepath, variant=variant)
    pairs = keys(found_method)
    found_params = Dict{Set{Tuple{Set{String}, String}}, Dict{String, Any}}()
    if !isempty(found_method)
        for pair in pairs
            for line_number in found_method[pair][selected_method]
                retrieved = Dict(zip(header, parseline(filepath, line_number)))
                assoc_pair = Set([(Set([retrieved["species1"]]), retrieved["site1"]), (Set([retrieved["species2"]]), retrieved["site2"])])
                found_params[assoc_pair] = retrieved
            end
        end
        # Check if reverse pair exists; if not, make create reverse pair equal to original pair
        for assoc_pair in keys(found_params)
            if length(assoc_pair) == 2
                pair_tuples = collect(assoc_pair)
                reverse_assoc_pair = Set([(pair_tuples[1][1], pair_tuples[2][2]), (pair_tuples[2][1], pair_tuples[1][2])])
                if !haskey(found_params, reverse_assoc_pair)
                    found_params[reverse_assoc_pair] =
                        found_params[assoc_pair]
                end
            end
        end
    end
    return found_params
end

function retrieveparams(components::Array{String, 1}, selected_method;
        customdatabase_like = "None", variant_like = "None", redirect_like = "None",
        customdatabase_unlike = "None", variant_unlike = "None", redirect_unlike = "None",
        customdatabase_assoc = "None", variant_assoc = "None", redirect_assoc = "None")
    params_like = retrieveparams_like(components, selected_method; customdatabase_filepath = customdatabase_like, variant = variant_like, redirect = redirect_like)
    params_unlike = retrieveparams_unlike(components, selected_method; customdatabase_filepath = customdatabase_unlike, variant = variant_unlike, redirect = redirect_unlike)
    params_assoc = retrieveparams_assoc(components, selected_method; customdatabase_filepath = customdatabase_assoc, variant = variant_assoc, redirect = redirect_assoc)
    return [params_like, params_unlike, params_assoc]
end

function filterparams(raw_params, like_params::T; unlike_params::T=Array{String,1}([]), assoc_params::T=Array{String,1}([])) where T<:Array{String,1}
    # Filters the raw parametrs from retrieveparameters into selected headers
    # Returns dictionaries where the keys are the header of the selected columns
    # One dictionary for each like, unlike, and assoc
    raw_params_like = raw_params[1]
    raw_params_unlike = raw_params[2]
    raw_params_assoc = raw_params[3]
    components = keys(raw_params_like)
    pairs = keys(raw_params_unlike)
    assoc_pairs = keys(raw_params_assoc)
    like_params_dict = Dict{String, Dict{Set{String}, Float64}}()
    unlike_params_dict = Dict{String, Dict{Set{String}, Float64}}()
    assoc_params_dict = Dict{String, Dict{Set{Tuple{Set{String}, String}}, Float64}}()
    for like_param in like_params
        like_params_dict[like_param] = Dict()
        for component in components
            param_value = raw_params_like[component][like_param]
            if !ismissing(param_value)
                push!(like_params_dict[like_param], component => param_value)
            end
        end
    end
    for unlike_param in unlike_params
        unlike_params_dict[unlike_param] = Dict()
        for pair in pairs
            if haskey(raw_params_unlike[pair], unlike_param)
                param_value = raw_params_unlike[pair][unlike_param]
                if !ismissing(param_value)
                    push!(unlike_params_dict[unlike_param], pair => param_value)
                end
            end
        end
    end
    for assoc_param in assoc_params
        assoc_params_dict[assoc_param] = Dict()
        for assoc_pair in assoc_pairs
            param_value = raw_params_assoc[assoc_pair][assoc_param]
            if !ismissing(param_value)
                push!(assoc_params_dict[assoc_param], assoc_pair => param_value)
            end
        end
    end
    return like_params_dict, unlike_params_dict, assoc_params_dict
end


text:
module TestZeroOneRollCall
    using Base.Test
    import RollCallDataIO

    encoding = RollCallDataIO.RollCallEncoding()
    votes = [1 2 3; 4 5 6; 7 8 9; 0 0 0]
    roll_call = RollCallDataIO.RollCall(votes, encoding)

    @test RollCallDataIO.transcode(1, encoding) === 1.0
    @test RollCallDataIO.transcode(2, encoding) === 1.0
    @test RollCallDataIO.transcode(3, encoding) === 1.0
    @test RollCallDataIO.transcode(4, encoding) === 0.0
    @test RollCallDataIO.transcode(5, encoding) === 0.0
    @test RollCallDataIO.transcode(6, encoding) === 0.0
    @test isnan(RollCallDataIO.transcode(7, encoding))
    @test isnan(RollCallDataIO.transcode(8, encoding))
    @test isnan(RollCallDataIO.transcode(9, encoding))
    @test isnan(RollCallDataIO.transcode(0, encoding))

    new_roll_call = convert(RollCallDataIO.ZeroOneRollCall, roll_call)
    @test isa(new_roll_call, RollCallDataIO.ZeroOneRollCall)

    new_votes = new_roll_call.votes
    @test size(new_votes) === (4, 3)

    @test new_votes[1, 1] === 1.0
    @test new_votes[1, 2] === 1.0
    @test new_votes[1, 3] === 1.0

    @test new_votes[2, 1] === 0.0
    @test new_votes[2, 2] === 0.0
    @test new_votes[2, 3] === 0.0

    @test isnan(new_votes[3, 1])
    @test isnan(new_votes[3, 2])
    @test isnan(new_votes[3, 3])

    @test isnan(new_votes[4, 1])
    @test isnan(new_votes[4, 2])
    @test isnan(new_votes[4, 3])
end

ext: jl

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
using Base.Test
using TestSetExtensions

test_path = splitdir(realpath(@__FILE__))[1]
# test_path = pwd()
modules_path = joinpath(splitdir(test_path)[1], "src/")
push!(LOAD_PATH, modules_path)

import FoodRecom

my_tests = ["GetData_test.jl"]

println("Running tests:")
for my_test in my_tests
  include(my_test)
end


text:
module TestBenchSmoke

using Test
using SyncBarriersBenchmarks: clear, setup_smoke

function test_bench_smoke()
    try
        local suite
        @test (suite = setup_smoke()) isa Any
        @test run(suite) isa Any
    finally
        clear()
    end
end

end  # module


text:
# Method when using Boykov-Kolmogorov as a subroutine
# Kishimoto algorithm
function kishimoto{T<:AbstractFloat}(
  flow_graph::DiGraph,                       # the input graph
  source::Int,                               # the source vertex
  target::Int,                               # the target vertex
  capacity_matrix::AbstractArray{T, 2},      # edge flow capacities
  flow_algorithm::BoykovKolmogorovAlgorithm, # keyword argument for algorithm
  routes::Int                                # keyword argument for routes
  )
  # Initialisation
  flow, F, labels = maximum_flow(flow_graph, source, target,
         capacity_matrix, algorithm = flow_algorithm)
  restriction = flow / routes
  flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
                       algorithm = flow_algorithm, restriction = restriction)

  # Loop condition : approximatively not equal is enforced by floating precision
  i = 1
  while flow < routes * restriction && flow ≉ routes * restriction
    restriction = (flow - i * restriction) / (routes - i)
    i += 1
    flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
                         algorithm = flow_algorithm, restriction = restriction)
  end

  # End
  return flow, F, labels
end

"""
Computes the maximum multiroute flow (for an integer number of routes)
between the source and target vertexes in a flow graph using the
[Kishimoto algorithm](http://dx.doi.org/10.1109/ICCS.1992.255031).
Returns the value of the multiroute flow as well as the final flow matrix,
along with a multiroute cut if Boykov-Kolmogorov is used as a subroutine.
Use a default capacity of 1 when the capacity matrix isn\'t specified.
Requires arguments:
- flow_graph::DiGraph                    # the input graph
- source::Int                            # the source vertex
- target::Int                            # the target vertex
- capacity_matrix::AbstractArray{T, 2}   # edge flow capacities
- flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
- routes::Int                            # keyword argument for routes
"""

function kishimoto{T<:AbstractFloat}(
  flow_graph::DiGraph,                   # the input graph
  source::Int,                           # the source vertex
  target::Int,                           # the target vertex
  capacity_matrix::AbstractArray{T, 2},  # edge flow capacities
  flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
  routes::Int                            # keyword argument for routes
  )
  # Initialisation
  flow, F = maximum_flow(flow_graph, source, target,
         capacity_matrix, algorithm = flow_algorithm)
  restriction = flow / routes

  flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
          algorithm = flow_algorithm, restriction = restriction)

  # Loop condition : approximatively not equal is enforced by floating precision
  i = 1
  while flow < routes * restriction && flow ≉ routes * restriction
    restriction = (flow - i * restriction) / (routes - i)
    i += 1
    flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
              algorithm = flow_algorithm, restriction = restriction)
  end

  # End
  return flow, F
end


text:
"""
    struct UnstructuredGrid{Dc,Dp,Tp,Ti,O} <: Grid{Dc,Dp}
      node_coordinates::Vector{Point{Dp,Tp}}
      cell_node_ids::Table{Ti,Int32}
      reffes::Vector{<:LagrangianRefFE{Dc}}
      cell_types::Vector{Int8}
    end
"""
struct UnstructuredGrid{Dc,Dp,Tp,O,Tn} <: Grid{Dc,Dp}
  node_coordinates::Vector{Point{Dp,Tp}}
  cell_node_ids::Table{Int32,Vector{Int32},Vector{Int32}}
  reffes::Vector{LagrangianRefFE{Dc}}
  cell_types::Vector{Int8}
  orientation_style::O
  facet_normal::Tn
  cell_map
  @doc """
      function UnstructuredGrid(
        node_coordinates::Vector{Point{Dp,Tp}},
        cell_node_ids::Table{Ti},
        reffes::Vector{<:LagrangianRefFE{Dc}},
        cell_types::Vector,
        orientation_style::OrientationStyle=NonOriented()) where {Dc,Dp,Tp,Ti}
      end

  Low-level inner constructor.
  """
  function UnstructuredGrid(
    node_coordinates::Vector{Point{Dp,Tp}},
    cell_node_ids::Table{Ti},
    reffes::Vector{<:LagrangianRefFE{Dc}},
    cell_types::Vector,
    orientation_style::OrientationStyle=NonOriented(),
    facet_normal=nothing) where {Dc,Dp,Tp,Ti}

    cell_map = _compute_cell_map(node_coordinates,cell_node_ids,reffes,cell_types)
    B = typeof(orientation_style)
    Tn = typeof(facet_normal)
    new{Dc,Dp,Tp,B,Tn}(
      node_coordinates,
      cell_node_ids,
      reffes,
      cell_types,
      orientation_style,
      facet_normal,
      cell_map)
  end
end

function _compute_cell_map(node_coords,cell_node_ids,ctype_reffe,cell_ctype)
  cell_coords = lazy_map(Broadcasting(Reindex(node_coords)),cell_node_ids)
  ctype_shapefuns = map(get_shapefuns,ctype_reffe)
  cell_shapefuns = expand_cell_data(ctype_shapefuns,cell_ctype)
  default_cell_map = lazy_map(linear_combination,cell_coords,cell_shapefuns)
  ctype_poly = map(get_polytope,ctype_reffe)
  has_affine_map =
    all(map(is_first_order,ctype_reffe)) &&
    ( all(map(is_simplex,ctype_poly)) || all(map(p->num_dims(p)==1,ctype_poly)) )
  if has_affine_map
    ctype_q0 = map(p->zero(first(get_vertex_coordinates(p))),ctype_poly)
    cell_q0 = expand_cell_data(ctype_q0,cell_ctype)
    default_cell_grad = lazy_map(∇,default_cell_map)
    origins = lazy_map(evaluate,default_cell_map,cell_q0)
    gradients = lazy_map(evaluate,default_cell_grad,cell_q0)
    cell_map = lazy_map(Fields.affine_map,gradients,origins)
  else
    cell_map = default_cell_map
  end
  Fields.MemoArray(cell_map)
end

"""
    UnstructuredGrid(grid::Grid)
"""
function UnstructuredGrid(grid::Grid)
  @assert is_regular(grid) "UnstructuredGrid constructor only for regular grids"
  node_coordinates = collect1d(get_node_coordinates(grid))
  cell_node_ids = Table(get_cell_node_ids(grid))
  reffes = get_reffes(grid)
  cell_types = collect1d(get_cell_type(grid))
  orien = OrientationStyle(grid)
  UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types,orien)
end

function UnstructuredGrid(grid::UnstructuredGrid)
  grid
end

OrientationStyle(
  ::Type{<:UnstructuredGrid{Dc,Dp,Tp,B}}) where {Dc,Dp,Tp,B} = B()

get_reffes(g::UnstructuredGrid) = g.reffes

get_cell_type(g::UnstructuredGrid) = g.cell_types

get_node_coordinates(g::UnstructuredGrid) = g.node_coordinates

get_cell_node_ids(g::UnstructuredGrid) = g.cell_node_ids

get_cell_map(g::UnstructuredGrid) = g.cell_map

function get_facet_normal(g::UnstructuredGrid)
  @assert g.facet_normal != nothing "This Grid does not have information about normals."
  g.facet_normal
end

# From ReferenceFE

"""
    UnstructuredGrid(reffe::LagrangianRefFE)

Build a grid with a single cell that is the given reference FE itself
"""
function UnstructuredGrid(reffe::LagrangianRefFE)
  node_coordinates = get_node_coordinates(reffe)
  cell_node_ids = Table([collect(1:num_nodes(reffe)),])
  reffes = [reffe,]
  cell_types = [1,]
  UnstructuredGrid(node_coordinates,cell_node_ids,reffes,cell_types)
end

# From Polytope

function UnstructuredGrid(::Type{ReferenceFE{D}},p::Polytope{D}) where D
  order = 1
  reffe = LagrangianRefFE(Float64,p,order)
  UnstructuredGrid(reffe)
end

"""
    UnstructuredGrid(::Type{ReferenceFE{d}},p::Polytope) where d
"""
function UnstructuredGrid(::Type{ReferenceFE{d}},p::Polytope) where d
  node_coordinates = get_vertex_coordinates(p)
  cell_node_ids = Table(get_faces(p,d,0))
  reffaces = get_reffaces(Polytope{d},p)
  cell_type = get_face_type(p,d)
  order = 1
  reffes = map( (f)-> LagrangianRefFE(Float64,f,order), reffaces)
  UnstructuredGrid(
    node_coordinates,
    cell_node_ids,
    reffes,
    cell_type)
end

# From coordinates

"""
    UnstructuredGrid(x::AbstractArray{<:Point})
"""
function UnstructuredGrid(x::AbstractArray{<:Point})
  np = length(x)
  node_coords = collect1d(x)
  cell_node_ids = identity_table(Int32,Int32,np)
  cell_type = fill(1,np)
  order = 1
  reffes = [LagrangianRefFE(Float64,VERTEX,order),]
  UnstructuredGrid(
    node_coords,
    cell_node_ids,
    reffes,
    cell_type)
end

# IO

function to_dict(grid::UnstructuredGrid)
  dict = Dict{Symbol,Any}()
  x = get_node_coordinates(grid)
  dict[:node_coordinates] = reinterpret(eltype(eltype(x)),x)
  dict[:Dp] = num_point_dims(grid)
  dict[:cell_node_ids] = to_dict(get_cell_node_ids(grid))
  dict[:reffes] = map(to_dict, get_reffes(grid))
  dict[:cell_type] = get_cell_type(grid)
  dict[:orientation] = is_oriented(grid)
  dict
end

function from_dict(::Type{UnstructuredGrid},dict::Dict{Symbol,Any})
  x = collect1d(dict[:node_coordinates])
  T = eltype(x)
  Dp = dict[:Dp]
  node_coordinates::Vector{Point{Dp,T}} = reinterpret(Point{Dp,T},x)
  cell_node_ids = from_dict(Table{Int32,Vector{Int32},Vector{Int32}},dict[:cell_node_ids])
  reffes = [ from_dict(LagrangianRefFE,reffe) for reffe in dict[:reffes]]
  cell_type::Vector{Int8} = dict[:cell_type]
  O::Bool = dict[:orientation]
  UnstructuredGrid(
    node_coordinates,
    cell_node_ids,
    reffes,
    cell_type,
    O ? Oriented() : NonOriented())
end

function simplexify(grid::UnstructuredGrid)
  reffes = get_reffes(grid)
  @notimplementedif length(reffes) != 1
  reffe = first(reffes)
  order = 1
  @notimplementedif get_order(reffe) != order
  p = get_polytope(reffe)
  ltcell_to_lpoints, simplex = simplexify(p)
  cell_to_points = get_cell_node_ids(grid)
  tcell_to_points = _refine_grid_connectivity(cell_to_points, ltcell_to_lpoints)
  ctype_to_reffe = [LagrangianRefFE(Float64,simplex,order),]
  tcell_to_ctype = fill(Int8(1),length(tcell_to_points))
  point_to_coords = get_node_coordinates(grid)
  UnstructuredGrid(
    point_to_coords,
    tcell_to_points,
    ctype_to_reffe,
    tcell_to_ctype,
    Oriented())
end

function _refine_grid_connectivity(cell_to_points::Table, ltcell_to_lpoints)
  data, ptrs = _refine_grid_connectivity(
    cell_to_points.data,
    cell_to_points.ptrs,
    ltcell_to_lpoints)
  Table(data,ptrs)
end

function _refine_grid_connectivity(
  cell_to_points_data::AbstractVector{T},
  cell_to_points_ptrs::AbstractVector{P},
  ltcell_to_lpoints) where {T,P}

  nltcells = length(ltcell_to_lpoints)
  ncells = length(cell_to_points_ptrs) - 1
  ntcells = ncells * nltcells

  tcell_to_points_ptrs = zeros(P,ntcells+1)

  _refine_grid_connectivity_count!(
    tcell_to_points_ptrs,
    ncells,
    ltcell_to_lpoints)

  length_to_ptrs!(tcell_to_points_ptrs)

  ndata = tcell_to_points_ptrs[end]-1

  tcell_to_points_data = zeros(T,ndata)

  _refine_grid_connectivity!(
    tcell_to_points_data,
    cell_to_points_data,
    cell_to_points_ptrs,
    ltcell_to_lpoints )

  (tcell_to_points_data, tcell_to_points_ptrs)

end

function  _refine_grid_connectivity_count!(
    tcell_to_points_ptrs,
    ncells,
    ltcell_to_lpoints)

  tcell = 1

  for cell in 1:ncells
    for lpoints in ltcell_to_lpoints
      tcell_to_points_ptrs[tcell+1] = length(lpoints)
      tcell +=1
    end
  end

end

function _refine_grid_connectivity!(
    tcell_to_points_data,
    cell_to_points_data,
    cell_to_points_ptrs,
    ltcell_to_lpoints )

  ncells = length(cell_to_points_ptrs) - 1

  k = 1
  for cell in 1:ncells
    a = cell_to_points_ptrs[cell]-1
    for lpoints in ltcell_to_lpoints
      for lpoint in lpoints
        point = cell_to_points_data[a+lpoint]
        tcell_to_points_data[k] = point
        k += 1
      end
    end
  end

end


text:
# automatically generated by the FlatBuffers compiler, do not modify

NamespaceA.NamespaceB.eval(quote

FlatBuffers.@STRUCT struct StructInNestedNS
    a::Int32
    b::Int32
end
FlatBuffers.@ALIGN(StructInNestedNS, 4)

StructInNestedNS(buf::AbstractVector{UInt8}) = FlatBuffers.read(StructInNestedNS, buf)
StructInNestedNS(io::IO) = FlatBuffers.deserialize(io, StructInNestedNS)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
# SPDX-License-Identifier: MIT

const Abstract3DTensor{T} = AbstractArray{T,3} where T

#======================================================================
Forward recursion
======================================================================#

function αrecursion!(α::AbstractMatrix{T}, π::AbstractVector{T},
                     Aᵀ::AbstractMatrix{T}, lhs::AbstractMatrix) where T
    N = size(lhs, 2)
    buffer = similar(lhs[:, 1])
    @views elmul!(α[:,1], π, lhs[:,1])
    @views for n in 2:N
        matmul!(buffer, Aᵀ, α[:,n-1])
        elmul!(α[:,n], buffer, lhs[:,n])
    end
    α
end

function αrecursion!(α::Abstract3DTensor{T}, π::AbstractVector{T},
                     Aᵀ::AbstractMatrix{T}, lhs::Abstract3DTensor{T}) where T
    N = size(lhs, 2)
    buffer = similar(α[:,1,:])
    @views elmul!(α[:,1,:], π, lhs[:,1,:])
    @views for n in 2:N
        matmul!(buffer, Aᵀ, α[:,n-1,:])
        elmul!(α[:,n,:], buffer, lhs[:,n,:])
    end
    α
end

#======================================================================
Backward recursion
======================================================================#

function βrecursion!(β::AbstractMatrix{T}, ω::AbstractVector{T},
                     A::AbstractMatrix{T}, lhs::AbstractMatrix{T}) where T
    N = size(lhs, 2)
    buffer = similar(lhs[:, 1])
    β[:, end] = ω
    @views for n in N-1:-1:1
        elmul!(buffer, β[:,n+1], lhs[:,n+1])
        matmul!(β[:,n], A, buffer)
    end
    β
end

function βrecursion!(β::Abstract3DTensor{T}, ω::AbstractVector{T},
                     A::AbstractMatrix{T}, lhs::Abstract3DTensor{T}) where T
    N = size(lhs, 2)
    buffer = fill!(similar(β[:,1,:]), one(T))
    @views elmul!(β[:,end,:], ω, buffer)
    @views for n in N-1:-1:1
        elmul!(buffer, β[:,n+1,:], lhs[:,n+1,:])
        matmul!(β[:,n,:], A, buffer)
    end
    β
end

#======================================================================
Generic forward-backward algorithm
======================================================================#

function αβrecursion(π::AbstractVector{T}, ω::AbstractVector{T},
                     A::AbstractMatrix{T}, Aᵀ::AbstractMatrix,
                     lhs::AbstractMatrix{T}) where T
    S, N = length(π), size(lhs, 2)
    α = similar(lhs, T, S, N)
    β = similar(lhs, T, S, N)

    αrecursion!(α, π, Aᵀ, lhs)
    βrecursion!(β, ω, A, lhs)

    γ = similar(lhs, T, S, N)
    elmul!(γ, α, β)
    sums = sum(γ, dims = 1)
    eldiv!(γ, γ, sums)

    γ, minimum(sums)
end

function αβrecursion(π::AbstractVector{T}, ω::AbstractVector{T},
                     A::AbstractMatrix{T}, Aᵀ::AbstractMatrix,
                     lhs::Abstract3DTensor{T}) where T
    S, N, B = size(lhs)
    α = similar(lhs, T, S, N, B)
    β = similar(lhs, T, S, N, B)

    αrecursion!(α, π, Aᵀ, lhs)
    βrecursion!(β, ω, A, lhs)

    γ = similar(lhs, T, S, N, B)
    elmul!(γ, α, β)
    sums = sum(γ, dims = 1)
    eldiv!(γ, γ, sums)

    γ, dropdims(minimum(sums, dims = (1, 2)), dims = (1, 2))
end

#======================================================================
Specialized algorithms
======================================================================#

_convert(T, ttl::Number) = convert(T, ttl)
_convert(T, ttl::AbstractVector) = copyto!(similar(ttl, T), ttl)

function stateposteriors(in_cfsm, in_lhs::AbstractArray{T}) where T <: Real
    SF = LogSemifield{T}
    cfsm = convert(CompiledFSM{SF}, in_cfsm)
    lhs = copyto!(similar(in_lhs, SF), in_lhs)
    γ, ttl = αβrecursion(cfsm.π, cfsm.ω, cfsm.A, cfsm.Aᵀ, lhs)
    out = copyto!(similar(in_lhs), γ)
    exp.(out), _convert(T, ttl)
end

function bestpath(in_cfsm, in_lhs::AbstractArray{T}) where T <: Real
    SF = TropicalSemifield{T}
    cfsm = convert(CompiledFSM{SF}, in_cfsm)
    lhs = copyto!(similar(in_lhs, SF), in_lhs)
    γ, _ = αβrecursion(cfsm.π, cfsm.ω, cfsm.A, cfsm.Aᵀ, lhs)
    out = copyto!(similar(in_lhs), γ)
    idxs = mapslices(argmax, out, dims = 1)
    dropdims(idxs, dims = 1), dropdims(lhs[idxs], dims = 1)
end

max_stars_repo_name: MartinKocour/MarkovModels

max_issues_repo_name: MartinKocour/MarkovModels

max_forks_repo_name: MartinKocour/MarkovModels


text:
if Base.libllvm_version >= v"7.0"
    include(joinpath("gcn", "math.jl"))
end
include(joinpath("gcn", "indexing.jl"))
include(joinpath("gcn", "assertion.jl"))
include(joinpath("gcn", "synchronization.jl"))
include(joinpath("gcn", "memory_static.jl"))
include(joinpath("gcn", "memory_dynamic.jl"))
include(joinpath("gcn", "hostcall.jl"))
include(joinpath("gcn", "output.jl"))


text:
# This example shows the T-intersection scenario -- nominal and failure scenario
using AdversarialDriving
using POMDPs, POMDPPolicies, POMDPSimulators

sut_agent = BlinkerVehicleAgent(up_left(id=1, s=25., v=15.), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(left_straight(id=2, s=0., v=15.0), TIDM(Tint_TIDM_template))
adv2 = BlinkerVehicleAgent(left_turnright(id=3, s=10., v=15.0), TIDM(Tint_TIDM_template))
adv3 = BlinkerVehicleAgent(right_straight(id=4, s=40., v=20.0), TIDM(Tint_TIDM_template))
adv4 = BlinkerVehicleAgent(right_turnleft(id=5, s=30., v=20.0), TIDM(Tint_TIDM_template))

mdp = AdversarialDrivingMDP(sut_agent, [adv1, adv2, adv3, adv4], Tint_roadway, 0.1)

null_action = actions(mdp)[1]
blinker_action = actions(mdp)[7]

# Nominal Behavior
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> null_action))
scenes_to_gif(state_hist(hist), mdp.roadway, "Tint_nominal.gif")

# Behavior with noise
blinker_pol = FunctionPolicy((s) -> !blinker(get_by_id(s, 2)) ? blinker_action : null_action)
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, blinker_pol)
scenes_to_gif(state_hist(hist), mdp.roadway, "Tint_failure.gif")


text:
using StanSample, MCMCChains

df = filter(row -> row[:age] >= 18,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "Howell1.csv"), DataFrame))

max_height = maximum(df[!, :height])
min_height = minimum(df[!, :height])

m4_2s = "
data {
  int N;
  real h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=$(min_height),upper=$(max_height)> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 0.1);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

m_4_2s = SampleModel("m4.2s", m4_2s);

m4_2_data = Dict("N" => size(df, 1), "h" => df[!, :height]);

rc = stan_sample(m_4_2s, data=m4_2_data);

if success(rc)
  chn = read_samples(m_4_2s; output_format=:mcmcchains)
  chn = replacenames(chn, Dict("mu" => "μ", "sigma" => "σ"))
  describe(chn)
end


text:
struct OceanModel{PS, P, T} <: AbstractOceanModel
    param_set::PS
    problem::P
    ρₒ::T
    cʰ::T
    cᶻ::T
    add_fast_substeps::T
    numImplSteps::T
    ivdc_dt::T
    αᵀ::T
    νʰ::T
    νᶻ::T
    κʰ::T
    κᶻ::T
    κᶜ::T
    fₒ::T
    β::T
    function OceanModel{FT}(
        param_set,
        problem;
        ρₒ = FT(1000),  # kg/m^3
        cʰ = FT(0),     # m/s
        cᶻ = FT(0),     # m/s
        add_fast_substeps = 0,
        numImplSteps = 0,
        ivdc_dt = FT(1),
        αᵀ = FT(2e-4),  # 1/K
        νʰ = FT(5e3),   # m^2/s
        νᶻ = FT(5e-3),  # m^2/s
        κʰ = FT(1e3),   # m^2/s
        κᶻ = FT(1e-4),  # vertical diffusivity, m^2/s
        κᶜ = FT(1e-4),  # convective adjustment vertical diffusivity, m^2/s
        fₒ = FT(1e-4),  # Hz
        β = FT(1e-11),  # Hz/m
    ) where {FT <: AbstractFloat}
        return new{typeof(param_set), typeof(problem), FT}(
            param_set,
            problem,
            ρₒ,
            cʰ,
            cᶻ,
            add_fast_substeps,
            numImplSteps,
            ivdc_dt,
            αᵀ,
            νʰ,
            νᶻ,
            κʰ,
            κᶻ,
            κᶜ,
            fₒ,
            β,
        )
    end
end

function calculate_dt(grid, ::OceanModel, _...)
    #=
      minΔx = min_node_distance(grid, HorizontalDirection())
      minΔz = min_node_distance(grid, VerticalDirection())

      CFL_gravity = minΔx / model.cʰ
      CFL_diffusive = minΔz^2 / (1000 * model.κᶻ)
      CFL_viscous = minΔz^2 / model.νᶻ

      dt = 1 // 2 * minimum([CFL_gravity, CFL_diffusive, CFL_viscous])
    =#
    # FT = eltype(grid)
    # dt = FT(1)

    return nothing
end

"""
    OceanDGModel()

helper function to add required filtering
not used in the Driver+Config setup
"""
function OceanDGModel(
    bl::OceanModel,
    grid,
    numfluxnondiff,
    numfluxdiff,
    gradnumflux;
    kwargs...,
)
    vert_filter = CutoffFilter(grid, polynomialorder(grid) - 1)
    exp_filter = ExponentialFilter(grid, 1, 8)

    flowintegral_dg = DGModel(
        FlowIntegralModel(bl),
        grid,
        numfluxnondiff,
        numfluxdiff,
        gradnumflux,
    )

    tendency_dg = DGModel(
        TendencyIntegralModel(bl),
        grid,
        numfluxnondiff,
        numfluxdiff,
        gradnumflux,
    )

    conti3d_dg = DGModel(
        Continuity3dModel(bl),
        grid,
        numfluxnondiff,
        numfluxdiff,
        gradnumflux,
    )
    FT = eltype(grid)
    conti3d_Q = init_ode_state(conti3d_dg, FT(0); init_on_cpu = true)

    ivdc_dg = DGModel(
        IVDCModel(bl),
        grid,
        numfluxnondiff,
        numfluxdiff,
        gradnumflux;
        direction = VerticalDirection(),
    )
    ivdc_Q = init_ode_state(ivdc_dg, FT(0); init_on_cpu = true) # Not sure this is needed since we set values later,
    # but we'll do it just in case!

    ivdc_RHS = init_ode_state(ivdc_dg, FT(0); init_on_cpu = true) # Not sure this is needed since we set values later,
    # but we'll do it just in case!

    ivdc_bgm_solver = BatchedGeneralizedMinimalResidual(
        ivdc_dg,
        ivdc_Q;
        max_subspace_size = 10,
    )

    modeldata = (
        dg_2D = kwargs[1][1],
        Q_2D = kwargs[1][2],
        vert_filter = vert_filter,
        exp_filter = exp_filter,
        flowintegral_dg = flowintegral_dg,
        tendency_dg = tendency_dg,
        conti3d_dg = conti3d_dg,
        conti3d_Q = conti3d_Q,
        ivdc_dg = ivdc_dg,
        ivdc_Q = ivdc_Q,
        ivdc_RHS = ivdc_RHS,
        ivdc_bgm_solver = ivdc_bgm_solver,
    )

    return DGModel(
        bl,
        grid,
        numfluxnondiff,
        numfluxdiff,
        gradnumflux;
        kwargs...,
        modeldata = modeldata,
    )
end

function vars_state(m::OceanModel, ::Prognostic, T)
    @vars begin
        u::SVector{2, T}
        η::T
        θ::T
    end
end

function init_state_prognostic!(m::OceanModel, Q::Vars, A::Vars, coords, t)
    return ocean_init_state!(m, m.problem, Q, A, coords, t)
end

function vars_state(m::OceanModel, ::Auxiliary, T)
    @vars begin
        w::T
        pkin::T         # kinematic pressure: ∫(-g αᵀ θ)
        wz0::T          # w at z=0
        u_d::SVector{2, T}  # velocity deviation from vertical mean
        ΔGu::SVector{2, T}
        y::T     # y-coordinate of the box
    end
end

function init_state_auxiliary!(
    m::OceanModel,
    state_aux::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        m,
        (m, A, tmp, geom) -> ocean_init_aux!(m, m.problem, A, geom),
        state_aux,
        grid,
        direction,
    )
end

function vars_state(m::OceanModel, ::Gradient, T)
    @vars begin
        u::SVector{2, T}
        ud::SVector{2, T}
        θ::T
    end
end

@inline function compute_gradient_argument!(
    m::OceanModel,
    G::Vars,
    Q::Vars,
    A,
    t,
)
    G.u = Q.u
    G.ud = A.u_d
    G.θ = Q.θ

    return nothing
end

function vars_state(m::OceanModel, ::GradientFlux, T)
    @vars begin
        ν∇u::SMatrix{3, 2, T, 6}
        κ∇θ::SVector{3, T}
    end
end

@inline function compute_gradient_flux!(
    m::OceanModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    ν = viscosity_tensor(m)
    #   D.ν∇u = ν * G.u
    D.ν∇u = @SMatrix [
        m.νʰ * G.ud[1, 1] m.νʰ * G.ud[1, 2]
        m.νʰ * G.ud[2, 1] m.νʰ * G.ud[2, 2]
        m.νᶻ * G.u[3, 1] m.νᶻ * G.u[3, 2]
    ]

    κ = diffusivity_tensor(m, G.θ[3])
    D.κ∇θ = κ * G.θ

    return nothing
end

@inline viscosity_tensor(m::OceanModel) = Diagonal(@SVector [m.νʰ, m.νʰ, m.νᶻ])

@inline function diffusivity_tensor(m::OceanModel, ∂θ∂z)

    if m.numImplSteps > 0
        κ = (@SVector [m.κʰ, m.κʰ, m.κᶻ * 0.5])
    else
        ∂θ∂z < 0 ? κ = (@SVector [m.κʰ, m.κʰ, m.κᶜ]) : κ = (@SVector [
            m.κʰ,
            m.κʰ,
            m.κᶻ,
        ])
    end

    return Diagonal(κ)
end

"""
    vars_integral(::OceanModel)

location to store integrands for bottom up integrals
∇hu = the horizontal divegence of u, e.g. dw/dz
"""
function vars_state(m::OceanModel, ::UpwardIntegrals, T)
    @vars begin
        ∇hu::T
        buoy::T
        #       ∫u::SVector{2, T}
    end
end

"""
    integral_load_auxiliary_state!(::OceanModel)

copy w to var_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case OceanModel
I -> array of integrand variables
Q -> array of state variables
A -> array of aux variables
"""
@inline function integral_load_auxiliary_state!(
    m::OceanModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∇hu = A.w # borrow the w value from A...
    I.buoy = grav(m.param_set) * m.αᵀ * Q.θ # buoyancy to integrate vertically from top (=reverse)
    #   I.∫u = Q.u

    return nothing
end

"""
    integral_set_auxiliary_state!(::OceanModel)

copy integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case OceanModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function integral_set_auxiliary_state!(m::OceanModel, A::Vars, I::Vars)
    A.w = I.∇hu
    A.pkin = -I.buoy
    #   A.∫u = I.∫u

    return nothing
end

"""
    vars_reverse_integral(::OceanModel)

location to store integrands for top down integrals
αᵀθ = density perturbation
"""
function vars_state(m::OceanModel, ::DownwardIntegrals, T)
    @vars begin
        buoy::T
    end
end

"""
    reverse_integral_load_auxiliary_state!(::OceanModel)

copy αᵀθ to var_reverse_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case OceanModel
I -> array of integrand variables
A -> array of aux variables
"""
@inline function reverse_integral_load_auxiliary_state!(
    m::OceanModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.buoy = A.pkin

    return nothing
end

"""
    reverse_integral_set_auxiliary_state!(::OceanModel)

copy reverse integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case OceanModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function reverse_integral_set_auxiliary_state!(
    m::OceanModel,
    A::Vars,
    I::Vars,
)
    A.pkin = I.buoy

    return nothing
end

@inline function flux_first_order!(
    m::OceanModel,
    F::Grad,
    Q::Vars,
    A::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        u = Q.u # Horizontal components of velocity
        θ = Q.θ
        w = A.w   # vertical velocity
        pkin = A.pkin
        v = @SVector [u[1], u[2], w]
        Iʰ = @SMatrix [
            1 -0
            -0 1
            -0 -0
        ]

        # ∇h • (g η)
        #- jmc: put back this term to check
        #       η = Q.η
        #       F.u += grav(m.param_set) * η * Iʰ

        # ∇ • (u θ)
        F.θ += v * θ

        # ∇h • pkin
        F.u += pkin * Iʰ

        # ∇h • (v ⊗ u)
        # F.u += v * u'
    end

    return nothing
end

@inline function flux_second_order!(
    m::OceanModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    HD::Vars,
    A::Vars,
    t::Real,
)
    # horizontal viscosity done in horizontal model
    #   F.u -= @SVector([0, 0, 1]) * D.ν∇u[3, :]'
    #- jmc: put back this term to check
    F.u -= D.ν∇u

    F.θ -= D.κ∇θ

    return nothing
end

@inline function source!(
    m::OceanModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        u = Q.u    # Horizontal components of velocity
        ud = A.u_d # Horizontal velocity deviation from vertical mean

        # f × u
        f = coriolis_force(m, A.y)
        # S.u -= @SVector [-f * u[2], f * u[1]]
        S.u -= @SVector [-f * ud[2], f * ud[1]]

        #- borotropic tendency adjustment
        S.u += A.ΔGu

        # switch this to S.η if you comment out the fast mode in MultistateMultirateRungeKutta
        S.η += A.wz0
    end

    return nothing
end

@inline coriolis_force(m::OceanModel, y) = m.fₒ + m.β * y

function update_auxiliary_state!(
    dg::DGModel,
    m::OceanModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    A = dg.state_auxiliary
    MD = dg.modeldata

    # `update_auxiliary_state!` gets called twice, once for the real elements
    # and once for the ghost elements. Only apply the filters to the real elems.
    if elems == dg.grid.topology.realelems
        # required to ensure that after integration velocity field is divergence free
        vert_filter = MD.vert_filter
        # Q[1] = u[1] = u, Q[2] = u[2] = v
        apply!(Q, (1, 2), dg.grid, vert_filter, direction = VerticalDirection())

        exp_filter = MD.exp_filter
        # Q[4] = θ
        apply!(Q, (4,), dg.grid, exp_filter, direction = VerticalDirection())
    end

    #----------
    # Compute Divergence of Horizontal Flow field using "conti3d_dg" DGmodel

    conti3d_dg = dg.modeldata.conti3d_dg
    # ct3d_Q = dg.modeldata.conti3d_Q
    # ct3d_dQ = similar(ct3d_Q)
    # fill!(ct3d_dQ, 0)
    #- Instead, use directly conti3d_Q to store dQ (since we will not update the state)
    ct3d_dQ = dg.modeldata.conti3d_Q

    ct3d_bl = conti3d_dg.balance_law

    # call "conti3d_dg" DGmodel
    # note: with "increment = false", just return tendency (no state update)
    p = nothing
    conti3d_dg(ct3d_dQ, Q, p, t; increment = false)

    # Copy from ct3d_dQ.θ which is realy ∇h•u into A.w (which will be integrated)
    function f!(::OceanModel, dQ, A, t)
        @inbounds begin
            A.w = dQ.θ
        end
    end
    update_auxiliary_state!(f!, dg, m, ct3d_dQ, t, elems)
    #----------

    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dg)

    # compute integrals for w and pkin
    indefinite_stack_integral!(dg, m, Q, A, t, elems) # bottom -> top
    reverse_indefinite_stack_integral!(dg, m, Q, A, t, elems) # top -> bottom

    # copy down wz0
    # We are unable to use vars (ie A.w) for this because this operation will
    # return a SubArray, and adapt (used for broadcasting along reshaped arrays)
    # has a limited recursion depth for the types allowed.
    nb_aux_m = number_states(m, Auxiliary())
    data_m = reshape(A.data, Nq^2, Nqk, nb_aux_m, nelemv, nelemh)

    # project w(z=0) down the stack
    index_w = varsindex(vars_state(m, Auxiliary(), FT), :w)
    index_wz0 = varsindex(vars_state(m, Auxiliary(), FT), :wz0)
    flat_wz0 = @view data_m[:, end:end, index_w, end:end, 1:nrealelemh]
    boxy_wz0 = @view data_m[:, :, index_wz0, :, 1:nrealelemh]
    boxy_wz0 .= flat_wz0

    # Compute Horizontal Flow deviation from vertical mean

    flowintegral_dg = dg.modeldata.flowintegral_dg
    flowint = flowintegral_dg.balance_law
    update_auxiliary_state!(flowintegral_dg, flowint, Q, 0, elems)

    ## get top value (=integral over full depth)
    nb_aux_flw = number_states(flowint, Auxiliary())
    data_flw = reshape(
        flowintegral_dg.state_auxiliary.data,
        Nq^2,
        Nqk,
        nb_aux_flw,
        nelemv,
        nelemh,
    )
    index_∫u = varsindex(vars_state(flowint, Auxiliary(), FT), :∫u)

    flat_∫u = @view data_flw[:, end:end, index_∫u, end:end, 1:nrealelemh]

    ## make a copy of horizontal velocity
    A.u_d .= Q.u

    ## and remove vertical mean velocity
    index_ud = varsindex(vars_state(m, Auxiliary(), FT), :u_d)
    boxy_ud = @view data_m[:, :, index_ud, :, 1:nrealelemh]
    boxy_ud .-= flat_∫u / m.problem.H

    return true
end

@inline wavespeed(m::OceanModel, n⁻, _...) =
    abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)

# We want not have jump penalties on η (since not a flux variable)
#@inline function update_penalty!(
#   ::Union{RusanovNumericalFlux, CentralNumericalFluxFirstOrder},
function update_penalty!(
    ::RusanovNumericalFlux,
    ::OceanModel,
    n⁻,
    λ,
    ΔQ::Vars,
    Q⁻,
    A⁻,
    Q⁺,
    A⁺,
    t,
)
    ΔQ.η = -0

    return nothing
end

@inline function boundary_state!(nf, ocean::OceanModel, args...)
    boundary_conditions = ocean.problem.boundary_conditions
    return ocean_boundary_state!(nf, boundary_conditions, ocean, args...)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# For loading / saving julia variables
#ToDo: Update for JLD2 or BSON +++

using FileIO, JLD

"""
## Save a variable so that the name is known

Extension for saving Julia objects must be JLD
"""
function save(fPath, saveS)
    JLD.save(fPath,  "saveS", saveS)
   # jldopen(fPath, "w") do file
   #    write(file, "saveS", saveS)
   # end
   return nothing
end

"""
## Load a variable saved with save
# Types of variables loaded must be defined before load can be called
"""
function load(fPath)
    saveS = JLD.load(fPath, "saveS")
   # saveS = jldopen(fPath, "r") do file
   #    read(file, "saveS")
   # end
   return saveS
end

# ------------


text:
using Documenter

makedocs(
    format = :html,
    sitename = "StatsMakie",
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Tutorial" => "manual/tutorial.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/JuliaPlots/StatsMakie.jl.git",
    target = "build",
)


text:
__precompile__()

module ULID

export ulid

# Crockford's base-32 encoding
const _ENCODING = "0123456789ABCDEFGHJKMNPQRSTVWXYZ"
const _BASE = 32

@inline function prng()
    return rand(RandomDevice(), UInt16) / 0xFFFF
end

function encodetime(now, n)
    s = Vector{Char}(n)
    @inbounds for i = n:-1:1
        m = now % _BASE
        s[i] = _ENCODING[m+1]
        now = (now - m) ÷ _BASE
    end
    return String(s)
end

function encoderandom(n)
    s = Vector{Char}(n)
    @inbounds for i = n:-1:1
        r = floor(Int, _BASE * prng()) + 1
        s[i] = _ENCODING[r]
    end
    return String(s)
end

"""
    ulid()

Generate a Universally Unique Lexicographically Sortable Identifier
(ULID) as a string.
"""
function ulid()
    now = trunc(Int, Base.Dates.datetime2unix(Base.Dates.now()))
    return encodetime(now, 10) * encoderandom(16)
end

end # module


text:
@testset "Power Models ACP tests" begin

@testset "case14" begin
    println("==================================")
    println("ACPPowerModel case14.m")
    println("==================================")

    result = run_ots("data/pglib_opf_case14_ieee.m", ACPPowerModel, juniper_strong_no_restart)

    status = result["status"]
    @test status == :Optimal || status == :LocalOptimal

    juniper_val = result["objective"]

    println("")
    println("Solution by Juniper")
    println("obj: ", juniper_val)

    @test isapprox(juniper_val, 6291.28, atol=1e0)
end

end


text:
"""
    angle(v::AbstractVector{<:Real}) -> Float

Angle the vector `v` makes with the x axis, throws error if vector is not of size 2
"""
function angle(v::AbstractVector{<:Real})
    if length(v)!=2 error("Vector must be of length 2") end
    atan(v[2],v[1])
end

angle(v::SVector{2, <:Real})  = atan(v[2],v[1])

"""
    sqrt(z, branchcutangle::Real=π) -> typeof(z)

Return \\sqrt{x}, analytically continued from real line with the branch cut from 0 in a straight line along `branchcutangle`
"""
function sqrt(z, branchcutangle::Real)

    if branchcutangle > 0
        if angle(z) > branchcutangle return -sqrt(z) end
    else
        if angle(z) < branchcutangle return -sqrt(z) end
    end

    return sqrt(z)
end

"""
    round(num::Real, base::Real) -> Real

Nearest multiple of `base` to `num`, result will be promoted to num and base
"""
round(num::Real, base::Real) = round(num/base)*base

"""
    ceil(num::Real, base::Real) -> Real

Smallest multiple of `base` greater than `num`, result will be promoted to num and base
"""
ceil(num::Real, base::Real) = ceil(num/base)*base

"""
    ceil(num::Real, base::Real) -> Real

Largest multiple of `base` smaller than `num`, result will be promoted to num and base
"""
floor(num::Real, base::Real) = floor(num/base)*base


text:
#
export semq
#
# SEM global --> local operator
#
function semq(E,n,bc)

Q = spzeros(E*n,E*(n-1)+1);

Id= sparse(I,n,n);
i = 1;
j = 1;
for e=1:E
    Q[i:(i+n-1),j:(j+n-1)] = Id;
    i += n;
    j += n-1;
end

if(bc)
    Q[end,1] = 1;
    Q=Q[:,1:end-1];
end

return Q
end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    s1 = sum([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
    s2 = fsum([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
    a = [1, 2, 3, 4]
    a_sum = sum(a)
    @assert(s1 == 0.9999999999999999)
    @assert(s2 == 1.0)
    @assert(a_sum == 10)
    @assert(sin(0) == 0)
    @assert(cos(0) == 1)
    @assert(tan(0) == 0)
    @assert(sin(deg2rad(30)) == 0.49999999999999994)
    @assert(cos(deg2rad(30)) == 0.8660254037844387)
    @assert(tan(deg2rad(30)) == (√(3) / 3))
    @assert(round(12.556, digits = 2) == 12.56)
end


text:
# This file constains the Clock tools for time synchronization of DsSimulator.

import Base: iterate, length

"""
    $TYPEDEF

`Clock` generates time instants to sample the systems. A `Clock` is an iterable

# Fields
    $TYPEDFIELDS

# Example
```jldoctest
julia> clk = Clock(0 : 2)
Clock(gen:0:2, paused:false)

julia> for t in clk
       @show t
       end
t = 0
t = 1
t = 2
```
"""
mutable struct Clock{T, CB}
    "Internal generator. Can be any iterable"
    generator::T
    "If true, clock iteration is halted"
    paused::Bool
    "Callback set. See [`Callback`](@ref)"
    callbacks::CB
    "Name of clock"
    name::Symbol
    "Identiy number of clock"
    uuid::UUID
    function Clock(generator::T; callbacks::CB=nothing, name=Symbol()) where {T, CB}
        if hasmethod(iterate, Tuple{T})
            new{T, CB}(generator, false, callbacks, name, uuid4())
        else
            error("$gen is not iterable")
        end
    end
end

# Backward compatibility
Clock(ti::Real, dt::Real, tf::Real; kwargs...) = Clock(ti:dt:tf; kwargs...)

show(io::IO, clk::Clock) = print(io, "Clock(gen:$(clk.generator), paused:$(clk.paused))")

##### Iteration
function iterate(clk::Clock, state...)
    clk.paused && (@warn "Clock is paused"; return nothing)
    iter = iterate(clk.generator, state...)
    applycallbacks(clk)
    iter
end

"""
    $SIGNATURES

Pauses `clk`. When paused, the iteration of clock is halted.
# Example
```jldoctest
julia> clk = Clock(0 : 5)
Clock(gen:0:5, paused:false)

julia> iterate(clk)
(0, 0)

julia> iterate(clk, 1)
(2, 2)

julia> pause!(clk)
Clock(gen:0:5, paused:true)

julia> iterate(clk, 1)
┌ Warning: Clock is paused
└ @ Causal ~/.julia/dev/Causal/src/components/sources/clock.jl:41
```
"""
pause!(clk::Clock) = (clk.paused = true; clk)

length(clk::Clock) = length(clk.generator)


text:
import ..UncertainDatasets:
	AbstractUncertainValueDataset,
	UncertainValueDataset

import ..UncertainDatasets:ConstrainedUncertainValueDataset

"""
	constrain(udata::Vector{<:AbstractUncertainValue},
		s::SamplingConstraint) -> ConstrainedUncertainValueDataset

Return a vector of uncertain value by applying the constraint `s` to each
uncertain value in `udata`.
"""
function constrain(udata::Vector{<:AbstractUncertainValue}, constraint::SamplingConstraint)
	[constrain(uval, constraint) for uval in udata]
end

"""
	constrain(udata::Vector{<:AbstractUncertainValue},
		constraints::Vector{T}) where {T<:SamplingConstraint} -> Vector{<:AbstractUncertainValue}

Return a vector of uncertain values by applying a different sampling constraint to each uncertain
value in `udata`.
"""
function constrain(udata::Vector{<:AbstractUncertainValue},
		constraints::Vector{T}) where {T <: SamplingConstraint}

	if length(udata) != length(constraints)
		error("Number of sampling constraints must match length of dataset.")
	end

	return [constrain(udata[i], constraints[i]) for i in 1:length(udata)]
end

"""
	constrain(udata::AbstractUncertainValueDataset,
		s::SamplingConstraint) -> ConstrainedUncertainValueDataset

Return a uncertain dataset by applying the constraint `s` to each
uncertain value in `udata`.
"""
function constrain(udata::DT, constraint::SamplingConstraint) where {
		DT <: AbstractUncertainValueDataset}

	ConstrainedUncertainValueDataset([constrain(uval, constraint) for uval in udata])
end

"""
	constrain(udata::AbstractUncertainValueDataset,
		constraints::Vector{T}) where {T<:SamplingConstraint} -> ConstrainedUncertainValueDataset

Return a uncertain dataset by applying a different sampling constraint to each uncertain
value in `udata`.
"""
function constrain(udata::DT, constraints::Vector{T}) where {
		T<:SamplingConstraint,
		DT <: AbstractUncertainValueDataset}

	if length(udata) != length(constraints)
		error("Number of sampling constraints must match length of dataset.")
	end

	n_vals = length(udata)

	ConstrainedUncertainValueDataset([constrain(udata[i], constraints[i]) for i in 1:n_vals])
end

# """
# 	constrain(udata::ConstrainedUncertainValueDataset,
# 		constraint::SamplingConstraint) -> ConstrainedUncertainValueDataset

# Return a uncertain dataset by applying the `constraint` to each
# uncertain value in `udata`.
# """
# constrain(udata::ConstrainedUncertainValueDataset, constraint::SamplingConstraint) =
# 	ConstrainedUncertainValueDataset([constrain(uval, constraint) for uval in udata])

# """
# 	constrain(udata::ConstrainedUncertainValueDataset,
# 		constraints::Vector{T}) where {T<:SamplingConstraint} ->  ConstrainedUncertainValueDataset

# Return a uncertain dataset by applying a different sampling constraint to each uncertain
# value in `udata`.
# """
# function constrain(udata::ConstrainedUncertainValueDataset,
# 		constraints::Vector{T}) where {T<:SamplingConstraint}
# 	if length(udata) != length(constraints)
# 		error("Number of sampling constraints must match length of dataset.")
# 	end

# 	n_vals = length(udata)

# 	ConstrainedUncertainValueDataset([constrain(udata[i], constraints[i]) for i in 1:n_vals])
# end

export constrain, verify_constraints


text:
@testset "102.binary-tree-level-order-traversal.jl" begin
    @test level_order(TreeNode{Int}([3, 9, 20, nothing, nothing, 15, 7])) ==
          [[3], [9, 20], [15, 7]]
end


text:
function compute_frequencies(Z::Matrix, W::Vector, q)
    N, M = size(Z)
    length(W) == M || error("incompatible length of W")
    Pij = zeros(q-1, q-1, N, N)
    for i in 1:N
        for j in i:N
            for s in 1:M
                coli = Z[i,s]
                colj = Z[j,s]
                (coli == q || colj == q) && continue
                Pij[colj,coli,j,i] += W[s]
                Pij[coli,colj,i,j] = Pij[colj,coli,j,i]
            end
        end
    end
    Pi = zeros(q-1,N)
    for i in 1:N
        for s in 1:M
            coli = Z[i,s]
            (coli == q) && continue
            Pi[coli,i] += W[s]
        end
    end
    return Pij,Pi
end

compute_frequencies(Z,W) = compute_frequencies(Z,W,maximum(Z))
compute_frequencies(Z) = compute_frequencies(Z,ones(size(Z,2))/size(Z,2),maximum(Z))

function add_pseudocount(Pi_true::Array{Float64,2}, Pij_true::Array{Float64,4},pc)
    0 <= pc <= 1 || error("pseudocount $pc should be ∈ [0,1]")
    s,N = size(Pi_true)
    q = s+1
    Pipc  = (1.0-pc)*Pi_true  .+ pc/q
    Pijpc = (1.0-pc)*Pij_true .+ pc/q^2

    for i in 1:N
        Pijpc[:,:,i,i] .= (1-pc) * Pij_true[:,:,i,i]
        for a in 1:s
            Pijpc[a,a,i,i] = Pipc[a,i]
        end
    end
    Pijpc,Pipc
end

function Cij(Pi,Pij)
    q,N = size(Pi)
    cij = Matrix{Float64}(undef, N*q,N*q)
    mi = Vector{Float64}(undef,N*q)
    for i in 1:N
        rowblk = ((i-1)*q+1):i*q
        vpi = @view Pi[:,i]
        mi[rowblk].= vpi
        for j in 1:N
            vpj = @view Pi[:,j]
            vpij = @view Pij[:,:,i,j]
            colblk = ((j-1)*q+1):j*q
            cij[rowblk,colblk] .= vpij - vpi*vpj'
        end
    end
    return cij,mi
end

"""
    covariance_matrix(Z, W; pc::Real=0)
Compute the covariance matrix from numerical alignment `Z` and weights `W`. The output is a `N(q-1) × N(q-1)`
(skipping color `q`) symmetric matrix.
## Keywords arguments:
*  `pc` in [0,1]: pseudocount [default = `0`]
end
"""
function covariance_matrix(Z, W; pc::Real=0)
    0<=pc<=1 || error("pc = $pc should be ∈ [0,1]")
    Pij,Pi = compute_frequencies(Z,W)
    if pc > 0
        Pijpc,Pipc=add_pseudocount(Pi,Pij,pc)
        return Cij(Pipc,Pijpc)

    end
    return Cij(Pi,Pij)
end

covariance_matrix(Z; kwds...) = covariance_matrix(Z,ones(size(Z,2))/size(Z,2); kwds...)

max_stars_repo_name: pagnani/CorrDCA

max_issues_repo_name: pagnani/CorrDCA

max_forks_repo_name: pagnani/CorrDCA


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

using Markdown
using InteractiveUtils

# ╔═╡ acd7da78-1d69-11eb-2517-e1bafdc14f4a
using Pkg, DrWatson

# ╔═╡ acd82a8c-1d69-11eb-10c9-0902b833d79e
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots, LaTeXStrings, Measures
end

# ╔═╡ 2e6c02cc-1d69-11eb-283a-1bcf9e0044bf
md"## Listing 1.10a"

# ╔═╡ acd8c1cc-1d69-11eb-3164-5b9e66fc3664
begin
	f(x,y) = x^2 + y^2
	f0(x) = f(x,0)
	f2(x) = f(x,2)
end

# ╔═╡ ace462fa-1d69-11eb-2b31-c1c9efe2376f
begin
	x = -5:0.1:5
	y = -5:0.1:5
	z = [ f(x,y) for y in y, x in x ]
	p2 = surface(x, y, z, c=cgrad([:blue, :red]), legend=:none, ylabel="y", zlabel="f(x,y)")
end

# ╔═╡ ad012874-1d69-11eb-24a0-e936c7c266d5
md"## End of listing 1.10a"

# ╔═╡ Cell order:
# ╟─2e6c02cc-1d69-11eb-283a-1bcf9e0044bf
# ╠═acd7da78-1d69-11eb-2517-e1bafdc14f4a
# ╠═acd82a8c-1d69-11eb-10c9-0902b833d79e
# ╠═acd8c1cc-1d69-11eb-3164-5b9e66fc3664
# ╠═ace462fa-1d69-11eb-2b31-c1c9efe2376f
# ╟─ad012874-1d69-11eb-24a0-e936c7c266d5


text:
mutable struct Script
    instructions::Vector{Any}
    Script(instructions::Nothing) = new(Union{UInt8, Vector{UInt8}}[])
    Script(instructions) = new(instructions)
end

function show(io::IO, z::Script)
    for instruction in z.instructions
        if typeof(instruction) <: Integer
            if haskey(OP_CODE_NAMES, instruction)
                print(io, "\n", OP_CODE_NAMES[Int(instruction)])
            else
                print(io, "\n", string("OP_CODE_", Int(instruction)))
            end
        elseif typeof(instruction) <: Vector{UInt8}
            print(io, "\n", bytes2hex(instruction))
        else
            print(io, "\n", instruction)
        end
    end
end

"""
    scriptparse(::GenericIOBuffer) -> Script

Returns a Script object from an IOBuffer
"""
function scriptparse(s::IOBuffer)
    length_ = read_varint(s)
    instructions = []
    count = 0
    while count < length_
        current = UInt8[]
        readbytes!(s, current, 1)
        count += 1
        current_byte = current[1]
        if current_byte >= 1 && current_byte <= 75
            n = current_byte
            push!(instructions, read(s, n))
            count += n
        elseif current_byte == 76
            # op_pushdata1
            n = read(s, 1)[1]
            push!(instructions, read(s, n))
            count += n + 1
        elseif current_byte == 77
            # op_pushdata2
            n = reinterpret(Int16, read(s, 2))[1]
            push!(instructions, read(s, n))
            count += n + 2
        else
            # op_code
            push!(instructions, current_byte)
        end
    end
    if count != length_
        error("Error: parsing Script failed")
    end
    return Script(instructions)
end

function rawserialize(s::Script)
    result = UInt8[]
    for instruction in s.instructions
        if typeof(instruction) == UInt8
            push!(result, instruction)
        else
            length_ = length(instruction)
            if length_ < 0x4b
                append!(result, UInt8(length_))
            elseif length_ > 0x4b && length_ < 0x100
                append!(result, 0x4c)
                append!(result, UInt8(length_))
            elseif length_ >= 0x100 && length_ <= 0x208
                append!(result, 0x4d)
                result += int2bytes(length_, 2)
            else
                error("too long an instruction")
            end
            append!(result, instruction)
        end
    end
    return result
end

function serialize(s::Script)
    result = rawserialize(s)
    total = length(result)
    prepend!(result, encode_varint(total))
    return result
end

"""
    evaluate(s::Script, z::Integer) -> Bool

Evaluate if Script is valid given the transaction signature hash
"""
function evaluate(s::Script, z::Integer, witness::Union{Script, Nothing}=nothing)
    instructions = copy(s.instructions)
    stack = Vector{UInt8}[]
    altstack = Vector{UInt8}[]
    while length(instructions) > 0
        instruction = popfirst!(instructions)
        if typeof(instruction) <: Integer
            operation = OP_CODE_FUNCTIONS[instruction]
            function badop(instruction::Integer)
                println("bad op: ", OP_CODE_NAMES[instruction])
            end
            if instruction in (99, 100)
                # op_if/op_notif require the  array
                if !operation(stack, instructions)
                    badop(instruction)
                    return false
                end
            elseif instruction in (107, 108)
                # op_toaltstack/op_fromaltstack require the altstack
                if !operation(stack, altstack)
                    badop(instruction)
                    return false
                end
            elseif instruction in (172, 173, 174, 175)
                if !operation(stack, z)
                    badop(instruction)
                    return false
                end
            elseif !operation(stack)
                badop(instruction)
                return false
            end
        else
            push!(stack, instruction)
            # p2sh rule. if the next three instructions are:
            # OP_HASH160 <20 byte hash> OP_EQUAL this is the RedeemScript
            # OP_HASH160 == 0xa9 && OP_EQUAL == 0x87
            if length(instructions) == 3 && instructions[1] == 0xa9 &&
               typeof(instructions[2]) == Vector{UInt8} && length(instructions[2]) == 20 &&
               instructions[3] == 0x87
                println(" ---- ==== !!!! P2SH Script Found !!!! ==== ---- ")
                redeem_script = encode_varint(length(instruction))
                append!(redeem_script, instruction)
                # we execute the next three op codes
                pop!(instructions)
                h160 = pop!(instructions)
                pop!(instructions)
                if !op_hash160(stack)
                    return false
                end
                push!(stack, h160)
                if !op_equal(stack)
                    return false
                end
                # final result should be a 1
                if !op_verify(stack)
                    println("bad p2sh h160")
                    return false
                end
                # hashes match! now add the RedeemScript
                stream = IOBuffer(redeem_script)
                append!(instructions, scriptparse(stream).instructions)
            end

            if witness != nothing
                # witness program version 0 rule. if stack instructions are:
                # 0 <20 byte hash> this is p2wpkh
                if length(stack) == 2 && stack[1] == [0x00] && length(stack[2]) == 20
                    println(" ---- ==== !!!! P2WPKH Script Found !!!! ==== ---- ")
                    h160 = pop!(stack)
                    pop!(stack)
                    append!(instructions, witness.instructions)
                    append!(instructions, p2pkh_script(h160).instructions)
                end
                # witness program version 0 rule. if stack instructions are:
                # 0 <32 byte hash> this is p2wsh
                if length(stack) == 2 && stack[1] == [0x00] && length(stack[2]) == 32
                    println(" ---- ==== !!!! P2WSH Script Found !!!! ==== ---- ")
                    h256 = pop!(stack)
                    pop!(stack)
                    append!(instructions, witness.instructions[1:end-1])
                    witness_script = witness.instructions[end]
                    if h256 != sha256(witness_script)
                        print("bad sha256")
                        return false
                    end
                    # hashes match! now add the Witness Script
                    stream = IOBuffer(append!(encode_varint(length(witness_script)), witness_script))
                    witness_script_instructions = scriptparse(stream).instructions
                    append!(instructions, witness_script_instructions)
                end
            end
        end
    end
    if length(stack) == 0
        return false
    end
    if pop!(stack) == Vector{UInt8}[]
        return false
    end
    return true
end

"""
Takes a hash160 && returns the p2pkh scriptPubKey
"""
function p2pkh_script(h160::Vector{UInt8})
    script = Union{UInt8, Vector{UInt8}}[]
    pushfirst!(script, 0x76, 0xa9)
    push!(script, h160, 0x88, 0xac)
    return Script(script)
end

"""
Takes a hash160 && returns the p2sh scriptPubKey
"""
function p2sh_script(h160::Vector{UInt8})
    script = Union{UInt8, Vector{UInt8}}[]
    pushfirst!(script, 0xa9)
    push!(script, h160, 0x87)
    return Script(script)
end

"""
Takes a hash160 && returns the p2wpkh ScriptPubKey
"""
function p2wpkh_script(h160::Vector{UInt8})
    return Script([0x00, h160])
end

"""
Takes a hash160 && returns the p2wsh ScriptPubKey
"""
function p2wsh_script(hash256::Vector{UInt8})
    return Script([0x00, h256])
end

function scripttype(script::Script)
    if is_p2pkh(script)
        return "P2PKH"
    elseif is_p2sh(script)
        return "P2SH"
    elseif is_p2wsh(script)
        return "P2WSH"
    elseif is_p2wpkh(script)
        return "P2WPKH"
    else
        return error("Unknown Script type")
    end
end

"""
Returns whether this follows the
OP_DUP OP_HASH160 <20 byte hash> OP_EQUALVERIFY OP_CHECKSIG pattern.
"""
function is_p2pkh(script::Script)
    return length(script.instructions) == 5 &&
        script.instructions[1] == 0x76 &&
        script.instructions[2] == 0xa9 &&
        typeof(script.instructions[3]) == Vector{UInt8} &&
        length(script.instructions[3]) == 20 &&
        script.instructions[4] == 0x88 &&
        script.instructions[5] == 0xac
end

"""
Returns whether this follows the
OP_HASH160 <20 byte hash> OP_EQUAL pattern.
"""
function is_p2sh(script::Script)
    return length(script.instructions) == 3 &&
           script.instructions[1] == 0xa9 &&
           typeof(script.instructions[2]) == Vector{UInt8} &&
           length(script.instructions[2]) == 20 &&
           script.instructions[3] == 0x87
end

function is_p2wpkh(script::Script)
    length(script.instructions) == 2 &&
    script.instructions[1] == 0x00 &&
    typeof(script.instructions[2]) == Vector{UInt8} &&
    length(script.instructions[2]) == 20
end

"""
Returns whether this follows the
OP_0 <20 byte hash> pattern.
"""
function is_p2wsh(script::Script)
    length(script.instructions) == 2 &&
    script.instructions[1] == 0x00 &&
    typeof(script.instructions[2]) == Vector{UInt8} &&
    length(script.instructions[2]) == 32

end

const H160_INDEX = Dict([
    ("P2PKH", 3),
    ("P2SH", 2)
])

"""
Returns the address corresponding to the script
"""
function script2address(script::Script, testnet::Bool)
    type = scripttype(script)
    h160 = script.instructions[H160_INDEX[type]]
    return h160_2_address(h160, testnet, type)
end

@deprecate scriptevaluate(s::Script, z::Integer) evaluate(s::Script, z::Integer)
@deprecate scriptserialize(s::Script) serialize(s::Script)


text:
function psse_bus_results_compare(file_name, pf_results)
    pf_result_bus = CSV.read(file_name, DataFrame)

    v_diff = Float64[]
    angle_diff = Float64[]
    number = Int[]
    for (ix, n) in enumerate(eachrow(pf_result_bus))
        push!(v_diff, pf_results["bus_results"][ix, "Vm"] - n."Voltage (pu)")
        push!(
            angle_diff,
            pf_results["bus_results"][ix, "θ"] - (n."Angle (deg)" * (π / 180)),
        )
        push!(number, pf_results["bus_results"][ix, "bus_number"])
    end
    return v_diff, angle_diff, number
end

function psse_gen_results_compare(file_name, system::PSY.System)
    base_power = get_base_power(system)
    pf_result_gen = CSV.read(file_name, DataFrame)
    p_diff = Float64[]
    q_diff = Float64[]
    names = String[]
    for row in eachrow(pf_result_gen)
        if ismissing(row."Bus  Number")
            continue
        end
        name = "generator-" * "$(row."Bus  Number")" * "-" * row."Id"
        gen = get_component(ThermalStandard, system, name)
        ap = get_active_power(gen)
        rp = get_reactive_power(gen)
        push!(p_diff, ap - row."PGen (MW)" / base_power)
        push!(q_diff, rp - row."QGen (Mvar)" / base_power)
        push!(names, name)
    end
    return p_diff, q_diff, names
end


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
"""
    module HANTS

A Julia module of the Harmonic ANalysis of Time Series (HANTS)

Wout Verhoef
NLR, Remote Sensing Dept.
June 1998

Convert to MATLAB:
Mohammad Abouali (2011) (BSD 2-Clause)

https://mabouali.wordpress.com/projects/harmonic-analysis-of-time-series-hants/

Convert to Julia:
Shen Ruoque (2019) (MIT License)

https://github.com/shenrq/HANTS.jl

## Modified

Use long name instead of acronym.

Add support to multidimensional arrays.
"""
module HANTS
using LinearAlgebra

export hants, reconstruct

isinvalid(x) = ismissing(x) || isnothing(x)

"""
    hants(y, fet, dod, δ; nbase, nfreq, validrange, tseries, outlier)

Apply the HANTS process on the series `y`.

## Arguments:

- `y`   : array of input sample values (e.g. NDVI values)
- `fet` : fit error tolerance (points deviating more than fet from curve fit are rejected)
- `dod` : degree of overdeterminedness (iteration stops if number of points reaches the
minimum required for curve fitting, plus dod). This is a safety measure
- `δ`   : (\\delta<tab>) small positive number (e.g. 0.1) to suppress high amplitudes
- `nbase` : length of the base period, measured in virtual samples
(days, decades, months, etc.)
- `nfreq` : number of frequencies to be considered above the zero frequency.
Set to 3 by default
- `validrange` : tuple of valid range minimum and maximum
(values outside the valid range are rejeced right away)
Set to the extremas by default.
- `tseries` : array same size as `y` of time sample indicators (indicates virtual sample
number relative to the base period); numbers in array `tseries` maybe greater
than nbase. Set to `1:length(tseries)` by default
- `outlier` : a symbol `:Hi` or `:Lo` indicating rejection of high or low outliers.
Set to `nothing` by default

## Outputs:

- `A` : returned array of amplitudes, first element is the average of the curve
- `φ` : (\\varphi<tab>) returned array of phases, first element is zero
- `yrec` : array holding reconstructed time series
"""
function hants(
    y::AbstractVector{T}, fet, dod, δ;
    nbase=length(y), nfreq=3, validrange=extrema(y[.!ianan.(y)]),
    tseries=1:length(y), outlier=nothing
) where {T<:AbstractFloat}

    ny = length(y)
    nr = min(2nfreq+1, ny)
    matirx = Matrix{T}(undef, 2nfreq+1, ny)
    A = Vector{T}(undef, nfreq+1)
    φ = Vector{T}(undef, nfreq+1)
    local zr, yrec

    soutlier = outlier in [:Hi, :High] ? -1 : outlier in [:Lo, :Low] ? 1 : 0

    low, high = Float64.(validrange)

    noutmax = ny - nr - dod
    matirx[1, :] .= 1

    ang = 2 * (0:nbase-1) / nbase
    cs = cospi.(ang); sn = sinpi.(ang)
    @inbounds @simd for i = 1:nfreq
        j = 1:ny
        index = @. 1 + mod(i * (tseries[j] - 1), nbase)
        matirx[2i  , j] = cs[index]
        matirx[2i+1, j] = sn[index]
    end

    p = low .< y .< high
    nout = ny - sum(p)

    if nout > noutmax
        error("nout > noutmax, get nout $nout and noutmax $noutmax")
        return
    end

    ready = false; nloop = 0; nloopmax = ny

    while (!ready) && (nloop < nloopmax)
        nloop += 1
        za = matirx * (p .* y)

        arr = matirx * diagm(0=>p) * matirx'
        arr += δ * I
        arr[1, 1] -= δ
        zr = arr \ za

        yrec = matirx' * zr
        diffvec = soutlier * (yrec - y)
        err = p .* diffvec

        rankvec = sortperm(err)

        maxerr = diffvec[rankvec[ny]]
        ready = maxerr ≤ fet || nout == noutmax

        if !ready
            i = ny; j = rankvec[i]
            while p[j] * diffvec[j] > maxerr / 2 && nout < noutmax
                p[j] = 0; nout += 1; i -= 1; j = rank(i)
            end
        end
    end

    A[1] = zr[1]
    φ[1] = 0

    push!(zr, 0)

    i = 2:2:nr
    ifr = 2:1:nr÷2+1
    ra = zr[i]; rb = zr[i.+1]
    A[ifr] = @. √(ra^2 + rb^2)
    phase = atand.(rb, ra)
    replace!(x -> x < 0 ? x + 360 : x, phase)
    φ[ifr] = phase

    A, φ, yrec
end

hants(
    y::AbstractVector{Union{Missing, Nothing, T}}, args...; kargs...
) where {T<:AbstractFloat} = hants(
    convert(Vector{T}, replace(x -> isinvalid(x) ? NaN : x, y)), args...; kargs...
)

hants(
    y::AbstractVector{Union{Missing, T}}, args...; kargs...
) where {T<:AbstractFloat} = hants(
    convert(Vector{T}, replace(x -> isinvalid(x) ? NaN : x, y)), args...; kargs...
)

hants(
    y::AbstractVector{Union{Nothing, T}}, args...; kargs...
) where {T<:AbstractFloat} = hants(
    convert(Vector{T}, replace(x -> isinvalid(x) ? NaN : x, y)), args...; kargs...
)

hants(
    y::AbstractVector{<:Union{Missing, Nothing, Integer}}, args...; kargs...
) = hants(
    convert(Vector{Union{Missing, Nothing, Float64}}, y), args...; kargs...
)

"""
    hants(arr, fet, dod, δ; dims, nbase, nfreq, validrange, tseries, outlier)

Apply the HANTS process on the multidimensional array `arr`
along the given dimension `dims`.
"""
function hants(
    arr::AbstractArray{T,N}, fet, dod, δ;
    dims::Integer, nbase=size(arr)[dims], nfreq=3,
    validrange=extrema(arr[.!isinvalid.(arr)]),
    tseries=1:size(arr)[dims], outlier=nothing
) where{T<:AbstractFloat, N}
    if dims > N error("dims must less than N, get dims $dims and N $N.") end

    arrsize = size(arr)
    Asize = copy(collect(arrsize))
    Asize[dims] = nfreq + 1

    ny = arrsize[dims]
    pdims = (dims, setdiff(1:ndims(arr), dims)...)
    arr_rec = Array{T, N}(undef, arrsize...)
    A = Array{T, N}(undef, Asize...)
    φ = Array{T, N}(undef, Asize...)
    isize = prod(arrsize)÷ny
    arrv = reshape(permutedims(arr, pdims), ny, isize)
    arr_recv = reshape(PermutedDimsArray(arr_rec, pdims), ny, isize)
    Av = reshape(PermutedDimsArray(A, pdims), nfreq+1, isize)
    φv = reshape(PermutedDimsArray(φ, pdims), nfreq+1, isize)

    @inbounds @simd for i = 1:isize
        Av[:, i], φv[:, i], arr_recv[:, i] = hants(
            arrv[:, i], fet, dod, δ; nbase=nbase, nfreq=nfreq,
            validrange=validrange, tseries=tseries, outlier=outlier
        )
    end

    A, φ, arr_rec
end

hants(
    arr::AbstractArray{Union{Missing, Nothing, T},N}, args...; kargs...
) where {T<:AbstractFloat, N} = hants(
    convert(Array{T,N}, replace(x -> isinvalid(x) ? NaN : x, arr)), args...; kargs...
)

hants(
    arr::AbstractArray{Union{Missing, T},N}, args...; kargs...
) where {T<:AbstractFloat, N} = hants(
    convert(Array{T,N}, replace(x -> isinvalid(x) ? NaN : x, arr)), args...; kargs...
)

hants(
    arr::AbstractArray{Union{Nothing, T},N}, args...; kargs...
) where {T<:AbstractFloat, N} = hants(
    convert(Array{T,N}, replace(x -> isinvalid(x) ? NaN : x, arr)), args...; kargs...
)

hants(
    arr::AbstractArray{Union{Missing, Nothing, Integer},N}, args...; kargs...
) where {N} = hants(
    convert(Array{Union{Missing, Nothing, Float64},N}, arr), args...; kargs...
)

"""
    reconstruct(A, φ, nbase)

Comput reconstructed series.

## Arguments:

- `A` : array of amplitudes, first element is the average of the curve
- `φ` : (\\varphi<tab>) array of phases, first element is zero
- `nbase` : length of the base period, measured in virtual samples
(days, decades, months, etc.)

## Outputs:

- `y` : reconstructed array of sample values
"""
function reconstruct(
    A::AbstractVector{T}, φ, nbase::Integer
) where {T<:AbstractFloat}
    nfreq = length(A)
    y = zeros(T, nbase)
    a_coef = @. A * cosd(φ)
    b_coef = @. A * sind(φ)
    @inbounds for i = 1:nfreq
        tt = @. (i - 1) * 2 * (0:nbase-1) / nbase
        @. y += a_coef[i] * cospi(tt) + b_coef[i] * sinpi(tt)
    end
    y
end

"""
    reconstruct(A, φ, nbase; dims::Integer)

Comput reconstructed multidimensional array.
"""
function reconstruct(
    A::AbstractArray{T,N}, φ, nbase::Integer; dims::Integer
) where{T<:AbstractFloat, N}
    if dims > N error("dims must less than N") end

    Asize = collect(size(A))
    arrsize = copy(Asize)
    arrsize[dims] = nbase

    nA = Asize[dims]
    pdims = (dims, setdiff(1:ndims(A), dims)...)
    arr = Array{T}(undef, arrsize...)
    isize = prod(Asize)÷nA
    arrv = reshape(PermutedDimsArray(arr, pdims), nbase, isize)
    Av = reshape(permutedims(A, pdims), nA, isize)
    φv = reshape(permutedims(φ, pdims), nA, isize)

    @inbounds @simd for i = 1:isize
        arrv[:, i] = reconstruct(Av[:, i], φv[:, i], nbase)
    end

    arr
end

end # module HANTS

max_stars_repo_name: shenrq/JuliaHANTS


text:
function iterative_first_improvement(tuning_run::Run,
                                     channel::RemoteChannel;
                                     cutoff::Integer = 10_000)
    function iterate(tuning_run::Run)
        return first_improvement(tuning_run, cutoff = cutoff)
    end

    technique(tuning_run,
              channel,
              iterate,
              name = "Iterative First Improvement")
end


text:
using NeuralVerification
using NeuralPriorityOptimizer
using NeuralPriorityOptimizer: cell_to_all_subcells
using LinearAlgebra
using LazySets

include(string(@__DIR__, "/utils.jl"))

function generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)
    # Read in the network
    network = read_nnet(network_file)

    # Define the coefficients for a linear objective
    coefficients = [-0.74; -0.44]
    full_input_region = Hyperrectangle(low=lbs, high=ubs)
    cells = cell_to_all_subcells(full_input_region, cells_per_dim)

    # Parameters for the solvers
    params = PriorityOptimizerParameters(stop_gap=1e-4, stop_frequency=10, max_steps=20000000) # for the priority solver

    # Now, time the solvers for each cell and record the optimal values
    times_max = zeros(cells_per_dim...)
    lower_bounds_max = zeros(cells_per_dim...)
    upper_bounds_max = zeros(cells_per_dim...)

    times_min = zeros(cells_per_dim...)
    lower_bounds_min = zeros(cells_per_dim...)
    upper_bounds_min = zeros(cells_per_dim...)

    i_linear = 1 # number the queries starting at 1
    for i in CartesianIndices(cells)
        cell = cells[i]
        # Run the priority optimizer then the mip splitting optimizer
        time_max = @elapsed x_star, lower_bound_max, upper_bound_max, steps_max = optimize_linear(network, cell, coefficients, params; maximize=true)
        times_max[i], lower_bounds_max[i], upper_bounds_max[i] = time_max, lower_bound_max, upper_bound_max

        time_min = @elapsed x_star, lower_bound_min, upper_bound_min, steps_min = optimize_linear(network, cell, coefficients, params; maximize=false)
        times_min[i], lower_bounds_min[i], upper_bounds_min[i] = time_min, lower_bound_min, upper_bound_min

	sample_time = @elapsed sample_lower, sample_upper = sample_based_bounds(network, cell, coefficients, 5000)

        println("Cell ", i)
        println("min and max solve steps: ", [steps_min, steps_max])
        println("lower, upper max: ", [lower_bound_max, upper_bound_max])
        println("lower, upper min: ", [lower_bound_min, upper_bound_min])
	println("so our control interval is: ", [lower_bound_min, upper_bound_max])
	println("compared to from sampling: ", [sample_lower, sample_upper])
        println("Time min: ", time_min, "  time max: ", time_max, " time sample: ", sample_time)

        epsilon = 0.2 # amount off the true maximum or minimum to point the assertion
        # Write out the query, it should cycle between:
        # (i) property has -0.74y0 - 0.44y1 <= minimum - epsilon which will be UNSAT (HOLDS that it's > minimum-epsilon)
        # (ii) property has -0.74y0 - 0.44y1 <= minimum + epsilon which will be SAT (VIOLATED that it's > minimum + epsilon)
        # (iii) property has -0.74y0 - 0.44y1 >= maximum + epsilon which will be UNSAT (HOLDS that it's < maximum + epsilon)
        # (iv) property has -0.74y0 - 0.44y1 >= maximum - epsilon which will be SAT (VIOLATED that it's < maximum - epsilon)
        output_file = string(query_base_name, "_", string(i_linear), ".vnnlib")
	# it shouldn't matter upper vs. lower bound here, since they should all be within 1e-4 of each other
	@assert abs(upper_bound_max - lower_bound_max) <= 1e-4 "gap should be <= 1e-4"
	@assert abs(upper_bound_min - lower_bound_min) <= 1e-4 "gap should be <= 1e-4"
        if i_linear % 4 == 0
            write_query(network, cell, "<=", lower_bound_min - epsilon, coefficients, output_file)
        elseif i_linear % 4 == 1
            write_query(network, cell, "<=", lower_bound_min + epsilon, coefficients, output_file)
        elseif i_linear % 4 == 2
            write_query(network, cell, ">=", upper_bound_max + epsilon, coefficients, output_file)
        elseif i_linear % 4 == 3
            write_query(network, cell, ">=", upper_bound_max - epsilon, coefficients, output_file)
        end
        i_linear = i_linear + 1
    end

    return times_max, lower_bounds_max, upper_bounds_max, times_min, lower_bounds_min, upper_bounds_min
end

function write_query(network, cell, symbol, value, coefficients, output_file)
    open(output_file, "w") do io
        print_var_definition(io, network)
        print_cell(io, cell)
        print_output_constraint(io, coefficients, symbol, value)
    end
end

function print_var_definition(io, network)
    # Define the input variables
    for i = 1:size(network.layers[1].weights, 2)
        println(io, string("(declare-const X_", i-1, " Real)"))
    end

    # Define the output variables
    for i = 1:size(network.layers[end].weights, 1)
        println(io, string("(declare-const Y_", i-1, " Real)"))
    end
end

function print_cell(io, cell)
    lbs, ubs = low(cell), high(cell)
    for (i, (lb, ub)) in enumerate(zip(lbs, ubs))
        println(io, string("(assert (>= X_", i-1, " ", lb, "))"))
        println(io, string("(assert (<= X_", i-1, " ", ub, "))"))
    end
end

# For now for simplicity we'll assume two-dimensional coefficients
function print_output_constraint(io, coefficients, symbol, value)
    println(io, string("(assert (", symbol, " (+ (* ", coefficients[1], " Y_0) (* ", coefficients[2], " Y_1)) ", value, "))"))
end

# Small is 0.03, medium is 0.06, large is 0.12
# we will do a Small, Medium, and Large set of queries
# at the negative edge of the state space and at the center of the state space
# I expect the center of the state space to be slightly more challenging

# Use the full best network for all queries
network_file = string(@__DIR__, "/networks/GANControl/full_mlp_best_conv.nnet")

# Small queries at the negative edge of the state space
# formed by splitting the state space into 5 cells along each axis
# leading to 25 queries with width [full latent, full latent, 0.03, 0.03]
query_base_name = "queries/negative_edge_small"
lbs = [-0.8, -0.8, -1.0, -1.0]
ubs = [0.8, 0.8, -0.85, -0.85]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Small Edge Region -----------\n")
#generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Medium queries at the negative edge of the state space
# formed by splitting the state space into 5 cells along each axis
# leading to 25 queries with width [full latent, full latent, 0.06, 0.06]
query_base_name = "queries/negative_edge_medium"
lbs = [-0.8, -0.8, -1.0, -1.0]
ubs = [0.8, 0.8, -0.7, -0.7]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Medium Edge Region -----------\n")
#generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Medium queries at the negative edge of the state space
# formed by splitting the state space into 5 cells along each axis
# leading to 25 queries with width [full latent, full latent, 0.12, 0.12]
query_base_name = "queries/negative_edge_large"
lbs = [-0.8, -0.8, -1.0, -1.0]
ubs = [0.8, 0.8, -0.4, -0.4]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Large Edge Region -----------\n")
#generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Small queries at the center of the state space
# formed by splitting the state space into 3 cells along each axis
# leading to 100 queries with width [full latent, full latent, 0.03, 0.03]
query_base_name = "queries/center_small"
lbs = [-0.8, -0.8, -0.075, -0.075]
ubs = [0.8, 0.8, 0.075, 0.075]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Small Center Region -----------\n")
#generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Medium queries at the center of the state space
# formed by splitting the state space into 3 cells along each axis
# leading to 100 queries with width [full latent, full latent, 0.03, 0.03]
query_base_name = "queries/center_medium"
lbs = [-0.8, -0.8, -0.15, -0.15]
ubs = [0.8, 0.8, 0.15, 0.15]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Medium Center Region -----------\n")
#generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Test small upper right corner of state space
query_base_name = "queries/positive_edge_small"
lbs = [-0.8, -0.8, 0.85, 0.85]
ubs = [0.8, 0.8, 1.0, 1.0]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Small Positive Edge Region -----------\n")
generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)

# Large queries at the center of the state space
# formed by splitting the state space into 3 cells along each axis
# leading to 100 queries with width [full latent, full latent, 0.03, 0.03]
query_base_name = "queries/center_large"
lbs = [-0.8, -0.8, -0.3, -0.3]
ubs = [0.8, 0.8, 0.3, 0.3]
cells_per_dim = [1, 1, 5, 5]
println("\n\n------------- Starting Large Center Region -----------\n")
generate_queries_for_region(lbs, ubs, cells_per_dim, network_file, query_base_name)


text:
using AIBECS, Test
using SparseArrays, LinearAlgebra
using Unitful
using Unitful: °, m, km
using WorldOceanAtlasTools
using DiffEqBase
using ForwardDiff, DualNumbers
using DataFrames
using Distributions
using Plots

# For CI, make sure the downloads do not hang
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

test_setup_only = [:TwoBoxModel, :Archer_etal_2000, :OCIM1, :OCIM0, :OCCA]
# Using `include` evaluates at global scope,
# so `Circulation` must be changed at the global scope too.
# This is why there is an `eval` in the for loop(s) below
@testset "test setup.jl only" for C in test_setup_only
    @testset "$C" begin
        eval(:(Circulation = $C))
        include("setup.jl") #TODO change back to just setup.jl
    end
end

test_plots = [:OCIM2]
@testset "Test setup, plots, and sources" for C in test_plots
    @testset "$C" begin
        eval(:(Circulation = $C))
        include("setup.jl")    #TODO change back to just setup.jl (can't remember why...)
        include("plots.jl")    #TODO change back to just setup.jl
        include("sources.jl")  #TODO change back to just setup.jl
    end
end

test_everything = [:Primeau_2x2x2]
@testset "test everything" for C in test_everything
    @testset "$C" begin
        eval(:(Circulation = $C))
        include("setup.jl")
        include("particles.jl")
        include("parameters.jl")
        include("bgc_functions.jl")
        include("gridtools.jl")
        include("cost_functions.jl")
        include("solvers.jl")
        include("time_steps.jl")
        include("derivatives.jl")
    end
end


text:
using PSFModels: Gaussian, Moffat, AiryDisk

# TODO - figure out why using `ones(1, 1)` fails only in test
@testset "image injection - $(typeof(kernel))" for kernel in [[0 0 0; 0 1 0; 0 0 0],
                                          Gaussian(1e-2),
                                          Moffat(1e-3),
                                          AiryDisk(1e-2)]
    frame = zeros(3, 3)
    @test inject(frame, kernel, (1, 3))[3, 1] ≈ 1
    @test inject(frame, kernel, (2, 2); A=2)[2, 2] ≈ 2
    @test inject(frame, kernel, (10, 10)) ≈ zeros(3, 3) atol = sqrt(eps(eltype(frame)))

    cube = zeros(10, 3, 3)
    @test inject(cube, kernel, (2, 2))[:, 2, 2] ≈ ones(10)
    @test inject(cube, 90ones(10), kernel, (2, 2))[:, 2, 2] ≈ ones(10)
end

@testset "cube generator - $(typeof(kernel))" for kernel in [[0 0 0; 0 1 0; 0 0 0],
                                          Gaussian(1e-2),
                                          Moffat(1e-3),
                                          AiryDisk(1e-2)]
    cube = zeros(10, 3, 3)
    angles = 90ones(10)
    gen = CubeGenerator(cube, angles, kernel)

    c = gen((2, 2))
    c2 = gen(Float32, (2, 2))
    c3 = gen(Polar(0, 0))
    @test c ≈ c2 ≈ c3
    @test eltype(c2) === Float32

    q = gen(ones(size(cube)), (2, 2))
    @test all(q .≈ c .+ 1)

    flat = gen(zero(flatten(cube)), (2, 2))
    @test flat ≈ flatten(c)
end

@testset "cube generator AnnulusView - $(typeof(kernel))" for kernel in [[0 0 0; 0 1 0; 0 0 0], Gaussian(1e-2)]
    cube = AnnulusView(zeros(10, 21, 21))
    angles = 90ones(10)
    gen = CubeGenerator(cube, angles, kernel)

    c = gen((11, 11))
    c2 = gen(Float32, (11, 11))
    c3 = gen(Polar(0, 0))
    @test c ≈ c2 ≈ c3
    @test eltype(c2) === Float32

    q = gen(ones(size(cube)), (11, 11))
    @test all(q .≈ c .+ 1)
end

@testset "inject av - $(typeof(kernel))" for kernel in [[0 0 0; 0 1 0; 0 0 0], Gaussian(1e-2)]
    cube = AnnulusView(zeros(10, 21, 21))
    angles = 90ones(10)

    c = inject(cube, angles, kernel, (11, 11))
    c2 = inject(cube, angles, kernel, Polar(0, 0))
    @test c ≈ c2
    @test c isa AnnulusView && c2 isa AnnulusView

    gen = CubeGenerator(cube, angles, kernel)
    flat = gen(cube(), (11, 11))
    @test flat ≈ c()
end

@testset "inject mav - $(typeof(kernel))" for kernel in [[0 0 0; 0 1 0; 0 0 0], Gaussian(1e-2)]
    cube = MultiAnnulusView(zeros(10, 21, 21), 2)
    angles = 90ones(10)

    c = inject(cube, angles, kernel, (11, 11))
    c2 = inject(cube, angles, kernel, Polar(0, 0))
    @test c ≈ c2
    @test c isa MultiAnnulusView && c2 isa MultiAnnulusView
end


text:
# Implements a non-linear fit

using Polynomials
using NeXLUncertainties: LabeledValues, Label
using NeXLSpectrum
using Gadfly

abstract type FitFunction end

struct CoeffN <: Label
    n::Int
end
struct FWHMLabel <: Label end
struct PeakI <: Label
    xray::CharXRay
end

struct PeakShapeFF <: FitFunction
    material::Material
    order::Int

    PeakShapeFF(mat::Material, cxrs::AbstractVector{CharXRay}, order::Int) = new(mat, order)
end

function fitfunc(ps::PeakShapeFF, args::LabeledValues)::Function
    det = BasicEDS(
        length(chs),
        channelcount::Int,
        LinearEnergyScale(args[CoeffN(0)], args[CoeffN(1)]),
        MnKaResolution(args[FWHMLabel()]),
        10
    )

    resp = detectorresponse(det, SDDEfficiency(AP33Model()))

    @show mnka
    imaxs = map(ps.xrays) do xray
        args[PeakI(xray)]*weight(xray)
    end
    @show imaxs
    exs = map(xray->energy(xray), ps.xrays)
    @show exs
    return ch -> begin
        e = efunc(ch)
        sum(i->imaxs[i]*profile(e, exs[i], mnka), eachindex(imaxs))
    end
end

function plotit()
    chs = 1:2048

    psff = PeakShapeFF(characteristic(n"Fe",alltransitions), 1)

    cxrs = characteristic(n"Fe", alltransitions)
    args = LabeledValues(
        map(cxr->PeakI(cxr)=> 1.0, cxrs)...,
        CoeffN(0) => 0.0,
        CoeffN(1) => 10.0,
        FWHMLabel() =>132.0
        )

    ff = fitfunc(psff, args)
    plot(x=chs, y = ff.(chs), Geom.line)
end