text:
export CoffeeBreak, Coffee, CoffeeLounge
export init, enjoy, destroy
export add_coffee_break, check_coffee_break, setup
export update_statistics, get_statistics, save_statistics, shutdown

@compat abstract type Coffee end
function init(::Coffee, ::Net) end
# The first parameter will be a CoffeeLounge, we put Any here because
# Julia do not have forward declaration
function enjoy(::Any, ::Coffee, ::Net, ::SolverState) end
function destroy(::Coffee, ::Net) end

type CoffeeBreak
  coffee        :: Coffee
  every_n_iter  :: Int
  every_n_epoch :: Int
end

############################################################
# Coffee break utilities
############################################################
#-- This function is to be called by the end-user
function setup_coffee_lounge(solver::Solver; save_into::AbstractString="", every_n_iter::Int=1, file_exists=:merge)
  solver.coffee_lounge.filename=save_into
  solver.coffee_lounge.save_every_n_iter=every_n_iter
  solver.coffee_lounge.file_exists=file_exists
end

function add_coffee_break(solver::Solver, coffee::Coffee; kw...)
  add_coffee_break(solver.coffee_lounge, coffee; kw...)
end

################################################################################
# Coffee Lounge
################################################################################
using HDF5, JLD

const StatisticsValue = AbstractFloat
const StatisticsRecords = Dict{Int, StatisticsValue}
type CoffeeLounge
  filename          :: AbstractString
  save_every_n_iter :: Int
  file_exists       :: Symbol # :overwrite, :panic, :merge

  statistics        :: Dict{AbstractString, StatisticsRecords}
  coffee_breaks     :: Vector{CoffeeBreak}

  stats_modified    :: Bool
  last_epoch        :: Int
  curr_iter         :: Int

  CoffeeLounge(;filename="", save_every_n_iter=1, file_exists=:merge) = begin
    lounge = new(filename, save_every_n_iter, file_exists)
    lounge.statistics = Dict{AbstractString, StatisticsRecords}()
    lounge.coffee_breaks = CoffeeBreak[]
    return lounge
  end
end

function setup(lounge::CoffeeLounge, state::SolverState, net::Net)
  if !isempty(lounge.filename)
    directory = dirname(lounge.filename)
    if !isempty(directory) && !isdir(directory)
      mkdir_p(directory)
    end

    if isfile(lounge.filename)
      if lounge.file_exists == :overwrite
        @warn("Overwriting existing coffee lounge statistics in $(lounge.filename)")
      elseif lounge.file_exists == :merge
        @info("Merging existing coffee lounge statistics in $(lounge.filename)")
        lounge.statistics = jldopen(lounge.filename) do file
          read(file, "statistics")
        end
      elseif lounge.file_exists == :panic
        error("File already exists for coffee lounge statistics store: $(lounge.filename)")
      else
        error("Coffee lounge: unknown behavior on file exists: $(lounge.file_exists)")
      end
    end
  end

  lounge.stats_modified = false
  lounge.last_epoch = get_epoch(net)
  lounge.curr_iter = state.iter

  for cb in lounge.coffee_breaks
    init(cb.coffee, net)
  end
end

function update_statistics(dummy, key::AbstractString, val::StatisticsValue)
  # dummy function used when you do not want to record statistics
end

function update_statistics(lounge::CoffeeLounge, key::AbstractString, val::StatisticsValue)
  dict = get(lounge.statistics, key, StatisticsRecords())
  dict[lounge.curr_iter] = val
  lounge.statistics[key] = dict
  lounge.stats_modified = true
end

function get_statistics(lounge::CoffeeLounge, key::AbstractString)
  if haskey(lounge.statistics, key)
    lounge.statistics[key]
  else
    StatisticsRecords()
  end
end

function save_statistics(lounge::CoffeeLounge)
  if !isempty(lounge.filename) && lounge.stats_modified
    jldopen(lounge.filename, "w") do file
      write(file, "statistics", lounge.statistics)
    end
    lounge.stats_modified = false
  end
end

function shutdown(lounge::CoffeeLounge, net::Net)
  for cb in lounge.coffee_breaks
    destroy(cb.coffee, net)
  end

  save_statistics(lounge)
end

function add_coffee_break(lounge::CoffeeLounge, coffee::Coffee; every_n_iter::Int=0, every_n_epoch::Int=0)
  cb = CoffeeBreak(coffee, every_n_iter, every_n_epoch)
  push!(lounge.coffee_breaks, cb)
end

function check_coffee_break(lounge::CoffeeLounge, state::SolverState, net::Net)
  lounge.curr_iter = state.iter

  for cb in lounge.coffee_breaks
    if cb.every_n_iter > 0
      if state.iter % cb.every_n_iter == 0
        enjoy(lounge, cb.coffee, net, state)
      end
    elseif cb.every_n_epoch > 0
      epoch = get_epoch(net)
      if epoch != lounge.last_epoch && epoch % cb.every_n_epoch == 0
        enjoy(lounge, cb.coffee, net, state)
      end
    end
  end

  if state.iter % lounge.save_every_n_iter == 0
    save_statistics(lounge)
  end

  lounge.last_epoch = get_epoch(net)
end

include("coffee/training-summary.jl")
include("coffee/validation-performance.jl")
include("coffee/snapshot.jl")


text:
using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
len = length
<<<<<<< HEAD
filename = "./methane.nc"
=======
filename = "./cridun.nc"
>>>>>>> 93d403fe11aa311eacdc73eb249f740f6c2cb4bf

data = ncdata.get(filename,grp=0)
specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])
rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]
rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i])
    end

  end
end

smiles =readtable("carbons.csv")
smiles = Set(smiles[:species])

t = 216

links = filter(i -> flux[i][t]>0 , 1:len(flux))
<<<<<<< HEAD
tflux = [log10(flux[i][t]) for i in links]

# log after combining

weight =  tflux - minimum(tflux)
=======
tflux = [flux[i][t] for i in links] #[log10(flux[i][t]) for i in links]
weight =  tflux #- minimum(tflux)
>>>>>>> 93d403fe11aa311eacdc73eb249f740f6c2cb4bf
#weight = 1 - weight/maximum(weight)

#1+normalise(tflux)

newflux = [flux[i][t] for i in 1:len(flux)]
counter = 0
for i in links
  newflux[i] = weight[counter+=1]
end

edge = filter(i -> newflux[i[3]]>0 , edges)
source = [i[1] for i in edge]
target=[i[2] for i in edge]
weighted = [newflux[i[3]] for i in edge]
#grouping  = [reactiontypes[i[3]] for i in edge]

val = [i>0. for i in Array(specs[t,:])]
values = Set([string(i) for i in names(specs)[val]])

novalues = !val
novalues = [string(i) for i in names(specs)[novalues]]

     @rput sourceed
     @rput target
     @rput weighted

     R"library(igraph)"
     R"library(rgexf)"

     R"el <- structure(list(V1 = source, V2 = target, weight = weighted), .Names = c('V1',
     'V2', 'weight'), class = 'data.frame', row.names = c(1:$(length(edge))
     ))"

     R"g <- graph.data.frame(el)"

#c only
     R"v = V(g)$name"
     @rget v
     v=Set(v)

     # to remove carbons comment below
<<<<<<< HEAD
     #smiles = Set()

     smiles = intersect(smiles,novalues)

     diff = [i for i in setdiff(v,smiles)]
     @rput diff
     R"g = delete.vertices(g,diff)"

=======

     #comment if using smiles strings, else  keep commented
     #smiles = intersect(smiles,values)
     #diff = [i for i in setdiff(v,smiles)]

    diff = [i for i in intersect(v,Set(vcat(["NO","NO2","NO3","OH","O3","O","HONO","HO2","H2O2","HSO3","H2","HNO3","HO2NO2","HSO3","N2O5","SO2","SO3","SA"],novalues)))]

     @rput diff
     R"g = delete.vertices(g,diff)"

     R"g = simplify(g, edge.attr.comb='sum')"

      R"weights = log10(E(g)$weight)"

      R"weights = (abs(min(weights))+weights)"
      R"E(g)$weight = 1e-10 +weights" #"/max(weights)"

     R"v = V(g)$name"
     @rget v

     v= [specs[Symbol(i)][t] for i in v ]

     @rput v

     R"v = log10(v)"
     R"v= (v+abs(min(v)))"
     R"v=1e-10 + v/max(v)"

     R"V(g)$conc = v"

>>>>>>> 93d403fe11aa311eacdc73eb249f740f6c2cb4bf
R"""g1.gexf <- igraph.to.gexf(g)
f <- file('togephi.gexf')
writeLines(g1.gexf$graph, con = f)
close(f)
"""

'''
strength in

library(Hmisc)
library(ggplot2)
require(reshape2)

Flux = strength(g,vids=V(g),mode=c("all"))
Flux = Flux_in/max(Flux_in)

Concentration = V(g)$conc

EigenCentrality = eigen_centrality(g, directed=TRUE, weights=E(g)$weight)$vector
EigenCentrality = EigenCentrality/max(EigenCentrality)

PageRank= page_rank(g, algo = c('prpack', 'arpack', 'power'), vids = V(g),
     directed = TRUE, damping = 0.85, personalized = NULL, weights = NULL,
     options = NULL)$vector
PageRank=PageRank/max(PageRank)

Betweenness=betweenness(net, directed=T)
Betweenness = Betweenness/max(Betweenness)

Hubs = hub_score(net)$vector
Authorities=authority_score(net, weights=NA)$vector

df = data.frame(Flux,EigenCentrality,PageRank,Betweenness,Authorities)

df1 <- melt(df ,  id.vars = 'Flux', variable.name = 'series')
df1$Value = df1$value

ggplot(df1, aes(Flux,Value)) + geom_point(aes(colour = series))+
geom_smooth(method='lm')

with(df1, qplot(Flux, Value, colour = series, shape = series, ) +
  geom_smooth(method='lm')    )

  df = data.frame(Concentration,EigenCentrality,PageRank,Betweenness,Authorities)

  df1 <- melt(df ,  id.vars = 'Concentration', variable.name = 'series')
  df1$Value = df1$value

  with(df1, qplot(Concentration, Value, colour = series, shape = series, )

  ggsave('plot.svg', plot = last_plot(), device = svg, path = NULL,
      scale = 1, width = 6, height = 2,
      dpi = 300, limitsize = FALSE,)

'''

print("fini")

ext: jl

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
using Base.Test, Whittaker

