text:
using NOAAData
using Base.Test

# write your own tests here
TOKEN = "QGdQVDottDzkuRlaZbtAtRRRWyNrBZiE"
stationid = "GHCND:USW00094728"
noaa = NOAA(TOKEN)

startdate = Date(2015, 1, 1)
enddate = Date(2015, 12, 31)

results = get(GHCND(), noaa, startdate, enddate, stationid)
@test results.count == 3873

df = DataFrame(results)
@test size(df) == (365, 14)
@test typeof(df) == DataFrames.DataFrame
@test df[:TMAX][1] == 3.9

it = NDSparse(results)
@test length(it) == 365
@test typeof(it) == IndexedTables.NDSparse{NamedTuples._NT_AWND_FMTM_PGTM_PRCP_SNOW_SNWD_TMAX_TMIN_WDF2_WDF5_WSF2_WSF5_WT{Float64,DateTime,DateTime,Float64,Float64,
                    Float64,Float64,Float64,Float64,Float64,Float64,Float64,String},Tuple{Date},IndexedTables.Columns{NamedTuples._NT_DATE{Date},
                    NamedTuples._NT_DATE{Array{Date,1}}},IndexedTables.Columns{NamedTuples._NT_AWND_FMTM_PGTM_PRCP_SNOW_SNWD_TMAX_TMIN_WDF2_WDF5_WSF2_WSF5_WT{
                        Float64,DateTime,DateTime,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,String},
                    NamedTuples._NT_AWND_FMTM_PGTM_PRCP_SNOW_SNWD_TMAX_TMIN_WDF2_WDF5_WSF2_WSF5_WT{Array{Float64,1},Array{DateTime,1},Array{DateTime,1},Array{Float64,1},
                    Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{String,1}}}}
@test it.data.columns[:TMAX][1] == 3.9

# aggregation with indexed tables
# max temp
maxtemp = IndexedTables.convertdim(map(IndexedTables.pick(:TMAX), it), 1, Dates.month, agg=max)
@test length(maxtemp) == 12
@test maxtemp[2] == 6.1

# mean high temp
meantemp = IndexedTables.convertdim(map(IndexedTables.pick(:TMAX), it), 1, Dates.month, vecagg=mean)
@test length(meantemp) == 12
@test meantemp[5] == 25.761290322580642

# let's try a new data type
startdate = Date(2011, 1, 1)
results = get(GSOM(), noaa, startdate, enddate, stationid)

@test results.count == 1739

df = DataFrame(results)
@test size(df) == (60, 39)

it = NDSparse(results)
@test length(it.data) == 60

# Test errors
startdate = Date(1985, 1, 1)
@test_throws ErrorException get(GSOM(), noaa, startdate, enddate, stationid)
@test_throws ErrorException get(GHCND(), noaa, startdate, enddate, stationid)


text:
using NashGames
using Documenter

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

makedocs(;
  modules = [NashGames],
  authors = "Tangi Migot tangi.migot@gmail.com",
  repo = "https://github.com/JuliaOptimizationVariationalAnalysis/NashGames.jl/blob/{commit}{path}#{line}",
  sitename = "NashGames.jl",
  format = Documenter.HTML(;
    prettyurls = get(ENV, "CI", "false") == "true",
    canonical = "https://JuliaOptimizationVariationalAnalysis.github.io/NashGames.jl",
    assets = String[],
  ),
  pages = ["Home" => "index.md"],
)

deploydocs(;
  repo = "github.com/JuliaOptimizationVariationalAnalysis/NashGames.jl",
  devbranch = "main",
)

ext: jl


text:
using Hadamard
using Base.Test

H8 = [   1     1     1     1     1     1     1     1
         1     1     1     1    -1    -1    -1    -1
         1     1    -1    -1    -1    -1     1     1
         1     1    -1    -1     1     1    -1    -1
         1    -1    -1     1     1    -1    -1     1
         1    -1    -1     1    -1     1     1    -1
         1    -1     1    -1    -1     1    -1     1
         1    -1     1    -1     1    -1     1    -1  ]

H8n = [  1     1     1     1     1     1     1     1
         1    -1     1    -1     1    -1     1    -1
         1     1    -1    -1     1     1    -1    -1
         1    -1    -1     1     1    -1    -1     1
         1     1     1     1    -1    -1    -1    -1
         1    -1     1    -1    -1     1    -1     1
         1     1    -1    -1    -1    -1     1     1
         1    -1    -1     1    -1     1     1    -1  ]

H8d = [  1     1     1     1     1     1     1     1
         1     1     1     1    -1    -1    -1    -1
         1     1    -1    -1     1     1    -1    -1
         1     1    -1    -1    -1    -1     1     1
         1    -1     1    -1     1    -1     1    -1
         1    -1     1    -1    -1     1    -1     1
         1    -1    -1     1     1    -1    -1     1
         1    -1    -1     1    -1     1     1    -1  ]

ieye(n) = eye(Int8, n)

@test ifwht(eye(8),1) == H8
@test ifwht(eye(8),2)' == H8
@test ifwht_natural(eye(8),1) == H8n
@test ifwht_natural(eye(8),2)' == H8n
@test ifwht_dyadic(eye(8),1) == H8d
@test ifwht_dyadic(eye(8),2)' == H8d

@test ifwht(ieye(8),1) == H8
@test ifwht(ieye(8),2)' == H8
@test ifwht_natural(ieye(8),1) == H8n
@test ifwht_natural(ieye(8),2)' == H8n
@test ifwht_dyadic(ieye(8),1) == H8d
@test ifwht_dyadic(ieye(8),2)' == H8d

H32 = [  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1
 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1
 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 -1 -1 -1 -1
 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1
 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1
 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1
 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1
 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1
 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1
 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 1 1 -1 -1
 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1
 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1
 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1
 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1
 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1
 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1
 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1
 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 1 1 -1
 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1
 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1
 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1
 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1
 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1
 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1
 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1
 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1
 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1
 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 ]

@test ifwht(eye(32),1) == H32
@test ifwht(eye(32),2)' == H32

@test ifwht(ieye(32),1) == H32
@test ifwht(ieye(32),2)' == H32

X = reshape(sin([1:1024*32;]), 1024,32);
norminf(A) = maximum(abs(A))

