text:
using KernelAbstractions.Extras.LoopInfo: @unroll

#####
##### Kernels that ensure 'no-flux' from second- and fourth-order diffusion operators.
##### Note that flux divergence associated with a flux boundary condition is added
##### in a separate step.
#####

@kernel function _fill_west_halo!(c, ::FBC, H, N)
    j, k = @index(Global, NTuple)

    @unroll for i in 1:H
        @inbounds c[1-i, j, k] = c[i, j, k]
    end
end

@kernel function _fill_south_halo!(c, ::FBC, H, N)
    i, k = @index(Global, NTuple)

    @unroll for j in 1:H
        @inbounds c[i, 1-j, k] = c[i, j, k]
    end
end

@kernel function _fill_bottom_halo!(c, ::FBC, H, N)
    i, j = @index(Global, NTuple)

    @unroll for k in 1:H
        @inbounds c[i, j, 1-k] = c[i, j, k]
    end
end

@kernel function _fill_east_halo!(c, ::FBC, H, N)
    j, k = @index(Global, NTuple)

    @unroll for i in 1:H
        @inbounds c[N+i, j, k] = c[N+1-i, j, k]
    end
end

@kernel function _fill_north_halo!(c, ::FBC, H, N)
    i, k = @index(Global, NTuple)

    @unroll for j in 1:H
        @inbounds c[i, N+j, k] = c[i, N+1-j, k]
    end
end

@kernel function _fill_top_halo!(c, ::FBC, H, N)
    i, j = @index(Global, NTuple)

    @unroll for k in 1:H
        @inbounds c[i, j, N+k] = c[i, j, N+1-k]
    end
end

#####
##### Kernel launchers for flux boundary conditions
#####

  fill_west_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :yz, _fill_west_halo!,   c, bc, grid.Hx, grid.Nx, dependencies=dep)
  fill_east_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :yz, _fill_east_halo!,   c, bc, grid.Hx, grid.Nx, dependencies=dep)
 fill_south_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :xz, _fill_south_halo!,  c, bc, grid.Hy, grid.Ny, dependencies=dep)
 fill_north_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :xz, _fill_north_halo!,  c, bc, grid.Hy, grid.Ny, dependencies=dep)
fill_bottom_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :xy, _fill_bottom_halo!, c, bc, grid.Hz, grid.Nz, dependencies=dep)
   fill_top_halo!(c, bc::FBC, arch, dep, grid, args...) = launch!(arch, grid, :xy, _fill_top_halo!,    c, bc, grid.Hz, grid.Nz, dependencies=dep)


text:
function gatingGraph(sample::EzXML.Node; channelMap::Dict=Dict(), transform::Function=x->asinh(x/250))

 	############# store strategy in graph
	name = findfirst("../SampleNode",sample)["name"]
	graph = MetaDiGraph{Int64,Bool}()
	set_props!(graph,Dict(:sample=>name))

	############################################## population names
	populations = findall("..//Population",sample)
	if length(populations)==0

		@warn("no populations found")
		return graph
	end

	for population ∈ populations
		Gate,name = findfirst("Gate",population), population["name"]

		id = Gate["gating:id"]
		parent_id = haskey(Gate,"gating:parent_id") ? Gate["gating:parent_id"] : nothing

		#################### iterate through disconnected polygons
		for gate ∈ eachelement(Gate)

			channels = map( dimension -> dimension["data-type:name"] ∈ keys(channelMap) ? channelMap[dimension["data-type:name"]]
				: throw("""$(dimension["data-type:name"]) not found in channels $(keys(channelMap))"""),
				findall("gating:dimension/data-type:fcs-dimension",gate) )

			@assert(gate.name ∈ ["PolygonGate","RectangleGate"], "$(gate.name) not supported for population label $name in sample")
			@assert(length(channels) ≤ 2, "length($channels)>2 in $(gate.name) not supported for population label $name in sample")

			if gate.name == "PolygonGate"

				vertices = map( coordinate-> parse(Float32,coordinate["data-type:value"]),
					findall("gating:vertex/gating:coordinate",gate) )
				polygon = map( (x,y)->SVector(transform(x),transform(y)), @view(vertices[1:2:end]), @view(vertices[2:2:end]) )

			elseif gate.name == "RectangleGate"

				minima = map( dimension-> haskey(dimension,"gating:min") ? parse(Float32,dimension["gating:min"]) : -Inf,
					findall("gating:dimension",gate) )
				maxima = map( dimension-> haskey(dimension,"gating:max") ? parse(Float32,dimension["gating:max"]) : Inf,
					findall("gating:dimension",gate) )

				polygon = map( (x,y)->SVector(transform(x),transform(y)),
					[first(minima),first(minima),first(maxima),first(maxima)],
					length(channels) == 1 ? [ -Inf, Inf, Inf, -Inf] : [ last(minima), last(maxima), last(maxima), last(minima)] )
			end
			push!(polygon,first(polygon))

			MetaGraphs.add_vertex!(graph) ################## store gate as vertex
			set_indexing_prop!( graph, MetaGraphs.nv(graph), :id,id )
			set_props!( graph, MetaGraphs.nv(graph), Dict(:channels=>channels,:polygon=>polygon,:name=>name) )

			################# connect parent gates to children
			(parent_id !== nothing) && MetaGraphs.add_edge!(graph,graph[parent_id,:id],MetaGraphs.nv(graph))
		end
    end

    return graph
end

function gatingGraph(sample::Nothing; kwargs...)
   return MetaDiGraph{Int64,Bool}()
end

function gate(graph::MetaDiGraph,idx::Integer;prefix::String="__gate__")
	channels, gateName = get_prop(graph,idx,:channels), prefix*get_prop(graph,idx,:name)

	return [first(channels),last(channels)] =>
		ByRow( (x,y)->inpolygon( SVector(x,y), get_prop(graph,idx,:polygon);

			in=true, on=false, out=false)
	) => gateName
end

function gate!(data::DataFrame,gating::MetaDiGraph,idx::Integer;prefix::String="__gate__")
	name = prefix*get_prop(gating,idx,:name)
	transform!( data, gate(gating,idx;prefix=prefix) )

	for parent ∈ MetaGraphs.inneighbors(gating,idx) ######### todo @gszep does not support multiple parents
		parent = prefix*get_prop(gating,parent,:name)
		transform!(data, [name,parent] => ByRow((x,y)->x&y) => name)
	end

	for child ∈ MetaGraphs.outneighbors(gating,idx)
		gate!(data,gating,child;prefix=prefix)
	end
end

function gate(data::DataFrame,gating::MetaDiGraph;prefix::String="__gate__")

	roots = [ idx for idx ∈ 1:MetaGraphs.nv(gating) if MetaGraphs.indegree(gating,idx)==0 ]
	names = unique([ prefix*get_prop(gating,idx,:name) for idx ∈ 1:MetaGraphs.nv(gating)]) # todo @gszep why wrap in unique()? possible bug?

	if length(roots) == 0
		return DataFrame("Unlabelled"=>fill(true,size(data,1)))
	end

	for idx ∈ roots
		gate!(data,gating,idx;prefix=prefix)
	end

	labels = select(data,names)
	select!(data,Not(names))

	rename!(labels,Dict([ name=>replace(name,prefix=>"") for name ∈ names ]))
	return labels
end

ext: jl


text:
__precompile__()

module StokesDiffEq

using Reexport
@reexport using DiffEqBase
import DiffEqBase: solve
import Base: size
using Parameters, VectorizedRoutines.Matlab

include("problems.jl")
include("solutions.jl")
include("stokes_solve.jl")
include("premade_problems.jl")

export solve, StokesProblem, StokesSolution, FDMMesh

export prob_stokes_homogenous, prob_stokes_dirichletzero

end # module


text:
export id_to_index,
    add_node!,
    delete_node!,
    update_node_name!,
    update_node_position!,
    update_node_params!,
    add_edge!,
    delete_edge!

id_to_index(id::Int, ng::NodeGraph) = findnext(n -> n.id == id, ng.nodes, 1)

function add_node!(ng::NodeGraph, type::Symbol, newnode_id::Int, new_pos::NTuple{2,Float64})
    new_params = create(Val{type}())
    push!(ng.nodes, Node(newnode_id, "new_$type", new_pos..., new_params))
end

function delete_node!(ng::NodeGraph, node_ind::Int)
    removed = popat!(ng.nodes, node_ind)
    filter!(e -> e.from != removed.id && e.to != removed.id, ng.edges)
end

function update_node_name!(ng::NodeGraph, node_ind::Int, new_name::String)
    ng.nodes[node_ind].name = new_name
end

function update_node_position!(ng::NodeGraph, node_ind::Int, x::Float64, y::Float64)
    ng.nodes[node_ind].x = x
    ng.nodes[node_ind].y = y
end

function update_node_params!(ng::NodeGraph, node_ind::Int, key::Symbol, value)
    parsed_value = JSON3.read(
        value,
        Dict{Symbol,typeof(getproperty(ng.nodes[node_ind].params, key))}
    )[:value]
    setproperty!(ng.nodes[node_ind].params, key, parsed_value)
end

function add_edge!(ng::NodeGraph, edge::Edge)
    edge in ng.edges && return false

    from_ind = id_to_index(edge.from, ng)
    to_ind = id_to_index(edge.to, ng)

    if validate_edge(ng.nodes[from_ind].params, ng.nodes[to_ind].params, edge)
        push!(ng.edges, edge)
        return true
    end
    return false
end

function delete_edge!(ng::NodeGraph, edge::Edge)
    filter!(e -> e != edge, ng.edges)
end


text:
using PyFortran90Namelists
using Test

@testset "PyFortran90Namelists.jl" begin
    include("tokenizer.jl")
end


text:
"""
`parseargs(args...)` parses arguments of `@container`, `@session` and
`@style`. Arguments of type Pair are returned in a dictionnary
Dict{String,String}. Other arguments are returned as strings.

##arguments
- `args...` : arguments as provided by the macro

"""
function parseargs(args...)
  objects = []
  opts = Dict{String,String}()
  for a in args
    if isa(a, Expr) && a.head == :(=>) && isa(a.args[1],Symbol)
      opts[string(a.args[1])] = tostring(a.args[2])
    else
      push!(objects,a)
    end
  end
  (objects, opts)
end

function tostring(a)
  isa(a, String) && return a
  # isa(a, Symbol) && return string(a)
  try
    sa = eval(current_module(), a)
    isa(sa, String) && return sa
    throw()
  catch
    error("cannot parse argument '$a' to a String")
  end
end

"""
`unfold(e)` transforms expressions `:(sym1.sym2.sym3)` into a tuple of symbols
`[:sym1, :sym2, :sym3]`.
"""
unfold(ex::Any) = Symbol[]
unfold(ex::Symbol) = Symbol[ex]
unfold(ex::String) = Symbol[Symbol(ex)]
unfold(ex::QuoteNode) = Symbol[ex.value]
function unfold(ex::Expr)
  if ex.head == :.
    lexpr = unfold(ex.args[1])
    rexpr = unfold(ex.args[2])
    length(lexpr)!=0 || length(rexpr)!=0 || return Symbol[]
    vcat(lexpr,rexpr)
  elseif ex.head == :quote
    Symbol[ex.args[1]]
  else
    Symbol[]
  end
end

# ex = :(abcd.qsdf.sd)
# res = unfold(:(abcd.qsdf.sd))
# unfold(:( aze-45 ))
# unfold(:( aze=45 ))
# unfold(:( aze => abcd.sdf ))
# unfold(:(abcd.yo))

function sendcurrent(command::String, args::Dict=Dict())
  currentSession == nothing && error("[send] no active session")

  send(currentSession,
       currentSession.active_container.nid,
       command, args)
end

function send(session::Session, nid::Int,
              command::String, args::Dict=Dict())
  msg = Dict{Symbol, Any}(:nid     => nid,
                          :command => command,
                          :args    => args)

  put!(session.channel, JSON.json(msg))
end

function htmlpath(path::String)
  npath = normpath(path)
  npath = replace(npath, "\\", "/")
  "file:///" * npath
end

function limited_println(ln::String)
  sl = chomp(ln)
  sl = length(sl)>200 ? sl[1:200]*"..." : sl
  println(sl)
end


text:
abstract type AbstractScenarioProblems{S <: AbstractScenario} end

abstract type AbstractBlockStructure{N} <: AbstractStochasticStructure{N} end

# Getters #
# ========================== #
function scenarioproblems(structure::AbstractBlockStructure{N}, s::Integer = 2) where N
    s == 1 && error("Stage 1 does not have scenario problems.")
    N == 2 && (s == 2 || error("Stage $s not available in two-stage model."))
    1 < s <= N || error("Stage $s not in range 2 to $N.")
    return structure.scenarioproblems[s-1]
end
function decision(structure::AbstractBlockStructure, index::MOI.VariableIndex)
    return decision(structure.decisions[1], index)
end
function decisions(structure::AbstractBlockStructure)
    return structure.decisions[1]
end
function all_decisions(structure::AbstractBlockStructure)
    return structure.decisions[1].undecided
end
function num_decisions(structure::AbstractBlockStructure{N}, s::Integer = 1) where N
    1 <= s < N || error("Stage $s not in range 1 to $(N - 1).")
    return num_decisions(structure.decisions[s])