@testset "whittaker2" begin
  @testset "out-of-place" begin
    tol = 0.0001

    y = [1.0, 1.2, 0, 1.6, 1.8]
    w = [1.0, 1.0, 0, 1.0, 1.0]
    lambda = 0.0001
    z = whittaker2(y, w, lambda)
    @test_approx_eq_eps z [0.999997, 1.20001, 1.4305, 1.64575, 1.80004] tol

    y = [1., 2, 0, 4, 5]
    w = [1, 1, 0, 1, 1]
    lambda = 0.0001
    z = whittaker2(y, w, lambda)
    @test_approx_eq_eps z [0.999992, 2.00002, 3.07626, 4.11438, 5.00005] tol
  end

  @testset "in-place" begin
    tol = 0.0001

    y = [1.0, 1.2, 0, 1.6, 1.8]
    w = [1.0, 1.0, 0, 1.0, 1.0]
    lambda = 0.0001
    whittaker2!(y, w, lambda)
    @test_approx_eq_eps y [0.999997, 1.20001, 1.4305, 1.64575, 1.80004] tol

    y = [1., 2, 0, 4, 5]
    w = [1, 1, 0, 1, 1]
    lambda = 0.0001
    whittaker2!(y, w, lambda)
    @test_approx_eq_eps y [0.999992, 2.00002, 3.07626, 4.11438, 5.00005] tol
  end
end


text:
function MLJModelInterface.fit(model::EvoTypes, verbosity::Int, A, y)
    fitresult, cache = init_evotree(model, A.matrix, y, verbosity=verbosity)
    grow_evotree!(fitresult, cache, verbosity=verbosity)
    report = (feature_importances = importance(fitresult, A.names), )
    return fitresult, cache, report
end

function okay_to_continue(new, old)
    new.nrounds - old.nrounds >= 0 &&
    new.loss == old.loss &&
    new.λ == old.λ &&
    new.γ == old.γ &&
    new.max_depth  == old.max_depth &&
    new.min_weight == old.min_weight &&
    new.rowsample ==  old.rowsample &&
    new.colsample ==  old.colsample &&
    new.nbins ==  old.nbins &&
    new.α ==  old.α &&
    new.metric ==  old.metric
end

# Generate names to be used by feature_importances in the report
MLJModelInterface.reformat(::EvoTypes, X, y) = ((matrix = MLJModelInterface.matrix(X), names = [name for name ∈ schema(X).names]), y)
MLJModelInterface.reformat(::EvoTypes, X) = ((matrix = MLJModelInterface.matrix(X), names = [name for name ∈ schema(X).names]),)
MLJModelInterface.reformat(::EvoTypes, X::AbstractMatrix, y) = ((matrix = X, names = ["feat_$i" for i in 1:size(X, 2)]), y)
MLJModelInterface.reformat(::EvoTypes, X::AbstractMatrix) = ((matrix = X, names = ["feat_$i" for i in 1:size(X, 2)]),)
MLJModelInterface.selectrows(::EvoTypes, I, A, y) = ((matrix = view(A.matrix, I, :), names = A.names), view(y, I))
MLJModelInterface.selectrows(::EvoTypes, I, A) = ((matrix = view(A.matrix, I, :), names = A.names),)

# For EarlyStopping.jl support
MLJModelInterface.iteration_parameter(::Type{<:EvoTypes}) = :nrounds

function MLJModelInterface.update(model::EvoTypes, verbosity::Integer, fitresult, cache, A, y)

    if okay_to_continue(model, cache.params)
        grow_evotree!(fitresult, cache, verbosity=verbosity)
    else
        fitresult, cache = init_evotree(model, A.matrix, y, verbosity=verbosity)
        grow_evotree!(fitresult, cache, verbosity=verbosity)
    end

    report = (feature_importances = importance(fitresult, A.names),)

    return fitresult, cache, report
end

function predict(::EvoTreeRegressor, fitresult, A)
    pred = predict(fitresult, A.matrix)
    return pred
end

function predict(::EvoTreeClassifier, fitresult, A)
    pred = predict(fitresult, A.matrix)
    return MLJModelInterface.UnivariateFinite(fitresult.levels, pred, pool=missing)
end

function predict(::EvoTreeCount, fitresult, A)
    λ = predict(fitresult, A.matrix)
    return [Distributions.Poisson(λᵢ) for λᵢ ∈ λ]
end

function predict(::EvoTreeGaussian, fitresult, A)
    pred = predict(fitresult, A.matrix)
    return [Distributions.Normal(pred[i,1], pred[i,2]) for i in 1:size(pred, 1)]
end

# Metadata
const EvoTreeRegressor_desc = "Regression models with various underlying methods: least square, quantile, logistic."
const EvoTreeClassifier_desc = "Multi-classification with softmax and cross-entropy loss."
const EvoTreeCount_desc = "Poisson regression fitting λ with max likelihood."
const EvoTreeGaussian_desc = "Gaussian maximum likelihood of μ and σ²."

MLJModelInterface.metadata_pkg.((EvoTreeRegressor, EvoTreeClassifier, EvoTreeCount, EvoTreeGaussian),
    name="EvoTrees",
    uuid="f6006082-12f8-11e9-0c9c-0d5d367ab1e5",
    url="https://github.com/Evovest/EvoTrees.jl",
    julia=true,
    license="Apache",
    is_wrapper=false)

MLJModelInterface.metadata_model(EvoTreeRegressor,
    input=Union{MLJModelInterface.Table(MLJModelInterface.Continuous),AbstractMatrix{MLJModelInterface.Continuous}},
    target=AbstractVector{<:MLJModelInterface.Continuous},
    weights=false,
    path="EvoTrees.EvoTreeRegressor",
    descr=EvoTreeRegressor_desc)

MLJModelInterface.metadata_model(EvoTreeClassifier,
        input=Union{MLJModelInterface.Table(MLJModelInterface.Continuous),AbstractMatrix{MLJModelInterface.Continuous}},
        target=AbstractVector{<:MLJModelInterface.Finite},
        weights=false,
        path="EvoTrees.EvoTreeClassifier",
        descr=EvoTreeClassifier_desc)

MLJModelInterface.metadata_model(EvoTreeCount,
    input=Union{MLJModelInterface.Table(MLJModelInterface.Continuous),AbstractMatrix{MLJModelInterface.Continuous}},
    target=AbstractVector{<:MLJModelInterface.Count},
    weights=false,
    path="EvoTrees.EvoTreeCount",
    descr=EvoTreeCount_desc)

MLJModelInterface.metadata_model(EvoTreeGaussian,
    input=Union{MLJModelInterface.Table(MLJModelInterface.Continuous),AbstractMatrix{MLJModelInterface.Continuous}},
    target=AbstractVector{<:MLJModelInterface.Continuous},
    weights=false,
    path="EvoTrees.EvoTreeGaussian",
    descr=EvoTreeGaussian_desc)

# shared metadata
# MLJModelInterface.package_name(::Type{<:EvoTypes}) = "EvoTrees"
# MLJModelInterface.package_uuid(::Type{<:EvoTypes}) = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
# MLJModelInterface.package_url(::Type{<:EvoTypes}) = "https://github.com/Evovest/EvoTrees.jl"
# MLJModelInterface.is_pure_julia(::Type{<:EvoTypes}) = true
#
# MLJModelInterface.load_path(::Type{<:EvoTreeRegressor}) = "EvoTrees.EvoTreeRegressor"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeRegressor}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeRegressor}) = AbstractVector{<:MLJModelInterface.Continuous}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeCount}) = "EvoTrees.EvoTreeCount"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeCount}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeCount}) = AbstractVector{<:MLJModelInterface.Count}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeClassifier}) = "EvoTrees.EvoTreeClassifier"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeClassifier}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeClassifier}) = AbstractVector{<:MLJModelInterface.Finite}
#
# MLJModelInterface.load_path(::Type{<:EvoTreeGaussian}) = "EvoTrees.EvoTreeGaussian"
# MLJModelInterface.input_scitype(::Type{<:EvoTreeGaussian}) = MLJModelInterface.Table(MLJModelInterface.Continuous)
# MLJModelInterface.target_scitype(::Type{<:EvoTreeGaussian}) = AbstractVector{<:MLJModelInterface.Continuous}

# function MLJ.clean!(model::EvoTreeRegressor)
#     warning = ""
#     if model.nrounds < 1
#         warning *= "Need nrounds ≥ 1. Resetting nrounds=1. "
#         model.nrounds = 1
#     end
#     if model.λ < 0
#         warning *= "Need λ ≥ 0. Resetting λ=0. "
#         model.λ = 0.0
#     end
#     if model.γ < 0
#         warning *= "Need γ ≥ 0. Resetting γ=0. "
#         model.γ = 0.0
#     end
#     if model.η <= 0
#         warning *= "Need η > 0. Resetting η=0.001. "
#         model.η = 0.001
#     end
#     if model.max_depth < 1
#         warning *= "Need max_depth ≥ 0. Resetting max_depth=0. "
#         model.max_depth = 1
#     end
#     if model.min_weight < 0
#         warning *= "Need min_weight ≥ 0. Resetting min_weight=0. "
#         model.min_weight = 0.0
#     end
#     if model.rowsample < 0
#         warning *= "Need rowsample ≥ 0. Resetting rowsample=0. "
#         model.rowsample = 0.0
#     end
#     if model.rowsample > 1
#         warning *= "Need rowsample <= 1. Resetting rowsample=1. "
#         model.rowsample = 1.0
#     end
#     if model.colsample < 0
#         warning *= "Need colsample ≥ 0. Resetting colsample=0. "
#         model.colsample = 0.0
#     end
#     if model.colsample > 1
#         warning *= "Need colsample <= 1. Resetting colsample=1. "
#         model.colsample = 1.0
#     end
#     if model.nbins > 250
#         warning *= "Need nbins <= 250. Resetting nbins=250. "
#         model.nbins = 250
#     end
#     return warning
# end


text:
@testset "OpenCL.Minver" begin
    @testset "OpenCL.Minver platform" begin
        for platform in cl.platforms()

            version = cl.opencl_version(platform)

            v11 = cl.min_v11(platform)
            v12 = cl.min_v12(platform)
            v20 = cl.min_v20(platform)
            v21 = cl.min_v21(platform)
            v22 = cl.min_v22(platform)

            @test v11 == (version >= v"1.1")
            @test v12 == (version >= v"1.2")
            @test v20 == (version >= v"2.0")
            @test v21 == (version >= v"2.1")
            @test v22 == (version >= v"2.2")
        end
    end

    @testset "OpenCL.Minver device" begin
        for platform in cl.platforms()
            for device in cl.devices(platform)
                version = cl.opencl_version(device)

                v11 = cl.min_v11(device)
                v12 = cl.min_v12(device)
                v20 = cl.min_v20(device)
                v21 = cl.min_v21(device)
                v22 = cl.min_v22(device)

                @test v11 == (version >= v"1.1")
                @test v12 == (version >= v"1.2")
                @test v20 == (version >= v"2.0")
                @test v21 == (version >= v"2.1")
                @test v22 == (version >= v"2.2")
            end
        end
    end