for f in (:fwht, :fwht_natural, :fwht_dyadic)
    fi = Symbol(string("i", f))
    @eval begin
        @test norminf(X - $fi($f(X))) < 1e-14
        @test norminf(X - $fi($f(X,1),1)) < 1e-14
        @test norminf(X - $fi($f(X,2),2)) < 1e-14
        @test norminf($f($f(X,1),2) - $f(X)) < 1e-14
        @test norminf($f(X',1)' - $f(X,2)) < 1e-14
    end
end

@test hadamard(8) == H8n
@test ifwht_natural(eye(32), 1) == hadamard(32)
@test ifwht_natural(eye(2), 1) == hadamard(2)

print("Checking unitarity of hadamard(n): ")
for i = 4:4:1000
    H = try
        convert(Matrix{Int}, hadamard(i))
    catch
        Int[]
    end
    if !isempty(H)
        print(i, ", ")
        @test norm(H'*H - size(H,1)*eye(H), 1) == 0
    end
end
println(".\nSuccess!")


text:
include("containers.jl")
include("math.jl")

# Store main timer for global timing of functions
const main_timer = TimerOutput()

# Always call timer() to hide implementation details
timer() = main_timer

# Initialize top-level parameters structure for program-wide parameters
const parameters = Dict{Symbol,Any}()

# Parse parameters file into global dict
function parse_parameters_file(filename)
  parameters[:default] = parsefile(filename)
  parameters[:default]["parameters_file"] = filename
end

# Return parameter by name, optionally taking a default value and a range of valid values.
#
# If no default value is specified, the parameter is required and the program
# stops if the parameter was not found. The range of valid parameters is used
# to restrict parameters to sane values.
function parameter(name, default=nothing; valid=nothing)
  if haskey(parameters[:default], name)
    # If parameter exists, use its value
    value = parameters[:default][name]
  else
    # Otherwise check whether a default is given and abort if not
    if default === nothing
      error("requested parameter '$name' does not exist and no default value was provided")
    else
      value = default
    end
  end

  # If a range of valid values has been specified, check parameter value against it
  if valid !== nothing
    if !(value in valid)
      error("'$value' is not a valid value for parameter '$name' (valid: $valid)")
    end
  end

  return value
end

"""
    examples_dir()

Return the directory where the example files provided with Trixi.jl are located. If Trixi is
installed as a regular package (with `]add Trixi`), these files are read-only and should *not* be
modified. To find out which files are available, use, e.g., `readdir`:

# Examples
```julia
julia> readdir(examples_dir())
4-element Array{String,1}:
 "2d"
 "3d"
 "README.md"
 "paper-self-gravitating-gas-dynamics"
```
"""
examples_dir() = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

"""
    get_examples()

Return a list of all example parameter files that are provided by Trixi.
"""
function get_examples()
  examples = String[]
  for (root, dirs, files) in walkdir(examples_dir())
    for f in files
      if endswith(f, ".toml")
        push!(examples, joinpath(root, f))
      end
    end
  end

  return examples
end

"""
    default_example()

Return the path to an example parameter file that can be used to quickly see Trixi in action.
"""
default_example() = joinpath(examples_dir(), "2d", "parameters.toml")

"""
    setparameter(name::String, value)

Set parameter with the specified `name` to the specified `value`.
"""
function setparameter(name::String, value)
  parameters[:default][name] = value
end

# Return true if parameter exists.
parameter_exists(name::String) = haskey(parameters[:default], name)

# Print informative message at startup
function print_startup_message()
  s = """

    ████████╗██████╗ ██╗██╗  ██╗██╗
    ╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
       ██║   ██████╔╝██║ ╚███╔╝ ██║
       ██║   ██╔══██╗██║ ██╔██╗ ██║
       ██║   ██║  ██║██║██╔╝ ██╗██║
       ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝
    """
  println(s)
end

"""
    get_name(x)

Returns a name of `x` ready for pretty printing.
By default, return `string(y)` if `x isa Val{y}` and return `string(x)` otherwise.

# Examples

```jldoctest
julia> Trixi.get_name("test")
"test"

julia> Trixi.get_name(Val(:test))
"test"
```
"""
get_name(x) = string(x)
get_name(::Val{x}) where x = string(x)


text:
"""

The exported names are
$(EXPORTS)
"""
module MultiField

using DocStringExtensions
using Gridap.Helpers
using Gridap.Algebra
using Gridap.Arrays
using Gridap.FESpaces
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.CellData
using Gridap.Fields

using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values
using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis
using Gridap.CellData: CellFieldAt
import Gridap.Fields: gradient, DIV, ∇∇

using ForwardDiff
using FillArrays
using SparseArrays
using LinearAlgebra
using BlockArrays

export num_fields
export compute_field_offsets
export restrict_to_field
export MultiFieldCellField
export MultiFieldFESpace
export MultiFieldFEFunction
export MultiFieldStyle
export ConsecutiveMultiFieldStyle

include("MultiFieldCellFields.jl")

include("MultiFieldFESpaces.jl")

include("MultiFieldFEFunctions.jl")

include("MultiFieldFEAutodiff.jl")

end # module


text:
using Flux

struct Xconv
    n_p::Int
    k::Int
    mlp1
    mlp2_dense
    mlp2_conv1
    mlp2_conv2
    conv
end

function Xconv(c1::Int,
    c2::Int,
    dim::Int,
    k::Int,
    n_p::Int,
    c_del::Int,
    dm::Int)
    Xconv(n_p,
        k,
        Chain(Conv((1, 1), dim=>c_del, elu), Conv((1, 1), c_del=>c_del, elu)),
        Conv((1, k), dim=>k*k, elu),
        DepthwiseConv((1, k), k=>k*k, elu),
        DepthwiseConv((1, k), k=>k*k),
        Chain(Conv((1, k), (c1+c_del)=>dm*(c1+c_del)), Conv((1, 1), dm*(c1+c_del)=>c2))
    )
end

function (m::Xconv)(p, P, F=nothing)
    # P : (n_p, k, 3, batch)
    # p : (n_p, 1, 3, batch)
    # F : (n_p, 1, c1, batch) or nothing

    p_shifted = P .- p

    F_del = m.mlp1(p_shifted) # F_del : (n_p, k, c_del, batch)
    F_star = F != nothing ? cat(F_del, F, dims=3) : F_del # F_star : (n_p, k, c_del + c1, batch)
    X = reshape(
        m.mlp2_conv2(reshape(
            m.mlp2_conv1(reshape(
                m.mlp2_dense(p_shifted),
                (m.n_p, m.k, m.k, :)
            )),
            (m.n_p, m.k, m.k, :)
        )),
        (m.n_p, m.k, m.k, :)
    )

    F_X = cat([cat(
        [reshape(X[i,:,:,j]*F_star[i,:,:,j], (1, size(X)[2], size(F_star)[3])) for j in 1:size(X)[end]]...,
        dims=4
    ) for i in 1:size(X)[1]]..., dims=1);

    F_p = m.conv(F_X)
    return F_p
end

Flux.@functor Xconv


text:
# test idea of Liouville evolution
using LinearAlgebra
using QuantumInformation

H_ctrl = [pi * sx, pi * sy]
H0 = sz

function to_superoperator(H)
    D = size(H)[1]
    kron(I(D), H) - kron(H', I(D))
end

H0 = to_superoperator(sz)
H_ctrls = to_superoperator.(H_ctrl)
duration = 1
n_slices = 10
drive = rand(2, n_slices)

function prop(x)
    Htotal = H0
    for k = 1:2
        Htotal = Htotal + x[k] * H_ctrls[k]
    end
    U = exp(-im * Htotal * dt)
end

p = []
for s = 1:n_slices
    p_n = prop(drive[:, s])

    append!(p, [p_n])
end

psi = [1; 0.0 + 0.0im]

rho0 = psi * psi'
# stack it into vector
rho0 = reshape(rho0, 4)

o = [rho0]
for s = 1:n_slices
    o_n = p[s] * o[s]
    append!(o, [o_n])
end

obs = to_superoperator(sz)
z = []
for s = 1:n_slices
    # temp = obs * o[s]
    kk = reshape(o[s], (2, 2))
    out = real(tr(sz * kk))
    append!(z, out)
end
using Plots
plot(z)

##

using FastClosures

prob = StateTransferProblem(
    B = [Sx, Sy],
    A = Sz,
    X_init = ρinit,
    X_target = ρfin,
    duration = 5.0,
    n_timeslices = 25,
    n_controls = 2,
    initial_guess = rand(2, 25),
)

ensemble_problem =
    ClosedEnsembleProblem(prob, 5, A_gens, B_gens, X_init_gens, X_target_gens, ones(5) / 5)

# sol = GRAPE(ens, inplace=true)
# @test sol.result[1].minimum - C1(ρfin, ρfin) < tol * 10

function ensemble_GRAPE!_unchanged(ensemble_problem, optim_options)
    # we expand the ensemble problem into an array of single problems that we can solve
    ensemble_problem_array = QuOptimalControl.init_ensemble(ensemble_problem)

    # to initialise the holding arrays we define a problem
    problem = ensemble_problem_array[1]
    # initialise holding arrays for inplace operations, these will be modified
    state_store, costate_store, generators, propagators, ofom, gradient = init_GRAPE(
        problem.X_init,
        problem.n_timeslices,
        ensemble_problem.n_ensemble,
        problem.A,
        problem.n_controls,
    )

    evolve_store = similar(state_store[1])
    weights = ensemble_problem.weights
    n_ensemble = ensemble_problem.n_ensemble

    function _to_optim!(F, G, x)
        fom = 0.0
        for k = 1:n_ensemble
            curr_prob = ensemble_problem_array[k]

            fom += @views QuOptimalControl._fom_and_gradient_GRAPE!(
                curr_prob.A,
                curr_prob.B,
                x,
                curr_prob.n_timeslices,
                curr_prob.duration,
                curr_prob.n_controls,
                gradient[k, :, :],
                state_store[:, k],
                costate_store[:, k],
                generators[:, k],
                propagators[:, k],
                curr_prob.X_init,
                curr_prob.X_target,
                evolve_store,
                curr_prob,
            ) * weights[k]
        end

        if G !== nothing
            @views G .= sum(gradient .* weights, dims = 1)[1, :, :]
        end
        if F !== nothing
            return fom
        end
    end
    return _to_optim!

    # init = problem.initial_guess

    # res = Optim.optimize(Optim.only_fg!(_to_optim!), init, Optim.LBFGS(), optim_options)
    # solres = SolutionResult([res], [res.minimum], [res.minimizer], problem)
    # return solres
end

prob = ensemble_problem

@code_warntype ensemble_GRAPE!_unchanged(ensemble_problem, 1.0)

topt = ensemble_GRAPE!_unchanged(ensemble_problem, 1.0)

@code_warntype topt(nothing, nothing, rand(2, 25))

@benchmark topt($nothing, $nothing, $rand(2, 25))

function ensemble_GRAPE!closure(ensemble_problem, optim_options)
    # we expand the ensemble problem into an array of single problems that we can solve
    ensemble_problem_array = QuOptimalControl.init_ensemble(ensemble_problem)

    # to initialise the holding arrays we define a problem
    problem = ensemble_problem_array[1]
    # initialise holding arrays for inplace operations, these will be modified
    state_store, costate_store, generators, propagators, ofom, gradient = init_GRAPE(
        problem.X_init,
        problem.n_timeslices,
        ensemble_problem.n_ensemble,
        problem.A,
        problem.n_controls,
    )

    evolve_store = similar(state_store[1])
    weights = ensemble_problem.weights
    n_ensemble = ensemble_problem.n_ensemble

    function _to_optim!(F, G, x)
        fom = 0.0
        for k = 1:n_ensemble
            curr_prob = ensemble_problem_array[k]

            fom += @views QuOptimalControl._fom_and_gradient_GRAPE!(
                curr_prob.A,
                curr_prob.B,
                x,
                curr_prob.n_timeslices,
                curr_prob.duration,
                curr_prob.n_controls,
                gradient[k, :, :],
                state_store[:, k],
                costate_store[:, k],
                generators[:, k],
                propagators[:, k],
                curr_prob.X_init,
                curr_prob.X_target,
                evolve_store,
                curr_prob,
            ) * weights[k]
        end

        if G !== nothing
            @views G .= sum(gradient .* weights, dims = 1)[1, :, :]
        end
        if F !== nothing
            return fom
        end
    end
    return _to_optim!

    # init = problem.initial_guess

    # res = Optim.optimize(Optim.only_fg!(_to_optim!), init, Optim.LBFGS(), optim_options)
    # solres = SolutionResult([res], [res.minimum], [res.minimizer], problem)
    # return solres
end

@code_warntype ensemble_GRAPE!closure(ensemble_problem, 1.0)

toptclosure = ensemble_GRAPE!closure(ensemble_problem, 1.0)

@code_warntype toptclosure(nothing, nothing, rand(2, 25))

function GRAPEfix!(problem, optim_options)

    # initialise holding arrays for inplace operations, these will be modified
    state_store, costate_store, generators, propagators, fom, gradient =
        init_GRAPE(problem.X_init, problem.n_timeslices, 1, problem.A, problem.n_controls)

    evolve_store = similar(state_store[1])
    k = 1
    function _to_optim!(F, G, x)
        fom = 0.0
        fom += @views _fom_and_gradient_GRAPE!(
            problem.A,
            problem.B,
            x,
            problem.n_timeslices,
            problem.duration,
            problem.n_controls,
            gradient[k, :, :],
            state_store[:, k],
            costate_store[:, k],
            generators[:, k],
            propagators[:, k],
            problem.X_init,
            problem.X_target,
            evolve_store,
            problem,
        )
        if G !== nothing
            @views G .= gradient
        end
        if F !== nothing
            return fom
        end
    end

    init = problem.initial_guess

    res = Optim.optimize(Optim.only_fg!(_to_optim!), init, Optim.LBFGS(), optim_options)
    solres = SolutionResult([res], [res.minimum], [res.minimizer], problem)
    return solres
end

using Optim

GRAPEfix!(prob, 1.0)

problem = prob
# initialise holding arrays for inplace operations, these will be modified
state_store, costate_store, generators, propagators, fom, gradient =
    init_GRAPE(problem.X_init, problem.n_timeslices, 1, problem.A, problem.n_controls)

evolve_store = similar(state_store[1])
k = 1
function _to_optim!(F, G, x)
    fom = 0.0
    fom += @views QuOptimalControl._fom_and_gradient_GRAPE!(
        problem.A,
        problem.B,
        x,
        problem.n_timeslices,
        problem.duration,
        problem.n_controls,
        gradient[k, :, :],
        state_store[:, k],
        costate_store[:, k],
        generators[:, k],
        propagators[:, k],
        problem.X_init,
        problem.X_target,
        evolve_store,
        problem,
    )
    if G !== nothing
        @views G .= gradient[1, :, :]
    end
    if F !== nothing
        return fom
    end
end

init = problem.initial_guess

res = Optim.optimize(Optim.only_fg!(_to_optim!), init, Optim.LBFGS())
solres = SolutionResult([res], [res.minimum], [res.minimizer], problem)


text:
mutable struct LangData
	langs::Dict{Symbol,Dict{Symbol,Dict}}
end
struct SubLangData
	ref::Dict
end
function Base.getindex(l::LangData,key::Symbol)
	return SubLangData(l.langs[lsetting[:lang]][key])
end
function Base.getindex(l::SubLangData,keys...)
	d=l
	for key in keys
		if haskey(d.ref,key)
			@inbounds r=d.ref[key]
			if isa(r,Dict)
				d=SubLangData(r)
			else
				return r
			end
		else
			return "$key"
		end
	end
	return d
end
function load_langdata_fromdir!(ld::LangData,lang::Symbol,dirpath::String)
	ld.langs[lang]=Dict{Symbol,Dict}()
	li=readdir(dirpath;sort=false)
	for i in li
		ld.langs[lang][Symbol(i)]=read_code(joinpath(dirpath,i))
	end
end


text:
function jump_randsdp(solver, seed, n, m, verbose = false)

    A, b, C = randsdp_data(seed, m, n)

    model = Model(with_optimizer(solver))
    @variable(model, X[1:n, 1:n], PSD)
    @objective(model, Min, sum(C[i, j] * X[i, j] for j in 1:n, i in 1:n))
    @constraint(model, ctr[k in 1:m], sum(A[k][i, j] * X[i, j] for j in 1:n, i in 1:n) == b[k])
    # @constraint(model, bla, sum(C[i, j] * X[i, j] for i in 1:n, j in 1:n)<=0.1)

    teste = @time optimize!(model)

    XX = value.(X)

    verbose && randsdp_eval(A,b,C,n,m,XX)

    objval = objective_value(model)
    stime = MOI.get(model, MOI.SolveTime())

    # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
    # @show fieldnames(tp)
    rank = -1
    try
        @show rank = model.moi_backend.optimizer.model.optimizer.sol.final_rank
    catch
    end
    status = 0
    if JuMP.termination_status(model) == MOI.OPTIMAL
        status = 1
    end
    return (objval, stime, rank, status, -1.0, -1.0)
end


text:
round(x::ArbFloat{P}, ::RoundingMode{:Up}) where {P} = ceil(x)
round(x::ArbFloat{P}, ::RoundingMode{:Down}) where {P} = floor(x)
round(x::ArbFloat{P}, ::RoundingMode{:ToZero}) where {P} = signbit(x) ? ceil(x) : floor(x)
round(x::ArbFloat{P}, ::RoundingMode{:FromZero}) where {P} = signbit(x) ? floor(x) : ceil(x)

round(x::ArbReal{P}, ::RoundingMode{:Up}) where {P} = ArbReal(round(ArbFloat(x), RoundUp))
round(x::ArbReal{P}, ::RoundingMode{:Down}) where {P} = ArbReal(round(ArbFloat(x), RoundDown))
round(x::ArbReal{P}, ::RoundingMode{:ToZero}) where {P} = ArbReal(round(ArbFloat(x), RoundToZero))
round(x::ArbReal{P}, ::RoundingMode{:FromZero}) where {P} = ArbReal(round(ArbFloat(x), RoundFromZero))

round(x::ArbComplex{P}, ::RoundingMode{:Up}) where {P} =
    ArbComplex(ArbReal(round(ArbFloat(real(x)), RoundUp)),
               ArbReal(round(ArbFloat(imag(x)), RoundUp)) )
round(x::ArbComplex{P}, ::RoundingMode{:Down}) where {P} =
    ArbComplex(ArbReal(round(ArbFloat(real(x)), RoundDown)),
               ArbReal(round(ArbFloat(imag(x)), RoundDown)) )
round(x::ArbComplex{P}, ::RoundingMode{:ToZero}) where {P} =
    ArbComplex(ArbReal(round(ArbFloat(real(x)), RoundToZero)),
               ArbReal(round(ArbFloat(imag(x)), RoundToZero)) )
round(x::ArbComplex{P}, ::RoundingMode{:FromZero}) where {P} =
    ArbComplex(ArbReal(round(ArbFloat(real(x)), RoundFromZero)),
               ArbReal(round(ArbFloat(imag(x)), RoundFromZero)) )

function round(x::ArbFloat{P}, ::RoundingMode{:Nearest}) where {P}
    absx = abs(x)
    fl = floor(absx)
    cl = fl + 1
    xfl = absx - fl
    xcl = cl - absx
    if xfl < xcl
        copysign(fl, x)
    elseif xfl > xcl
        copysign(cl, x)
    elseif iseven(BigInt(fl))
        copysign(fl, x)
    else
        copysign(cl, x)
    end
end

round(x::ArbReal{P}, ::RoundingMode{:Nearest}) where {P} =
    ArbReal(round(ArbFloat(x), RoundNearest))
round(x::ArbComplex{P}, ::RoundingMode{:Nearest}) where {P} =
    ArbComplex(ArbReal(round(ArbFloat(real(x)), RoundNearest)),
               ArbReal(round(ArbFloat(imag(x)), RoundNearest)) )

#=
round(x::ArbFloat{P}) where {P} = round(x, RoundNearest)
round(x::ArbReal{P}) where {P} = round(x, RoundNearest)
round(x::ArbComplex{P}) where {P} = round(x, RoundNearest)

round!(x::ArbFloat{P}) where {P} = round!(x, RoundNearest)
round!(x::ArbReal{P}) where {P} = round!(x, RoundNearest)
round!(x::ArbComplex{P}) where {P} = round!(x, RoundNearest)
=#

round(x::ArbFloat; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round(x,RoundNearest) : round(x,RoundNearest,sigdigits=sigdigits,base=base)
round(x::ArbReal; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round(x,RoundNearest) : round(x,RoundNearest,sigdigits=sigdigits,base=base)
round(x::ArbComplex; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round(x,RoundNearest) : round(x,RoundNearest,sigdigits=sigdigits,base=base)

round!(x::ArbFloat; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round!(x,RoundNearest) : round!(x,RoundNearest,sigdigits=sigdigits,base=base)
round!(x::ArbReal; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round!(x,RoundNearest) : round!(x,RoundNearest,sigdigits=sigdigits,base=base)
round!(x::ArbComplex; sigdigits::Integer=0, base::Integer=10) = sigdigits==0 ? round!(x,RoundNearest) : round!(x,RoundNearest,sigdigits=sigdigits,base=base)

#=
function round(x::ArbComplex{P}, roundingmode::RoundingMode) where {P}
    z = ArbComplex{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(acb_set_round), Cint, (Ref{ArbComplex}, Ref{ArbComplex}, Clong, Cint), z, x, P, rounding)
    return z
end
=#

# int arf_set_round(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)
function roundbits(x::ArbFloat{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    sigbits >= P && return ArbFloat{sigbits}(x)
    z = ArbFloat{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, sigbits, rounding)
    return z
end

roundbits(x::ArbFloat{P}, sigbits::Int) where {P} = roundbits(x, RoundNearest, sigbits)

function roundbits!(x::ArbFloat{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    z = roundbits(x, roundingmode, sigbits)
    return ArbFloat{sigbits+extrabits()}(z)
end

roundbits!(x::ArbFloat{P}, sigbits::Int) where {P} = roundbits!(x, RoundNearest, sigbits)

rounddigits(x::ArbFloat{P}, roundingmode::RoundingMode, digits::Int) where {P} =
    roundbits(x, roundingmode, bits4digits(digits))

rounddigits(x::ArbFloat{P}, digits::Int) where {P} = roundbits(x, RoundNearest, bits4digits(digits))

rounddigits!(x::ArbFloat{P}, roundingmode::RoundingMode, digits::Int) where {P} =
    roundbits!(x, roundingmode, bits4digits(digits))

rounddigits!(x::ArbFloat{P}, digits::Int) where {P} = roundbits!(x, RoundNearest, bits4digits(digits))

function roundbits(x::ArbReal{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    sigbits >= P && return ArbReal{sigbits}(x)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_round), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, sigbits)
    return z
end

function roundbits!(x::ArbReal{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    z = roundbits(x, roundingmode, sigbits)
    P2 = sigbits + extrabits()
    return ArbReal{P2}(z)
end

roundbits(x::ArbReal{P}, sigbits::Int) where {P} = roundbits(x, RoundNearest, sigbits)

rounddigits(x::ArbReal{P}, roundingmode::RoundingMode, digits::Int) where {P} =
    roundbits(x, roundingmode, bits4digits(digits))

rounddigits(x::ArbReal{P}, digits::Int) where {P} = roundbits(x, RoundNearest, bits4digits(digits))

rounddigits!(x::ArbReal{P}, roundingmode::RoundingMode, digits::Int) where {P} =
    roundbits!(x, roundingmode, bits4digits(digits))

rounddigits!(x::ArbReal{P}, digits::Int) where {P} = roundbits!(x, RoundNearest, bits4digits(digits))

function roundbits(x::ArbComplex{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    sigbits >= P && return ArbComplex{sigbits}(x)
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_round), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, sigbits)
    return z
end

function roundbits!(x::ArbComplex{P}, roundingmode::RoundingMode, sigbits::Int) where {P}
    z = roundbits(x, roundingmode, sigbits)
    P2 = sigbits + extrabits()
    return ArbComplex{P2}(z)
end

roundbits(x::ArbComplex{P}, sigbits::Int) where {P} = roundbits(x, RoundNearest, sigbits)

rounddigits(x::ArbComplex{P}, roundingmode::RoundingMode, digits::Int) where {P} =
    roundbits(x, roundingmode, bits4digits(digits))

rounddigits(x::ArbComplex{P}, digits::Int) where {P} = roundbits(x, RoundNearest, bits4digits(digits))

function round(x::ArbFloat{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   r = round!(x, roundingmode, sigdigits=sigdigits, base=base)
   s = string(r)
   result = ArbFloat{P}(s)
   return result
end

function round!(x::ArbFloat{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   if base==10
       return rounddigits!(x, roundingmode, sigdigits)
   elseif base == 2
       return roundbits!(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

function round(x::ArbReal{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   if base==10
       return rounddigits(x, roundingmode, sigdigits)
   elseif base==2
       return roundbits(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

function round!(x::ArbReal{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   if base==10
       return rounddigits!(x, roundingmode, sigdigits)
   elseif base == 2
       return roundbits!(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

function round(x::ArbComplex{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   if base==10
       return rounddigits(x, roundingmode, sigdigits)
   elseif base==2
       return roundbits(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

function round!(x::ArbComplex{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P}
   if base==10
       return rounddigits!(x, roundingmode, sigdigits)
   elseif base == 2
       return roundbits!(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

function roundfrac(x::ArbFloat{P}, roundingmode::RoundingMode, sigdigits::Integer, base::Integer) where {P}
   if base==10
       return rounddigits(x, roundingmode, sigdigits)
   elseif base == 2
       return roundbits(x, roundingmode, sigdigits)
   else
       throw(ErrorException("unsupported base ($base)"))
   end
end

#=
function round(x::ArbFloat{P}, bitprecision::Int) where {P}
    z = ArbFloat{P}()
    rounding = match_rounding_mode(RoundNearest)
    rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

# void arb_set_round(arb_t y, const arb_t x, slong prec)
function round(x::ArbReal{P}, bitprecision::Int) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_round), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, bitprecision)
    return z
end

#void acb_set_round(acb_t z, const acb_t x, slong prec)
function round(x::ArbComplex{P}, bitprecision::Int) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_round), Cint, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, bitprecision)
    return z
end
=#

#=
for A in (:ArbFloat, :ArbReal, :ArbComplex)
  @eval begin

    round(::Type{T}, x::$A{P}, rounding::RoundingMode=RoundNearest) where {P,T} =
       T(round(x, digits=P, base=rounding))

    function round(x::$A{P}, rounding::RoundingMode=RoundNearest;
                   sigdigits::Integer = 0, digits::Integer = 0, base = 10) where {P}
        sigdigits = max(sigdigits, digits)
        if base == 10
            round(x, digits=digits4bits(sigdigits), base=rounding)
        elseif base ==  2
            round(x, digits=sigdigits, base=rounding)
        else
            throw(DomainError("base ($base) not supported"))
        end
     end
  end

end
=#

#=
  round([T,] x, [r::RoundingMode])
  round(x, [r::RoundingMode]; digits::Integer=0, base = 10)
  round(x, [r::RoundingMode]; sigdigits::Integer, base = 10)

  Rounds the number x.

  Without keyword arguments, x is rounded to an integer value, returning a value of type T, or of the same type of x
  if no T is provided. An InexactError will be thrown if the value is not representable by T, similar to convert.

  If the digits keyword argument is provided, it rounds to the specified number of digits after the decimal place (or
  before if negative), in base base.

  If the sigdigits keyword argument is provided, it rounds to the specified number of significant digits, in base
  base.

  The RoundingMode r controls the direction of the rounding; the default is RoundNearest, which rounds to the nearest
  integer, with ties (fractional values of 0.5) being rounded to the nearest even integer. Note that round may give
  incorrect results if the global rounding mode is changed (see rounding).
=#


text:
function pascaltriangle!(Msize::Int64, Np::Int64, matp::Matrix{Int64})

    Lm = Msize+1
    Ln = Np + 1

    matp .= 0 # matp = zeros(Int,Lm,Ln)
    matp[:,1] .= 1

    for jj = 2:Lm
        for ii = 2:Ln
            matp[jj,ii] = sum(matp[jj-1,1:ii])
        end
    end

    # return matp

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
### routines for checking proximity or intervals for floating-point numbers.

function isnumericallyclose(x::T, y::T, tol = eps(T)*2) where T
    if abs(x-y) < tol
        return true
    end

    return false
end

function isnumericallyin(x::T, a::T, b::T, tol = eps(x)*2)::Bool where T
    if x > a-tol && x < b+tol
        return true
    end

    return false
end

# remove points that are too close to the entry before it.
# also sorts the collect A.
function makenotclose!(A, tol)::Nothing
    remove_list = Vector{Int}(undef,0)
    sort!(A)

    for i = 2:length(A)
        if isnumericallyclose(A[i-1], A[i], tol)
            push!(remove_list, i)
        end
    end

    deleteat!(A,remove_list)

    return nothing
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
#######################################################################
# Data types for the analysis of lower bounds
#######################################################################

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Analysis
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#--------------------------------------------------
# Simple search of text for `<:` and `>:`
#--------------------------------------------------

@enum ConstrKind subtc suptc

# Regex expressions corresponding to
# subtype and supertype constraints
const CONSTRAINT_PATTERN = Dict(
    subtc => "<:",
    suptc => ">:"
)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Gathering statistics
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#--------------------------------------------------
# Cumulative lower-bounds stat
#--------------------------------------------------

# Information about textual and parsing-based
# appearance of lower-bound constraints
struct TxtConstrStat
    subConsr :: UInt     # number of `<:` in text
    supConsr :: UInt     # number of `>:` in text
end
TxtConstrStat() = TxtConstrStat(0, 0)

# Frequencies of lower-bound values
LBValsFreq = Multiset

# Information about lower-bound constraints
struct LBStat
    lbs         :: UInt     # number of lower bounds
    lbsUnique   :: UInt     # number of unique lower-bound values
    lbsFreq     :: LBValsFreq   # frequencies of lower-bound values
end
LBStat() = LBStat(0, 0, LBValsFreq())
LBStat(lbsFreq :: LBValsFreq) =
    LBStat(length(lbsFreq), lengthUnique(lbsFreq), lbsFreq)

#--------------------------------------------------
# File stat
#--------------------------------------------------

# Information about lower bounds in a file:
# - `txtStat`   textual bounds info
# - `err`       possibly exception during processing
# - `lbStat`    possibly proper lb-stat if `txtStat` is non-vacuous
struct FileLBInfo
    txtStat :: TxtConstrStat
    err     :: Union{Exception, Nothing}
    lbStat  :: Union{LBStat, Nothing}
end
FileLBInfo(txtStat :: TxtConstrStat) =
    FileLBInfo(txtStat, nothing, nothing)
FileLBInfo(txtStat :: TxtConstrStat, err :: Exception) =
    FileLBInfo(txtStat, err, nothing)
FileLBInfo(txtStat :: TxtConstrStat, lbStat :: LBStat) =
    FileLBInfo(txtStat, nothing, lbStat)

#--------------------------------------------------
# Package stat
#--------------------------------------------------

# Files statistics (fileName => statistics)
FilesLBInfo = Dict{String, FileLBInfo}

# Single package statistics
mutable struct PackageStat
    name             :: String
    hasSrc           :: Bool
    totalFilesNum    :: UInt # number of source Julia files
    failedFilesNum   :: UInt # number of files that failed to process
    nonVacFilesNum   :: UInt # number of files with lower bounds
    filesInfo        :: FilesLBInfo # fileName => FileLBInfo
    lbStat           :: LBStat      # package cumulative statistics
end
PackageStat(name :: String, hasSrc :: Bool) =
    PackageStat(name, hasSrc, 0, 0, 0, FilesLBInfo(), LBStat())

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Base functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#--------------------------------------------------
# Equality
#--------------------------------------------------

Base.:(==)(v1 :: TxtConstrStat, v2 :: TxtConstrStat) = structEqual(v1, v2)

Base.:(==)(v1 :: LBStat, v2 :: LBStat) = structEqual(v1, v2)

Base.:(==)(v1 :: FileLBInfo, v2 :: FileLBInfo) = structEqual(v1, v2)

Base.:(==)(v1 :: PackageStat, v2 :: PackageStat) = structEqual(v1, v2)

#--------------------------------------------------
# Show
#--------------------------------------------------

Base.show(io :: IO, un :: UInt) = print(io, string(un, base=10))

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
abstract type AbstractGrid{T, N} <: AbstractArray{T, N} end

"""
    struct Grid{T, N, S <: AbstractCoordinateSystem, AT} <: AbstractGrid{T, N}

Collection of `N` axes that define the dimensions of the grid needed to calculate
[`ElectricPotential`](@ref), [`ElectricField`](@ref) or [`WeightingPotential`](@ref).

## Parametric types
* `T`: Tick type (element type) of the axes.
* `N`: Dimension of the grid.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `axes::AT`: Tuple of length `N` containing `DiscreteAxis` for each dimension of the grid.

See also [`DiscreteAxis`](@ref).
"""
struct Grid{T, N, S <: AbstractCoordinateSystem, AT} <: AbstractGrid{T, N}
    axes::AT
end

const CartesianGrid{T, N} = Grid{T, N, Cartesian}
const CartesianGrid1D{T} = CartesianGrid{T, 1}
const CartesianGrid2D{T} = CartesianGrid{T, 2}
const CartesianGrid3D{T} = CartesianGrid{T, 3}
const CylindricalGrid{T} = Grid{T, 3, Cylindrical}
#const RadialGrid{T} = Grid{T, 1, Radial}
#const PolarGrid{T} = Grid{T, 2, Polar}
#const SphericalGrid{T} = Grid{T, 3, Spherical}

CylindricalGrid{T}(a) where {T} = Grid{T, 3, Cylindrical, typeof(a)}(a)
CartesianGrid3D{T}(a) where {T} = Grid{T, 3, Cartesian, typeof(a)}(a)

@inline size(grid::Grid{T, N, S}) where {T, N, S} = size.(grid.axes, 1)
@inline length(grid::Grid{T, N, S}) where {T, N, S} = prod(size(grid))
@inline getindex(grid::Grid{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = broadcast(getindex, grid.axes, I)
@inline getindex(grid::Grid{T, N, S}, i::Int) where {T, N, S} = getproperty(grid, :axes)[i]
@inline getindex(grid::Grid{T, N, S}, s::Symbol) where {T, N, S} = getindex(grid, Val{s}())

@inline getproperty(grid::Grid{T, N, S}, s::Symbol) where {T, N, S} = getproperty(grid, Val{s}())
@inline getproperty(grid::Grid{T}, ::Val{:axes}) where {T} = getfield(grid, :axes)

@inline getproperty(grid::CylindricalGrid{T}, ::Val{:axes}) where {T} = getfield(grid, :axes)
@inline getproperty(grid::CylindricalGrid{T}, ::Val{:r}) where {T} = @inbounds grid.axes[1]
@inline getproperty(grid::CylindricalGrid{T}, ::Val{:φ}) where {T} = @inbounds grid.axes[2]
@inline getproperty(grid::CylindricalGrid{T}, ::Val{:z}) where {T} = @inbounds grid.axes[3]
@inline getproperty(grid::CartesianGrid3D{T}, ::Val{:x}) where {T} = @inbounds grid.axes[1]
@inline getproperty(grid::CartesianGrid3D{T}, ::Val{:y}) where {T} = @inbounds grid.axes[2]
@inline getproperty(grid::CartesianGrid3D{T}, ::Val{:z}) where {T} = @inbounds grid.axes[3]

@inline getindex(grid::CylindricalGrid{T}, ::Val{:r}) where {T} = @inbounds grid.axes[1]
@inline getindex(grid::CylindricalGrid{T}, ::Val{:φ}) where {T} = @inbounds grid.axes[2]
@inline getindex(grid::CylindricalGrid{T}, ::Val{:z}) where {T} = @inbounds grid.axes[3]
@inline getindex(grid::CartesianGrid3D{T}, ::Val{:x}) where {T} = @inbounds grid.axes[1]
@inline getindex(grid::CartesianGrid3D{T}, ::Val{:y}) where {T} = @inbounds grid.axes[2]
@inline getindex(grid::CartesianGrid3D{T}, ::Val{:z}) where {T} = @inbounds grid.axes[3]

@inline GridPoint(grid::Grid{T, 3, Cylindrical}, inds::NTuple{3, Int}) where {T} =
    CylindricalPoint{T}(broadcast(i -> grid.axes[i].ticks[inds[i]], (1, 2, 3)))
@inline GridPoint(grid::Grid{T, 3, Cartesian}, inds::NTuple{3, Int}) where {T} =
    CartesianPoint{T}(broadcast(i -> grid.axes[i].ticks[inds[i]], (1, 2, 3)))

function sizeof(grid::Grid{T, N, S}) where {T, N, S}
    return sum( sizeof.(grid.axes) )
end

function print(io::IO, grid::Grid{T, N, S}) where {T, N, S}
    print(io, "Grid{$T, $N, $S}", grid.axes)
end
function println(io::IO, grid::Grid{T, N, S}) where {T, N, S}
    println(" Grid{$T, $N, $S}")
    for (i, ax) in enumerate(grid.axes)
        println(io, "  Axis $(i): ", ax)
    end
end
show(io::IO, grid::Grid{T, N, S}) where {T, N, S} = print(io, grid)
show(io::IO, ::MIME"text/plain", grid::Grid{T, N, S}) where {T, N, S} = show(io, grid)

function check_grid(grid::CylindricalGrid{T})::Nothing where {T}
    nr::Int, nφ::Int, nz::Int = size(grid)
    @assert iseven(nz) "GridError: Field simulation algorithm in cylindrical coordinates needs an even number of grid points in z. This is not the case. #z-ticks = $(nz)."
    @assert (iseven(nφ) || (nφ == 1)) "GridError: Field simulation algorithm in cylindrical coordinates needs an even number of grid points in φ or just one point (2D). This is not the case. #φ-ticks = $(nφ)."
    return nothing
end

function check_grid(grid::CartesianGrid3D{T})::Nothing where {T}
    nx::Int, ny::Int, nz::Int = size(grid)
    @assert iseven(nx) "GridError: Field simulation algorithm in cartesian coordinates needs an even number of grid points in x. This is not the case. #x-ticks = $(nx)."
    return nothing
end

function get_coordinate_system(grid::Grid{T, N, S}) where {T, N, S}
    return S
end
function get_number_of_dimensions(grid::Grid{T, N, S}) where {T, N, S}
    return N
end
function Base.eltype(grid::Grid{T, N, S})::DataType where {T, N, S}
    return T
end

function get_boundary_types(grid::Grid{T, N, S}) where {T, N, S}
   return get_boundary_types.(grid.axes)
end

# Tuples with ticks to sample with differently spaced ticks
const CartesianTicksTuple{T} = NamedTuple{(:x,:y,:z), NTuple{3,Vector{T}}}
const CylindricalTicksTuple{T} = NamedTuple{(:r,:φ,:z), NTuple{3,Vector{T}}}

TicksTuple(grid::CartesianGrid3D{T}) where {T} = (x = grid.axes[1].ticks, y = grid.axes[2].ticks, z = grid.axes[3].ticks)
TicksTuple(grid::CylindricalGrid{T}) where {T} = (r = grid.axes[1].ticks, φ = grid.axes[2].ticks, z = grid.axes[3].ticks)

function Grid(nt::NamedTuple)
    if nt.coordtype == "cylindrical"
        axr::DiscreteAxis = DiscreteAxis(nt.axes.r, unit=u"m")
        axφ::DiscreteAxis = DiscreteAxis(nt.axes.phi, unit=u"rad")
        axz::DiscreteAxis = DiscreteAxis(nt.axes.z, unit=u"m")
        T = typeof(axr.ticks[1])
        return CylindricalGrid{T}( (axr, axφ, axz) )
    elseif nt.coordtype == "cartesian"
        axx::DiscreteAxis = DiscreteAxis(nt.axes.x, unit=u"m")
        axy::DiscreteAxis = DiscreteAxis(nt.axes.y, unit=u"m")
        axz = DiscreteAxis(nt.axes.z, unit=u"m")
        T = typeof(axx.ticks[1])
        return CartesianGrid3D{T}( (axx, axy, axz) )
    else
        error("`coordtype` = $(nt.coordtype) is not valid.")
    end
end

Base.convert(T::Type{Grid}, x::NamedTuple) = T(x)

function NamedTuple(grid::CylindricalGrid{T}) where {T}
    axr::DiscreteAxis{T} = grid.axes[1]
    axφ::DiscreteAxis{T} = grid.axes[2]
    axz::DiscreteAxis{T} = grid.axes[3]
    return (
        coordtype = "cylindrical",
        ndims = 3,
        axes = (
            r   = NamedTuple(axr, unit=u"m"),
            phi = NamedTuple(axφ, unit=u"rad"),
            z   = NamedTuple(axz, unit=u"m"),
        )
    )
end
function NamedTuple(grid::CartesianGrid3D{T}) where {T}
    axx::DiscreteAxis{T} = grid.axes[1]
    axy::DiscreteAxis{T} = grid.axes[2]
    axz::DiscreteAxis{T} = grid.axes[3]
    return (
        coordtype = "cartesian",
        ndims = 3,
        axes = (
            x = NamedTuple(axx, unit=u"m"),
            y = NamedTuple(axy, unit=u"m"),
            z = NamedTuple(axz, unit=u"m"),
        )
    )
end

Base.convert(T::Type{NamedTuple}, x::Grid) = T(x)

function find_closest_gridpoint(pt::CylindricalPoint{T}, grid::CylindricalGrid{T})::NTuple{3,Int} where {T <: SSDFloat}
    return (searchsortednearest(grid.axes[1].ticks, pt.r), searchsortednearest(grid.axes[2].ticks, pt.φ), searchsortednearest(grid.axes[3].ticks, pt.z))
end
function find_closest_gridpoint(pt::CartesianPoint{T}, grid::CylindricalGrid{T})::NTuple{3,Int} where {T <: SSDFloat}
    find_closest_gridpoint(CylindricalPoint(pt),grid)
end

function find_closest_gridpoint(pt::CartesianPoint{T}, grid::CartesianGrid3D{T})::NTuple{3,Int} where {T <: SSDFloat}
    @inbounds return (searchsortednearest(grid.axes[1].ticks, pt.x), searchsortednearest(grid.axes[2].ticks, pt.y), searchsortednearest(grid.axes[3].ticks, pt.z))
end
function find_closest_gridpoint(pt::CylindricalPoint{T}, grid::CartesianGrid3D{T})::NTuple{3,Int} where {T <: SSDFloat}
    find_closest_gridpoint(CartesianPoint(pt),grid)
end

multiplicity(g::Grid) = prod(multiplicities(g))

function multiplicities(g::CylindricalGrid{T}) where {T}
    mr = one(T)
    mφ = T(2π) / width(g.axes[2].interval)
    mz = multiplicity(g.axes[3], Cartesian)
    mr, mφ, mz
end

multiplicities(g::CartesianGrid3D) = broadcast(ax -> multiplicity(ax, Cartesian), g.axes)

function voxel_widths(grid::CartesianGrid3D{T}, i1::Int, i2::Int, i3::Int) where {T}
    wx::T = grid[1].ticks[i1 + 1] - grid[1].ticks[i1]
    wy::T = grid[2].ticks[i2 + 1] - grid[2].ticks[i2]
    wz::T = grid[3].ticks[i3 + 1] - grid[3].ticks[i3]
    wx, wy, wz
end
function voxel_widths(grid::CylindricalGrid{T}, i1::Int, i2::Int, i3::Int) where {T}
    wr::T =  grid[1].ticks[i1 + 1] - grid[1].ticks[i1]
    wφ::T = (grid[2].ticks[i2 + 1] - grid[2].ticks[i2]) * (grid[1].ticks[i1 + 1] + grid[1].ticks[i1])/2
    wz::T =  grid[3].ticks[i3 + 1] - grid[3].ticks[i3]
    wr, wφ, wz
end

voxel_volume(grid::CylindricalGrid{T}, i1::Int, i2::Int, i3::Int, w1::T, w2::T, w3::T) where {T} =
    (grid[2].ticks[i2 + 1] - grid[2].ticks[i2]) * w3 * (grid[1].ticks[i1 + 1]^2 - grid[1].ticks[i1]^2) / 2
voxel_volume(grid::CartesianGrid3D{T}, i1::Int, i2::Int, i3::Int, w1::T, w2::T, w3::T) where {T} =
    w1 * w2 * w3

function get_extended_midpoints_grid(grid::Grid{T,3}) where {T}
    ticks = broadcast(i -> midpoints(get_extended_ticks(grid.axes[i])), (1,2,3))
    axes = broadcast(i -> typeof(grid.axes[i])(grid.axes[i].interval, ticks[i]) , (1,2,3))
    typeof(grid)(axes)
end


text:
#Get maximum time parameters
function get_t_max( p::Param )

    #Get Variables
    fdata   = "./data"
    fname   = fdata * p.folder *  "$(p.fname).h5"
    fid         = h5open( fname , "r")

    #Max iteration and time
    max_iter = read_attribute(fid, "max_iter")
    max_time = read_attribute(fid, "max_time")

    close(fid)
    return max_iter , max_time

end

# Get displacement at a given iteration
function get_fields( p::Param , itt::Int)

    #Get Variables
    fdata   = "./data"
    fname   = fdata * p.folder *  "$(p.fname).h5"
    fid         = h5open( fname , "r")

    max_iter = read_attribute(fid, "max_iter")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Discretize space
    dx  = (xmax - xmin)/(Nx-1)
    x   = [xmin + (i-1)*dx for i in 1:Nx]

    dy  = (ymax - ymin)/(Ny-1)
    y   = [ymin + (i-1)*dy for i in 1:Ny]

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        x       = [xmin + (i-1)*dx for i in 1:Nx]

        dy      = (ymax - ymin)/(Ny)
        y       = [ymin + (i-1)*dy for i in 1:Ny]
    end

    #Check it time iteration is OK
    if( itt > max_iter)
        println("ERROR: out of bounds (max iteration is" , max_iter,")")
        return
    end

    #Get Displacement
    g1      = open_group(fid, "ψ")

    dname   = string("iteration_" , itt)
    dset    = open_dataset(g1, dname )
    ttt     = read_attribute(dset,"time")

    ψ = reshape(Array(dset) , Nx,Ny )

    #Get Derivative
    g1      = open_group(fid, "dψ")

    dname   = string("iteration_" , itt)
    dset    = open_dataset(g1, dname )

    dψ = reshape(Array(dset),Nx,Ny)

    close(fid)

    #Return all data
    return ttt , x , y , ψ , dψ

end

function get_field_slice(p::Param , coordinate::Symbol , itt::Int , slice_index::Int)

   if( coordinate == :x)

        #Check bounds
        if(slice_index > p.xnodes)
            println("ERROR: Index out of bounds. Max index is " , p.xnodes)
            return
        end

        ttt , x , y , ψ , dψ = get_fields( p , itt)

        return ttt , x , ψ[slice_index,:] , dψ[slice_index,:]

    elseif( coordinate == :y)

    #Check bounds
        if(slice_index > p.ynodes)
            println("ERROR: Index out of bounds. Max index is " , p.ynodes)
            return
        end

        ttt , x , y , ψ , dψ = get_fields( p , itt)

        return ttt , y , ψ[slice_index,:] , dψ[:,slice_index]

    else
        println("ERROR: " , coordinate , "is not a valid argument. Valid arguments are :x and :y")
        return
    end

end

function get_field_point(p::Param, xi_read::Int , yi_read::Int)

    #Check bounds
    if(xi_read > p.xnodes || yi_read > p.ynodes)
        println("ERROR: Index out of bounds. Max index is " , p.xnodes)
        return
    end

    #Max iteration time
    max_sim_iter , _ = get_t_max(p)

    #Storage vectors
    pos_vector = Float64[]
    dpos_vector = Float64[]
    time_vector = Float64[]

    #ReadOnce
    t, x , y , ψ , dψ = get_fields(p,1)

    println(x[xi_read] , "  " , y[yi_read] )

    #Cycle
    for i in 1:max_sim_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(p,i)

        push!(time_vector, t)
        push!(pos_vector, ψ[xi_read , yi_read] )
        push!(dpos_vector, dψ[xi_read , yi_read] )

    end

    #return vectors
    return time_vector , pos_vector , dpos_vector

end

function get_field_point(p::Param, x_read::Float64 , y_read::Float64)

    #Get associated index
    xi_read =floor(Int , p.xnodes * (x_read - p.xmin) / (p.xmax - p.xmin))
    yi_read =floor(Int , p.ynodes * (y_read - p.ymin) / (p.ymax - p.ymin))

    #Check bounds
    if(xi_read > p.xnodes || yi_read > p.ynodes)
        println("ERROR: Index out of bounds. Max index is " , p.xnodes)
        return
    end

    #Max iteration time
    max_sim_iter , _ = get_t_max(p)

    #Storage vectors
    pos_vector = Float64[]
    dpos_vector = Float64[]
    time_vector = Float64[]

    #ReadOnce
    t, x , y , ψ , dψ = get_fields(p,1)

    println(x[xi_read] , "  " , y[yi_read] )

    #Cycle
    for i in 1:max_sim_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(p,i)

        push!(time_vector, t)
        push!(pos_vector, ψ[xi_read , yi_read] )
        push!(dpos_vector, dψ[xi_read , yi_read] )

    end

    #return vectors
    return time_vector , pos_vector , dpos_vector

end


text:
# some functions for dealing with strings
# there ought to be a regular expression way to doing this
function findWord(str::AbstractString, substr::AbstractString)
# returns the entire word containing a substring

substr_len = length(substr)
end_index = 0  # location of end of substring
start_index = 0  # starting index of word
for i=1:(length(str) - substr_len + 1)
  tmp = str[i:(i + substr_len - 1)]  # get the current substring

  if tmp == substr
    end_index = i + substr_len - 1
    start_index = getWordStart(str, end_index, ' ') + 1
    end_index = getWordStart(str, end_index, '/') - 1  # remove libmds.so from end
  end
end

return str[start_index:end_index]

end

function getWordStart(str::AbstractString, end_index::Integer, delim::Char)
# look backwards in the string for a space

  start_index = 1  # return 1 if no space found
  for i=end_index:(-1):1
    char_i = str[i]
#    println("index = ", i, " , ", "char = ", char_i, " , val = ", Int32(char_i))
    if char_i == delim
      start_index = i
      break
    end

  end

  return start_index
end


text:
## FrameGeom.jl -- for geometrical description of Frame contents

export phi, psi, ramachandranangles, omega
export kappa, alpha, hsspangles

const EXTRA_ANGLE_DOCS = """
Returns `nothing` if cannot be defined.

Use the chain dictionary `chaindict` as reference. This can be supplied if
previously calculated, so there'd be no need to calculate for each index.
Defaults to the builtin `getchaindict` function.

**Note**: It is recommended to use _only_ when all indices `i` can be found in
the chain dictionary. Furthermore, it is also recommended to make sure chains
are those of amino acid residues.
"""

"""
    phi(frame, i[, chaindict]) ∈ [-pi, pi], Nothing

Get the phi (``ϕ``) angle of the residue in the frame pointed by index `i`.
$EXTRA_ANGLE_DOCS
"""
function phi(fr::StructureFrame, i, chaindict=getchaindict(fr))
    ## get chain of index i; check if previous is a protein residue
    ## in the same chain (using chaindict)
    chain_ = _getchain(i, chaindict)
    ## totally skip if
    isnothing(chain_) && return nothing

    j = i-1     ## index of previous residue; check if in same chain
    j ∈ chaindict[chain_] || return nothing

    (; res_list, at_pos) = fr
    resi = res_list[i]; resj = res_list[j]

    ## carbonyl carbon of prev residue
    Cprevid = getatom(resj, :C)
    ## N, CA, C of pointed residue
    Nid, CAid, Cid = getatom(resi, :N), getatom(resi, :CA), getatom(resi, :C)

    return broadcast((Cprevid, Nid, CAid, Cid)) do index
        getindex(at_pos, index)
    end |> q->dihedral(q...)
end

"""
    psi(frame, i[, chaindict]) ∈ [-pi, pi], Nothing

Get the psi (``ψ``) angle of the residue in the frame pointed by index `i`.
$EXTRA_ANGLE_DOCS
"""
function psi(fr::StructureFrame, i, chaindict=getchaindict(fr))
    ## get chain of index i; check if previous is a protein residue
    ## in the same chain (using chaindict)
    chain_ = _getchain(i, chaindict)
    ## totally skip if
    isnothing(chain_) && return nothing

    j = i+1     ## index of next residue; check if in same chain
    j ∈ chaindict[chain_] || return nothing

    (; res_list, at_pos) = fr
    resi = res_list[i]; resj = res_list[j]

    ## N, CA, C of pointed residue
    Nid, CAid, Cid = getatom(resi, :N), getatom(resi, :CA), getatom(resi, :C)
    ## N of next residue
    Nnextid = getatom(resj, :N)

    return broadcast((Nid, CAid, Cid, Nnextid)) do index
        getindex(at_pos, index)
    end |> q->dihedral(q...)
end

"""
    omega(frame, i[, chaindict]) ∈ [-pi, pi], Nothing

Get the omega (``ω``) angle of the residue in the frame pointed by index `i`.
$EXTRA_ANGLE_DOCS
"""
function omega(fr::StructureFrame, i, chaindict=getchaindict(fr))
    ## get chain of index i; check if previous is a protein residue
    ## in the same chain (using chaindict)
    chain_ = _getchain(i, chaindict)
    ## totally skip if
    isnothing(chain_) && return nothing

    j = i-1     ## index of previous residue; check if in same chain
    j ∈ chaindict[chain_] || return nothing

    (; res_list, at_pos) = fr
    resi = res_list[i]; resj = res_list[j]

    CAprevid, Cprevid = getatom(resj, :CA), getatom(resj, :C)
    Nid, CAid = getatom(resi, :N), getatom(resi, :CA)

    return broadcast((CAprevid, Cprevid, Nid, CAid)) do index
        getindex(at_pos, index)
    end |> q->dihedral(q...)
end

"""
    kappa(frame, i[, chaindict]) ∈ [0, pi], Nothing

Get the kappa (``κ``) angle of the residue in the fram pointed by index `i`.
$EXTRA_ANGLE_DOCS
"""
function kappa(fr::StructureFrame, i, chaindict=getchaindict(fr))
    ## TODO: get definition from DSSP or HSSP references
    chain_ = _getchain(i, chaindict)
    ## totally skip if
    isnothing(chain_) && return nothing

    (i-2 ∉ chaindict[chain_] || i+2 ∉ chaindict[chain_]) && return nothing

    (; res_list, at_pos) = fr

    idcs = (i-2, i, i+2)

    return map(idcs) do index
        res = res_list[index]
        getatom(res, :CA) |> m->at_pos[m]
    end |> q->anglespan(q...)
end

"""
    alpha(frame, i[, chaindict]) ∈ [-pi, pi], Nothing
    chirality(frame, i, chaindict) ∈ (-1, 0, 1), Nothing

Get the alpha (``α``) angle of the residue in the frame pointed by index `i`.
Chirality (in HSSP) is defined as the sign of alpha.
$EXTRA_ANGLE_DOCS
"""
function alpha(fr::StructureFrame, i, chaindict=getchaindict(fr))
    ## get chain of residue pointed to by i
    chain_ = _getchain(i, chaindict)
    ## totally skip if
    isnothing(chain_) && return nothing

    (i-1 ∉ chaindict[chain_] || i+2 ∉ chaindict[chain_]) && return nothing

    (; res_list, at_pos) = fr

    idcs = ntuple(j->i+j-2, 4)

    return map(idcs) do index
        res = res_list[index]
        getatom(res, :CA) |> m->at_pos[m]
    end |> q->dihedral(q...)
end
## what is chirality when alpha is ±π
@doc (@doc alpha) function chirality(fr, i, chaindict=getchaindict(fr))
    α = alpha(fr, i, chaindict)
    isnothing(α) && return nothing
    return sign(α)
end

"""
    ramachandranangles(frame, i[, chaindict]) -> Tuple
    ProtStructs.ϕψ(frame, i[, chaindict])

Get the Ramachandran angles of the residue in the frame pointed by index `i`.
Uses `phi` and `psi` functions as defined. Returns a tuple based on that order.
Shorthand above is provided, but not exported.

$EXTRA_ANGLE_DOCS
"""
function ramachandranangles(fr::StructureFrame, i, chaindict=getchaindict(fr))
    return phi(fr, i, chaindict), psi(fr, i, chaindict)
end
@doc (@doc ramachandranangles) const ϕψ = ramachandranangles

"""
    hsspangles(frame, i[, chaindict]) -> Tuple
    ProtStructs.κα(frame, i[, chaindict])

Get `kappa` and `alpha` angles for the residue in the frame pointed by index `i`.
These angles are used by HSSP in searching for homologous structures.
Shorthand above is provided, but not exported.

$EXTRA_ANGLE_DOCS
"""
function hsspangles(fr::StructureFrame, i, chaindict=getchaindict(fr))
    return kappa(fr, i, chaindict), alpha(fr, i, chaindict)
end
@doc (@doc hsspangles) const κα = hsspangles

@noinline function _getchain(i, cdict)
    ## does a linear search
    ## TODO: better search methods (max 20 elements ig)
    for (chain, range) in cdict
        i ∈ range && return chain
    end
    return nothing
end


text:
# Matern 1/2 ARD covariance Function

@doc """
# Description
Constructor for the ARD Matern 1/2 kernel (covariance)

k(x,x') = σ²exp(-d/L^2), where d = |x-x'| and L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Mat12Ard <: StationaryARD
    ℓ2::Vector{Float64}      # Log of length scale
    σ2::Float64              # Log of signal std
    Mat12Ard(ll::Vector{Float64}, lσ::Float64) = new(2.0*exp(ll),exp(2.0*lσ))
end

function set_params!(mat::Mat12Ard, hyp::Vector{Float64})
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat12 kernel only has $(num_params(mat)) parameters"))
    d = length(mat.ℓ2)
    mat.ℓ2  = exp(2.0*hyp[1:d])
    mat.σ2 = exp(2.0*hyp[d+1])
end

get_params(mat::Mat12Ard) = [log(mat.ℓ2)/2.0; log(mat.σ2)/2.0]
get_param_names(mat::Mat12Ard) = [get_param_names(mat.ℓ2, :ll); :lσ]
num_params(mat::Mat12Ard) = length(mat.ℓ2) + 1

metric(mat::Mat12Ard) = WeightedEuclidean(1.0./(mat.ℓ2))
cov(mat::Mat12Ard, r::Float64) = mat.σ2*exp(-r)

function grad_kern(mat::Mat12Ard, x::Vector{Float64}, y::Vector{Float64})
    r = distance(mat, x, y)
    exp_r = exp(-r)
    wdiff = ((x-y).^2)./mat.ℓ2

    g1 = (mat.σ2*wdiff*exp_r)/r
    g2 = 2.0*mat.σ2*exp_r

    return [g1; g2]
end

function grad_stack!(stack::AbstractArray, mat::Mat12Ard, X::Matrix{Float64}, data::StationaryARDData)
    d = size(X,1)
    R = distance(mat, X, data)
    stack[:,:,d+1] = mat.σ2 * exp(-R)
    ck = view(stack, :, :, d+1)
    broadcast!(*, view(stack, :, :, 1:d), data.dist_stack, reshape(1.0./mat.ℓ2, (1,1,d)), ck)
    broadcast!(/, view(stack, :, :, 1:d), view(stack, :, :, 1:d), R)
    stack[:,:, d+1] = 2.0 * ck
    return stack
end


text:
const DIR = [SVector( 0, 1), # north
             SVector( 1, 0), # east
             SVector( 0,-1), # south
             SVector(-1, 0)] # west

const PING = 5
const ENGAGE = 6

const ACT_LETTER = ['N', 'E', 'S', 'W', 'P', 'G']

const END_KILL = SubState([-1,-1], [-1,-1], -1, false)

# values for aware
# const ISLEFT = 1
# const ISRIGHT = 2

# """
# Return direction index from `from` to `to`
# """
# function dirto(from::Pos, to::Pos)
#
# end
#
# function leftright(from::Pos, dir::Int, to::Pos)
#
# end
# =#


text:
# f(A,x,b)=A*x+b
# Note that the bias (a column vector) gets expanded here to match the dimension of A*x

FAXplusBias(A,X,b)=(A*X+b*ones(1,size(X,2)),nothing)

function FAXplusBias_inplace(value,aux,A,X,b)
    copy!(value,A*X+b*ones(1,size(X,2)))
end

function DAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X,b)
    if derivativeIDX==1
        axpy!(1.0,grad_c*X',grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,A'*grad_c,grad_n)
    elseif derivativeIDX==3
        axpy!(1.0,sum(grad_c,2),grad_n)
    end
end

if PROC=="GPU"
#    function FAXplusBias(A::CudaArray,X::CudaArray,b::CudaArray)
#        ons=CudaArray(Float64,(1,size(X,2))); fill!(ons,1.0);
#        #return (Fxpy(FAX(A,X)[1],FAX(b,ons)[1])[1],nothing)
#        #return (Fxpy(A*X,b*ons)[1],nothing)
#        return (A*X+b*ons,nothing)
#    end

    function FAXplusBias_inplace(value,aux,A::CudaArray{Float64},X::CudaArray{Float64},b::CudaArray{Float64})
        ons=CudaArray(Float64,(1,size(X,2))); fill!(ons,1.0)
#        Fxpy_inplace(value,aux,FAX(A,X)[1],FAX(b,ons)[1])
        gemm!('N','N',1.0,b,ons,0.0,value)
        gemm!('N','N',1.0,A,X,1.0,value)
        free(ons)
    end

    function DAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float64},X::CudaArray{Float64},b::CudaArray{Float64})
        if derivativeIDX==1
            gemm!('N','T',1.0,grad_c,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,grad_c,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float64,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,grad_c,ons,1.0,grad_n)
            free(ons)
        end
    end

    function FAXplusBias_inplace(value,aux,A::CudaArray{Float32},X::CudaArray{Float32},b::CudaArray{Float32})
        ons=CudaArray(Float32,(1,size(X,2))); fill!(ons,Float32(1.0))
#        Fxpy_inplace(value,aux,FAX(A,X)[1],FAX(b,ons)[1])
        gemm!('N','N',Float32(1.0),b,ons,Float32(0.0),value)
        gemm!('N','N',Float32(1.0),A,X,Float32(1.0),value)
        free(ons)
    end

    function DAXplusBias(derivativeIDX,f_c,faux_c,grad_c::CudaArray{Float32},grad_n::CudaArray{Float32},A::CudaArray{Float32},X::CudaArray{Float32},b::CudaArray{Float32})
        if derivativeIDX==1
            gemm!('N','T',Float32(1.0),grad_c,X,Float32(1.0),grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',Float32(1.0),A,grad_c,Float32(1.0),grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float32,(size(X,2),1)); fill!(ons,Float32(1.0))
            gemm!('N','N',Float32(1.0),grad_c,ons,Float32(1.0),grad_n)
            free(ons)
        end
    end

end

Derivative[FAXplusBias]=DAXplusBias
Inplace[FAXplusBias]=FAXplusBias_inplace

AXplusBias(A,X,b)=ADnode(FAXplusBias,[A X b])
export AXplusBias

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function jointplot(; x::AbstractVector = Int[], y::AbstractVector = Int[], smoothing::Int = 4)

    @assert sum(Missings.ismissing.(x)) == 0 "x has missing data"
    @assert sum(Missings.ismissing.(y)) == 0 "y has missing data"

    #Calculate kde
    a = kde((collect(skipmissing(x)),collect(skipmissing(y))))

    #Parse data from grid to 1D array
    x_mod = []
    y_mod = []
    d_mod = []

    for _x in enumerate(a.x)
        for _y in enumerate(a.y)
            push!(x_mod, _x[2])
            push!(y_mod, _y[2])
            push!(d_mod, a.density[_x[1], _y[1]])
        end
    end

    min_d = minimum(d_mod)
    max_d = maximum(d_mod)
    intervals = (max_d - min_d)/10
    center =  (max_d - min_d)/2

    #Plot
    a = heatmap(x = x_mod, y = y_mod, color = d_mod)
    a.width = a.height = 500
    a.scales[1]._type = a.scales[2]._type = "linear"
    a.marks[1].properties.enter.width = a.marks[1].properties.enter.height = VegaValueRef(value = smoothing)

    a.scales[2].zero = false

    a.scales[1].nice = a.scales[2].nice = true
    a.scales[3].domain =  collect(min_d:intervals:max_d)
    a.legends[1].values = [min_d, center, max_d]

    return a

end


text:
using Documenter
using Intervals

# Including these modules results in printing without the module names
# https://github.com/JuliaDocs/Documenter.jl/issues/944
using Dates
using TimeZones

makedocs(;
    modules=[Intervals],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/invenia/Intervals.jl/blob/{commit}{path}#L{line}",
    sitename="Intervals.jl",
    authors="Invenia Technical Computing Corporation",
    checkdocs=:exports,
    strict=true,
)

deploydocs(;
    repo="github.com/invenia/Intervals.jl",
)


text:
using WignerFamilies
using PyPlot
using OffsetArrays
using BenchmarkTools

w = WignerF(Float64, 4000, 4000, 0, 0)
w3j = get_wigner_array(w)

function test(w_, w3j_)
    for i in 1:2000
        wigner3j_f!(w_, w3j_)
    end
end

@btime test($w, $w3j)

##

function test2(w_, w3j_)
    for i in 1:2000
        WignerFamilies.classical_wigner3j_m0!(w_, w3j_)
    end
end

@btime test2($w, $w3j)

##


text:
mutable struct Nnet
    y::Array{Int64,1}
    X::Array{Float64,2}
    N::Int64
    D::Int64
    hidden_layers::Int64
    nodes::Int64
    weights::Int64
    λ::Array{Float64,2}
    μ::Array{Float64,1}
    σ::Array{Float64,1}
    prior::Distributions.Normal{Float64}
    iter::Int64
end

"Construct Neural Network"
function Nnet(y::Array{Int64, 1},
              X::Array{Float64, 2};
              hidden_layers::Int64=1,
              nodes::Int64=2,
              prior::Distributions.Normal{Float64}=Normal())

    N = length(y)
    D = size(X)[2]

    # still need to account for more than one hidden layer
    weights = D * nodes + nodes + # input to layer 1
              # layer n to layer n + 1
              (hidden_layers - 1) * ((nodes + 1) * nodes)  +
              nodes + 1 # final layer to f(X)

    # initialize variational parameters
    λ = rand(Normal(), weights, 2)
    truncate!(λ, weights)

    μ = λ[:, 1]
    σ = softplus.(λ[:, 2])

    Nnet(y, X, N, D,
         hidden_layers, nodes, weights,
         λ, μ, σ, prior, 0)
end

# estimate of evidence lower bound
function ℒ(nn::Nnet, μ, σ)
    # extract elements from Nnet
    y = nn.y
    X = nn.X
    prior = nn.prior
    hidden_layer = nn.hidden_layers
    nodes = nn.nodes

    # draw values from q
    ϵ = rand(Normal(), length(μ))
    z = μ + ϵ .* σ

    # initialize vectors for calculating ELBO
    # note that the type is important for AD
    log_lik = zeros(eltype(μ), 1)
    log_prior = zeros(eltype(μ), length(μ))
    entropy = zeros(eltype(μ), length(μ))
    f = zeros(eltype(μ), length(μ))

    # get number of covariates
    N, D = size(X)

    # store layer calculations
    layer = zeros(eltype(μ), nodes)
    layer_p1 = zeros(eltype(μ), nodes)
    out = zeros(eltype(μ), N)

    for n in 1:N
        out[n] = fwd_prop(X[n, :],
                          nodes, hidden_layer,
                          z,
                          layer, layer_p1)

        log_lik += logpdf(Bernoulli(out[n]), y[n])
    end

    # construct ELBO estimate
    for i in 1:length(μ)
        # first calculate log prior and entropy
        log_prior[i] = logpdf(nn.prior, z[i])
        entropy[i] = logpdf(Normal(μ[i], σ[i]), z[i])
        f[i] = log_lik[1] + log_prior[i] - entropy[i]
    end

    return f
end

function fit(nn::Nnet,
             opt::GradDescent.Optimizer;
             tol=0.005,
             maxiter=1000,
             msg=100,
             elboiter=10)
    N = length(nn.y)
    D = nn.weights

    mean_δ = 1.0
    i = 1

    while mean_δ > tol && i < maxiter
        nn.μ = nn.λ[:, 1]
        nn.σ = softplus.(nn.λ[:, 2])

        g = ForwardDiff.jacobian(ϕ -> ℒ(nn, ϕ[1:D], ϕ[(D+1):2D]),
                                 vcat(nn.μ, nn.σ))
        ∇ℒ = hcat(diag(g[1:D, 1:D]),
                  diag(g[1:D, (D+1):2D]))

        δ = update(opt, ∇ℒ)
        nn.λ += δ
        truncate!(nn.λ, D)

        mean_δ = mean(δ .^ 2)

        if (i % msg == 0)
            println(i, " ", mean_δ)
        end

        i += 1
    end

    # calculate elbo
    nn.μ = nn.λ[:, 1]
    nn.σ = softplus.(nn.λ[:, 2])
    elbo = 0

    for i in 1:elboiter
        elbo += mean(ℒ(nn, nn.μ, nn.σ))
    end

    elbo = elbo / elboiter

    return elbo
end

function predict(nn::Nnet,
                 X::Array{Float64, 2})
    # extract elements from Nnet
    hidden_layer = nn.hidden_layers
    nodes = nn.nodes

    # get number of covariates
    N, D = size(X)

    # store layer calculations
    layer = zeros(nodes)
    layer_p1 = zeros(nodes)
    out = zeros(N)

    # for right now, fix z at mean
    z = nn.μ

    for n in 1:N
        out[n] = fwd_prop(X[n, :],
                          nodes, hidden_layer,
                          z,
                          layer, layer_p1)
    end

    return out
end


text:
using Graphs

typealias TermVertex KeyVertex{Term}

import Base: ==, isequal

isequal(tv1::TermVertex, tv2::TermVertex) = isequal(tv1.key, tv2.key)
==(tv1::TermVertex, tv2::TermVertex) = isequal(tv1, tv2)

type GOGraph
    ids::Vector{TermID}
    upward::IncidenceList{TermVertex,Edge{TermVertex}}
    downward::IncidenceList{TermVertex,Edge{TermVertex}}

    function GOGraph(filepath::String; mode::Int=2)
        ids, upward, downward = gograph_builder(filepath, mode)
        new(ids, upward, downward)
    end
end

type ISAVisitor <: AbstractGraphVisitor
    target::TermVertex
    found::Bool

    function ISAVisitor(target::TermVertex)
        new(target, false)
    end
end

function Graphs.discover_vertex!(visitor::ISAVisitor, v::TermVertex)
    # if discovered the target, immediately finish traversing
    if v == visitor.target
        visitor.found = true
        return false
    end
    true
end

# check whether term1 *is a* term2 or not
function is_a(go::GOGraph, term1::Term, term2::Term)
    id = binary_search(go.ids, term1.id)
    v = TermVertex(id, term1)
    pid = binary_search(go.ids, term2.id)
    pv = TermVertex(pid, term2)
    visitor = ISAVisitor(pv)
    traverse_graph(go.upward, DepthFirst(), v, visitor)
    visitor.found
end

function binary_search(ids::Vector{TermID}, id::TermID)
    isempty(ids) && return 0
    l = 1
    n = endof(ids)
    while n >= l
        m = div(n - l, 2) + l
        if id < ids[m]
            n = m - 1
        elseif id > ids[m]
            l = m + 1
        else
            return m
        end
    end
    0
end

function search_vertex(go::GOGraph, term::Term)
    id = binary_search(go.ids, term.id)
    id == 0 && error("term $term does not exist in the ontology graph")
    TermVertex(id, term)
end

function gograph_builder(filepath::String, mode::Int)
    parser = OBOParser(filepath, mode=mode)
    terms = Term[]

    # accumulate all terms
    for term in eachterm(parser)
        push!(terms, term)
    end

    sort!(terms)
    upward = inclist(TermVertex)
    downward = inclist(TermVertex)
    ids = TermID[]

    # add term vertices
    for (i, term) in enumerate(terms)
        v = KeyVertex{Term}(i, term)
        add_vertex!(upward, term)
        add_vertex!(downward, term)
        push!(ids, term.id)
    end

    # add edges between two terms
    for (i, term) in enumerate(terms)
        v = TermVertex(i, term)
        for pid in term.isa
            ip = binary_search(ids, pid)
            @assert ip != 0
            parent = terms[ip]
            vp = TermVertex(ip, parent)
            add_edge!(upward, v, vp)
            add_edge!(downward, vp, v)
        end
    end

    ids, upward, downward
end

function parents(go::GOGraph, term::Term)
    v = search_vertex(go, term)
    Term[p.key for p in out_neighbors(v, go.upward)]
end

function children(go::GOGraph, term::Term)
    v = search_vertex(go, term)
    Term[p.key for p in out_neighbors(v, go.downward)]
end


text:
using Weber
using Base.Test
include("find_timing.jl")

import Weber: prepare!, handle

struct TestPrepareMoment <: Weber.SimpleMoment
  event::Symbol
end
function handle(exp::Weber.Experiment,queue::Weber.MomentQueue,
                moment::TestPrepareMoment,x)
  Weber.dequeue!(queue)
  true
end

prepare!(moment::TestPrepareMoment) = record(moment.event)

ks,vs,_ = find_timing() do
  addtrial(TestPrepareMoment(:a),moment(0.5s),moment(),moment(),
           TestPrepareMoment(:b),moment(record,:b_post),moment(0.5s))
  addtrial(moment(),TestPrepareMoment(:c),moment(record,:c_post),
           timeout(() -> nothing,iskeydown,0.5s),TestPrepareMoment(:d))
  addtrial(moment(100ms),TestPrepareMoment(:e) >> TestPrepareMoment(:f))
end
prepare_timing = Dict(k => v for (k,v) in zip(ks,vs))

expanding_prepare,_,_ = find_timing() do
  @addtrials let test = 1
    addtrial(TestPrepareMoment(:a))
  end
end

struct TestPrepareError <: Weber.SimpleMoment end
function handle(exp::Weber.Experiment,queue::Weber.MomentQueue,
                moment::TestPrepareError,x)
  Weber.dequeue!(queue)
  true
end
struct TestPrepareException <: Exception end
prepare!(moment::TestPrepareError,time::Float64) =
  isinf(time) ? throw(TestPrepareException()) : record(:success)

function cause_prepare_error1()
  find_timing() do
    addtrial(moment(0.5s),timeout(() -> nothing,iskeydown,0.5s),TestPrepareError())
  end
end

prepare_noerror,_,_ = find_timing() do
  addtrial(moment(0.5s),TestPrepareError())
end

@testset "Moment Preparation" begin
  @test :a in keys(prepare_timing)
  @test :d in keys(prepare_timing)
  @test :e in keys(prepare_timing)
  @test :f in keys(prepare_timing)
  @test :a in expanding_prepare
  @test abs(prepare_timing[:b_post] - prepare_timing[:b] - 0.5) < 0.25
  @test abs(prepare_timing[:c_post] - prepare_timing[:c] - 0.5) < 0.25
  @test_throws TestPrepareException cause_prepare_error1()
  @test prepare_noerror == [:success]
end


text:
# паковка и распаковка массивов в байтовые массивы и обратно

"""
pack: vector of numbers => byte vector
"""
pack_vec(vect::Vector{T}) where T <: Number = reinterpret(UInt8, vect) |> collect

"""
pack: tuple of vectors => byte vector
* сейчас сделано с промежуточной аллокацией памяти
"""
# This works because, unlike for any other type,
# tuples are covariant in their subtypes
pack_vec(vecs::Tuple{Vararg{Vector}}) = zip(vecs...) |> collect |> pack_vec

"""
pack: vector of tuples => byte vector
"""
@generated function pack_vec(vect::Vector{T}) where T <: Tuple
    types = fieldtypes(T)
    exprs = Expr[]
    k = 0
    for T in types
        k += 1
        ex = quote
            let
                p::Ptr{$T} = pointer(result, ind)
                unsafe_store!(p, x[$k], 1)
                ind += sizeof($T)
            end
        end
        push!(exprs, ex)
    end
    loop_unroll = Expr(:block, exprs...)

    out_expr = quote
        elsize = reduce(+, sizeof.($types))
        len = length(vect) * elsize
        result = Vector{UInt8}(undef, len)

        GC.@preserve result begin
            ind::Int = 1
            @inbounds for x in vect
                $loop_unroll
            end
        end
        return result
    end

    return out_expr
end

# inner function to get type inside generated function
_gettype(::Type{Type{T}}) where {T} = T

"""
unpack: byte vector => vector of numbers
"""
unpack_vec(vect::Vector{UInt8}, type::Type) = reinterpret(type, vect) |> collect

"""
unpack: byte vector => vector of tuples
"""
@generated function unpack_vec(vect::Vector{UInt8}, types...)
    exprs = Expr[]
    for T in types
        ex = quote
            let
                p::Ptr{_gettype($T)} = pointer(vect, ind)
                x = unsafe_load(p)
                ind += sizeof(_gettype($T))
                x
            end
        end
        push!(exprs, ex)
    end
    loop_unroll = :(tuple($(exprs...)))

    out_expr = quote
        elsize = reduce(+, sizeof.(types))
        len = sizeof(vect) ÷ elsize
        result = Vector{Tuple{types...}}(undef, len)

        GC.@preserve vect begin
            ind::Int = 1
            @inbounds for i in 1:len
                result[i] = $loop_unroll
            end
        end
        return result
    end

    return out_expr
end

"""
unpack: byte vector => tuple of vectors
* такого пока что нет - надо преобразовывать отдельно в StructArray
"""


text:
module MultiIndices

import Base: sub2ind, size, show, getindex

# * Abstract spaces

abstract type Space{N,O} end

dims(mi::Space{N,O}) where {N,O} = N

order(mi::Space{N,O}) where {N,O} = O
function order(𝔄::Space{N𝔄,O𝔄},𝔅::Space{N𝔅,O𝔅}) where {N𝔄,O𝔄,N𝔅,O𝔅}
    mO𝔄 = maximum(O𝔄)
    (O𝔄..., (o+mO𝔄 for o in O𝔅)...)
end

Base.getindex(mi::Space, args...) = sub2ind(mi, args...)

function sub2ind(mi::Space, args...)
    s = size(mi)
    O = order(mi)
    sub2ind(tuple((s[o] for o in O)...),
            (args[o] for o in O)...)
end

# ** Sum space

struct SumSpace{S<:Space{N,O} where {N,O}} <: Space{1,(1,)}
    𝔖::Vector{S}
end

dims(mi::SumSpace{S}) where S <: Space{N,O} where {N,O} = 1
order(mi::SumSpace{S}) where S <: Space{N,O} where {N,O} = (1,)

function size(mi::SumSpace)
    sizes = map(mi.𝔖) do s
        prod(size(s))
    end
    sum(sizes)
end

getindex(mi::MultiIndices.SumSpace, i::Integer, args...) = mi.𝔖[i][args...] # This does not work as intended, returns inner dimension index only
getindex(mi::MultiIndices.SumSpace, i::Integer) = mi.𝔖[i]

"""
    𝔄 ⊕ 𝔅

Construct the sum space of which `𝔄` and `𝔅` are terms.
"""
⊕(𝔄::S, 𝔅::S) where S <: Space =
    SumSpace([𝔄, 𝔅])

"""
    𝔖 ⊕ 𝔄

Construct the sum space of which `𝔖` and `𝔄` are terms.
"""
⊕(𝔖::SumSpace{S}, 𝔄::S) where S <: Space =
    SumSpace([𝔖.𝔖..., 𝔄])

# ** Product space

struct ProductSpace{N,O} <: Space{2,(1,2)} where {N,O}
    𝔄::Space
    𝔅::Space
end
function ProductSpace(𝔄::Space,𝔅::Space)
    N = dims(𝔄) + dims(𝔅)
    O = order(𝔄, 𝔅)
    ProductSpace{N,O}(𝔄, 𝔅)
end

dims(mi::ProductSpace{N,O}) where {N,O} = N
order(mi::ProductSpace{N,O}) where {N,O} = O

function size(mi::ProductSpace)
    (size(mi.𝔄)...,size(mi.𝔅)...)
end

"""
    𝔄 ⊗ 𝔅

Construct the product space of which `𝔄` and `𝔅` are factors. `𝔄` will
be the "fast dimension".
"""
⊗(𝔄::Space, 𝔅::Space) = ProductSpace(𝔄, 𝔅)

# ** Copy space

struct CopySpace{N,O} <: Space{2,(1,2)} where {N,O}
    𝔖::Space
    n::Integer
end
function CopySpace(𝔖::Space, n::Integer)
    N = dims(𝔖) + 1
    O = order(𝔖)
    O = (O..., length(O)+1)
    CopySpace{N,O}(𝔖, n)
end
dims(mi::CopySpace{N,O}) where {N,O} = N
order(mi::CopySpace{N,O}) where {N,O} = O
size(mi::CopySpace) = (size(mi.𝔖)..., mi.n)

"""
    𝔖 ⊗ n

Construct the product space from `N` copies of `𝔖`. `𝔖` will
be the "fast dimension".
"""
⊗(𝔖::Space, n::Integer) = CopySpace(𝔖, n)

# * Concrete spaces
# ** Cartesian coordinate systems

struct Cartesian{N,O} <: Space{N,O}
    dims::Tuple
end

size(mi::Cartesian) = mi.dims

Cartesian(i::Integer) = Cartesian{1,(1,)}((i,))
Cartesian{N,O}(args...) where {N,O} = Cartesian{N,O}(args)

CartesianXY = Cartesian{2, (1,2)}
CartesianYX = Cartesian{2, (2,1)}

CartesianXYZ = Cartesian{3, (1,2,3)}
CartesianXZY = Cartesian{3, (1,3,2)}

CartesianYXZ = Cartesian{3, (2,1,3)}
CartesianYZX = Cartesian{3, (2,3,1)}

CartesianZXY = Cartesian{3, (3,1,2)}
CartesianZYX = Cartesian{3, (3,2,1)}

# ** Polar coordinates

struct Polar{O} <: Space{2,O}
    nr::Integer
    nϕ::Integer
end
size(mi::Polar) = (mi.nr, mi.nϕ)

# ** Spherical coordinates

struct Spherical{O} <: Space{3,O}
    nr::Integer
    nθ::Integer
    nϕ::Integer
end
size(mi::Spherical) = (mi.nr, mi.nθ, mi.nϕ)

# *** Spherical harmonics

# For expansions in spherical harmonics, given an ℓ, m ∈
# [-ℓ,ℓ]. However, it may be desired to truncate the expansion not
# only in ℓ, but also in m. The "fast dimension" should be m.

struct SphericalHarmonics <: Space{2,(2,1)}
    nℓ::Integer
    nm::Integer
end
SphericalHarmonics(nℓ) = SphericalHarmonics(nℓ, 2nℓ-1)

# It would be nice if ℓ was zero-based, and m from -ℓ–ℓ. Maybe
# "custom" indexing can be implemented like this?
function (mi::SphericalHarmonics)(ℓ::Integer, m::Integer)
    ℓ ∉ 0:mi.nℓ-1 && error("Invalid value of ℓ = $(ℓ) ∉ [0,$(mi.nℓ-1)]")
    abs(m) > ℓ && error("Invalid value of m = $(m) ∉ [-$(ℓ),$(ℓ)]")
    2abs(m) > mi.nm && error("Invalid value of |m| = $(abs(m)) > $((mi.nm-1)/1) (max m)")

end

# This function actually returns the *maximum* nm, that is nm for the
# largest ℓ. The expansion is not "rectangular", but rather "triangular":
#=

              ℓ

    0   1    2    3    4   ...

                      -4
                 -3   -3
            -2   -2   -2
       -1   -1   -1   -1
m   0   0    0    0    0
        1    1    1    1
             2    2    2
                  3    3
                       4
=#
size(mi::SphericalHarmonics) = (mi.nℓ, mi.nm)

struct RotationSpace <: Space{1, (1,)}
    ℓ::Integer
    mmax::Integer
    RotationSpace(ℓ::Integer, mmax::Integer) = new(ℓ, min(ℓ, mmax))
end
RotationSpace(ℓ::Integer) = RotationSpace(ℓ, ℓ)
size(mi::RotationSpace) = 2mi.mmax + 1

ℓs = "spdfghiklmnoqrtuvwxyz"
ℓm(ℓ::Integer, m::Integer) = (ℓ+1, ℓ+m+1)
ℓm(ℓ::Char, m::Integer) = ℓm(search(ℓs, ℓ)-1, m)

function show(io::IO, mi::RotationSpace)
    write(io, "Rot{ℓ = ")
    if mi.ℓ + 1 < length(ℓs)
        write(io, "$(ℓs[mi.ℓ+1])")
    else
        write(io, "$(mi.ℓ)")
    end
    if mi.mmax < mi.ℓ
        write(io, ", mmax=$(mi.mmax)")
    end
    write(io, "}")
end

# This coordinate system is actually proper 3d, but in cylindrically
# symmetric problems (e.g. linear polarization), the inital m quantum
# number is conserved, and the computational effort is 2d.

struct SphericalHarmonics2d <: Space{1, (1,)}
    nℓ::Integer
    m₀::Integer
end
SphericalHarmonics2d(nℓ) = SphericalHarmonics2d(nℓ,0)

# Spin angular momentum; given a spin magnitude s, the projection can
# take values in the interval mₛ ∈ [-s,s].

# ** Spin

struct Spin{S} <: Space{1, (1,)} end
size(mi::Spin{S}) where S = (Integer(2S+1),)

# * Exports

export Space,
    size, getindex, sub2ind, dims, order,
    ⊕, ⊗,
    Cartesian, CartesianXY, CartesianYX,
    CartesianXYZ, CartesianXZY, CartesianYXZ,
    CartesianYZX, CartesianZXY, CartesianZYX,
    Polar, Spherical,
    SphericalHarmonics, SphericalHarmonics2d,
    RotationSpace, ℓm,
    Spin

# * End

end # module


text:
using ExaTron
using Random
using LinearAlgebra
using SparseArrays
using CUDA
using CUDA.CUSPARSE

using Profile
using PProf

function build_problem(; n=10, gpu=true)
    Random.seed!(1)
    m = 0
    P = sparse(Diagonal(rand(n)) + 2.0 * sparse(I, n, n))
    q = randn(n)
    u =   1. * rand(n)
    l = - 100. * rand(n)
    Iz, Jz, vals = findnz(P)

    if gpu
        P = CuSparseMatrixCSC(P)
        q = CuArray(q)
        u = CuArray(u)
        l = CuArray(l)
    end

    eval_f(x) = 0.5 * x' * P * x + q' * x

    function eval_g(x, g)
        fill!(g, 0)
        mul!(g, P, x)
        g .+= q
    end

    function eval_h(x, mode, rows, cols, obj_factor, lambda, values)
        if mode == :Structure
            @inbounds for i in 1:nnz(P)
                rows[i] = Iz[i]
                cols[i] = Jz[i]
            end
        else
            copy!(values, vals)
        end
    end

    return ExaTron.createProblem(n, l, u, nnz(P), eval_f, eval_g, eval_h)
end

prob = build_problem(; n=1000)
prob.x .= 0.5 .* (prob.x_l .+ prob.x_u)
@profile ExaTron.solveProblem(prob)
pprof()

# CUDA.@profile ExaTron.solveProblem(prob)


text:
using Base.Test
using DataFrames

let
	#
	# NA's
	#

	@assert isna(NAtype())
	@assert isna(NA)

	#
	# DataVector's
	#

	dv = DataArray([1, 2, 3], falses(3))
	@assert isequal(dv.data, [1, 2, 3])
	@assert isequal(dv.na, falses(3))
	@assert isequal(dv, DataArray([1, 2, 3], [false, false, false]))
	@assert isequal(dv, DataArray([1, 2, 3]))

	dv = DataArray(trues(3), falses(3))
	@assert isequal(dv.data, [true, true, true])
	@assert isequal(dv.na, falses(3))
	@assert isequal(dv, DataArray(trues(3)))

	dv = DataArray([1, 2, 3], falses(3))
	@assert isequal(dv, DataArray(1:3))
	@assert isequal(dv, DataArray(DataArray([1, 2, 3])))

	dv = DataArray(Int, 3)
	@assert isequal(eltype(dv), Int)
	@assert isequal(dv.na, trues(3))

	# dv = DataArray(3)
	# @assert isequal(eltype(dv), Float64)
	# @assert isequal(dv.na, trues(3))

	# dv = DataArray()
	# @assert isequal(eltype(dv), Float64)
	# @assert isequal(dv.na, trues(0))

	@assert isequal(datazeros(3), DataArray(zeros(3)))
	@assert isequal(datazeros(Int, 3), DataArray(zeros(Int, 3)))
	@assert isequal(dataones(3), DataArray(ones(3)))
	@assert isequal(dataones(Int, 3), DataArray(ones(Int, 3)))
	@assert isequal(datafalses(3), DataArray(falses(3)))
	@assert isequal(datatrues(3), DataArray(trues(3)))

	dv = DataVector[1, 2, NA]
	@assert dv[1] == 1
	@assert dv[2] == 2
	@assert isna(dv[3])
	@assert isequal(eltype(dv), Int)

	#
	# PooledDataArray's
	#

	pdv = PooledDataArray([1, 2, 3], falses(3))
	@assert all(pdv .== [1, 2, 3])
	@assert all(isna(pdv) .== falses(3))

	@assert isequal(pdv, PooledDataArray([1, 2, 3], [false, false, false]))
	@assert isequal(pdv, PooledDataArray([1, 2, 3]))

	pdv = PooledDataArray(trues(3), falses(3))
	@assert all(pdv .== [true, true, true])
	@assert all(isna(pdv) .== falses(3))
	@assert isequal(pdv, PooledDataArray(trues(3)))

	pdv = PooledDataArray([1, 2, 3], falses(3))
	@assert isequal(pdv, PooledDataArray(1:3))
	@assert isequal(pdv, PooledDataArray(PooledDataArray([1, 2, 3])))

	pdv = PooledDataArray(Int, 3)
	@assert isequal(eltype(pdv), Int)
	@assert all(isna(pdv) .== trues(3))

	# pdv = PooledDataArray(3)
	# @assert isequal(eltype(pdv), Float64)
	# @assert all(isna(pdv) .== trues(3))

	# pdv = PooledDataArray()
	# @assert isequal(eltype(pdv), Float64)
	# @assert all(isna(pdv) .== trues(0))

	@assert isequal(pdatazeros(3), PooledDataArray(zeros(3)))
	@assert isequal(pdatazeros(Int, 3), PooledDataArray(zeros(Int, 3)))
	@assert isequal(pdataones(3), PooledDataArray(ones(3)))
	@assert isequal(pdataones(Int, 3), PooledDataArray(ones(Int, 3)))
	@assert isequal(pdatafalses(3), PooledDataArray(falses(3)))
	@assert isequal(pdatatrues(3), PooledDataArray(trues(3)))

	pdv = PooledDataVector[1, 2, NA]
	@assert pdv[1] == 1
	@assert pdv[2] == 2
	@assert isna(pdv[3])
	@assert isequal(eltype(pdv), Int)

	#
	# DataMatrix
	#

	dm = DataArray([1 2; 3 4], falses(2, 2))
	@assert isequal(dm.data, [1 2; 3 4])
	@assert isequal(dm.na, falses(2, 2))

	@assert isequal(dm, DataArray([1 2; 3 4], [false false; false false]))
	@assert isequal(dm, DataArray([1 2; 3 4]))

	dm = DataArray(trues(2, 2), falses(2, 2))
	@assert isequal(dm.data, trues(2, 2))
	@assert isequal(dm.na, falses(2, 2))

	@assert isequal(dm, DataArray(trues(2, 2)))

	#DataMatrix(dvzeros(3), dvzeros(3))
	#DataMatrix(1:3, 1:3)

	@assert isequal(DataArray([1 2; 3 4]), DataArray(DataArray([1 2; 3 4])))

	dm = DataArray(Int, 2, 2)
	@assert isequal(eltype(dm), Int)
	@assert isequal(dm.na, trues(2, 2))

	# dm = DataArray(2, 2)
	# @assert isequal(eltype(dm), Float64)
	# @assert isequal(dm.na, trues(2, 2))

	# dm = DataArray(Int)
	# @assert isequal(eltype(dm), Int)
	# @assert isequal(dm.na, trues(0, 0))

	# dm = DataArray()
	# @assert isequal(eltype(dm), Float64)
	# @assert isequal(dm.na, trues(0, 0))

	@assert isequal(datazeros(2, 2), DataArray(zeros(2, 2)))
	@assert isequal(datazeros(Int, 2, 2), DataArray(zeros(Int, 2, 2)))

	@assert isequal(dataones(2, 2), DataArray(ones(2, 2)))
	@assert isequal(dataones(Int, 2, 2), DataArray(ones(Int, 2, 2)))

	@assert isequal(datafalses(2, 2), DataArray(falses(2, 2)))
	@assert isequal(datatrues(2, 2), DataArray(trues(2, 2)))

	@assert isequal(dataeye(3, 2), DataArray(eye(3, 2)))
	@assert isequal(dataeye(2), DataArray(eye(2)))
	@assert isequal(datadiagm(Float64[pi, pi]), DataArray(diagm(Float64[pi, pi])))

	#
	# DataFrame
	#

	df = DataFrame()
	@assert isequal(df.columns, {})
	# TODO: Get this to work
	#@assert isequal(df.colindex, Index())

	df = DataFrame({datazeros(3), dataones(3)}, Index(["x1", "x2"]))
	@assert nrow(df) == 3
	@assert ncol(df) == 2

	# TODO: Make isequal fail if colnames don't match
	@assert isequal(df, DataFrame({datazeros(3), dataones(3)}))
	@assert isequal(df, DataFrame(quote x1 = [0.0, 0.0, 0.0]; x2 = [1.0, 1.0, 1.0] end))

	@assert isequal(df, DataFrame([0.0 1.0; 0.0 1.0; 0.0 1.0], ["x1", "x2"]))
	@assert isequal(df, DataFrame([0.0 1.0; 0.0 1.0; 0.0 1.0]))
	@assert isequal(df, DataFrame(datazeros(3), dataones(3)))

	# TODO: Fill these in
	# From (Associative): ???
	# From (Vector, Vector, Groupings): ???

	@assert isequal(df, DataFrame({"x1" => [0.0, 0.0, 0.0], "x2" => [1.0, 1.0, 1.0]}))
	@assert isequal(df, DataFrame({"x1" => [0.0, 0.0, 0.0], "x2" => [1.0, 1.0, 1.0], "x3" => [2.0, 2.0, 2.0]}, ["x1", "x2"]))

	df = DataFrame(Int, 2, 2)
	@assert size(df) == (2, 2)
	@assert all(coltypes(df) .== {Int, Int})
	@assert all(isna(df))

	df = DataFrame(2, 2)
	@assert size(df) == (2, 2)
	@assert all(coltypes(df) .== {Float64, Float64})
	@assert all(isna(df))

	df = DataFrame({Int, Float64}, ["x1", "x2"], 2)
	@assert size(df) == (2, 2)
	@assert all(coltypes(df) .== {Int, Float64})
	@assert all(isna(df))

	@assert isequal(df, DataFrame({Int, Float64}, 2))
end