end
function scenario(structure::AbstractBlockStructure, i::Integer, s::Integer = 2)
    scenario(scenarioproblems(structure, s), i)
end
function scenarios(structure::AbstractBlockStructure, s::Integer = 2)
    scenarios(scenarioproblems(structure, s))
end
function expected(structure::AbstractBlockStructure, s::Integer = 2)
    return expected(scenarioproblems(structure, s))
end
function scenario_type(structure::AbstractBlockStructure, s::Integer = 2)
    return scenario_type(scenarioproblems(structure, s))
end
function probability(structure::AbstractBlockStructure, i::Integer, s::Integer = 2)
    return probability(scenarioproblems(structure, s), i)
end
function stage_probability(structure::AbstractBlockStructure, s::Integer = 2)
    return probability(scenarioproblems(structure, s))
end
function subproblem(structure::AbstractBlockStructure, i::Integer, s::Integer = 2)
    return subproblem(scenarioproblems(structure, s), i)
end
function subproblems(structure::AbstractBlockStructure, s::Integer = 2)
    return subproblems(scenarioproblems(structure, s))
end
function num_subproblems(structure::AbstractBlockStructure, s::Integer = 2)
    return num_subproblems(scenarioproblems(structure, s))
end
function num_scenarios(structure::AbstractBlockStructure, s::Integer = 2)
    return num_scenarios(scenarioproblems(structure, s))
end
deferred(structure::AbstractBlockStructure{N}) where N = deferred(structure, Val(N))
deferred(structure::AbstractBlockStructure, ::Val{1}) = deferred_first_stage(structure)
function deferred(structure::AbstractBlockStructure, ::Val{N}) where N
    return deferred_stage(structure, N) || deferred(structure, Val(N-1))
end
deferred_first_stage(structure::AbstractBlockStructure) = false
function deferred_stage(structure::AbstractBlockStructure{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == 1 && return deferred_first_stage(structure)
    num_subproblems(structure, s) < num_scenarios(structure, s)
end
function distributed(structure::AbstractBlockStructure{N}, s) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == 1 && return false
    return distributed(scenarioproblems(structure, s))
end

# ========================== #

# Setters
# ========================== #
function add_scenario!(structure::AbstractBlockStructure, scenario::AbstractScenario, stage::Integer = 2)
    add_scenario!(scenarioproblems(structure, stage), scenario)
    return nothing
end
function add_worker_scenario!(structure::AbstractBlockStructure, scenario::AbstractScenario, w::Integer, stage::Integer = 2)
    add_scenario!(scenario(structure, stage), scenario, w)
    return nothing
end
function add_scenario!(scenariogenerator::Function, structure::AbstractBlockStructure, stage::Integer = 2)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage))
    return nothing
end
function add_worker_scenario!(scenariogenerator::Function, structure::AbstractBlockStructure, w::Integer, stage::Integer = 2)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage), w)
    return nothing
end
function add_scenarios!(structure::AbstractBlockStructure, scenarios::Vector{<:AbstractScenario}, stage::Integer = 2)
    add_scenarios!(scenarioproblems(structure, stage), scenarios)
    return nothing
end
function add_worker_scenarios!(structure::AbstractBlockStructure, scenarios::Vector{<:AbstractScenario}, w::Integer, stage::Integer = 2)
    add_scenarios!(scenarioproblems(structure, stage), scenarios, w)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, structure::AbstractBlockStructure, n::Integer, stage::Integer = 2)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n)
    return nothing
end
function add_worker_scenarios!(scenariogenerator::Function, structure::AbstractBlockStructure, n::Integer, w::Integer, stage::Integer = 2)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n, w)
    return nothing
end
function sample!(structure::AbstractBlockStructure, sampler::AbstractSampler, n::Integer, stage::Integer = 2)
    sample!(scenarioproblems(structure, stage), sampler, n)
    return nothing
end
# ========================== #

# Includes
# ========================== #
include("scenarioproblems.jl")
include("vertical.jl")
include("horizontal.jl")


text:
using BinaryBuilder

name = "nlminb"
version = v"0.1.0"