end


text:
@defcomp impactsealevelrise begin
    regions = Index()

    wetval = Variable(index=[time,regions]) # Valuation of wetlands
    wetlandloss = Variable(index=[time,regions]) # Amount of wetland lost
    cumwetlandloss = Variable(index=[time,regions]) # Cumulative wetland loss over time
    wetlandgrowth = Variable(index=[time,regions])
    wetcost = Variable(index=[time,regions]) # The cost of wetland actually lost in this year

    dryval = Variable(index=[time,regions]) # Valuation of dryland
    landloss = Variable(index=[time,regions]) # Dryland lost in this year
    cumlandloss = Variable(index=[time,regions]) # Cumulative dryland loss over time
    drycost = Variable(index=[time,regions]) # The cost of dryland actually lost in this year

    npprotcost = Variable(index=[time,regions])
    npwetcost = Variable(index=[time,regions])
    npdrycost = Variable(index=[time,regions])

    protlev = Variable(index=[time,regions]) # Share of coastline protected in this time period
    protcost = Variable(index=[time,regions]) # The cost of protection in this time period

    enter = Variable(index=[time,regions]) # Migration in (in people)
    leave = Variable(index=[time,regions]) # Migration out (in people)
    entercost = Variable(index=[time,regions]) # Cost to the economy of immigration
    leavecost = Variable(index=[time,regions]) # Cost to the economy of emigration

    imigrate = Variable(index=[regions,regions]) # Migration matrix

    incdens     = Parameter(default = 0.000635) # Normalization income density
    emcst       = Parameter(default = 3) # emigration loss benchmark value
    immcst      = Parameter(default = 0.4) # immigration loss benchmark
    dvydl       = Parameter(default = 1) # income density elasticity of dryland value
    wvel        = Parameter(default = 1.16) # wetland value income elasticity
    wvbm        = Parameter(default = 0.00588) # wetland value benchmark value
    slrwvpopdens0 = Parameter(default = 27.5937717888728) # wetland value income normalization value
    wvpdl       = Parameter(default = 0.47)
    wvsl        = Parameter(default = -0.11)
    dvbm        = Parameter(default = 0.004) # dryland value benchmark value
    slrwvypc0   = Parameter(default = 25000)

    pc          = Parameter(index=[regions])
    slrprtp     = Parameter(index=[regions])
    wmbm        = Parameter(index=[regions])
    dlbm        = Parameter(index=[regions]) # dryland loss benchmark value
    drylandlossparam = Parameter(index=[regions])
    wlbm        = Parameter(index=[regions]) # wetland loss benchmark per slr change
    coastpd     = Parameter(index=[regions]) # coastal population density
    wetmax      = Parameter(index=[regions])
    wetland90   = Parameter(index=[regions])
    maxlandloss = Parameter(index=[regions])

    sea = Parameter(index=[time]) # sea-level rise in m since pre-industrial

    migrate = Parameter(index=[regions,regions])

    income = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])
    area = Parameter(index=[time,regions]) # Land area km^2

    function run_timestep(p, v, d, t)

        if is_first(t)
            for r1 in d.regions
                for r2 in d.regions
                    immsumm = 0
                    for i in d.regions
                        immsumm += p.migrate[i, r1]
                    end
                    v.imigrate[r1, r2] = p.migrate[r2, r1] / immsumm
                end

                t0 = TimestepIndex(1)
                v.landloss[t0, r1] = 0.0
                v.cumlandloss[t0, r1] = 0.0
                v.cumwetlandloss[t0, r1] = 0.0
                v.wetlandgrowth[t0, r1] = 0.0
            end
        else
            # slr in m/year
            ds = p.sea[t] - p.sea[t - 1]

            for r in d.regions
                ypc = p.income[t, r] / p.population[t, r] * 1000.0 # income per capita
                ypcprev = p.income[t - 1, r] / p.population[t - 1, r] * 1000.0
                ypcgrowth = ypc / ypcprev - 1.0

                if gettime(t) == 1951
                    ypcgrowth = 0
                end

                # Needs to be in $bn per km^2
                # Income is in billion, area is in km^2
                incomedens = p.income[t, r] / p.area[t, r]

                incomedensprev = p.income[t - 1, r] / p.area[t - 1, r]

                incomedensgrowth = incomedens / incomedensprev - 1.0

                # In population/km^2
                ## population is in million, area is in km^2
                popdens = p.population[t, r] / p.area[t, r] * 1000000.0
                popdensprev = p.population[t - 1, r] / p.area[t - 1, r] * 1000000.0
                popdensgrowth = popdens / popdensprev - 1.0

                # Unit of dryval is $bn/km^2
                v.dryval[t, r] = p.dvbm * (incomedens / p.incdens)^p.dvydl

                # Unit of wetval is $bn/km^2
                v.wetval[t, r] = p.wvbm *
                    (ypc / p.slrwvypc0)^p.wvel *
                    (popdens / p.slrwvpopdens0)^p.wvpdl *
                    ((p.wetland90[r] - v.cumwetlandloss[t - 1, r]) / p.wetland90[r])^p.wvsl

                potCumLandloss = min(p.maxlandloss[r], p.dlbm[r] * p.sea[t]^p.drylandlossparam[r])

                potLandloss = potCumLandloss - v.cumlandloss[t - 1, r]

                ## If sea levels fall, no protection is build
                if ds < 0
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 0
                # If the discount rate < -100% people will not build protection
                elseif (1.0 + p.slrprtp[r] + ypcgrowth) < 0.0
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 0
                # Dryland value is worthless
                elseif (1.0 + p.dvydl * incomedensgrowth) < 0.0
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 0
                # Is protecting the coast infinitly expensive?
                elseif (1.0 / (1.0 + p.slrprtp[r] + ypcgrowth)) >= 1
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 0
                # Is dryland infinitly valuable?
                elseif ((1.0 + p.dvydl * incomedensgrowth) / (1.0 + p.slrprtp[r] + ypcgrowth)) >= 1.0
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 1
                # Is wetland infinitly valuable?
                elseif ((1.0 + p.wvel * ypcgrowth + p.wvpdl * popdensgrowth + p.wvsl * v.wetlandgrowth[t - 1, r]) / (1.0 + p.slrprtp[r] + ypcgrowth)) >= 1.0
                    v.npprotcost[t, r] = 0
                    v.npwetcost[t, r] = 0
                    v.npdrycost[t, r] = 0
                    v.protlev[t, r] = 0
                else
                    # NPV of protecting the whole coast
                    # pc is in $bn/m
                    v.npprotcost[t, r] = p.pc[r] * ds * (1.0 + p.slrprtp[r] + ypcgrowth) / (p.slrprtp[r] + ypcgrowth)

                    # NPV of wetland
                    if (1.0 + p.wvel * ypcgrowth + p.wvpdl * popdensgrowth + p.wvsl * v.wetlandgrowth[t - 1, r]) < 0.0
                        v.npwetcost[t, r] = 0
                    else
                        v.npwetcost[t, r] = p.wmbm[r] * ds * v.wetval[t, r] * (1.0 + p.slrprtp[r] + ypcgrowth) / (p.slrprtp[r] + ypcgrowth - p.wvel * ypcgrowth - p.wvpdl * popdensgrowth - p.wvsl * v.wetlandgrowth[t - 1, r])
                    end

                    # NPV of dryland
                    if (1.0 + p.dvydl * incomedensgrowth) < 0.0
                        v.npdrycost[t, r] = 0
                    else
                        v.npdrycost[t, r] = potLandloss * v.dryval[t, r] * (1 + p.slrprtp[r] + ypcgrowth) / (p.slrprtp[r] + ypcgrowth - p.dvydl * incomedensgrowth)
                    end

                    # Calculate protection level
                    v.protlev[t, r] = max(0.0, 1.0 - 0.5 * (v.npprotcost[t, r] + v.npwetcost[t, r]) / v.npdrycost[t, r])

                    if v.protlev[t, r] > 1
                        error("protlevel >1 should not happen")
                    end
                end

                # Calculate actual wetland loss and cost
                v.wetlandloss[t, r] = min(
                    p.wlbm[r] * ds + v.protlev[t, r] * p.wmbm[r] * ds,
                    p.wetmax[r] - v.cumwetlandloss[t - 1, r])

                v.cumwetlandloss[t, r] = v.cumwetlandloss[t - 1, r] + v.wetlandloss[t, r]

                # Calculate wetland growth
                v.wetlandgrowth[t, r] = (p.wetland90[r] - v.cumwetlandloss[t, r]) / (p.wetland90[r] - v.cumwetlandloss[t - 1, r]) - 1.0

                v.wetcost[t, r] = v.wetval[t, r] * v.wetlandloss[t, r]

                v.landloss[t, r] = (1.0 - v.protlev[t, r]) * potLandloss

                v.cumlandloss[t, r] = v.cumlandloss[t - 1, r] + v.landloss[t, r]
                v.drycost[t, r] = v.dryval[t, r] * v.landloss[t, r]

                v.protcost[t, r] = v.protlev[t, r] * p.pc[r] * ds

                if v.landloss[t, r] < 0
                    v.leave[t, r] = 0
                else
                    v.leave[t, r] = p.coastpd[r] * popdens * v.landloss[t, r]
                end

                v.leavecost[t, r] = p.emcst * ypc * v.leave[t, r] / 1000000000
            end

            for destination in d.regions
                enter = 0.0
                for source in d.regions
                    enter += v.leave[t, source] * v.imigrate[source, destination]
                end
                v.enter[t, destination] = enter
            end

            for r in d.regions
                ypc = p.income[t, r] / p.population[t, r] * 1000.0
                v.entercost[t, r] = p.immcst * ypc * v.enter[t, r] / 1000000000
            end
        end
    end
end


text:
using KernelAbstractions: @kernel, @index

const R_Earth = 6371.0e3    # Mean radius of the Earth [m] https://en.wikipedia.org/wiki/Earth

struct LatitudeLongitudeGrid{FT, TX, TY, TZ, M, MY, FX, FY, FZ, VX, VY, VZ, Arch} <: AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ, Arch}
    architecture :: Arch
    Nx :: Int
    Ny :: Int
    Nz :: Int
    Hx :: Int
    Hy :: Int
    Hz :: Int
    Lx :: FT
    Ly :: FT
    Lz :: FT
    # All directions can be either regular (FX, FY, FZ) <: Number
    # or stretched (FX, FY, FZ) <: AbstractVector
    Δλᶠᵃᵃ :: FX
    Δλᶜᵃᵃ :: FX
    λᶠᵃᵃ  :: VX
    λᶜᵃᵃ  :: VX
    Δφᵃᶠᵃ :: FY
    Δφᵃᶜᵃ :: FY
    φᵃᶠᵃ  :: VY
    φᵃᶜᵃ  :: VY
    Δzᵃᵃᶠ :: FZ
    Δzᵃᵃᶜ :: FZ
    zᵃᵃᶠ  :: VZ
    zᵃᵃᶜ  :: VZ
    # Precomputed metrics M <: Nothing means metrics will be computed on the fly
    Δxᶠᶜᵃ :: M
    Δxᶜᶠᵃ :: M
    Δxᶠᶠᵃ :: M
    Δxᶜᶜᵃ :: M
    Δyᶠᶜᵃ :: MY
    Δyᶜᶠᵃ :: MY
    Azᶠᶜᵃ :: M
    Azᶜᶠᵃ :: M
    Azᶠᶠᵃ :: M
    Azᶜᶜᵃ :: M
    # Spherical radius
    radius :: FT

    function LatitudeLongitudeGrid{TX, TY, TZ}(architecture::Arch,
                                               Nλ, Nφ, Nz,
                                               Hλ, Hφ, Hz,
                                               Lλ::FT, Lφ::FT, Lz::FT,
                                               Δλᶠᵃᵃ :: FX, Δλᶜᵃᵃ :: FX,
                                                λᶠᵃᵃ :: VX,  λᶜᵃᵃ :: VX,
                                               Δφᵃᶠᵃ :: FY, Δφᵃᶜᵃ :: FY,
                                                φᵃᶠᵃ :: VY,  φᵃᶜᵃ :: VY,
                                               Δzᵃᵃᶠ :: FZ, Δzᵃᵃᶜ :: FZ,
                                                zᵃᵃᶠ :: VZ,  zᵃᵃᶜ :: VZ,
                                               Δxᶠᶜ::M, Δxᶜᶠ::M,
                                               Δxᶠᶠ::M, Δxᶜᶜ::M,
                                               Δyᶠᶜ::MY, Δyᶜᶠ::MY,
                                               Azᶠᶜ::M, Azᶜᶠ::M, Azᶠᶠ::M, Azᶜᶜ::M,
                                               radius::FT) where {Arch, FT,
                                                                  TX, TY, TZ,
                                                                  FX, FY, FZ,
                                                                  VX, VY, VZ,
                                                                  M, MY}

        return new{FT, TX, TY, TZ, M, MY, FX, FY, FZ, VX, VY, VZ, Arch}(architecture,
                                                                        Nλ, Nφ, Nz,
                                                                        Hλ, Hφ, Hz,
                                                                        Lλ, Lφ, Lz,
                                                                        Δλᶠᵃᵃ, Δλᶜᵃᵃ, λᶠᵃᵃ, λᶜᵃᵃ,
                                                                        Δφᵃᶠᵃ, Δφᵃᶜᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
                                                                        Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ,
                                                                        Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ,
                                                                        Δyᶠᶜ, Δyᶜᶠ,
                                                                        Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ, radius)
    end
end

const XRegLatLonGrid = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Number}
const YRegLatLonGrid = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any,    <:Number}
const ZRegLatLonGrid = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any,    <:Any, <:Number}
const HRegLatLonGrid = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Number, <:Number}

regular_dimensions(::ZRegLatLonGrid) = tuple(3)

"""
    LatitudeLongitudeGrid([architecture = CPU(), FT = Float64];
                          size,
                          longitude,
                          latitude,
                          z,
                          radius = R_Earth,
                          topology = nothing,
                          precompute_metrics = true,
                          halo = nothing)

Creates a `LatitudeLongitudeGrid` with coordinates `(λ, φ, z)` denoting longitude, latitude,
and vertical coordinate respectively.

Positional arguments
====================

- `architecture`: Specifies whether arrays of coordinates and spacings are stored
                  on the CPU or GPU. Default: `CPU()`.

- `FT` : Floating point data type. Default: `Float64`.

Keyword arguments
=================

- `size` (required): A 3-tuple prescribing the number of grid points each direction.

- `longitude`, `latitude`, `z` (required): Each is either a
                                           (i) 2-tuple that specify the end points of the domain,
                                           (ii) one-dimensional array specifying the cell interface locations or
                                           (iii) a single-argument function that takes an index and returns
                                                 cell interface location.
  **Note**: the latitude and longitude coordinates extents are expected in degrees.

- `radius`: The radius of the sphere the grid lives on. By default is equal to the radius of Earth.

- `topology`: Tuple of topologies (`Flat`, `Bounded`, `Periodic`) for each direction. The vertical
              `topology[3]` must be `Bounded`, while the latitude-longitude topologies can be
              `Bounded`, `Periodic`, or `Flat`.

- `precompute_metrics`: Boolean specifying whether to precompute horizontal spacings and areas.
                        Default: `true`. When `false`, horizontal spacings and areas are computed
                        on-the-fly during a simulation.

- `halo`: A 3-tuple of integers specifying the size of the halo region of cells surrounding
          the physical interior. The default is 3 halo cells in every direction.
"""
function LatitudeLongitudeGrid(architecture::AbstractArchitecture = CPU(),
                               FT::DataType = Float64;
                               size,
                               longitude,
                               latitude,
                               z,
                               radius = R_Earth,
                               topology = nothing,
                               precompute_metrics = true,
                               halo = nothing)

    Nλ, Nφ, Nz, Hλ, Hφ, Hz, latitude, longitude, topology =
        validate_lat_lon_grid_args(latitude, longitude, size, halo, topology)

    # Calculate all direction (which might be stretched)
    # A direction is regular if the domain passed is a Tuple{<:Real, <:Real},
    # it is stretched if being passed is a function or vector (as for the VerticallyStretchedRectilinearGrid)
    TX, TY, TZ = topology

    Lλ, λᶠᵃᵃ, λᶜᵃᵃ, Δλᶠᵃᵃ, Δλᶜᵃᵃ = generate_coordinate(FT, TX, Nλ, Hλ, longitude, architecture)
    Lφ, φᵃᶠᵃ, φᵃᶜᵃ, Δφᵃᶠᵃ, Δφᵃᶜᵃ = generate_coordinate(FT, TY, Nφ, Hφ, latitude,  architecture)
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, TZ, Nz, Hz, z,         architecture)

    preliminary_grid = LatitudeLongitudeGrid{TX, TY, TZ}(architecture,
                                                         Nλ, Nφ, Nz,
                                                         Hλ, Hφ, Hz,
                                                         Lλ, Lφ, Lz,
                                                         Δλᶠᵃᵃ, Δλᶜᵃᵃ, λᶠᵃᵃ, λᶜᵃᵃ,
                                                         Δφᵃᶠᵃ, Δφᵃᶜᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
                                                         Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ,
                                                         (nothing for i=1:10)..., FT(radius))

    return !precompute_metrics ? preliminary_grid : with_precomputed_metrics(preliminary_grid)
end

# architecture = CPU() default, assuming that a DataType positional arg
# is specifying the floating point type.
LatitudeLongitudeGrid(FT::DataType; kwargs...) = LatitudeLongitudeGrid(CPU(), FT; kwargs...)

""" Return a reproduction of `grid` with precomputed metric terms. """
function with_precomputed_metrics(grid)
    arch = architecture(grid)
    Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Δyᶠᶜ, Δyᶜᶠ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ = allocate_metrics(grid)
    wait(device_event(arch))

    precompute_curvilinear_metrics!(grid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ )
    wait(device_event(arch))

    Δyᶠᶜ, Δyᶜᶠ = precompute_Δy_metrics(grid, Δyᶠᶜ, Δyᶜᶠ)

    Nλ, Nφ, Nz = size(grid)
    Hλ, Hφ, Hz = halo_size(grid)
    TX, TY, TZ = topology(grid)

    return LatitudeLongitudeGrid{TX, TY, TZ}(architecture(grid),
                                             Nλ, Nφ, Nz,
                                             Hλ, Hφ, Hz,
                                             grid.Lx, grid.Ly, grid.Lz,
                                             grid.Δλᶠᵃᵃ, grid.Δλᶜᵃᵃ, grid.λᶠᵃᵃ, grid.λᶜᵃᵃ,
                                             grid.Δφᵃᶠᵃ, grid.Δφᵃᶜᵃ, grid.φᵃᶠᵃ, grid.φᵃᶜᵃ,
                                             grid.Δzᵃᵃᶠ, grid.Δzᵃᵃᶜ, grid.zᵃᵃᶠ, grid.zᵃᵃᶜ,
                                             Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ,
                                             Δyᶠᶜ, Δyᶜᶠ,
                                             Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ, grid.radius)
end

function validate_lat_lon_grid_args(latitude, longitude, size, halo, topology)

    λ₁, λ₂ = get_domain_extent(longitude, size[1])
    @assert λ₁ < λ₂ && λ₂ - λ₁ ≤ 360

    φ₁, φ₂ = get_domain_extent(latitude, size[2])
    @assert -90 <= φ₁ < φ₂ <= 90

    (φ₁ == -90 || φ₂ == 90) &&
        @warn "Are you sure you want to use a latitude-longitude grid with a grid point at the pole?"

    Lλ = λ₂ - λ₁

    if !isnothing(topology)
        TX, TY, TZ = topology
        TZ === Bounded || throw(ArgumentError("z topology must be Bounded"))
    else
        TX = Lλ == 360 ? Periodic : Bounded
        TY = Bounded
        TZ = Bounded
    end

    Nλ, Nφ, Nz = N = validate_size(TX, TY, TZ, size)
    Hλ, Hφ, Hz = H = validate_halo(TX, TY, TZ, halo)

    return Nλ, Nφ, Nz, Hλ, Hφ, Hz, latitude, longitude, (TX, TY, TZ)
end

function Base.summary(grid::LatitudeLongitudeGrid)
    FT = eltype(grid)
    TX, TY, TZ = topology(grid)
    metric_computation = isnothing(grid.Δxᶠᶜᵃ) ? "without precomputed metrics" : "with precomputed metrics"

    return string(size_summary(size(grid)),
                  " LatitudeLongitudeGrid{$FT, $TX, $TY, $TZ} on ", summary(architecture(grid)),
                  " with ", size_summary(halo_size(grid)), " halo",
                  " and ", metric_computation)
end