# Collection of sources required to build NLopt
sources = [
    GitSource("https://github.com/geo-julia/nlminb.f.git",
              "6b2d460b866911e8bfca32b1217c08fb78c3111f"), # v0.1.0
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nlminb.f

## compile by hand
if [[ "${proc_family}" == "intel" ]]; then
    FLAGS="-mfpmath=sse -msse2 -mstackrealign"
fi
CFLAGS="-fPIC -DNDEBUG  -Iinclude -O2 -Wall -std=gnu99 ${FLAGS}"
FFLAGS="-fPIC -fno-optimize-sibling-calls -O2 ${FLAGS}"

gfortran ${FFLAGS} -c src/portsrc.f -o src/portsrc.o
gfortran ${FFLAGS} -c src/d1mach.f -o src/d1mach.o
gcc ${CFLAGS} -c src/port2.c -o src/port.o
gfortran -shared -static-libgcc -lm -o ${libdir}/libnlminb.${dlext} src/portsrc.o src/d1mach.o src/port.o -lopenblas
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
mutable struct PowerLawModel{T <: SSDFloat} <: AbstractTemperatureModel{T}

    #temperatures needed for the drift velocity rescaling
    temperature::T      # actual temperature of the crystal
    reftemperature::T   # temperature at which the reference values were measured

    #fit parameters for electrons (e) and holes (h) for <100> and <111> axes
    p1e100::T
    p1e111::T
    p1h100::T
    p1h111::T
end

function PowerLawModel(config_file::Dict; T::Type{<:AbstractFloat} = Float32)::PowerLawModel
    return PowerLawModel{T}(config_file::Dict)
end

function PowerLawModel{T}(config_file::Dict; temperature::Union{Missing, T} = missing)::PowerLawModel where {T <: SSDFloat}
    config_file["temperature_dependence"]["model"] != "PowerLaw" ? error() : nothing
    ismissing(temperature) ? temperature = config_file["temperature_dependence"]["temperature"] : nothing
    m = PowerLawModel{T}(
        temperature, #temperature = config_file["temperature_dependence"]["temperature"]
        config_file["drift"]["velocity"]["temperature"],
        config_file["temperature_dependence"]["parameters"]["e100"]["p1"],
        config_file["temperature_dependence"]["parameters"]["e111"]["p1"],
        config_file["temperature_dependence"]["parameters"]["h100"]["p1"],
        config_file["temperature_dependence"]["parameters"]["h111"]["p1"]
    )
end

function scale_to_given_temperature(m::PowerLawModel{T})::NTuple{4,T} where {T <: SSDFloat}

    function f(p1::T = 0.0)::T where {T <: SSDFloat}
        return (m.reftemperature/m.temperature)^1.5
    end

    scale_e100 = f(m.p1e100)
    scale_e111 = f(m.p1e111)
    scale_h100 = f(m.p1h100)
    scale_h111 = f(m.p1h111)

    return scale_e100, scale_e111, scale_h100, scale_h111
end


text:
using LowLevelFloatFunctions
using Test

sqrt2₆₄ = sqrt(2.0); sqrt2₃₂ = sqrt(2.0f0); sqrt2₁₆ = sqrt(Float16(2.0));

@testset "constants" begin
    @test precision(Float16) == 11
    @test precision(Float32) == 24
    @test precision(Float64) == 53

    @test sign_bits(Float16) == 1
    @test sign_bits(Float32) == 1
    @test sign_bits(Float64) == 1

    @test significand_bits(Float16) == 10
    @test significand_bits(Float32) == 23
    @test significand_bits(Float64) == 52

    @test exponent_bits(Float16) == 5
    @test exponent_bits(Float32) == 8
    @test exponent_bits(Float64) == 11

    @test exponent_max(Float16) == 15
    @test exponent_max(Float32) == 127
    @test exponent_max(Float64) == 1023

    @test exponent_bias(Float16) == 15
    @test exponent_bias(Float32) == 127
    @test exponent_bias(Float64) == 1023

    @test exponent_field_max(Float16) == 16
    @test exponent_field_max(Float32) == 128
    @test exponent_field_max(Float64) == 1024

    @test exponent_min(Float16) == -14
    @test exponent_min(Float32) == -126
    @test exponent_min(Float64) == -1022
end

@testset "value_extraction" begin
    @test sign(-sqrt2₆₄) === -1.0
    @test sign(sqrt2₃₂) === 1.0f0
    @test sign(-sqrt2₁₆) === Float16(-1.0)
    @test exponent(-sqrt2₆₄) == 0
    @test exponent(sqrt2₃₂) == 0
    @test exponent(-sqrt2₁₆) == 0
    @test significand(-sqrt2₆₄) === -1.4142135623730951
    @test significand(sqrt2₃₂) === 1.4142135f0
    @test significand(-sqrt2₁₆) === Float16(-1.414)
end

@testset "field_get" begin
    @test sign_field(-sqrt2₆₄) === 0x0000000000000001
    @test sign_field(sqrt2₃₂) === 0x00000000
    @test sign_field(-sqrt2₁₆) === 0x0001
    @test exponent_field(-sqrt2₆₄) === 0x00000000000003ff
    @test exponent_field(sqrt2₃₂) === 0x0000007f
    @test exponent_field(-sqrt2₁₆) === 0x000f
    @test significand_field(sqrt2₆₄) === 0x0006a09e667f3bcd
    @test significand_field(sqrt2₃₂) === 0x003504f3
    @test significand_field(sqrt2₁₆) === 0x01a8
end

@testset "field_set" begin
    @test sign_field(-sqrt2₆₄, 0%UInt64) === 1.4142135623730951
    @test exponent_field(sqrt2₆₄, exponent_field(sqrt2₆₄)+one(UInt64)) === 2.8284271247461903
    @test significand_field(sqrt2₃₂, significand_field(sqrt2₃₂) - one(UInt32)) === prevfloat(sqrt2₃₂)
end

@testset "characterization" begin
    @test sign_bits(Float64) == 1
    @test exponent_bits(Float32) == 8
    @test significand_bits(Float16) == 10
    @test exponent_field_max(Float64) === 0x0000000000000400
    @test exponent_max(Float64) == 1023
    @test exponent_min(Float64) == -1022
    @test exponent_bias(Float32) == 127
end

@testset "utilitiarian" begin
    @test bitwidth(Float64) == 64
    @test bitwidth(Float32) == 32
end


text:
const AVec = AbstractVector
const AMat = AbstractMatrix
const AArr = AbstractArray

"""
    iswet(grd)

Returns the 3D BitArray of wet boxes of the grid.
"""
iswet(g::T) where {T<:OceanGrid} = g.wet3D
export iswet

"""
    latvec(grd)

Returns the vector (wet points) of latitudes (units stripped).
"""
latvec(g::T) where {T<:OceanGrid} = ustrip.(g.lat_3D[iswet(g)])

"""
    lonvec(grd)

Returns the vector (wet points) of longitudes (units stripped).
"""
lonvec(g::T) where {T<:OceanGrid} = ustrip.(g.lon_3D[iswet(g)])

"""
    depthvec(grd)

Returns the vector (wet points) of depths (units stripped).
"""
depthvec(g::T) where {T<:OceanGrid} = ustrip.(g.depth_3D[iswet(g)])

"""
    topdepthvec(grd)

Returns the vector (wet points) of top depths (units stripped).
(top depths = depth of the top of the box.)
"""
topdepthvec(g::T) where {T<:OceanGrid} = ustrip.(g.depth_top_3D[iswet(g)])

"""
    bottomdepthvec(grd)

Returns the vector (wet points) of bottom depths (units stripped).
(bottom depths = depth of the bottom of the box.)
"""
bottomdepthvec(g::T) where {T<:OceanGrid} = 2depthvec(g) - topdepthvec(g)
function isseafloorvec(g::T) where {T<:OceanGrid}
    nlat, nlon, ndepth = size(g)
    iswetext = cat(iswet(g), falses(nlat,nlon,1), dims=3)
    isseafloor3D =  @views iswetext[:,:,1:ndepth] .& .!iswetext[:,:,2:ndepth+1]
    return isseafloor3D[iswet(g)]
end

"""
    latlondepthvecs(g)

Returns `latvec(g), lonvec(g), depthvec(g)`.
"""
latlondepthvecs(g::T) where {T<:OceanGrid} = latvec(g), lonvec(g), depthvec(g)

export latvec, lonvec, depthvec, latlondepthvecs, topdepthvec, bottomdepthvec
export isseafloorvec

#===================================
1D Vector <-> 3D array conversions
===================================#
"""
    rearrange_into_3Darray(x, wet3D::BitArray)

Returns a 3D array of `x` rearranged to the `true` entries of `wet3D`.
Entries where `wet3D` is `false` are filled with `NaN`s.
"""
function rearrange_into_3Darray(x, wet3D::BitArray)
    iwet = findall(wet3D)
    x3d = fill(NaN * unit(eltype(x)), size(wet3D))
    x3d[iwet] .= x
    return x3d
end
"""
    rearrange_into_3Darray(x, grd)

Returns a 3D array of `x` rearranged to the wet boxes of the grid.
"""
rearrange_into_3Darray(x, grd) = rearrange_into_3Darray(x, grd.wet3D)

"""
    vectorize(x3D, grd)

Returns the 1D vector corresponding to the 3D field of `x3D` at the wet boxes of `grd`.

This function essentially returns `x3D[iswet(grd)]` but with some additional checks.
"""
function vectorize(x3D::AArr{T,3} where T, grd)
    size(x3D) ≠ size(grd) && error("The size of your 3D array does not match the grid size")
    x3D[iswet(grd)]
end

export rearrange_into_3Darray, vectorize

#============
Interpolation
============#
"""
    interpolate(x, g)

Returns an "interpolation object" for fast interpolation.

It must be queried within the range of the grid in longitude...
"""
function Interpolations.interpolate(x, g::T) where {T<:OceanGrid}
    knots = (ustrip.(g.lat), periodic_longitude(ustrip.(g.lon)), ustrip.(g.depth))
    A = periodic_longitude(rearrange_into_3Darray(x, g))
    itp = interpolate(knots, A, Gridded(Constant()))
    return extrapolate(itp, (Flat(), Periodic(), Flat())) # extrapolate periodically
end
function Interpolations.interpolate(x, g::OceanGrid, lat, lon, depth; itp=interpolate(x,g))
    return itp(ustrip.(lat), ustrip.(lon), ustrip.(depth))
end
function Interpolations.interpolate(x, g::OceanGrid, lats::Vector, lons::Vector, depths::Vector; itp=interpolate(x,g))
    return [itp(y,x,z) for (y,x,z) in zip(ustrip.(lats), ustrip.(lons), ustrip.(depths))]
end
function Interpolations.interpolate(x, g::OceanGrid, obs; itp=interpolate(x,g))
    return interpolate(x, g, obs.lat, obs.lon, obs.depth; itp=itp)
end
Interpolations.interpolate(x, g::OceanGrid, ::Missing; itp=nothing) = missing
export interpolate

"""
    iswet(g, lat, lon, depth)
    iswet(g, lats, lons, depths)
    iswet(g, obs)

Returns a `BitArray` of the provided locations that are "in" a wet box of `g`.

Technically, this uses a nearest neighbour interpolation algorithm, so the bounding boxes
will not work perfectly.
If AIBECS will solely rely on this nearest neighbour interpolation,
then it might be a good idea to replace Interpolations.jl with NearestNeighbors.jl.
"""
function iswet(g, args...; itp=interpolate(1:count(iswet(g)),g))
    J = interpolate(1:count(iswet(g)), g, args...; itp=itp)
    return (!isnan).(J)
end
iswet(g, ::Missing; itp=nothing) = missing

"""
    interpolationmatrix(g, lats, lons, depths)
    interpolationmatrix(g, obs)

Returns the sparse matrix `M` such that `M*x` is the model vector `x` interpolated
onto the provided locations in the obs table.

Technically, this requires a linear interpolation, in the sense that
`interpolation(x) = M * x` for some `M`, which seems to require a nearest neighbor
interpolation.
If AIBECS will solely rely on this nearest neighbour interpolation,
then it might be a good idea to replace Interpolations.jl with NearestNeighbors.jl.
Also, the few observations that lie in model land (instead of water) are discarded.

In the future, a more generic approach will use an sparse-aware AD for any type of
interpolation, and use an interpolation that somehow does not discard observations
made at locations not inside an ocean grid cell.
"""
function interpolationmatrix(g, args...; itp=interpolate(1:count(iswet(g)),g))
    J = interpolate(1:count(iswet(g)), g, args...; itp=itp)
    iwet = findall((!isnan).(J))
    J = Int.(J[iwet])
    I = collect(1:length(iwet))
    return sparse(I, J, true, length(I), count(iswet(g)))
end
interpolationmatrix(g, ::Missing; itp=nothing) = missing
export interpolationmatrix

#=
    lats = range(ustrip(grd.lat[1]), ustrip(grd.lat[end]), length=length(grd.lat))
    lons = range(ustrip(grd.lon[1]), ustrip(grd.lon[1])+360, length=length(grd.lon)+1)
    stp = Interpolations.scale(itp, lats, lons, 1:ndepths) # re-scale to the actual domain
    etp = extrapolate(stp, (Line(), Periodic(), Line())) # periodic longitude

    xs = range(0, 2π, length=11)
    ys = range(-π/2, π/2, length=6)
    A = f.(xs,ys')
    itp = interpolate(A, BSpline(Linear())) # interpolate linearly between the data points
    stp = scale(itp, xs, ys) # re-scale to the actual domain
    etp = extrapolate(stp, (Periodic(), Throw())) # extrapolate periodically
=#

periodic_longitude(x3D::AArr{T,3}) where T = view(x3D,:,[1:size(x3D,2); 1],:)
periodic_longitude(v::AVec{T}) where T = [v; v[1]+360]

#===================================
Off-diaongal matrices
===================================#

"""
    buildIbelow(wet3D, iwet)

Builds the shifted-diagonal sparse matrix of the indices of below neighbours.

Ibelow[i,j] = 1 if the box represented by the linear index i
lies directly below the box represented by the linear index j.
"""
function buildIbelow(wet3D, iwet)
    nlat, nlon, ndepth = size(wet3D)
    n = nlon * nlat * (ndepth + 1)
    In = sparse(I, n, n)
    idx = zeros(Int, nlat, nlon, ndepth + 1)
    idx[:] = 1:n
    idx .= idx[:, :, [2:ndepth + 1; 1]]      # downward shift
    return In[idx[:], :][iwet, iwet]
end
"""
    buildIbelow(grd)

Builds the shifted-diagonal sparse matrix of the indices of below neighbours for `grid`.
See `buildIbelow(wet3D, iwet)`.
"""
buildIbelow(grd) = buildIbelow(grd.wet3D, indices_of_wet_boxes(grd))
export buildIbelow

"""
    buildIabove(wet3D, iwet)

Builds the shifted-diagonal sparse matrix of the indices of above neighbours.

Iabove[i,j] = 1 if the box represented by the linear index i
lies directly above the box represented by the linear index j.
"""
buildIabove(wet3D, iwet) = copy(transpose(buildIbelow(wet3D, iwet)))
"""
    buildIabove(grd)

Builds the shifted-diagonal sparse matrix of the indices of above neighbours for `grd`.
See `buildIabove(wet3D, iwet)`.
"""
buildIabove(grd) = buildIabove(grd.wet3D, indices_of_wet_boxes(grd))
export buildIabove

#===================================
Functions returning useful arrays,
vectors, and constants
===================================#

"""
    array_of_volumes(grd)

Returns the 3D array of volumes of grd boxes.
"""
array_of_volumes(grd) = grd.volume_3D
export array_of_volumes

"""
    indices_of_wet_boxes(wet3D::BitArray)

Returns the vector of the indices of wet grid boxes.
"""
indices_of_wet_boxes(wet3D::BitArray) = findall(vec(wet3D))
"""
    indices_of_wet_boxes(grd)

Returns the vector of the indices of wet grid boxes.
"""
indices_of_wet_boxes(grd) = indices_of_wet_boxes(grd.wet3D)
export indices_of_wet_boxes

"""
    number_of_wet_boxes(wet3D::BitArray)

Returns the number of wet grid boxes.
"""
number_of_wet_boxes(wet3D::BitArray) = sum(wet3D)
"""
    number_of_wet_boxes(grd)

Returns the number of wet grid boxes.
"""
number_of_wet_boxes(grd) = number_of_wet_boxes(grd.wet3D)
export number_of_wet_boxes

"""
    vector_of_volumes(grid)

Returns the vector of volumes of wet boxes.
"""
vector_of_volumes(grd) = array_of_volumes(grd)[indices_of_wet_boxes(grd.wet3D)]
volumevec(grd) = ustrip.(vector_of_volumes(grd))
export vector_of_volumes, volumevec

"""
    vector_of_depths(grd)

Returns the vector of depths of the center of wet boxes.
"""
vector_of_depths(grd) = grd.depth_3D[indices_of_wet_boxes(grd.wet3D)]
export vector_of_depths

"""
    vector_of_top_depths(grd)

Returns the vector of depths of the top of wet boxes.
"""
vector_of_top_depths(grd) = grd.depth_top_3D[indices_of_wet_boxes(grd.wet3D)]
export vector_of_top_depths

"""
    vector_of_bottom_depths(grd)

Returns the vector of depths of the bottom of wet boxes.
"""
vector_of_bottom_depths(grd) = 2vector_of_depths(grd) - vector_of_top_depths(grd)
export vector_of_bottom_depths

#===================================
Functions return weighted norms
===================================#

"""
    weighted_norm²(x, w)

Returns the square of the weighted norm of `x` using weights `w`.
"""
weighted_norm²(x, w) = transpose(x) * Diagonal(w) * x

"""
    weighted_norm(x, v)

Returns the weighted norm of `x` using weights `w`.
"""
weighted_norm(x, w) = sqrt(weighted_norm²(x, w))

"""
    weighted_mean(x, w)

Returns the weighted mean of `x` using weights `w`.
"""
weighted_mean(x, w) = transpose(w) * x / sum(w)

"""
    vmean(x, w, I)

Returns the weighted mean of `x` using weights `w`, but only over indices `I`.
This is useful to get the observed mean.
(Because there are some grid boxes without observations.)
"""
weighted_mean(x, w, I) = transpose(w[I]) * x[I] / sum(w[I])
export weighted_norm²

#===============
Slices and views
===============#
# depth helper functions
finddepthindex(depth::Quantity, grd) = findmin(abs.(grd.depth .- depth))[2]
convertdepth(depth::Real) = depth * u"m"
convertdepth(depth::Quantity{U, Unitful.𝐋, V}) where {U,V} = depth
convertdepth(x) = error("Not a valid depth")
finddepthindex(depths::Tuple, grd) = ((finddepthindex(d, grd) for d in depths)...,)
convertdepth(depths::Tuple) = ((convertdepth(d) for d in depths)...,)
# lon helper functions
findlonindex(lon::Quantity, grd) = findmin(abs.(grd.lon .- lon))[2]
convertlon(x::Real) = mod(x * °, 360°)
convertlon(x::Quantity) = unit(x) == ° ? mod(x, 360°) : error("Not a valid lon")
# lat helper functions
findlatindex(lat::Quantity, grd) = findmin(abs.(grd.lat .- lat))[2]
convertlat(y::Real) = y * °
convertlat(y::Quantity) = unit(y) == ° ? y : error("Not a valid lat")
export convertlat, convertlon
"""
    HorizontalSlice(x, grd; depth)

Returns a horizontal slice of tracer `x` at depth `depth`.
"""
function horizontalslice(x3D, grd; depth=nothing)
    isnothing(depth) && error("you must specify the depth, e.g., `depth=100`")
    depth = convertdepth(depth)
    iz = finddepthindex(depth, grd)
    return view(x3D, :, :, iz)
end
"""
    meridionalslice(x, grd; lon)

Returns a meridional slice of tracer `x` at longitude `lon`.
"""
function meridionalslice(x3D, grd; lon=nothing)
    isnothing(lon) && error("you must specify the longitude, e.g., `lon=100`")
    lon = convertlon(lon)
    ix = findlonindex(lon, grd)
    return PermutedDimsArray(view(x3D,:,ix,:), (2,1))
end
"""
    zonalslice(x, grd; lat)

Returns a zonal slice of tracer `x` at latitidue `lat`.
"""
function zonalslice(x3D, grd; lat=nothing)
    isnothing(lat) && error("you must specify the latitude, e.g., `lat=-30`")
    lat = convertlat(lat)
    iy = findlatindex(lat, grd)
    return PermutedDimsArray(view(x3D,iy,:,:), (2,1))
end

"""
    depthprofile(x, grd; lonlat)

Returns the profile of tracer `x` interpolated at `lonlat=(x,y)` coordinates.
"""
function depthprofile(x, grd; lonlat=nothing)
    isnothing(lonlat) && error("you must specify the coordinates, e.g., `lonlat=(-30,30)`")
    lon, lat = lonlat
    lon, lat = convertlon(lon), convertlat(lat)
    x3D = rearrange_into_3Darray(x, grd)
    iy, ix = findlatindex(lat, grd), findlonindex(lon, grd)
    return view(x3D,iy,ix,:)
end
export depthprofile

#=====================
Integrals and averages
=====================#
"""
    integral(x, grd; mask=iswet(grd), dims=(1,2,3))

Returns a "generic" integral of tracer `x` along `dims`.
"""
function integral(x3D, grd, mask3D=iswet(grd); dims=(1,2,3))
    w3D = (1 ∈ dims ? grd.δy_3D : 1) .* (2 ∈ dims ? grd.δx_3D : 1) .* (3 ∈ dims ? grd.δz_3D : 1)
    return dropdims(nansum(x3D .* w3D .* mask3D, dims=dims) ./ any(.!isnan.(mask3D), dims=dims), dims=dims)
end
nansum(x; kwargs...) = sum(x .* .!isnan.(x); kwargs...)
∫dxdydz(x3D, grd, mask3D=iswet(grd)) = integral(x3D, grd, mask3D)[1]
∫dx(x3D, grd, mask3D=iswet(grd))     = integral(x3D, grd, mask3D; dims=2)
∫dy(x3D, grd, mask3D=iswet(grd))     = integral(x3D, grd, mask3D; dims=1)
∫dz(x3D, grd, mask3D=iswet(grd))     = integral(x3D, grd, mask3D; dims=3)
∫dxdy(x3D, grd, mask3D=iswet(grd))   = integral(x3D, grd, mask3D; dims=(1,2))
∫dxdz(x3D, grd, mask3D=iswet(grd))   = integral(x3D, grd, mask3D; dims=(2,3))
∫dydz(x3D, grd, mask3D=iswet(grd))   = integral(x3D, grd, mask3D; dims=(1,3))

# integral convenience aliases
totalintegral = ∫dxdydz
∫dV = ∫dxdydz
∫d³r = ∫dxdydz
zonalintegral = ∫dx
meridionalintegral = ∫dy
verticalintegral = ∫dz
horizontalintegral = ∫dxdy
export ∫dxdydz, ∫dV, ∫d³r, ∫dx, ∫dy, ∫dz, ∫dxdy, ∫dydz, ∫dxdz

# Create mean/average functions from integral functions
for f in (:horizontal, :vertical, :meridional, :zonal, :total)
    fint = Symbol(string(f, "integral"))
    fmean = Symbol(string(f, "mean"))
    favg = Symbol(string(f, "average"))
    fstd = Symbol(string(f, "std"))
    @eval begin
        $fmean(x, grd, args...; kwargs...) = $fint(x, grd, args...; kwargs...) ./ $fint(1.0, grd, args...; kwargs...)
        $favg = $fmean
        $fstd(x, grd, args...; kwargs...) = sqrt.($favg(abs2.(x - paint3D($favg(x, grd, args...; kwargs...), grd)), grd, args...; kwargs...))
        export $favg
    end
end

"""
    paint3D(x, grd)

Returns a 3D array of the size of `grd` that copies the values of `x` along the dimensions it lacks.

This is useful, e.g., when `x` is a 1D profile, a 2D map, or a 2D zonal average and you want to "copy"
it onto the full 3D array.
"""
function paint3D(x, grd)
    length(unique(size(grd))) ≠ length(size(grd)) && @warn "paint3D might not work because of the grid size"
    dimsout = [findfirst(nx .== collect(size(grd))) for nx in size(x)]
    szin, szout = [1,1,1], collect(size(grd))
    szout[dimsout] .= 1
    szin[dimsout] .= size(grd)[dimsout]
    return repeat(reshape(x, Tuple(szin)), outer=szout)
end
paint3D(x::Number, grd) = x * iswet(grd)
paint3D(x::AArr{T,3} where T, grd) = x

# Allow all the integral/average functions to work when `x` is supplied as a vector
for f in (:∫dxdy,    :horizontalmean, :horizontalstd,
          :∫dz,      :verticalmean,   :verticalstd,
          :∫dy,      :meridionalmean, :meridionalstd,
          :∫dx,      :zonalmean,      :zonalstd,
          :∫dxdydz,  :totalmean,      :totalstd)
    @eval begin
        $f(x::AVec, grd, mask=1; kwargs...) = $f(rearrange_into_3Darray(x, grd), grd, rearrange_into_3Darray(mask, grd); kwargs...)
        export $f
    end
end

#================
Vertical sections
================#

# Functions extending 3D arrays and 1D longitudinal vectors
# to be used by interpolation functions.
# These assume that the longitude is in (0,360) and that
# the longitude is stored in the 2nd dimension (lat,lon,depth)
function lonextend(x3D::AArr{T,3}) where T
    x3Dext = Array{T,3}(undef, (size(x3D) .+ (0,2,0))...)
    x3Dext[:,2:end-1,:] .= x3D
    x3Dext[:,1,:] .= x3D[:,end,:]
    x3Dext[:,end,:] .= x3D[:,1,:]
    return x3Dext
end
function cyclicallon(lon::AVec{T}) where T<:Quantity
    lonext = Vector{T}(undef, length(lon) + 2)
    lonext[2:end-1] .= lon
    lonext[1] = lon[end] - 360°
    lonext[end] = lon[1] + 360°
    return lonext
end
cyclicallon(lon::AVec) = cyclicallon(convertlon.(lon))
# TODO make the function below more generic and be capable of taking in other
# representations than just OceanographyCruises.jl `CruiseTrack`
"""
    verticalsection(x3D, grd; ct=nothing)

Returns a tuple of `(d, depth, x2D)` where the 2D array is the vertical section that follows the cruise track `ct` and `d` is the distance.
This is made to work with the OceanographyCruises.jl package.
"""
function verticalsection(x3D, grd; ct=nothing)
    isnothing(ct) && error("you must include the cruise track with `ct=...`")
    # Create the interpolation object
    x3D = lonextend(x3D)
    knots = ustrip.((grd.lat, cyclicallon(grd.lon), grd.depth))
    itp = interpolate(knots, x3D, Gridded(Constant()))
    # convert lat and lon (TODO check if necessary)
    ctlats = [convertlat(st.lat) for st in ct.stations]
    ctlons = [convertlon(st.lon) for st in ct.stations]
    # compute distance
    distance = vcat(0.0u"km", cumsum(diff(ct)))
    # remove identical stations (lat,lon)
    idx = findall(diff(distance) .> 0u"km")
    push!(idx, length(ct))
    d = view(distance, idx)
    y = view(ctlats, idx)
    x = view(ctlons, idx)
    return d, grd.depth, [itp(ustrip.((lat, lon, d))...) for d in grd.depth, (lat,lon) in zip(y, x)]
end

function verticalsection2(x3D, grd; ct=nothing)
    isnothing(ct) && error("you must include the cruise track with `ct=...`")
    # convert lat and lon (TODO check if necessary)
    ctlats = [convertlat(st.lat) for st in ct.stations]
    ctlons = [convertlon(st.lon) for st in ct.stations]
    bl = baselon(ctlons)
    ctlons = shiftlon(ctlons, bl=bl)
    # compute distance
    distance = vcat(0.0u"km", cumsum(diff(ct)))
    # remove identical stations (lat,lon)
    idx = findall(diff(distance) .> 0u"km")
    push!(idx, length(ct))
    d = view(distance, idx)
    y = view(ctlats, idx)
    x = view(ctlons, idx)
	# intersections
    xe = lonedges(grd)
    ye = latedges(grd)
	dix = intersections(d, x, shiftlon(xe, bl=bl))
	diy = intersections(d, y, ye)
	ext_d = [d[1]; sort(vcat(dix, diy)); d[end]] # extended distance with intersections values
	# mid points (to get z value between intersections)
	mid_d = [0.5(ext_d[i] + ext_d[i+1]) for i in 1:length(ext_d)-1]
	itpx = LinearInterpolation(d, x)
    itpy = LinearInterpolation(d, y)
	x2 = itpx(mid_d)
    y2 = itpy(mid_d)
	ix = [searchsortedfirst(xe, shiftlon(x, bl=0°), lt=≤) - 1 for x in x2]
    iy = [searchsortedfirst(ye, y) - 1 for y in y2]
    return ext_d, depthedges(grd), [x3D[j,i,k] for k in 1:length(grd.depth), (j,i) in zip(iy,ix)]
end

# aux functions for vertical section delimited by edges of boxes
edges(x::AVec, dx::AVec) = [x .- 0.5 .* dx; x[end] + 0.5dx[end]]
edges(grd::OceanGrid) = lonedges(grd), latedges(grd), depthedges(grd)
lonedges(grd::OceanGrid) = edges(grd.lon, grd.δlon)
latedges(grd::OceanGrid) = edges(grd.lat, grd.δlat)
depthedges(grd::OceanGrid) = edges(grd.depth, grd.δdepth)

#======================================================
Functions for plotting scatter transect on top of section heatmap
======================================================#
"""
    intersections(x::Vector{T}, y, ylevs)

Returns an array of values in `x` coordinates of the intersections
of the linearly interpolated plot of `y` as a function of `x` with the `ylevs`.
"""
function intersections(x::AArr{T}, y, ylevs) where T
	out = T[]
	for i in 1:length(x) - 1
		mini_x = view(x, i:i+1)
		mini_y = view(y, i:i+1)
		idx = sortperm(mini_y)
		a, b = view(mini_y, idx)
        if a ≠ b
		    itp = LinearInterpolation(mini_y[idx], mini_x[idx], extrapolation_bc=Line())
		    for x_intersect in itp(ylevs[a .< ylevs .≤ b])
		    	push!(out, x_intersect)
		    end
        end
	end
	return out
end
shiftlon(x; bl=baselon(x)) = @. mod(mod(x - bl, 360°), 360°) + bl # enforces 0 ≤ m < 360 (https://github.com/JuliaLang/julia/issues/36310)
baselon(x) = (any(0° .≤ x .< 90°) && any(270° .≤ x .< 360°)) ? -180.0° : 0.0°

#=======================================================
Allow functions to work when `x` is supplied as a vector
=======================================================#

for f in (:horizontalslice, :meridionalslice, :zonalslice, :verticalsection, :verticalsection2)
    @eval begin
        $f(x::AVec, grd; kwargs...) = $f(rearrange_into_3Darray(x, grd), grd; kwargs...)
        export $f
    end
end

#=======================================================
Functions to grid/interpolate
=======================================================#

function lonextend(x2D::AArr{T,2}) where T
    x2Dext = Array{eltype(x2D),2}(undef, (size(x2D) .+ (0,2))...)
    x2Dext[:,2:end-1] .= x2D
    x2Dext[:,1] .= x2D[:,end]
    x2Dext[:,end] .= x2D[:,1]
    return x2Dext
end

"""
    regrid(x2D, lat, lon, grd)
    regrid(x3D, lat, lon, depth, grd)

Returns `x2D` (or `x3D`) interpolated onto `grd`.
"""
function regrid(x2D::AArr{T,2} where T, lat, lon, grd)
    size(x2D) ≠ (length(lat), length(lon)) && error("Dimensions of input and lat/lon don't match")
    # need to rearrange in case data was in (-180, 180)
    lon = convertlon.(lon)
    ilon = sortperm(lon)
    sort!(lon)
    x2D = view(x2D, :, ilon)
    knots = convertlat.(lat), cyclicallon(lon)
    x2D = lonextend(x2D)
    itp = LinearInterpolation(knots, x2D, extrapolation_bc = Flat())
    return [itp(y, x) for y in grd.lat, x in grd.lon]
end
function regrid(x3D::AArr{T,3} where T, lat, lon, depth, grd; interpolate_nans=false)
    size(x3D) ≠ (length(lat), length(lon), length(depth)) && error("Dimensions of input and lat/lon/depth don't match")
    if interpolate_nans
        any(isnan, x3D) && @warn "Be aware that NaNs will propagate with interpolation!"
    else
        x3D = inpaint(x3D)
        @warn "NaNs have been inpainted!"
    end
    x3D = lonextend(x3D)
    knots = convertlat.(lat), cyclicallon(lon), convertdepth.(depth)
    itp = LinearInterpolation(knots, x3D, extrapolation_bc = Flat())
    return [itp(y, x, z) for y in grd.lat, x in grd.lon, z in grd.depth]
end

#TODO Rename to rebin
function regrid(vs::AVec{T}, lats, lons, depths, grd) where T
    out = zeros(T, count(iswet(grd)))
    for (i, v) in zip(regrid_indices(lats, lons, depths, tree(grd)), vs)
        out[i] += v
    end
    return out
end
function regrid_indices(lats, lons, depths, tree)
    nind = knn(tree, [mod.(lons, 360)'; lats'; depths'], 1)[1]
    return [ind[1] for ind in nind]
end
tree(grd) = KDTree([lonvec(grd)'; latvec(grd)'; depthvec(grd)'])

export regrid

function regridandpaintsurface(x2D::AArr{T,2} where T, lat, lon, grd)
    x2D2 = regrid(x2D, lat, lon, grd)
    x3D = zeros(size(grd))
    x3D[:,:,1] .= x2D2
    return vectorize(x3D, grd)
end
export regridandpaintsurface

"""
    zcumsum(x, grd)

Returns the cumulative sum of `x` along the `z` dimension.
"""
function zcumsum(x, grd)
    x3D = rearrange_into_3Darray(x, grd)
    return vectorize(cumsum(x3D, dims=3), grd)
end
export zcumsum


text:
function updateμ!{T <: AbstractFloat}(m::GeneralizedLinearMixedModel{T})
    y, dist, link, η, μ = m.y, m.dist, m.link, m.η, m.μ
    wt, wrkresid, wrkwt, dres = m.wt, m.wrkresid, m.wrkwt, m.devresid

    priorwts = !isempty(wt)

    @inbounds Threads.@threads for i = eachindex(η)
        yi, ηi = y[i], η[i]
        μi = μ[i] = linkinv(link, ηi)
        dμdηi = mueta(link, ηi)
        vari = GLM.glmvar(dist, link, μi, ηi)
        wrkresid[i] = (yi - μi)/dμdηi
        wti = priorwts ? wt[i] : one(T)
        dres[i] = devresid(dist, yi, μi, wti)
        wrkwt[i] = wti * abs2(dμdηi) / vari
    end
    m
end

function updateweights{T <: AbstractFloat}(dres::Vector{T}, wrkwt::Vector{T}, wt::Vector{T})
    @inbounds @simd for i in eachindex(wt)
        wti = wt[i]
        dres[i] *= wti
        wrkwt[i] *= wti
    end
end

function updateμ!{T<:AbstractFloat, D<:Union{Bernoulli, Binomial}, L<:LogitLink}(m::GeneralizedLinearMixedModel{T,D,L})
    y, η, μ, wrkres, wrkwt, dres = m.y, m.η, m.μ, m.wrkresid, m.wrkwt, m.devresid

    @inbounds Threads.@threads for i in eachindex(η)
        ηi = clamp(η[i], -20.0, 20.0)
        ei = exp(-ηi)
        opei = 1 + ei
        μi = μ[i] = inv(opei)
        dμdη = wrkwt[i] = ei / abs2(opei)
        yi = y[i]
        wrkres[i] = (yi - μi) / dμdη
        dres[i] = -2 * (yi == 1 ? log(μi) : yi == 0 ? log1p(-μi) :
            (yi * (log(μi) - log(yi)) + (1 - yi) * (log1p(-μi) - log1p(-yi))))
    end

    if !isempty(m.wt)
        updateweights(dres, wrkwt, m.wt)
    end
end

function updateμ!{T<:AbstractFloat, D<:Poisson, L<:LogLink}(m::GeneralizedLinearMixedModel{T,D,L})
    y, η, μ, wrkres, wrkwt, dres = m.y, m.η, m.μ, m.wrkresid, m.wrkwt, m.devresid

    @inbounds Threads.@threads for i in eachindex(η)
        ηi = η[i]
        μi = μ[i] = exp(ηi)
        dμdη = wrkwt[i] = μi
        yi = y[i]
        wrkres[i] = (yi - μi) / dμdη
        dres[i] = 2 * (StatsFuns.xlogy(yi, yi / μi) - (yi - μi))
    end

    if !isempty(m.wt)
        updateweights(dres, wrkwt, m.wt)
    end
end

function wrkresp!{T <: AbstractFloat}(v::DenseVecOrMat{T}, m::GeneralizedLinearMixedModel{T})
    if isempty(m.offset)
        return broadcast!(+, v, m.η, m.wrkresid)
    end
    η, offset, wrkresid = m.η, m.offset, m.wrkresid
    @inbounds @simd for i in eachindex(η)
        v[i] = wrkresid[i] + η[i] - offset[i]
    end
    v
end


text:
#excitation transfers - we default to regular transfers when no better candidate is found
exci_transfer_left(v,A,B=A) = transfer_left(v,A,B)
exci_transfer_right(v,A,B=A) = transfer_right(v,A,B)
exci_transfer_left(v,A,B,C) = transfer_left(v,A,B,C)
exci_transfer_right(v,A,B,C) = transfer_right(v,A,B,C)
exci_transfer_left(v,A,B,C,D) = transfer_left(v,A,B,C,D)
exci_transfer_right(v,A,B,C,D) = transfer_right(v,A,B,C,D)

#transfer, but the upper A is an excited tensor
exci_transfer_left(v::MPSBondTensor, A::MPOTensor, Ab::MPSTensor) = @tensor t[-1 -2;-3] := v[1,2]*A[2,3,-2,-3]*conj(Ab[1,3,-1])
exci_transfer_right(v::MPSBondTensor, A::MPOTensor, Ab::MPSTensor) = @tensor t[-1 -2;-3] := A[-1,3,-2,1]*v[1,2]*conj(Ab[-3,3,2])

#transfer, but there is both a utility leg and an mpo leg that is passed through
exci_transfer_left(v::MPOTensor, A::MPSTensor, Ab::MPSTensor=A) = @tensor v[-1 -2;-3 -4] := v[1,-2,-3,2]*A[2,3,-4]*conj(Ab[1,3,-1])
exci_transfer_right(v::MPOTensor, A::MPSTensor, Ab::MPSTensor=A) = @tensor v[-1 -2;-3 -4] := A[-1,3,1]*v[1,-2,-3,2]*conj(Ab[-4,3,2])

#transfer, but the upper A is an excited tensor and there is an mpo leg being passed through
exci_transfer_left(v::MPSTensor, A::MPOTensor, Ab::MPSTensor) = @tensor t[-1 -2;-3 -4] := v[1,-2,2]*A[2,3,-3,-4]*conj(Ab[1,3,-1])
exci_transfer_right(v::MPSTensor, A::MPOTensor, Ab::MPSTensor) = @tensor t[-1 -2;-3 -4] := A[-1,3,-2,1]*v[1,-3,2]*conj(Ab[-4,3,2])

#mpo transfer, but with A an excitation-tensor
exci_transfer_left(v::MPSTensor,O::MPOTensor,A::MPOTensor,Ab::MPSTensor) = @tensor t[-1 -2;-3 -4] := v[4,5,1]*A[1,3,-3,-4]*O[5,2,-2,3]*conj(Ab[4,2,-1])
exci_transfer_right(v::MPSTensor,O::MPOTensor,A::MPOTensor,Ab::MPSTensor) = @tensor t[-1 -2;-3 -4] := A[-1,1,-2,5]*O[-3,3,4,1]*conj(Ab[-4,3,2])*v[5,4,2]

#mpo transfer, with an excitation leg
exci_transfer_left(v::MPOTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor=A) = @tensor v[-1 -2;-3 -4] := v[4,5,-3,1]*A[1,3,-4]*O[5,2,-2,3]*conj(Ab[4,2,-1])
exci_transfer_right(v::MPOTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor=A) = @tensor v[-1 -2;-3 -4] := A[-1,1,5]*O[-3,3,4,1]*conj(Ab[-4,3,2])*v[5,-2,4,2]

#A is an excitation tensor; with an excitation leg
exci_transfer_left(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::M,Ab::V=A) where V<:MPSTensor where M <:MPOTensor = exci_transfer_left(M,vec,ham,pos,A,Ab)
exci_transfer_right(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::M,Ab::V=A) where V<:MPSTensor where M <:MPOTensor = exci_transfer_right(M,vec,ham,pos,A,Ab)

#v has an extra excitation leg
exci_transfer_left(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::M,Ab::M=A) where V<:MPOTensor where M <:MPSTensor = exci_transfer_left(V,vec,ham,pos,A,Ab)
exci_transfer_right(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::M,Ab::M=A) where V<:MPOTensor where M <:MPSTensor = exci_transfer_right(V,vec,ham,pos,A,Ab)

function exci_transfer_left(RetType,vec,ham::MPOHamiltonian,pos,A,Ab=A)
    toreturn = Array{RetType,1}(undef,length(vec));
    assigned = [false for i in 1:ham.odim]

    for (j,k) in keys(ham,pos)
        if assigned[k]
            if isscal(ham,pos,j,k)
                toreturn[k]+=ham.Os[pos,j,k]*exci_transfer_left(vec[j],A,Ab)
            else
                toreturn[k]+=exci_transfer_left(vec[j],ham[pos,j,k],A,Ab)
            end
        else
            if isscal(ham,pos,j,k)
                toreturn[k]=ham.Os[pos,j,k]*exci_transfer_left(vec[j],A,Ab)
            else
                toreturn[k]=exci_transfer_left(vec[j],ham[pos,j,k],A,Ab)
            end
            assigned[k]=true
        end
    end

    for k in 1:ham.odim
        if !assigned[k]
            #prefereably this never happens, because it's a wasted step
            #it's also avoideable with a little bit more code
            toreturn[k]=exci_transfer_left(vec[1],ham[pos,1,k],A,Ab)
        end
    end

    return toreturn
end
function exci_transfer_right(RetType,vec,ham::MPOHamiltonian,pos,A,Ab=A)
    toreturn = Array{RetType,1}(undef,length(vec));
    assigned = [false for i in 1:ham.odim]

    for (j,k) in keys(ham,pos)
        if assigned[j]
            if isscal(ham,pos,j,k)
                toreturn[j]+=ham.Os[pos,j,k]*exci_transfer_right(vec[k],A,Ab)
            else
                toreturn[j]+=exci_transfer_right(vec[k],ham[pos,j,k],A,Ab)
            end

        else
            if isscal(ham,pos,j,k)
                toreturn[j]=ham.Os[pos,j,k]*exci_transfer_right(vec[k],A,Ab)
            else
                toreturn[j]=exci_transfer_right(vec[k],ham[pos,j,k],A,Ab)
            end
            assigned[j]=true
        end
    end

    for j in 1:ham.odim
        if !assigned[j]
            toreturn[j]=exci_transfer_right(vec[1],ham[pos,j,1],A,Ab)
        end
    end

    return toreturn
end

#=
original code :
(lBsEr,convhist)=linsolve(RecursiveVec(lBs...),RecursiveVec(lBs...),GMRES()) do y
    x=collect(y.vecs)
    tor=reduce((a,b)->transfer_left(a,ham,b,mpsright.AR[b],mpsleft.AL[b])*exp(-1im*p),1:len,init=x)

    if(trivial)
        for i in ids
            @tensor tor[i][-1,-2,-3,-4]-=tor[i][1,-2,2,-4]*r_RL(mpsleft)[2,1]*l_RL(mpsleft)[-1,-3]
        end
    end

    tor=x-tor
    return  RecursiveVec(tor...)
end
lBsE=lBsEr.vecs

but this can be made faster; using the fact that the hamiltonion is upper-triangular, which is what we do here
=#

function left_excitation_transfer_system(lBs,ham,mpsleft::InfiniteMPS,mpsright::InfiniteMPS,trivial,p)
    len = ham.period
    found=zero.(lBs)
    ids = collect(Iterators.filter(x->isid(ham,x),1:ham.odim));

    for i in 1:ham.odim

        #this operation can be sped up by at least a factor 2;  found mostly consists of zeros
        start = found
        for k in 1:len
            start = exci_transfer_left(start,ham,k,mpsright.AR[k],mpsleft.AL[k])*exp(conj(1im*p))

            if trivial
                for l in ids[2:end-1]
                    @tensor start[l][-1,-2,-3,-4]-=start[l][1,-2,-3,2]*r_RL(mpsright,k)[2,1]*l_RL(mpsright,k)[-1,-4]
                end
            end
        end

        #either the element ham_ii exists; in which case we have to solve a linear system
        #otherwise it's easy and we already know found[i]
        if reduce((a,b)->contains(ham,b,i,i),1:len,init=true)
            (found[i],convhist)=linsolve(lBs[i]+start[i],lBs[i]+start[i],GMRES()) do y
                x=reduce((a,b)->exci_transfer_left(a,ham[b,i,i],mpsright.AR[b],mpsleft.AL[b])*exp(conj(1im*p)),1:len,init=y)

                if trivial && i in ids
                    @tensor x[-1,-2,-3,-4]-=x[1,-2,-3,2]*r_RL(mpsleft)[2,1]*l_RL(mpsleft)[-1,-4]
                end

                return y-x
            end
            if convhist.converged<1
                @info "left $(i) excitation inversion failed normres $(convhist.normres)"
            end
        else
            found[i]=lBs[i]+start[i]
        end
    end
    return found
end

function right_excitation_transfer_system(rBs,ham,mpsleft,mpsright::InfiniteMPS,trivial,p)
    len = ham.period
    found=zero.(rBs)
    ids = collect(Iterators.filter(x->isid(ham,x),1:ham.odim));
    for i in ham.odim:-1:1

        start = found
        for k in len:-1:1
            start = exci_transfer_right(start,ham,k,mpsleft.AL[k],mpsright.AR[k])*exp(1im*p)

            if trivial
                for l in ids[2:end-1]
                    @tensor start[l][-1,-2,-3,-4]-=start[l][1,-2,-3,2]*l_LR(mpsright,k)[2,1]*r_LR(mpsright,k)[-1,-4]
                end
            end
        end

        if reduce((a,b)->contains(ham,b,i,i),1:len,init=true)
            (found[i],convhist)=linsolve(rBs[i]+start[i],rBs[i]+start[i],GMRES()) do y
                x=reduce((a,b)->exci_transfer_right(a,ham[b,i,i],mpsleft.AL[b],mpsright.AR[b])*exp(1im*p),len:-1:1,init=y)

                if trivial && i in ids
                    @tensor x[-1,-2,-3,-4]-=x[1,-2,-3,2]*l_LR(mpsright)[2,1]*r_LR(mpsright)[-1,-4]
                end

                return y-x
            end
            if convhist.converged<1
                @info "right $(i) excitation inversion failed normres $(convhist.normres)"
            end
        else
            found[i]=rBs[i]+start[i]
        end
    end
    return found
end


text:
"""
    function triggered(hits::Vector{T}) where {T<:TriggeredHit}

Return a `Vector` of triggered hits.
"""
triggered(hits::Vector{T}) where {T<:Union{Hit, CalibratedHit}} = filter(h->h.triggered, hits)

"""
    function nfoldhits(hits::Vector{T}, Δt, n) where {T<:AbstractDAQHit}

Create a `Vector` with hits contributing to `n`-fold coincidences within a time
window of Δt.
"""
function nfoldhits(hits::Vector{T}, Δt, n) where {T<:AbstractDAQHit}
    hit_map = DefaultDict{Integer}{Vector{T}}(() -> T[])
    for hit ∈ sort(hits)
        push!(hit_map[hit.dom_id], hit)
    end
    chits = Vector{T}()
    for (dom_id, dom_hits) ∈ hit_map
        bag = Vector{T}()
        push!(bag, dom_hits[1])
        t0 = dom_hits[1].t
        for hit in dom_hits[2:end]
            if hit.t - t0 > Δt
                if length(bag) >= n
                    append!(chits, bag)
                end
                bag = Vector{T}()
            end
            push!(bag, hit)
            t0 = hit.t
        end
    end
    return chits
end

"""
    function count_multiplicities(hits::Vector{T}, tmax=20) where {T<:AbstractHit}

Calculate the multiplicities for a given time window. Two arrays are
are returned, one contains the multiplicities, the second one the IDs
of the coincidence groups.
The hits should be sorted by time and then by dom_id.
"""
function count_multiplicities(hits::Vector{T}, tmax=20) where {T<:AbstractHit}
    n = length(hits)
    mtp = ones(Int32, n)
    cid = zeros(Int32, n)
    idx0 = 1
    _mtp = 1
    _cid = idx0
    t0 = hits[idx0].t
    dom_id = hits[idx0].dom_id
    for i in 2:n
        hit = hits[i]
        if hit.dom_id != dom_id
            mtp[idx0:i-1] .= _mtp
            cid[idx0:i-1] .= _cid
            dom_id = hit.dom_id
            t0 = hit.t
            _mtp = 1
            _cid += 1
            idx0 = i
            continue
        end
        Δt = hit.t - t0
        if Δt > tmax
            mtp[idx0:i] .= _mtp
            cid[idx0:i] .= _cid
            _mtp = 0
            _cid += 1
            idx0 = i
            t0 = hit.t
        end
        _mtp += 1
    end
    mtp[idx0:end] .= _mtp
    cid[idx0:end] .= _cid
    mtp, cid
end

"""
    function count_multiplicities!(hits::Vector{CalibratedHit}, tmax=20)

Counts the multiplicities and modifies the .multiplicity field of the hits.
Important: the hits have to be sorted by time and then by DOM ID first.
"""
function count_multiplicities!(hits::Vector{CalibratedHit}, tmax=20)
    _mtp = 0
    _cid = 0
    t0 = 0
    dom_id = 0
    hit_buffer = Vector{CalibratedHit}()

    function process_buffer()
        while !isempty(hit_buffer)
            _hit = pop!(hit_buffer)
            _hit.multiplicity.count = _mtp
            _hit.multiplicity.id = _cid
        end
    end

    function reset()
        _mtp = 1
        _cid += 1
    end

    for hit ∈ hits
        if length(hit_buffer) == 0
            reset()
            push!(hit_buffer, hit)
            t0 = hit.t
            dom_id = hit.dom_id
            continue
        end
        if hit.dom_id != dom_id
            process_buffer()
            push!(hit_buffer, hit)
            t0 = hit.t
            dom_id = hit.dom_id
            reset()
            continue
        end
        Δt = hit.t - t0
        if Δt > tmax
            process_buffer()
            push!(hit_buffer, hit)
            t0 = hit.t
            reset()
        else
            push!(hit_buffer, hit)
            _mtp += 1
        end
    end
    if length(hit_buffer) > 0
        process_buffer()
    end
    return
end

"""
    function domhits(hits::Vector{T}) where {T<:AbstractDAQHit}

Sort hits by DOM ID and put them into a dictionary.
"""
function domhits(hits::Vector{T}) where {T<:AbstractDAQHit}
    hit_map = DefaultDict{Integer}{Vector{T}}(() -> T[])
    for hit ∈ hits
        push!(hit_map[hit.dom_id], hit)
    end
    hit_map
end

"""
    function duhits(hits::Vector{T}) where {T<:CalibratedHit}

Sort hits by DU and put them into a dictionary.
"""
function duhits(hits::Vector{T}) where {T<:CalibratedHit}
    hit_map = DefaultDict{Integer}{Vector{T}}(() -> T[])
    for hit ∈ hits
        push!(hit_map[hit.du], hit)
    end
    hit_map
end

"""
    function totcut(hits::Vector{T}, tot) where {T<:AbstractDAQHit}

Return a vector of hits with ToT >= `tot`.
"""
function totcut(hits::Vector{T}, tot) where {T<:AbstractDAQHit}
    return filter(h->h.tot >= tot, hits)
end

"""
    function nphes(tot)

Returns the estimated number of photoelectrons for a given ToT.
"""
function nphes(tot)
    if tot <= 20
        return 0.0
    end
    if tot <= 26
        return 1.0
    end
    if tot < 170
        return 1.0 + (tot - 26)/(1/0.28)
    end
    return 40.0 + (255 - tot)*2.0
end

"""
    function slew(tot)

Return the time slewing for a ToT.
"""
function slew(tot)
    p₀ = 7.70824
    p₁ = 0.00879447
    p₂ = -0.0621101
    p₃ = -1.90226

    p₀ * ℯ^(p₁*√tot + p₂*tot) + p₃
end


text:
# first sorts x, then forces the entries in x such that f(x) has alternating sign.
function forcealternating(x_in::Vector{T}, f::Function) where T
    @assert length(x_in) > 1

    x = sort(x_in)

    i = 2
    while i <= length(x)

        fx_im1 = f(x[i-1])
        u::Bool = sign(fx_im1) > 0

        fx_i = f(x[i])
        v = sign(fx_i) > 0

        if !xor(u,v)
            # same sign for error evaluations of x[i-1] and x[i].
            # Remove the position with smaller absolute error.
            if abs(fx_i) > abs(fx_im1)
                deleteat!(x,i-1)

            else
                deleteat!(x,i)
            end
        else
            i += 1
        end
    end

    # the above loop sometimes miss the end points of the position.
    # check if endpoints have different error signs. If so, add them back.
    if xor( sign(f(x_in[1])) > 0, sign(f(x[1])) > 0 )
        insert!(x, 1, x_in[1])
    end

    if xor( sign(f(x_in[end])) > 0, sign(f(x[end])) > 0 )
        push!(x, x_in[end])
    end

    return x
end

# try ForwardDiff.derivative(f,x) someday.
function numericalchkextrema(f::Function, x::T, tol)::Bool where T
    return numericalchkextrema(f,x,convert(T,tol))
end

# TODO try passing df instead of f.
function numericalchkextrema(f::Function, x::T, tol::T)::Bool where T

    if isnumericallyclose(FiniteDiff.finite_difference_derivative(f,x), zero(T), tol)
        #δ = 1e-3
        df_δ = FiniteDiff.finite_difference_derivative(f,x+tol)
        df_mδ = FiniteDiff.finite_difference_derivative(f,x-tol)
        #return true
        return xor(df_δ > zero(T),df_mδ > zero(T))
    end

    return false
end

function returnNclosestzeroderivatives(a::Vector{T}, 𝑒::Function, N::Int, tol::T) where T
    @assert length(a) >= N

    abs_𝑑out_𝑑x = collect( abs(FiniteDiff.finite_difference_derivative(𝑒, a[i])) for i = 1:length(a) )

    sort_ind = sortperm(abs_𝑑out_𝑑x)

    out = a[sort_ind[1:N]]

    if out[end] > tol

        # if the smallest zero is larger than the input tolerance, then just return the original array unmodified.
        return a, false
    end

    return out, true
end

# array version of samplepoints(f::Function,...). F is f(k,L), k = 0,1,...,L
function samplepoints(F::Vector{T}, a::T, b::T)::Vector{T} where T
    @assert a < b

    two = one(T) + one(T)
    Δ = b-a
    return collect( (F[i]+one(T))/two*Δ+a for i = 1:length(F) )
end

# f generates N points between the interval [-1,1].
function samplepoints(f::Function, a::T, b::T, N::Int)::Vector{T} where T
    @assert a < b
    L = N-1

    two = one(T) + one(T)
    Δ = b-a
    return collect( (f(k,L)+one(T))/two*Δ+a for k = 0:L )
end

function packagefiltersolutionfast( 𝓧_new,
                                    h_history,
                                    wfunc,
                                    zfunc,
                                    dfunc,
                                    filter_type_num::Int)

    c,ν, discrepancy, xq, yq_bary, yq_cheby = BarycentrictoChebyshevfast(𝓧_new, wfunc, dfunc, filter_type_num)
    println("conversion discrepany between barycenter and Chebyshev is ", discrepancy)

    ######
    L = length(c)

    L_old = L

    # get impulse response.
    𝑏 = c
    if filter_type_num == 2
        𝑏 = converttotype2(c)
    elseif filter_type_num == 3
        𝑏 = converttotype3(c)
    elseif filter_type_num ==4
        𝑏 = converttotype4(c)
    end

    q = 1
    if filter_type_num == 1 || filter_type_num == 2
        q = 0
    end

    gethfunc = gethtype4
    if filter_type_num == 1
        L = length(c)-1

        gethfunc = gethtype1
    elseif filter_type_num == 2
        gethfunc = gethtype2
    elseif filter_type_num ==3
        gethfunc = gethtype3
    end

    h = gethfunc(𝑏)

    return h
end

function packagefiltersolution(𝓧_new,
                                h_history,
                                wfunc,
                                zfunc,
                                dfunc,
                                filter_type_num::Int,
                                fig_num::Int,
                                savefig_flag::Bool,
                                plot_flag::Bool)

    @time c,ν, discrepancy, xq, yq_bary, yq_cheby = BarycentrictoChebyshev(𝓧_new, wfunc, dfunc, filter_type_num)
    println("conversion discrepany between barycenter and Chebyshev is ", discrepancy)

    ######
    L = length(c)

    L_old = L

    # get impulse response.
    𝑏 = c
    if filter_type_num == 2
        𝑏 = converttotype2(c)
    elseif filter_type_num == 3
        𝑏 = converttotype3(c)
    elseif filter_type_num ==4
        𝑏 = converttotype4(c)
    end

    q = 1
    if filter_type_num == 1 || filter_type_num == 2
        q = 0
    end

    gethfunc = gethtype4
    if filter_type_num == 1
        L = length(c)-1

        gethfunc = gethtype1
    elseif filter_type_num == 2
        gethfunc = gethtype2
    elseif filter_type_num ==3
        gethfunc = gethtype3
    end

    h = gethfunc(𝑏)

    if !plot_flag
        return h, fig_num
    end

    ######### plot to interpret results.

    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(xq, yq_bary, "--", label = "yq_bary")
    # PyPlot.plot(xq, yq_cheby, label = "yq_cheby")
    # PyPlot.title("Barycenter vs. Chebyshev")
    # PyPlot.legend()
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_Barycenter_to_Chebyshev.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    ## plot level error history.

    # PyPlot.figure(fig_num)
    # fig_num += 1
    #
    # PyPlot.plot(collect(1:length(h_history)), abs.(h_history), ".-")
    # PyPlot.title("level error history")
    #
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_level_error_history.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    ## plot resultant filter's response.
    f(xx) = wfunc(xx)*dfunc(xx)
    L = length(𝓧_new)-2
    xq = convert(Vector{BigFloat},collect(LinRange(-1,1,5*L)))

    f_itp_xq, f_itp_𝓧_new, 𝑒_xq, 𝑒_𝓧_new = visualizefiltersolution(𝓧_new,xq,wfunc,dfunc,f)

    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(xq, f.(xq), "--", label = "f = w*d")
    # PyPlot.plot(xq, f_itp_xq, label = "solution")
    # PyPlot.plot(xq, dfunc.(xq), label = "d")
    # PyPlot.plot(𝓧_new, f_itp_𝓧_new,".", label = "reference positions")
    # PyPlot.title("solution filter's response")
    # PyPlot.legend()
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_solution.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(xq, 𝑒_xq, label = "error")
    # PyPlot.plot(𝓧_new, 𝑒_𝓧_new,".", label = "reference positions")
    # PyPlot.title("error")
    # PyPlot.legend()
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_solution_error.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    ###### analyze

    # prepare.
    N = length(h)

    a = 0.0 #-π
    b = 2*π # π
    ω_set = convert(Vector{BigFloat}, collect( LinRange(a,b,N*20) ) )
    ω_set_fft = collect( LinRange(a,b-2*π/N,N))

    # this is the verified one.
    y = Vector{Complex{BigFloat}}(undef,0)
    y2 = Vector{BigFloat}(undef,0)

    if filter_type_num == 1
        y = collect( evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
        y2 = collect( sum( c[k+1]*cos(k*ω_set[i]) for k = 0:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 2
        y = collect( cos(ω_set[i]/2)*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
        y2 = collect( sum( 𝑏[k]*cos((k-0.5)*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 3
       y = collect( sin(ω_set[i])*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
       y2 = collect( sum( 𝑏[k]*sin(k*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )

    elseif filter_type_num == 4
       y = collect( sin(ω_set[i]/2)*evalunivariateChebyshevpolynomial(c,cos(ω_set[i]))*exp( im*(q*π - (N-1)*ω_set[i])/2 ) for i = 1:length(ω_set) )
       y2 = collect( sum( 𝑏[k]*sin((k-0.5)*ω_set[i]) for k = 1:L ) for i = 1:length(ω_set) )
    end

    y = abs.(y)
    y2 = abs.(y2)
    discrepancy = LinearAlgebra.norm(y-y2)
    println("conversion discrepany between y and y2 is ", discrepancy)

    ##### Visualize impulse response.
    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(h, label = "h")
    # PyPlot.plot(h, ".")
    # title_string = Printf.@sprintf("check if the recovered h has the characteristics of type %d", filter_type_num)
    # PyPlot.title(title_string)
    # PyPlot.legend()
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_solution_impulse_rsp.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    ##### prepare fft and DFT from impulse response.

    impulse_rsp = fft(convert(Vector{Float64},h))
    mag_rsp_fft = abs.(impulse_rsp)
    phase_rsp_fft = angle.(impulse_rsp)

    impulse_rsp_G = collect( naiveimpulsersp(h,ω_set[i]) for i = 1:length(ω_set) )
    mag_rsp_G = abs.(impulse_rsp_G)
    phase_rsp_G = angle.(impulse_rsp_G)

    println("discrepancy between y and G is ", LinearAlgebra.norm(y-mag_rsp_G))

    ##### visualize.

    ## plot the polynomial specified by b, see if it is in agreement with
    #   the fft/DFT response.
    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(ω_set_fft, mag_rsp_fft, label = "fft: solution discrete")
    # PyPlot.plot(ω_set_fft, mag_rsp_fft, ".")
    #
    # PyPlot.plot(ω_set, mag_rsp_G, label = "G: solution continuous")
    #
    # y_d = abs.(wfunc.(cos.(ω_set)).*dfunc.(cos.(ω_set)))
    # PyPlot.plot(ω_set, y_d, label = "y_d: target")
    #
    # PyPlot.title("Magnitude response")
    # PyPlot.legend()
    #
    # if savefig_flag
    #     save_name = Printf.@sprintf("./outputs/type_%s_solution_magnitude_rsp.png",
    #                                     filter_type_num)
    #     PyPlot.savefig(save_name)
    #     sleep(save_delay)
    #     PyPlot.close(fig_num)
    # end

    return h, fig_num
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
# magma_cgesvd wrapper
function magmaCgesvd(jobu, jobvt, m, n, A, lda, s, U,
	ldu, VT, ldvt, work, lwork, rwork, info)
    ccall((:magma_cgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{ComplexF32},
			    Cint, PtrOrCuPtr{Float32}, PtrOrCuPtr{ComplexF32},
				Cint, PtrOrCuPtr{ComplexF32}, Cint,
				PtrOrCuPtr{ComplexF32}, Cint, PtrOrCuPtr{Float32},
				PtrOrCuPtr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, rwork, info)
end


text:
"Do one iteration and get the current value of a Xorshift RNG object."
xorshift_next

# We'll define a new type to support UInt32 type, so no need for rand(T, seed) now.

let doc_xorshift64(r) = """
```julia
$r <: AbstractXorshift64
$r([seed])
```

$r RNG. The `seed` will be automatically convert to an `UInt64` number.
"""
    @doc doc_xorshift64("Xorshift64") Xorshift64
    @doc doc_xorshift64("Xorshift64Star") Xorshift64Star
end

let doc_xorshift128(r) = """
```julia
$r <: AbstractXorshift128
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of two `UInt64`s, or an `Integer` which will be automatically convert to
an `UInt128` number.
"""
    @doc doc_xorshift128("Xorshift128") Xorshift128
    @doc doc_xorshift128("Xorshift128Star") Xorshift128Star
    @doc doc_xorshift128("Xorshift128Plus") Xorshift128Plus
end

let doc_xorshift1024(r) = """
```julia
$r <: AbstractXorshift1024
$r([seed...])
```

$r RNG. The `seed` can be a `Tuple` of 16 `UInt64`s, or several (no more than 16) `Integer`s which will all
be automatically converted to `UInt64` numbers.
"""
    @doc doc_xorshift1024("Xorshift1024") Xorshift1024
    @doc doc_xorshift1024("Xorshift1024Star") Xorshift1024Star
    @doc doc_xorshift1024("Xorshift1024Plus") Xorshift1024Plus
end

let doc_xoroshiro128(r) = """
```julia
$r <: AbstractXoroshiro128
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of two `UInt64`s, or an `Integer` which will be automatically convert to
an `UInt128` number.
"""
    @doc doc_xoroshiro128("Xoroshiro128") Xoroshiro128
    @doc doc_xoroshiro128("Xoroshiro128Star") Xoroshiro128Star
    @doc doc_xoroshiro128("Xoroshiro128Plus") Xoroshiro128Plus
end


text:
"""
    `getz(tply::AbstractArray{<:Real,1}, nply::AbstractArray{<:Integer,1})`

    `getz(lam::Laminate)`

Returns a laminate's z-coordinates (coordinates of top and bottom of laminas)
given the thickness of plies in each lamina and the number of plies in each
lamina.
"""
function getz(tply::AbstractArray{<:Real,1}, nply::AbstractArray{<:Integer,1})
    ttotal = 0.0
    tlam = zeros(eltype(tply), length(tply)+1)
    for i = 1:length(tply)
        ttotal += tply[i]*nply[i]
        tlam[i+1] = ttotal
    end
    tlam .-= ttotal/2.0
    return tlam
end
getz(lam::Laminate) = getz(lam.tply, lam.nply)

"""
    `getQ(e1::Real, e2::Real, g12::Real, nu12::Real, theta::Real=0.0)`

    `getQ(mat::Material, theta::Real)`

    `getQ(mat::Material, lam::Laminate)`

Returns Q matrix rotated `theta` degrees
"""
function getQ(e1::Real, e2::Real, g12::Real, nu12::Real, theta::Real=0.0)
    nud = 1.0-e2/e1*nu12^2.0
    q11 = e1/nud
    q22 = e2/nud
    q12 = nu12*e2/nud
    q66 = g12
    return rotQ(q11, q12, q22, q66, theta)
end
getQ(mat::Material, theta::Real=0.0) = getQ(mat.e1, mat.e2, mat.g12, mat.nu12, theta)
getQ(mat::Array{Material,1}, lam::Laminate) = getQ.(mat, lam.theta)

"""
    `rotQ(q::AbstractArray{<:Real,2}, theta::Real)`

    `rotQ(qxx::Real, q12::Real, q22::Real, q66::Real, theta::Real)`

Rotates Q matrix by theta degrees
"""
function rotQ!(q::AbstractArray{<:Real,2}, theta::Real)

    if theta != zero(eltype(theta))
        c = cosd(theta)
        s = sind(theta)
        c2 = c^2
        s2 = s^2
        cs = c*s

        # qbar = inv(tsigma)*q*teps

        # tsigma = [c^2.0 s^2.0 2*c*s; s^2.0 c^2.0 -2*c*s; -c*s c*s c^2-s^2]

        inv_tsig11 = c2
        inv_tsig12 = s2
        inv_tsig13 = -2*cs
        inv_tsig21 = s^2
        inv_tsig22 = c^2
        inv_tsig23 = 2*cs
        inv_tsig31 = cs
        inv_tsig32 = -cs
        inv_tsig33 = 2*c2-1

        # teps = LinearAlgebra.Diagonal([1.0,1.0,2.0])*tsigma*LinearAlgebra.Diagonal([1.0,1.0,0.5])

        teps11 = c2
        teps12 = s2
        teps13 = cs
        teps21 = s2
        teps22 = c2
        teps23 = -cs
        teps31 = -2*cs
        teps32 = 2*cs
        teps33 = c2-s2

        # inv(tsigma)*q

        inv_tsig_q11 = inv_tsig11*q[1,1] + inv_tsig12*q[2,1] + inv_tsig13*q[3,1]
        inv_tsig_q12 = inv_tsig11*q[1,2] + inv_tsig12*q[2,2] + inv_tsig13*q[3,2]
        inv_tsig_q13 = inv_tsig11*q[1,3] + inv_tsig12*q[2,3] + inv_tsig13*q[3,3]
        inv_tsig_q21 = inv_tsig21*q[1,1] + inv_tsig22*q[2,1] + inv_tsig23*q[3,1]
        inv_tsig_q22 = inv_tsig21*q[1,2] + inv_tsig22*q[2,2] + inv_tsig23*q[3,2]
        inv_tsig_q23 = inv_tsig21*q[1,3] + inv_tsig22*q[2,3] + inv_tsig23*q[3,3]
        inv_tsig_q31 = inv_tsig31*q[1,1] + inv_tsig32*q[2,1] + inv_tsig33*q[3,1]
        inv_tsig_q32 = inv_tsig31*q[1,2] + inv_tsig32*q[2,2] + inv_tsig33*q[3,2]
        inv_tsig_q33 = inv_tsig31*q[1,3] + inv_tsig32*q[2,3] + inv_tsig33*q[3,3]

        # inv(tsigma)*q*teps

        q[1,1] = inv_tsig_q11*teps11 + inv_tsig_q12*teps21 + inv_tsig_q13*teps31
        q[1,2] = inv_tsig_q11*teps12 + inv_tsig_q12*teps22 + inv_tsig_q13*teps32
        q[1,3] = inv_tsig_q11*teps13 + inv_tsig_q12*teps23 + inv_tsig_q13*teps33
        q[2,1] = inv_tsig_q21*teps11 + inv_tsig_q22*teps21 + inv_tsig_q23*teps31
        q[2,2] = inv_tsig_q21*teps12 + inv_tsig_q22*teps22 + inv_tsig_q23*teps32
        q[2,3] = inv_tsig_q21*teps13 + inv_tsig_q22*teps23 + inv_tsig_q23*teps33
        q[3,1] = inv_tsig_q31*teps11 + inv_tsig_q32*teps21 + inv_tsig_q33*teps31
        q[3,2] = inv_tsig_q31*teps12 + inv_tsig_q32*teps22 + inv_tsig_q33*teps32
        q[3,3] = inv_tsig_q31*teps13 + inv_tsig_q32*teps23 + inv_tsig_q33*teps33

    end

    return q
end

function rotQ(q::AbstractArray{<:Real,2}, theta::Real)
    qbar = copy(q)
    return rotQ!(qbar, theta)
end

function rotQ(q11::Real, q12::Real, q22::Real, q66::Real, theta::Real)
  q0 = zero(eltype(q11))
  q = [q11 q12 q0; q12 q22 q0; q0 q0 q66]
  return rotQ!(q, theta)
end

"""
    `getABD(matid::AbstractArray{<:Integer,1},
        nply::AbstractArray{<:Integer,1}, tply::AbstractArray{<:Real,1},
        theta::AbstractArray{<:Real,1},
        q::AbstractArray{<:AbstractArray{<:Real,2},1})`

    `getABD(lam::Laminate, q::AbstractArray{<:AbstractArray{<:Real,2},1})`

Returns A, B, and D matrices

# Arguments:
* `matid::AbstractArray{<:Integer,1}`: material id of each lamina
* `nply::AbstractArray{<:Integer,1}`: number of plies in each lamina
* `tply::AbstractArray{<:Real,1}`: thickness of a ply (m) in each lamina
* `theta::AbstractArray{<:Real,1}`: orientation (deg) of each lamina
* `q::AbstractArray{<:AbstractArray{<:Real,2}}`: Stiffness matrix of each lamina
"""
function getABD!(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2}, D::AbstractArray{<:Real,2},
    matid::AbstractArray{<:Integer,1}, nply::AbstractArray{<:Integer,1}, tply::AbstractArray{<:Real,1},
    theta::AbstractArray{<:Real,1}, q::AbstractArray{<:AbstractArray{<:Real,2},1})

    typeof(A)

    R = eltype(tply)

    nlam = length(nply)

    z = getz(tply, nply)

    qbar = zeros(R, 3, 3)

    # Loop through layers filling in ABD matrix
    for k = 1:nlam
        # Rotate material stiffness properties by specifed angle
        imat = matid[k]
        qbar .= q[imat]
        rotQ!(qbar, theta[k])

        # Find lumped stiffness properties for the laminate
        for i = 1:3
            for j = 1:3
                A[i,j] = A[i,j]+qbar[i,j]*(z[k+1]-z[k])
                B[i,j] = B[i,j]+qbar[i,j]/2.0*(z[k+1]^2.0-z[k]^2.0)
                D[i,j] = D[i,j]+qbar[i,j]/3.0*(z[k+1]^3.0-z[k]^3.0)
            end
        end
    end

    return A, B, D
end

function getABD(matid::AbstractArray{<:Integer,1}, nply::AbstractArray{<:Integer,1}, tply::AbstractArray{<:Real,1},
    theta::AbstractArray{<:Real,1}, q::AbstractArray{<:AbstractArray{<:Real,2},1})

    R = eltype(tply)

    A = zeros(R, 3, 3)
    B = zeros(R, 3, 3)
    D = zeros(R, 3, 3)

    return getABD!(A, B, D, matid, nply, tply, theta, q)
end

getABD(lam::Laminate, q::AbstractArray{<:AbstractArray{<:Real,2},1}) =
    getABD(lam.matid, lam.nply, lam.tply, lam.theta, q)

"""
    compliance(A::Array{<:Real,2}, B::Array{<:Real,2}, D::Array{<:Real,2})
Returns alpha, beta, and delta
"""
function compliance(A::AbstractArray{<:Real,2}, B::AbstractArray{<:Real,2},
    D::AbstractArray{<:Real,2})
    invA = inv(A)
    delta = inv(D-B*invA*B)
    beta = -invA*B*delta
    alfa = invA-beta*B*invA
    return alfa, beta, delta
end # compliance


text:
abstract type AbstractTargetingBehavior <: Any end

"Vessels target fishing behavior at random within the fishery domain."
struct RandomTargeting <: AbstractTargetingBehavior end

function target(rng::Random.AbstractRNG,
                t::RandomTargeting,
                Ω::DiscreteFisheryDomain)
    sample(rng, Ω)
end

function target(t::RandomTargeting,
                Ω::DiscreteFisheryDomain)
    target(Random.GLOBAL_RNG, t, Ω)
end

## FIXME: Should reimplement these (or at least FixedTargeting) as iterators
"""
Vessels target fixed locations in given order. This struct is mutable and has
state to track which survey station is next. State can be reset using the
`reset!` function.
"""
mutable struct FixedTargeting{T} <: AbstractTargetingBehavior
    locations::Vector{T}
    state::T

    FixedTargeting(locs::Vector{T}) where T<:Integer = new{T}(locs, one(T))
end

function target(t::FixedTargeting,
                domain::AbstractFisheryDomain)
    loc = t.locations[t.state]
    t.state += 1
    loc
end

## RNG is ignored, but allowed for compatibility with other targeting behaviors.
function target(rng::Random.AbstractRNG,
                t::FixedTargeting,
                domain::AbstractFisheryDomain)
    target(t, domain)
end

length(t::FixedTargeting) = length(t.locations)

"""
    reset!(FT::FixedTargeting{T}) where T
    reset!(Targeting::AbstractTargetingBehavior)

For `FixedTargeting` objects, resets `state` to one in place. For other
targeting types, does (and returns) nothing.
"""
reset!(FT::FixedTargeting{T}) where T = FT.state = one(T)
reset!(Targeting::AbstractTargetingBehavior) = nothing
reset!(target::AbstractTargetingBehavior, pop::PopState) = reset!(target)

abstract type AbstractPreferentialTargeting <: AbstractTargetingBehavior end

"""
Vessels target according to some function or vector of preferences.

Accepts either a matrix of preference weights (don't necessarily need to add
to one) the same dimensions as the DiscreteFisheryDomain, or a function that
accepts a single argument (typically a Tuple or Array with two elements) of a
location and returns a preference weight.
"""
struct PreferentialTargeting{T} <: AbstractPreferentialTargeting
    preference::T

    PreferentialTargeting(preference::W) where W<:Weights = new{W}(preference)
end

## Outer constructor fallback to construct Weights object for the inner constructor
PreferentialTargeting(preference) = PreferentialTargeting(weights(preference))

function PreferentialTargeting(pref::Tp, Ω::DiscreteFisheryDomain) where {Tp <: AbstractArray}
    all(size(pref) .== size(Ω)) ||
        throw(DimensionMismatch("Preference array and domain dimensions must match."))
    PreferentialTargeting(weights(pref))
end
function PreferentialTargeting(f::Function, Ω::DiscreteFisheryDomain)
    pref = f.(Ω.locs)
    PreferentialTargeting(weights(pref))
end

function target(rng::Random.AbstractRNG,
                t::AbstractPreferentialTargeting,
                Ω::DiscreteFisheryDomain)
    sample(rng, Ω, t.preference)
end
function target(t::AbstractPreferentialTargeting,
                Ω::DiscreteFisheryDomain)
    target(Random.GLOBAL_RNG, t, Ω)
end

mutable struct DynamicPreferentialTargeting{W, F} <: AbstractPreferentialTargeting
    preference::W
    pref_fn::F

    function DynamicPreferentialTargeting(preference::W, pref_fn::F) where
        {W<:Weights, F<:Function}
        new{W, F}(preference, pref_fn)
    end
end

function DynamicPreferentialTargeting(preference::AbstractArray,
                                      pref_fn::Function)
    DynamicPreferentialTargeting(weights(pref_fn(preference)), pref_fn)
end

function reset!(dynpref::DynamicPreferentialTargeting, pop::PopState)
    dynpref.preference = weights(dynpref.pref_fn(pop.P))
end


text:
## Need to work on item delegates: http://qt-project.org/doc/qt-4.8/demos-spreadsheet-spreadsheetdelegate-cpp.html

## WE follow R's qtdataframe and embed roles into the data frame via
## column names so columns like x__DisplayRole, x__BackgroundRole,
## ... x__DecorationRole will have info looked up in.  we have issues
## storing the PyObjects in the data frame (they go out of scope when
## the view is shown) so instead you store strings that get mapped to
## pyobjects. E.g, for BackgrounRole a color name, for Decoration, a
## filename of an icon, ..
##
## Otherwise, we use these defaults:

display_role(x::DataFrames.DataArray, row::Int) = x[row]
display_role{T <: Real}(x::DataFrames.DataArray{T}, row::Int) = x[row]

edit_role(x::DataFrames.DataArray, row::Int) = x[row]
edit_role{T <: Real}(x::DataFrames.DataArray{T}, row::Int) = x[row]

text_alignment_role(x::DataFrames.DataArray, row::Integer) = convert(Int, qt_enum("AlignLeft"))
text_alignment_role{T <: Real}(x::DataFrames.DataArray{T}, row::Integer) = convert(Int, qt_enum("AlignRight"))

background_role(x::DataFrames.DataArray, row::Integer) = Qt[:QBrush](Qt[:QColor](0,0,0,0))
foreground_role(x::DataFrames.DataArray, row::Integer) = Qt[:QBrush](Qt[:QColor]("black"))

tool_tip_role(x::DataFrames.DataArray, row::Integer) = nothing
whats_this_role(x::DataFrames.DataArray, row::Integer) = nothing

function make_role(r::String, value)
    ## some roles require us to do things
    if  r == "TextAlignmentRole"
        qt_enum(value)                  # "AlignRight"
    elseif r == "DecorationRole"
        Qt[:QIcon](value)                 # filename
    elseif r == "BackgroundRole"
        Qt[:QBrush](Qt[:QColor](value))     # "blue"
    elseif r == "ForegroundRole"
        Qt[:QBrush](Qt[:QColor](value))
    else
        value
    end
end

## A DataFrameModel allows the contents of a data frame to set the model items
## Changes to the data frame are propogated to the views, one need not update the model.
type DataFrameModel <: PySide.QtWidget
    w
    d
    function DataFrameModel(parent)
        qnew_class("DataFrameModel", "QtCore.QAbstractTableModel")
        m = qnew_class_instance("DataFrameModel")
        qinvoke(m, :setParent, project(parent))
        new(m, nothing)
    end
end

function DataFrameModel(d, parent; editable=false)
    m = DataFrameModel(parent)
    set_items(m, d; editable=editable)
end

function set_items(m::DataFrameModel, d::DataFrames.DataFrame; editable=false)
    m.d = d

    role_re  = r"__[a-zA-Z]+Role$"
    nms = filter(u -> !ismatch(role_re, string(u)), names(d))

    project(m)[:rowCount] = (index) -> size(d)[1]
    project(m)[:columnCount] = (index) -> length(nms)

    ## data property
## Qt::DisplayRole	0	The key data to be rendered in the form of text. (QString)
## Qt::DecorationRole	1	The data to be rendered as a decoration in the form of an icon. (QColor, QIcon or QPixmap)
## Qt::EditRole	2	The data in a form suitable for editing in an editor. (QString)
## Qt::ToolTipRole	3	The data displayed in the item's tooltip. (QString)
## Qt::StatusTipRole	4	The data displayed in the status bar. (QString)
## Qt::WhatsThisRole	5	The data displayed for the item in "What's This?" mode. (QString)
## Qt::SizeHintRole	13	The size hint for the item that will be supplied to views. (QSize)
## Qt::FontRole	6	The font used for items rendered with the default delegate. (QFont)
## Qt::TextAlignmentRole	7	The alignment of the text for items rendered with the default delegate. (Qt::AlignmentFlag)
## Qt::BackgroundRole	8	The background brush used for items rendered with the default delegate. (QBrush)
## Qt::BackgroundColorRole	8	This role is obsolete. Use BackgroundRole instead.
## Qt::ForegroundRole	9	The foreground brush (text color, typically) used for items rendered with the default delegate. (QBrush)
## Qt::TextColorRole	9	This role is obsolete. Use ForegroundRole instead.
## Qt::CheckStateRole	10	This role is used to obtain the checked state of an item. (Qt::CheckState)
## Qt::InitialSortOrderRole	14	This role is used to obtain the initial sort order of a header view section. (Qt::SortOrder). This role was introduced in Qt 4.8.

    function do_role(idx, role)
        ## role is Int64

        row = idx[:row]() + 1
        col = idx[:column]() + 1
        nm = nms[col]

        ## check the rol
        cnames = map(string, names(m.d))
        roles = ["DisplayRole", "EditRole", "TextAlignmentRole", "BackgroundRole", "ForegroundRole", "ToolTipRole", "WhatsThisRole"]
        function role_default(r, row, col)
            if r == "DisplayRole"
                display_role(m.d[:,col], row)
            elseif r == "EditRole"
                edit_role(m.d[:,col], row)
            elseif r == "TextAlignmentRole"
                text_alignment_role(m.d[:,col], row)
            elseif r == "BackgroundRole"
                background_role(m.d[:,col], row)
            elseif r == "ForegroundRole"
                foreground_role(d[:,col], row)
            elseif r == "ToolTipRole"
                tool_tip_role(m.d[:,col], row)
            elseif r == "WhatsThisRole"
                whats_this_role(m.d[:,col], row)
            else
                nothing
            end
        end
        for r in roles
            if role == convert(Int, qt_enum(r))
                role_name = string(nm) * "__" * r
                out = (any(cnames .== role_name) ? make_role(r, m.d[row, symbol(role_name)]) : role_default(r, row, nm))
                return(out)
            end
        end
        return(nothing)
    end
    m.w[:data] = do_role

    ##  Well behaved models will also implement headerData().
    function header_data(section::Int, orient, role)
        if orient.o ==  qt_enum("Horizontal").o #  match pointers
            ## column, section is column
            role == convert(Int, qt_enum("DisplayRole")) ?  nms[section + 1] : nothing
        else
             role == convert(Int, qt_enum("DisplayRole")) ?  string(section + 1) : nothing
        end
    end
    project(m)[:headerData] = header_data

    ## Editable models need to implement setData(), and implement flags() to return a value containing Qt::ItemIsEditable.
    if editable
        function flags(idx)
            if !idx[:isValid]() return(qt_enum("ItemIsEnabled")) end

            row = idx[:row]()
            col = idx[:column]() + 1
            qt_enum(["ItemIsSelectable", "ItemIsEditable", "ItemIsEnabled"], how="|")
        end
        project(m)[:flags] = flags

        function set_data(idx, value, role)
            if !idx[:isValid]() return end

            row = idx[:row]() + 1
            col = idx[:column]() + 1
            nm = nms[col]
            m.d[row, nm] = value

            project(m)[:dataChanged][:emit](idx, idx)
        end
        project(m)[:setData] = set_data

    end
    m
end

    ## return array of selected values in [row,col] format
function selected_indicies(view)
    sel = view[:selectionModel]()
    if !sel[:hasSelection]() return [] end

    idxs = sel[:selectedIndexes]()
    map(u -> [u[:row]() + 1, u[:column]() + 1], idxs)
end

## ## Test it out
## testing = false
## if testing
##     ## test roles
##     d = DataFrame(x=[randstring(10) for i in 1:10],
##                   x__ToolTipRole = [randstring(10) for i in 1:10],
##                   x__BackgroundRole= rep("yellow", 10),
##                   y = 1:10
##                   )
##     view = Qt.QTableView()
##     m = DataFrameModel(d, parent=view)
##     view[:setModel](m)
##     raise(view)

##     ## test scaling
##     n = 10^5
##     d = DataFrame(x = randn(n),
##                   y = 1:n)

##     view = Qt.QTableView()
##     m = DataFrameModel(d, parent=view)
##     view[:setModel](m)
##     raise(view)
## end


text:
using DSGE
using HDF5, Base.Test

path = dirname(@__FILE__)

# # Test hessian! in context of model
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)

# Setup paths

data = h5read("$path/../reference/hessian.h5","data")'
mode = h5read("$path/../reference/hessian.h5","paramsmode")
hessian_expected = h5read("$path/../reference/hessian.h5","hessian")

# Test subset of hessian elements.
para_free      = [!θ.fixed for θ in m.parameters]
para_free_inds = find(para_free)

max_free_ind = DSGE.n_hessian_test_params(m)
if max_free_ind < maximum(para_free_inds)
    para_free_inds = para_free_inds[1:max_free_ind]
end

hessian, _ = hessian!(m, mode, data; verbose=:none)

# The values here are liable to differ between machines, based on different
# architectures/multithreading. Hessian values tend to be large, so these larger tolerances
# are needed.
expect = hessian_expected[1:max_free_ind, 1:max_free_ind]
actual = hessian[1:max_free_ind, 1:max_free_ind]

@testset "Check Hessian calculation" begin
    @test @test_matrix_approx_eq_eps expect actual 0.1 3.0
end

m.testing = false

nothing