function Base.show(io::IO, grid::LatitudeLongitudeGrid)
    TX, TY, TZ = topology(grid)

    λ₁, λ₂ = domain(topology(grid, 1), grid.Nx, grid.λᶠᵃᵃ)
    φ₁, φ₂ = domain(topology(grid, 2), grid.Ny, grid.φᵃᶠᵃ)
    z₁, z₂ = domain(topology(grid, 3), grid.Nz, grid.zᵃᵃᶠ)

    x_summary = domain_summary(TX(), "λ", λ₁, λ₂)
    y_summary = domain_summary(TY(), "φ", φ₁, φ₂)
    z_summary = domain_summary(TZ(), "z", z₁, z₂)

    longest = max(length(x_summary), length(y_summary), length(z_summary))

    x_summary = "longitude: " * dimension_summary(TX(), "λ", λ₁, λ₂, grid.Δλᶜᵃᵃ, longest - length(x_summary))
    y_summary = "latitude:  " * dimension_summary(TY(), "φ", φ₁, φ₂, grid.Δφᵃᶜᵃ, longest - length(y_summary))
    z_summary = "z:         " * dimension_summary(TZ(), "z", z₁, z₂, grid.Δzᵃᵃᶜ, longest - length(z_summary))

    print(io, summary(grid), '\n',
          "├── ", x_summary, '\n',
          "├── ", y_summary, '\n',
          "└── ", z_summary)
end

# Node by node
@inline xnode(::Face,   i, grid::LatitudeLongitudeGrid) = @inbounds grid.λᶠᵃᵃ[i]
@inline xnode(::Center, i, grid::LatitudeLongitudeGrid) = @inbounds grid.λᶜᵃᵃ[i]
@inline ynode(::Face,   j, grid::LatitudeLongitudeGrid) = @inbounds grid.φᵃᶠᵃ[j]
@inline ynode(::Center, j, grid::LatitudeLongitudeGrid) = @inbounds grid.φᵃᶜᵃ[j]
@inline znode(::Face,   k, grid::LatitudeLongitudeGrid) = @inbounds grid.zᵃᵃᶠ[k]
@inline znode(::Center, k, grid::LatitudeLongitudeGrid) = @inbounds grid.zᵃᵃᶜ[k]

all_x_nodes(::Type{Face},   grid::LatitudeLongitudeGrid) = grid.λᶠᵃᵃ
all_x_nodes(::Type{Center}, grid::LatitudeLongitudeGrid) = grid.λᶜᵃᵃ
all_y_nodes(::Type{Face},   grid::LatitudeLongitudeGrid) = grid.φᵃᶠᵃ
all_y_nodes(::Type{Center}, grid::LatitudeLongitudeGrid) = grid.φᵃᶜᵃ
all_z_nodes(::Type{Face},   grid::LatitudeLongitudeGrid) = grid.zᵃᵃᶠ
all_z_nodes(::Type{Center}, grid::LatitudeLongitudeGrid) = grid.zᵃᵃᶜ

@inline x_domain(grid::LatitudeLongitudeGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TX, grid.Nx, grid.λᶠᵃᵃ)
@inline y_domain(grid::LatitudeLongitudeGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TY, grid.Ny, grid.φᵃᶠᵃ)
@inline z_domain(grid::LatitudeLongitudeGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TZ, grid.Nz, grid.zᵃᵃᶠ)

@inline cpu_face_constructor_x(grid::XRegLatLonGrid) = x_domain(grid)
@inline cpu_face_constructor_y(grid::YRegLatLonGrid) = y_domain(grid)
@inline cpu_face_constructor_z(grid::ZRegLatLonGrid) = z_domain(grid)

function with_halo(new_halo, old_grid::LatitudeLongitudeGrid)

    size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size     = pop_flat_elements(size, topo)
    new_halo = pop_flat_elements(new_halo, topo)

    new_grid = LatitudeLongitudeGrid(architecture(old_grid), eltype(old_grid);
                                     size = size, halo = new_halo,
                                     longitude = x, latitude = y, z = z,
                                     precompute_metrics = metrics_precomputed(old_grid))

    return new_grid
end

function on_architecture(new_arch::AbstractArchitecture, old_grid::LatitudeLongitudeGrid)
    old_properties = (old_grid.Δλᶠᵃᵃ, old_grid.Δλᶜᵃᵃ, old_grid.λᶠᵃᵃ,  old_grid.λᶜᵃᵃ,
                      old_grid.Δφᵃᶠᵃ, old_grid.Δφᵃᶜᵃ, old_grid.φᵃᶠᵃ,  old_grid.φᵃᶜᵃ,
                      old_grid.Δzᵃᵃᶠ, old_grid.Δzᵃᵃᶜ, old_grid.zᵃᵃᶠ,  old_grid.zᵃᵃᶜ,
                      old_grid.Δxᶠᶜᵃ, old_grid.Δxᶜᶠᵃ, old_grid.Δxᶠᶠᵃ, old_grid.Δxᶜᶜᵃ,
                      old_grid.Δyᶠᶜᵃ, old_grid.Δyᶜᶠᵃ,
                      old_grid.Azᶠᶜᵃ, old_grid.Azᶜᶠᵃ, old_grid.Azᶠᶠᵃ, old_grid.Azᶜᶜᵃ)

    new_properties = Tuple(arch_array(new_arch, p) for p in old_properties)

    TX, TY, TZ = topology(old_grid)

    return LatitudeLongitudeGrid{TX, TY, TZ}(new_arch,
                                             old_grid.Nx, old_grid.Ny, old_grid.Nz,
                                             old_grid.Hx, old_grid.Hy, old_grid.Hz,
                                             old_grid.Lx, old_grid.Ly, old_grid.Lz,
                                             new_properties...,
                                             old_grid.radius)
end

function Adapt.adapt_structure(to, grid::LatitudeLongitudeGrid)
    TX, TY, TZ = topology(grid)
    return LatitudeLongitudeGrid{TX, TY, TZ}(nothing,
                                             grid.Nx, grid.Ny, grid.Nz,
                                             grid.Hx, grid.Hy, grid.Hz,
                                             grid.Lx, grid.Ly, grid.Lz,
                                             Adapt.adapt(to, grid.Δλᶠᵃᵃ),
                                             Adapt.adapt(to, grid.Δλᶜᵃᵃ),
                                             Adapt.adapt(to, grid.λᶠᵃᵃ),
                                             Adapt.adapt(to, grid.λᶜᵃᵃ),
                                             Adapt.adapt(to, grid.Δφᵃᶠᵃ),
                                             Adapt.adapt(to, grid.Δφᵃᶜᵃ),
                                             Adapt.adapt(to, grid.φᵃᶠᵃ),
                                             Adapt.adapt(to, grid.φᵃᶜᵃ),
                                             Adapt.adapt(to, grid.Δzᵃᵃᶠ),
                                             Adapt.adapt(to, grid.Δzᵃᵃᶜ),
                                             Adapt.adapt(to, grid.zᵃᵃᶠ),
                                             Adapt.adapt(to, grid.zᵃᵃᶜ),
                                             Adapt.adapt(to, grid.Δxᶠᶜᵃ),
                                             Adapt.adapt(to, grid.Δxᶜᶠᵃ),
                                             Adapt.adapt(to, grid.Δxᶠᶠᵃ),
                                             Adapt.adapt(to, grid.Δxᶜᶜᵃ),
                                             Adapt.adapt(to, grid.Δyᶠᶜᵃ),
                                             Adapt.adapt(to, grid.Δyᶜᶠᵃ),
                                             Adapt.adapt(to, grid.Azᶠᶜᵃ),
                                             Adapt.adapt(to, grid.Azᶜᶠᵃ),
                                             Adapt.adapt(to, grid.Azᶠᶠᵃ),
                                             Adapt.adapt(to, grid.Azᶜᶜᵃ),
                                             grid.radius)
end

#####
##### On-the-fly computation of LatitudeLongitudeGrid metrics
#####

@inline hack_cosd(φ) = cos(π * φ / 180)
@inline hack_sind(φ) = sin(π * φ / 180)

@inline Δxᶠᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶜᵃ[j]) * deg2rad(grid.Δλᶠᵃᵃ[i])
@inline Δxᶜᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶠᵃ[j]) * deg2rad(grid.Δλᶜᵃᵃ[i])
@inline Δxᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶠᵃ[j]) * deg2rad(grid.Δλᶠᵃᵃ[i])
@inline Δxᶜᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶜᵃ[j]) * deg2rad(grid.Δλᶜᵃᵃ[i])
@inline Δyᶜᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶠᵃ[j])
@inline Δyᶠᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶜᵃ[j])
@inline Azᶠᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶜᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶜᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))

@inline Δxᶠᶜᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶜᵃ[j]) * deg2rad(grid.Δλᶠᵃᵃ)
@inline Δxᶜᶠᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶠᵃ[j]) * deg2rad(grid.Δλᶜᵃᵃ)
@inline Δxᶠᶠᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶠᵃ[j]) * deg2rad(grid.Δλᶠᵃᵃ)
@inline Δxᶜᶜᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius * hack_cosd(grid.φᵃᶜᵃ[j]) * deg2rad(grid.Δλᶜᵃᵃ)
@inline Δyᶜᶠᵃ(i, j, k, grid::YRegLatLonGrid) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶠᵃ)
@inline Δyᶠᶜᵃ(i, j, k, grid::YRegLatLonGrid) = @inbounds grid.radius * deg2rad(grid.Δφᵃᶜᵃ)
@inline Azᶠᶜᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶜᶠᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶠᶠᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶜᶜᵃ(i, j, k, grid::XRegLatLonGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))

#####
##### Utilities to precompute metrics
#####

@inline metrics_precomputed(::LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, Nothing}) = false
@inline metrics_precomputed(::LatitudeLongitudeGrid) = true

#####
##### Kernels that precompute the z- and x-metric
#####

@inline metric_worksize(grid::LatitudeLongitudeGrid)  = (length(grid.Δλᶜᵃᵃ), length(grid.φᵃᶜᵃ) - 1)
@inline metric_workgroup(grid::LatitudeLongitudeGrid) = (16, 16)

@inline metric_worksize(grid::XRegLatLonGrid)  =  length(grid.φᵃᶜᵃ) - 1
@inline metric_workgroup(grid::XRegLatLonGrid) =  16

function precompute_curvilinear_metrics!(grid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ)

    arch = grid.architecture

    workgroup, worksize  = metric_workgroup(grid), metric_worksize(grid)
    curvilinear_metrics! = precompute_metrics_kernel!(Architectures.device(arch), workgroup, worksize)

    event = curvilinear_metrics!(grid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ; dependencies=device_event(arch))
    wait(event)

    return nothing
end

@kernel function precompute_metrics_kernel!(grid::LatitudeLongitudeGrid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ)
    i, j = @index(Global, NTuple)

    # Manually offset x- and y-index
    i′ = i + grid.Δλᶜᵃᵃ.offsets[1]
    j′ = j + grid.φᵃᶜᵃ.offsets[1] + 1

    @inbounds begin
        Δxᶠᶜ[i′, j′] = Δxᶠᶜᵃ(i′, j′, 1, grid)
        Δxᶜᶠ[i′, j′] = Δxᶜᶠᵃ(i′, j′, 1, grid)
        Δxᶠᶠ[i′, j′] = Δxᶠᶠᵃ(i′, j′, 1, grid)
        Δxᶜᶜ[i′, j′] = Δxᶜᶜᵃ(i′, j′, 1, grid)
        Azᶠᶜ[i′, j′] = Azᶠᶜᵃ(i′, j′, 1, grid)
        Azᶜᶠ[i′, j′] = Azᶜᶠᵃ(i′, j′, 1, grid)
        Azᶠᶠ[i′, j′] = Azᶠᶠᵃ(i′, j′, 1, grid)
        Azᶜᶜ[i′, j′] = Azᶜᶜᵃ(i′, j′, 1, grid)
    end
end

@kernel function precompute_metrics_kernel!(grid::XRegLatLonGrid, Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ)
    j = @index(Global, Linear)

    # Manually offset y-index
    j′ = j + grid.φᵃᶜᵃ.offsets[1] + 1

    @inbounds begin
        Δxᶠᶜ[j′] = Δxᶠᶜᵃ(1, j′, 1, grid)
        Δxᶜᶠ[j′] = Δxᶜᶠᵃ(1, j′, 1, grid)
        Δxᶠᶠ[j′] = Δxᶠᶠᵃ(1, j′, 1, grid)
        Δxᶜᶜ[j′] = Δxᶜᶜᵃ(1, j′, 1, grid)
        Azᶠᶜ[j′] = Azᶠᶜᵃ(1, j′, 1, grid)
        Azᶜᶠ[j′] = Azᶜᶠᵃ(1, j′, 1, grid)
        Azᶠᶠ[j′] = Azᶠᶠᵃ(1, j′, 1, grid)
        Azᶜᶜ[j′] = Azᶜᶜᵃ(1, j′, 1, grid)
    end
end

#####
##### Kernels that precompute the y-metric
#####

function precompute_Δy_metrics(grid::LatitudeLongitudeGrid, Δyᶠᶜ, Δyᶜᶠ)
    arch = grid.architecture
    precompute_Δy! = precompute_Δy_kernel!(Architectures.device(arch), 16, length(grid.Δφᵃᶜᵃ) - 1)
    event = precompute_Δy!(grid, Δyᶠᶜ, Δyᶜᶠ; dependencies=device_event(arch))

    wait(event)
    return Δyᶠᶜ, Δyᶜᶠ
end

function  precompute_Δy_metrics(grid::YRegLatLonGrid, Δyᶠᶜ, Δyᶜᶠ)
    Δyᶜᶠ =  Δyᶜᶠᵃ(1, 1, 1, grid)
    Δyᶠᶜ =  Δyᶠᶜᵃ(1, 1, 1, grid)
    return Δyᶠᶜ, Δyᶜᶠ
end

@kernel function precompute_Δy_kernel!(grid, Δyᶠᶜ, Δyᶜᶠ)
    j  = @index(Global, Linear)

    # Manually offset y-index
    j′ = j + grid.Δφᵃᶜᵃ.offsets[1] + 1

    @inbounds begin
        Δyᶜᶠ[j′] = Δyᶜᶠᵃ(1, j′, 1, grid)
        Δyᶠᶜ[j′] = Δyᶜᶠᵃ(1, j′, 1, grid)
    end
end

#####
##### Metric memory allocation
#####

function allocate_metrics(grid::LatitudeLongitudeGrid)
    FT = eltype(grid)

    # preallocate quantities to ensure correct type and size
    grid_metrics = (:Δxᶠᶜ,
                    :Δxᶜᶠ,
                    :Δxᶠᶠ,
                    :Δxᶜᶜ,
                    :Azᶠᶜ,
                    :Azᶜᶠ,
                    :Azᶠᶠ,
                    :Azᶜᶜ)

    arch = grid.architecture

    if typeof(grid) <: XRegLatLonGrid
        offsets     = grid.φᵃᶜᵃ.offsets[1]
        metric_size = length(grid.φᵃᶜᵃ)
    else
        offsets     = (grid.Δλᶜᵃᵃ.offsets[1], grid.φᵃᶜᵃ.offsets[1])
        metric_size = (length(grid.Δλᶜᵃᵃ)   , length(grid.φᵃᶜᵃ))
    end

    for metric in grid_metrics
        parentM        = Symbol(metric, :_parent)
        @eval $parentM = zeros($FT, $metric_size...)
        @eval $metric  = OffsetArray(arch_array($arch, $parentM), $offsets...)
    end

    if typeof(grid) <: YRegLatLonGrid
        Δyᶠᶜ = FT(0.0)
        Δyᶜᶠ = FT(0.0)
    else
        parentC = zeros(FT, length(grid.Δφᵃᶜᵃ))
        parentF = zeros(FT, length(grid.Δφᵃᶜᵃ))
        Δyᶠᶜ    = OffsetArray(arch_array(arch, parentC), grid.Δφᵃᶜᵃ.offsets[1])
        Δyᶜᶠ    = OffsetArray(arch_array(arch, parentF), grid.Δφᵃᶜᵃ.offsets[1])
    end

    return Δxᶠᶜ, Δxᶜᶠ, Δxᶠᶠ, Δxᶜᶜ, Δyᶠᶜ, Δyᶜᶠ, Azᶠᶜ, Azᶜᶠ, Azᶠᶠ, Azᶜᶜ
end

#####
##### Get minima of grid
#####

function min_Δx(grid::LatitudeLongitudeGrid)
    topo = topology(grid)
    if topo[1] == Flat
        return Inf
    else
        ϕᵃᶜᵃ_max = maximum(abs, ynodes(Center, grid))
        return grid.radius * cosd(ϕᵃᶜᵃ_max) * deg2rad(min_number_or_array(grid.Δλᶜᵃᵃ))
    end
end

function min_Δy(grid::LatitudeLongitudeGrid)
    topo = topology(grid)
    if topo[2] == Flat
        return Inf
    else
        return grid.radius * deg2rad(min_number_or_array(grid.Δφᵃᶜᵃ))
    end
end

function min_Δz(grid::LatitudeLongitudeGrid)
    topo = topology(grid)
    if topo[3] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δzᵃᵃᶜ)
    end
end


text:
module ApproxFunSingularities
using Base, LinearAlgebra, Reexport, IntervalSets, DomainSets, Statistics

@reexport using ApproxFunBase
@reexport using ApproxFunOrthogonalPolynomials

import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSpace, VFun, RowVector,
                    UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                    HeavisideSpace, PointSpace,
                    IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
                    AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer,
                    AbstractTransformPlan, TransformPlan, ITransformPlan,
                    ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral,
                    ConcreteVolterra, Volterra, VolterraWrapper,
                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                    Conversion, defaultConversion, defaultcoefficients, default_Fun, Multiplication, Derivative, Integral, bandwidths,
                    ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                    DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral, IntegralWrapper,
                    ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven, Dirichlet, ConcreteDirichlet,
                    TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                    hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible,
                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform,
                    CanonicalTransformPlan, ICanonicalTransformPlan,
                    Integral,
                    domainspace, rangespace, boundary,
                    union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points,
                    rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp,
                    linesum, differentiate, integrate, linebilinearform, bilinearform,
                    UnsetNumber, coefficienttimes, subspace_coefficients, sumspacecoefficients, specialfunctionnormalizationpoint,
                    Segment, IntervalOrSegmentDomain, PiecewiseSegment, isambiguous, Vec, eps, isperiodic,
                    arclength, complexlength,
                    invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                    reverseorientation, checkpoints, evaluate, extrapolate, mul_coefficients, coefficients, isconvertible,
                    clenshaw, ClenshawPlan, sineshaw,
                    toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex,
                    SpaceOperator, ZeroOperator, InterlaceOperator,
                    interlace!, reverseeven!, negateeven!, cfstype, pad!, alternatesign!, mobius,
                    extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots,splitatroots,
                    chebmult_getindex, intpow, alternatingsum,
                    domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!, dotu,
                    components, promoterangespace, promotedomainspace,
                    block, blockstart, blockstop, blocklengths, isblockbanded, pointscompatible, affine_setdiff, complexroots,
                    ∞, gamma

import ApproxFunOrthogonalPolynomials: order

import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
            Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,
            dimension

import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                >=, /, ^, \, ∪, transpose, size, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
                similar, map, vcat, hcat, hvcat, show, summary, stride, sum, cumsum, sign, imag, conj, inv,
                complex, reverse, exp, sqrt, abs, abs2, sign, issubset, values, in, first, last, rand, intersect, setdiff,
                isless, union, angle, join, isnan, isapprox, isempty, sort, merge, promote_rule,
                minimum, maximum, extrema, argmax, argmin, findmax, findmin, isfinite,
                zeros, zero, one, promote_rule, repeat, length, resize!, isinf,
                getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
                split, exp, log

import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot, cross,
                qr, qr!, rank, isdiag, istril, istriu, issymmetric, ishermitian,
                Tridiagonal, diagm, diagm_container, factorize, nullspace,
                Hermitian, Symmetric, adjoint, transpose, char_uplo

include("divide_singularity.jl")
include("JacobiWeight.jl")
include("JacobiWeightOperators.jl")
include("JacobiWeightChebyshev.jl")
include("LogWeight.jl")
include("ExpWeight.jl")

/(c::Number,f::Fun{Ultraspherical{λ,DD,RR}}) where {λ,DD,RR} = c/Fun(f,Chebyshev(domain(f)))
/(c::Number,f::Fun{<:PolynomialSpace{<:IntervalOrSegment}}) = c/Fun(f,Chebyshev(domain(f)))
/(c::Number,f::Fun{<:Chebyshev}) = setdomain(c/setcanonicaldomain(f),domain(f))

scaleshiftdomain(f::Fun,sc,sh) = setdomain(f,sc .* domain(f) .+ sh)

function /(c::Number,f::Fun{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR}
    fc = setcanonicaldomain(f)
    d=domain(f)
    # if domain f is small then the pts get projected in
    tol = 200eps(promote_type(typeof(c),cfstype(f)))*norm(f.coefficients,1)

    # we prune out roots at the boundary first
    if ncoefficients(f) == 1
        return Fun(c/f.coefficients[1],space(f))
    elseif ncoefficients(f) == 2
        if isempty(roots(f))
            return \(Multiplication(f,space(f)),c;tolerance=0.05tol)
        elseif isapprox(fc.coefficients[1],fc.coefficients[2])
            # we check directly for f*(1+x)
            return Fun(JacobiWeight(-1,0,space(f)),[c/fc.coefficients[1]])
        elseif isapprox(fc.coefficients[1],-fc.coefficients[2])
            # we check directly for f*(1-x)
            return Fun(JacobiWeight(0,-1,space(f)),[c/fc.coefficients[1]])
        else
            # we need to split at the only root
            return c/splitatroots(f)
        end
    elseif abs(first(fc)) ≤ tol
        #left root
        g=divide_singularity((1,0),fc)
        p=c/g
        x = Fun(identity,domain(p))
        return scaleshiftdomain(p/(1+x),complexlength(d)/2,mean(d) )
    elseif abs(last(fc)) ≤ tol
        #right root
        g=divide_singularity((0,1),fc)
        p=c/g
        x=Fun(identity,domain(p))
        return scaleshiftdomain(p/(1-x),complexlength(d)/2,mean(d) )
    else
        r = roots(fc)

        if length(r) == 0
            return \(Multiplication(f,space(f)),c;tolerance=0.05tol)
        elseif abs(last(r)+1.0)≤tol  # double check
            #left root
            g=divide_singularity((1,0),fc)
            p=c/g
            x=Fun(identity,domain(p))
            return scaleshiftdomain(p/(1+x),complexlength(d)/2,mean(d) )
        elseif abs(last(r)-1.0)≤tol  # double check
            #right root
            g=divide_singularity((0,1),fc)
            p=c/g
            x=Fun(identity,domain(p))
            return scaleshiftdomain(p/(1-x),complexlength(d)/2,mean(d) )
        else
            # no roots on the boundary
            return c/splitatroots(f)
        end
    end
end

function ^(f::Fun{<:PolynomialSpace}, k::Real)
    T = cfstype(f)
    RT = real(T)
    # Need to think what to do if this is ever not the case..
    sp = space(f)
    fc = setcanonicaldomain(f) #Project to interval
    csp = space(fc)

    r = sort(roots(fc))
    #TODO divideatroots
    @assert length(r) <= 2

    if length(r) == 0
        setdomain(Fun((x->x^k) ∘ fc,csp),domain(f))  # using ∘ supports fast transforms for fc
    elseif length(r) == 1
        @assert isapprox(abs(r[1]),1)

        if isapprox(r[1], 1)
            Fun(JacobiWeight(zero(RT),k,sp),coefficients(divide_singularity(true,fc)^k,csp))
        else
            Fun(JacobiWeight(k,zero(RT),sp),coefficients(divide_singularity(false,fc)^k,csp))
        end
    else
        @assert isapprox(r[1],-1)
        @assert isapprox(r[2],1)

        Fun(JacobiWeight(k,k,sp),coefficients(divide_singularity(fc)^k,csp))
    end
end

# function log{MS<:MappedSpace}(f::Fun{MS})
#     g=log(Fun(f.coefficients,space(f).space))
#     Fun(g.coefficients,MappedSpace(domain(f),space(g)))
# end

# project first to [-1,1] to avoid issues with
# complex derivative
function log(f::Fun{<:PolynomialSpace{<:ChebyshevInterval}})
    r = sort(roots(f))
    #TODO divideatroots
    @assert length(r) <= 2

    if length(r) == 0
        cumsum(differentiate(f)/f)+log(first(f))
    elseif length(r) == 1
        @assert isapprox(abs(r[1]),1)

        if isapprox(r[1],1.)
            g=divide_singularity(true,f)
            lg=Fun(LogWeight(0.,1.,Chebyshev()),[1.])
            if isapprox(g,1.)  # this means log(g)~0
                lg
            else # log((1-x)) + log(g)
                lg⊕log(g)
            end
        else
            g=divide_singularity(false,f)
            lg=Fun(LogWeight(1.,0.,Chebyshev()),[1.])
            if isapprox(g,1.)  # this means log(g)~0
                lg
            else # log((1+x)) + log(g)
                lg⊕log(g)
            end
       end
    else
        @assert isapprox(r[1],-1)
        @assert isapprox(r[2],1)

        g=divide_singularity(f)
        lg=Fun(LogWeight(1.,1.,Chebyshev()),[1.])
        if isapprox(g,1.)  # this means log(g)~0
            lg
        else # log((1+x)) + log(g)
            lg⊕log(g)
        end
    end
end

# JacobiWeight explodes, we want to ensure the solution incorporates the fact
# that exp decays rapidly
exp(f::Fun{<:JacobiWeight}) = setdomain(exp(setdomain(f, ChebyshevInterval())), domain(f))
function exp(f::Fun{<:JacobiWeight{<:Any,<:ChebyshevInterval}})
    S=space(f)
    q=Fun(S.space,f.coefficients)
    if isapprox(S.α,0.) && isapprox(S.β,0.)
        exp(q)
    elseif S.β < 0 && isapprox(first(q),0.)
        # this case can remove the exponential decay
        exp(Fun(f,JacobiWeight(S.β+1,S.α,S.space)))
    elseif S.α < 0 && isapprox(last(q),0.)
        exp(Fun(f,JacobiWeight(S.β,S.α+1,S.space)))
    elseif S.β > 0 && isapproxinteger(S.β)
        exp(Fun(f,JacobiWeight(0.,S.α,S.space)))
    elseif S.α > 0 && isapproxinteger(S.α)
        exp(Fun(f,JacobiWeight(S.β,0.,S.space)))
    else
        #find normalization point
        xmax,opfxmax,opmax=specialfunctionnormalizationpoint(exp,real,f)

        if S.α < 0 && S.β < 0
            # provided both are negative, we get exponential decay on both ends
            @assert real(first(q)) < 0 && real(last(q)) < 0
            s=JacobiWeight(2.,2.,domain(f))
        elseif S.β < 0 && isapprox(S.α,0.)
            @assert real(first(q)) < 0
            s=JacobiWeight(2.,0.,domain(f))
        elseif S.α < 0 && isapprox(S.β,0.)
            @assert real(last(q)) < 0
            s=JacobiWeight(0.,2.,domain(f))
        else
            error("exponential of fractional power, not implemented")
        end

        D=Derivative(s)
        B=Evaluation(s,xmax)

        \([B,D-f'],Any[opfxmax,0.];tolerance=eps(cfstype(f))*opmax)
    end
end

## Root finding for JacobiWeight expansion

# Add endpoints for JacobiWeight
# TODO: what about cancellation?
function roots(f::Fun{S,T}) where {S<:JacobiWeight,T}
    sp=space(f)
    d=domain(sp)
    rts=roots(Fun(sp.space,f.coefficients))
    if sp.β > 0
        rts=[leftendpoint(d);rts]
    end
    if sp.α > 0
        rts=[rts;rightendpoint(d)]
    end
    rts
end

for Func in (:DefiniteIntegral,:DefiniteLineIntegral)
    @eval begin
        #TODO: this may be misleading
        $Func(d::IntervalOrSegment) = $Func(JacobiWeight(-0.5,-0.5,Chebyshev(d)))
        function $Func(α::Number,β::Number,d::IntervalOrSegment)
            @assert α == β
            @assert round(Int,α+.5) == α+.5
            @assert round(Int,α+.5) >= 0
            $Func(JacobiWeight(α,β,Ultraspherical(round(Int,α+.5),d)))
        end
        $Func(α::Number,β::Number) = $Func(α,β,ChebyshevInterval())
    end
end

## Integration
function integrate(f::Fun{LW}) where LW <: LaguerreWeight{LL} where LL <: Laguerre
    α = space(f).α;
    n = length(f.coefficients);
    if space(f).space.α != α
        throw(ArgumentError("`integrate` is applicable if only LaguerreWeight parameter and Laguerre parameter are equal."))
    else
        if n == 0
            Fun(0)
        else
            if !isinteger(α) || α == 2
                if n == 1
                    f̃ = Fun(f, JacobiWeight(α, 0, Chebyshev(Ray())));
                    g = integrate(f̃);
                    g = g - last(g)
                else
                    if f.coefficients[1] == 0
                        Fun(WeightedLaguerre(α + 1), f.coefficients[2:end] ./ (1:n-1))
                    else
                        f₀ = Fun(WeightedLaguerre(α), [f.coefficients[1]]);
                        f₀ = Fun(f₀, JacobiWeight(α, 0, Chebyshev(Ray())));
                        g₀ = integrate(f₀);
                        g₀ = g₀ - last(g₀);
                        g = Fun(WeightedLaguerre(α + 1), f.coefficients[2:end] ./ (1:n-1));
                        g₀ + g
                    end
                end
            else
                if n == 1
                    f̃ = Fun(f, Chebyshev(Ray()));
                    g = integrate(f̃);
                    g = g - last(g)
                else
                    if f.coefficients[1] == 0
                        Fun(WeightedLaguerre(α + 1), f.coefficients[2:end] ./ (1:n-1))
                    else
                        f₀ = Fun(WeightedLaguerre(α), [f.coefficients[1]]);
                        f₀ = Fun(f₀, Chebyshev(Ray()));
                        g₀ = integrate(f₀);
                        g₀ = g₀ - last(g₀);
                        g = Fun(WeightedLaguerre(α + 1), f.coefficients[2:end] ./ (1:n-1));
                        g₀ + g
                    end
                end
            end
        end
    end
end

end


text:
using DataFrames
using CSV
using Chain
using Antenna
using Interpolations

print("\ec")
extract = (filename) -> @chain CSV.File(filename)  DataFrame _[181:end, 2:end] Matrix circshift(_, (0, -180))
extract_theta_phi = (filename) -> @chain CSV.File(filename)  DataFrame _[181:end, 2:end] Matrix circshift(_, (0, -180))
# theta = @chain CSV.File("test/elem/mag_theta.csv")  DataFrame _[181:end,1] deg2rad.(_)
theta = collect(0.:180.) .|> deg2rad
phi = collect(-180.:180.) .|> deg2rad

mag_theta = extract("test/elem/mag_theta.csv")
mag_phi = extract("test/elem/mag_phi.csv")
ang_theta = extract("test/elem/ang_theta.csv")
ang_phi = extract("test/elem/ang_phi.csv")

rE_theta = [mag*exp(-1im*deg2rad(ang)) for (mag,ang) in zip(mag_theta, ang_theta)]
rE_phi = [mag*exp(-1im*deg2rad(ang)) for (mag,ang) in zip(mag_phi, ang_phi)]

anten_pattern(
    θ=LinearInterpolation((theta, phi), rE_theta),
    ϕ=LinearInterpolation((theta, phi), rE_phi)
)


text:
#=
The following code is modified from
https://github.com/billmclean/GaussQuadrature.jl with the original license:

> The MIT License (MIT)

> Copyright (c) 2013 billmclean

> Permission is hereby granted, free of charge, to any person obtaining a copy
> of this software and associated documentation files (the "Software"), to deal
> in the Software without restriction, including without limitation the rights
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> copies of the Software, and to permit persons to whom the Software is
> furnished to do so, subject to the following conditions:

> The above copyright notice and this permission notice shall be included in
> all copies or substantial portions of the Software.

> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> SOFTWARE.
=#

# October 2013 by Bill McLean, School of Maths and Stats,
# The University of New South Wales.
#
# Based on earlier Fortran codes
#
# gaussq.f original version 20 Jan 1975 from Stanford
# gaussq.f modified 21 Dec by Eric Grosse
# gaussquad.f95 Nov 2005 by Bill Mclean
#
# This module provides functions to compute the abscissae x[j] and
# weights w[j] for the classical Gauss quadrature rules, including
# the Radau and Lobatto variants.  Thus, the sum
#
#           n
#           ∑  w[j] f(x[j])
#          j=1
#
# approximates
#
#           hi
#          ∫  f(x) w(x) dx
#          lo
#
# where the weight function w(x) and interval lo < x < hi are as shown
# in the table below.
#
# Name                      Interval     Weight Function
#
# Legendre                 -1 < x < 1          1
# Chebyshev (first kind)   -1 < x < 1     1 / sqrt(1-x²)
# Chebyshev (second kind)  -1 < x < 1       sqrt(1-x²)
# Jacobi                   -1 < x < 1     (1-x)ᵅ (1+x)ᵝ
# Laguerre                  0 < x < ∞     xᵅ exp(-x)
# Hermite                  -∞ < x < ∞      exp(-x²)
#
# In addition to these classical rules, the module generates Gauss rules
# for logarithmic weights of the form
#
#    w(x) = x^ρ log(1/x)   for 0 < x < 1.
#
# For the Jacobi and Laguerre rules we require α > -1 and
# β > -1, so that the weight function is integrable.  Likewise, for
# log weight we require ρ > -1.
#
# Use the endpoint argument to include one or both of the end points
# of the interval of integration as an abscissa in the quadrature
# rule, as follows.
#
# endpoint = NeitherEndPoint()   Default      lo < x[j] < hi, j = 1:n.
# endpoint = LeftEndPoint()      Left Radau   lo = x[1] < x[j] < hi, j = 2:n.
# endpoint = RightEndPoint()     Right Radau  lo < x[j] < x[n] = hi, j = 1:n-1.
# endpoint = BothEndPoint()      Lobatto      lo = x[1] < x[j] < x[n] = hi, j = 2:n-1.
#
#
# The code uses the Golub and Welsch algorithm, in which the abscissae
# x[j] are the eigenvalues of a symmetric tridiagonal matrix whose
# entries depend on the coefficients in the 3-term recurrence relation
# for the othonormal polynomials generated by the weighted inner product.
#
# References:
#
#   1.  Golub, G. H., and Welsch, J. H., Calculation of Gaussian
#       quadrature rules, Mathematics of Computation 23 (April,
#       1969), pp. 221-230.
#   2.  Golub, G. H., Some modified matrix eigenvalue problems,
#       Siam Review 15 (april, 1973), pp. 318-334 (section 7).
#   3.  Stroud and Secrest, Gaussian Quadrature Formulas, Prentice-
#       Hall, Englewood Cliffs, N.J., 1966.

# Enumeration type used to specify which endpoints of the integration
# interval should be included amongst the quadrature points: neither,
# left, right, or both.

abstract type EndPoint end
struct NeitherEndPoint <: EndPoint end
struct LeftEndPoint <: EndPoint end
struct RightEndPoint <: EndPoint end
struct BothEndPoint <: EndPoint end

"""
    x, w = legendregauss(T, n, endpoint::EndPoint=Bennu.NeitherEndPoint())

Returns points `x` and weights `w` for the `n`-point Gauss-Legendre rule
for the interval `-1 < x < 1` with weight function `w(x) = 1`.

Use `endpoint=LeftEndPoint()`, `RightEndPoint() ` or `BothEndPoints()` for the
left Radau, right Radau, or Lobatto rules, respectively.
"""
function legendregauss(::Type{T}, n, endpoint=NeitherEndPoint()) where {T}
    @assert n ≥ 1
    a, b = legendrecoefficients(T, n)
    return gaussrule(-one(T), one(T), a, b, endpoint)
end

"""
    x, w = legendregauss(n, endpoint::EndPoint=Bennu.NeitherEndPoint())

Convenience function with type `T = Float64`:
"""
legendregauss(n, endpoint::EndPoint=NeitherEndPoint()) = legendregauss(Float64, n, endpoint)

"""
    x, w = legendregausslobatto(T, n)

Returns points `x` and weights `w` for the `n`-point Legendre-Gauss-Lobatto rule
for the interval `-1 ≤ x ≤ 1` with weight function `w(x) = 1`.
"""
function legendregausslobatto(::Type{T}, n) where {T}
    return legendregauss(T, n, BothEndPoint())
end

"""
    x, w = legendregausslobatto(n)

Convenience function with type `T = Float64`:
"""
legendregausslobatto(n) = legendregausslobatto(Float64, n)

function legendrecoefficients(::Type{T}, n) where {T}
    a = zeros(T, n)
    b = zeros(T, n + 1)
    b[1] = sqrt(convert(T, 2))
    for k in 2:(n + 1)
        b[k] = (k - 1) / sqrt(convert(T, (2k - 1) * (2k - 3)))
    end
    return a, b
end

"""
    x, w = gaussrule(lo, hi, a, b, endpoint, maxiterations=100)

Generates the points `x` and weights `w` for a Gauss rule with weight
function `w(x)` on the interval `lo < x < hi`.

The arrays `a` and `b` hold the coefficients (as given, for instance, by
`legendrecoefficients`) in the three-term recurrence relation for the monic
orthogonal polynomials `p(0,x)`, `p(1,x)`, `p(2,x)`, ... , that is,

    p(k, x) = (x-a[k]) p(k-1, x) - b[k]² p(k-2, x),    k ≥ 1,

where `p(0, x) = 1` and, by convention, `p(-1, x) = 0` with

              hi
    b[1]^2 = ∫  w(x) dx.
             lo

Thus, `p(k, x) = xᵏ + lower degree terms` and

     hi
    ∫  p(j, x) p(k, x) w(x) dx = 0 if j ≠ k.
    lo
"""
function gaussrule(lo, hi, a, b, endpoint::EndPoint, maxiterations=100)
    T = promote_type(typeof(lo), typeof(hi), eltype(a), eltype(b))
    n = length(a)
    @assert length(b) == n + 1
    if endpoint isa LeftEndPoint
        if n == 1
            a[1] = lo
        else
            a[n] = tridiagonalshiftsolve(n, lo, a, b) * b[n]^2 + lo
        end
    elseif endpoint isa RightEndPoint
        if n == 1
            a[1] = hi
        else
            a[n] = tridiagonalshiftsolve(n, hi, a, b) * b[n]^2 + hi
        end
    elseif endpoint isa BothEndPoint
        if n == 1
            error("Must have at least two points for both ends.")
        end
        g = tridiagonalshiftsolve(n, lo, a, b)
        t1 = (hi - lo) / (g - tridiagonalshiftsolve(n, hi, a, b))
        b[n] = sqrt(t1)
        a[n] = lo + g * t1
    end
    w = zero(a)
    tridiagonaleigenproblem!(a, b, w, maxiterations)
    for i in 1:n
        w[i] = (b[1] * w[i])^2
    end
    idx = sortperm(a)
    # Ensure end point values are exact.
    if endpoint isa LeftEndPoint || endpoint isa BothEndPoint
        a[idx[1]] = lo
    elseif endpoint isa RightEndPoint || endpoint isa BothEndPoint
        a[idx[n]] = hi
    end
    return a[idx], w[idx]
end

function tridiagonalshiftsolve(n, shift, a, b) where {T}
    #
    # Perform elimination to find the nth component s = delta[n]
    # of the solution to the nxn linear system
    #
    #     ( J_n - shift I_n ) delta = e_n,
    #
    # where J_n is the symmetric tridiagonal matrix with diagonal
    # entries a[i] and off-diagonal entries b[i], and e_n is the nth
    # standard basis vector.
    #
    t = a[1] - shift
    for i in 2:(n - 1)
        t = a[i] - shift - b[i]^2 / t
    end
    return one(t) / t
end

function tridiagonaleigenproblem!(d, e, z, maxiterations)
    #
    # Finds the eigenvalues and first components of the normalised
    # eigenvectors of a symmetric tridiagonal matrix by the implicit
    # QL method.
    #
    # d[i]   On entry, holds the ith diagonal entry of the matrix.
    #        On exit, holds the ith eigenvalue.
    #
    # e[i]   On entry, holds the [i,i-1] entry of the matrix for
    #        i = 2, 3, ..., n.  (The value of e[1] is not used.)
    #        On exit, e is overwritten.
    #
    # z[i]   On exit, holds the first component of the ith normalised
    #        eigenvector associated with d[i].
    #
    # maxiterations The maximum number of QL iterations.
    #
    # Martin and Wilkinson, Numer. Math. 12: 377-383 (1968).
    # Dubrulle, Numer. Math. 15: 450 (1970).
    # Handbook for Automatic Computation, Vol ii, Linear Algebra,
    #        pp. 241-248, 1971.
    #
    # This is a modified version of the Eispack routine imtql2.
    #
    T = promote_type(eltype(d), eltype(e), eltype(z))
    n = length(z)
    z[1] = one(T)
    z[2:n] .= zero(T)
    e[n + 1] = zero(T)

    if n == 1 # Nothing to do for a 1x1 matrix.
        return
    end
    for l in 1:n
        for j in 1:maxiterations
            # Look for small off-diagonal elements.
            m = n
            for i in l:(n - 1)
                if abs(e[i + 1]) <= eps(T) * (abs(d[i]) + abs(d[i + 1]))
                    m = i
                    break
                end
            end
            p = d[l]
            if m == l
                continue
            end
            if j == maxiterations
                msg = string("No convergence after ", j, " iterations",
                             " (try increasing maxiterations)")
                error(msg)
            end
            # Form shift
            g = (d[l + 1] - p) / (2 * e[l + 1])
            r = hypot(g, one(T))
            g = d[m] - p + e[l + 1] / (g + copysign(r, g))
            s = one(T)
            c = one(T)
            p = zero(T)
            for i in (m - 1):-1:l
                f = s * e[i + 1]
                b = c * e[i + 1]
                if abs(f) < abs(g)
                    s = f / g
                    r = hypot(s, one(T))
                    e[i + 2] = g * r
                    c = one(T) / r
                    s *= c
                else
                    c = g / f
                    r = hypot(c, one(T))
                    e[i + 2] = f * r
                    s = one(T) / r
                    c *= s
                end
                g = d[i + 1] - p
                r = (d[i] - g) * s + 2 * c * b
                p = s * r
                d[i + 1] = g + p
                g = c * r - b
                # Form first component of vector.
                f = z[i + 1]
                z[i + 1] = s * z[i] + c * f
                z[i] = c * z[i] - s * f
            end
            d[l] -= p
            e[l + 1] = g
            e[m + 1] = zero(T)
        end
    end
    return
end