text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Check that the QuadraticProgrammingProblem is valid.
"""
function validate(p::QuadraticProgrammingProblem)
  error_found = false
  if length(p.variable_lower_bound) != length(p.variable_upper_bound)
    @info "$(length(p.variable_lower_bound)) == length(p.variable_lower_bound)
      != length(p.variable_upper_bound) = $(length(p.variable_upper_bound))"
    error_found = true
  end
  if length(p.variable_lower_bound) != length(p.objective_vector)
    @info "$(length(p.variable_lower_bound)) == length(p.variable_lower_bound)
      != length(p.objective_vector)= $(length(p.objective_vector))"
    error_found = true
  end
  if size(p.constraint_matrix, 1) != length(p.right_hand_side)
    @info "$(size(p.constraint_matrix,1)) == size(p.constraint_matrix,1)
      != length(p.right_hand_side)= $(length(p.right_hand_side))"
    error_found = true
  end
  if size(p.constraint_matrix, 2) != length(p.objective_vector)
    @info "$(size(p.constraint_matrix,2)) == size(p.constraint_matrix,2)
      != length(p.objective_vector)= $(length(p.objective_vector))"
    error_found = true
  end
  if size(p.objective_matrix) !=
     (length(p.objective_vector), length(p.objective_vector))
    @info "$(size(p.objective_matrix)) == size(p.objective_matrix)
      is not square with length $(length(p.objective_vector))"
    error_found = true
  end
  if sum(p.variable_lower_bound .== Inf) > 0
    @info "sum(p.variable_lower_bound .== Inf) ==
            $(sum(p.variable_lower_bound .== Inf))"
    error_found = true
  end
  if sum(p.variable_upper_bound .== -Inf) > 0
    @info "sum(p.variable_upper_bound .== -Inf) ==
            $(sum(p.variable_lupper_bound .== -Inf))"
    error_found = true
  end
  if any(isnan, p.variable_lower_bound) || any(isnan, p.variable_upper_bound)
    @info "NaN found in variable bounds of QuadraticProgrammingProblem."
    error_found = true
  end
  if any(!isfinite, p.right_hand_side)
    @info "NaN or Inf found in right hand side of QuadraticProgrammingProblem."
    error_found = true
  end
  if any(!isfinite, p.objective_vector)
    @info "NaN or Inf found in objective vector of QuadraticProgrammingProblem."
    error_found = true
  end
  if any(!isfinite, nonzeros(p.constraint_matrix))
    @info "NaN or Inf found in constraint matrix of
            QuadraticProgrammingProblem."
    error_found = true
  end
  if any(!isfinite, nonzeros(p.objective_matrix))
    @info "NaN or Inf found in objective matrix of QuadraticProgrammingProblem."
    error_found = true
  end

  if error_found
    error("Error found when validating QuadraticProgrammingProblem. See @info
    for details.")
  end

  return true
end

"""
Returns the l2 norm of each row or column of a matrix. The method rescales
the sum-of-squares computation by the largest absolute value if nonzero in order
to avoid overflow.

# Arguments
- `matrix::SparseMatrixCSC{Float64, Int64}`: a sparse matrix.
- `dimension::Int64`: the dimension we want to compute the norm over. Must be
  1 or 2.

# Returns
An array with the l2 norm of a matrix over the given dimension.
"""
function l2_norm(matrix::SparseMatrixCSC{Float64,Int64}, dimension::Int64)
  scale_factor = vec(maximum(abs, matrix, dims = dimension))
  scale_factor[iszero.(scale_factor)] .= 1.0
  if dimension == 1
    scaled_matrix = matrix * Diagonal(1 ./ scale_factor)
    return scale_factor .*
           vec(sqrt.(sum(t -> t^2, scaled_matrix, dims = dimension)))
  end

  if dimension == 2
    scaled_matrix = Diagonal(1 ./ scale_factor) * matrix
    return scale_factor .*
           vec(sqrt.(sum(t -> t^2, scaled_matrix, dims = dimension)))
  end
end

"""
Removes the empty rows of a quadratic programming problem.

# Arguments
- `problem::QuadraticProgrammingProblem`: The input quadratic programming
  problem. This is modified to store the transformed problem.
"""
function remove_empty_rows(problem::QuadraticProgrammingProblem)
  seen_row = falses(size(problem.constraint_matrix, 1))
  for row in SparseArrays.rowvals(problem.constraint_matrix)
    seen_row[row] = true
  end
  empty_rows = findall(.!seen_row)

  for row in empty_rows
    if row > problem.num_equalities && problem.right_hand_side[row] > 0.0
      error("The problem is infeasible.")
    elseif row <= problem.num_equalities && problem.right_hand_side[row] != 0.0
      error("The problem is infeasible.")
    end
  end

  if !isempty(empty_rows)
    problem.constraint_matrix = problem.constraint_matrix[seen_row, :]
    problem.right_hand_side = problem.right_hand_side[seen_row]
    num_empty_equalities = sum(empty_rows .<= problem.num_equalities)
    problem.num_equalities -= num_empty_equalities
  end
  return empty_rows
end

"""
Removes the empty columns of a quadratic programming problem.

# Arguments
- `problem::QuadraticProgrammingProblem`: The input linear programming problem.
  This is modified to store the transformed problem. `objective_constant` is
  updated assuming that the eliminated variables are assigned to the best
  possible values.
"""
function remove_empty_columns(problem::QuadraticProgrammingProblem)

  # TODO: Adapt the implementation for quadratic objectives.
  @assert iszero(problem.objective_matrix)
  is_empty_column = [
    isempty(nzrange(problem.constraint_matrix, col)) for
    col in 1:size(problem.constraint_matrix, 2)
  ]
  empty_columns = findall(is_empty_column)
  if isempty(empty_columns)
    return empty_columns
  end

  for col in empty_columns
    objective_coef = problem.objective_vector[col]
    if objective_coef >= 0
      problem.objective_constant +=
        problem.variable_lower_bound[col] * objective_coef
    else
      problem.objective_constant +=
        problem.variable_upper_bound[col] * objective_coef
    end
  end
  is_non_empty = .!is_empty_column
  problem.constraint_matrix = problem.constraint_matrix[:, is_non_empty]
  problem.objective_vector = problem.objective_vector[is_non_empty]
  problem.variable_lower_bound = problem.variable_lower_bound[is_non_empty]
  problem.variable_upper_bound = problem.variable_upper_bound[is_non_empty]
  problem.objective_matrix =
    problem.objective_matrix[is_non_empty, is_non_empty]
  return empty_columns
end

"""
Modifies qp by transforming any finite variable bounds into linear constraints.
"""
function transform_bounds_into_linear_constraints(
  qp::FirstOrderLp.QuadraticProgrammingProblem,
)
  finite_lower_bound_indices = findall(isfinite.(qp.variable_lower_bound))
  finite_upper_bound_indices = findall(isfinite.(qp.variable_upper_bound))

  row_indices =
    1:(length(finite_lower_bound_indices)+length(finite_upper_bound_indices))
  column_indices = [finite_lower_bound_indices; finite_upper_bound_indices]
  nonzeros = [
    ones(length(finite_lower_bound_indices))
    -ones(length(finite_upper_bound_indices))
  ]

  identity_block = SparseArrays.sparse(
    row_indices,
    column_indices,
    nonzeros,
    length(row_indices),
    length(qp.variable_lower_bound),
  )
  qp.constraint_matrix = [qp.constraint_matrix; identity_block]
  qp.right_hand_side = [
    qp.right_hand_side
    qp.variable_lower_bound[finite_lower_bound_indices]
    -qp.variable_upper_bound[finite_upper_bound_indices]
  ]
  qp.variable_lower_bound .= -Inf
  qp.variable_upper_bound .= Inf
  return
end

struct PresolveInfo
  original_primal_size::Int64
  original_dual_size::Int64
  empty_rows::Vector{Int64}
  empty_columns::Vector{Int64}
  variable_lower_bound::Vector{Float64}
  variable_upper_bound::Vector{Float64}
end

"""
Presolve of a quadratic program. The quadratic program modified in place.
Returns a data structure with information allowing the presolve to be `undone`.
"""
function presolve(
  qp::QuadraticProgrammingProblem;
  verbosity::Int64 = 1,
  transform_bounds::Bool = false,
)
  saved_variable_lower_bound = copy(qp.variable_lower_bound)
  saved_variable_upper_bound = copy(qp.variable_upper_bound)

  original_dual_size, original_primal_size = size(qp.constraint_matrix)
  empty_rows = FirstOrderLp.remove_empty_rows(qp)
  # TODO: Remove this check after remove_empty_columns supports
  # quadratic objectives.
  if iszero(qp.objective_matrix)
    empty_columns = FirstOrderLp.remove_empty_columns(qp)
  else
    empty_columns = Array{Int64,1}()
  end

  if verbosity >= 1
    FirstOrderLp.check_for_singleton_constraints(qp)
  end
  # TODO: Write function check_for_singleton_variables(lp)

  if transform_bounds
    FirstOrderLp.transform_bounds_into_linear_constraints(qp)
  end

  return PresolveInfo(
    original_primal_size,
    original_dual_size,
    empty_rows,
    empty_columns,
    saved_variable_lower_bound,
    saved_variable_upper_bound,
  )
end

function check_for_singleton_constraints(qp::QuadraticProgrammingProblem)
  num_single = 0
  num_nz_by_row = zeros(Int64, size(qp.constraint_matrix, 1))
  for row_ind in SparseArrays.rowvals(qp.constraint_matrix)
    num_nz_by_row[row_ind] += 1
  end

  num_single = sum(num_nz_by_row .== 1)
  if num_single > 0
    println("$num_single constraints involving exactly a single variable")
    # TODO: Eliminate constraints involving a single variable by
    # replacing with bounds variables or by fixing the value.
  end
end

"""
Given a solution to the preprocessed problem this function recovers a solution
to the original problem.

# Arguments
- `solution::Vector{Float64}`: The solution after preprocessing.
- `empty_indicies::Vector{Int64}`: Indicies corresponding to portions of the
solution that were eliminated in preprocessing. When we recover the original
solution these indicies will have zero values.
- `original_size::Int64`: Size of the solution vector before preprocessing.
"""
function recover_original_solution(
  solution::Vector{Float64},
  empty_indicies::Vector{Int64},
  original_size::Int64,
)
  # construct a list of nonempty indicies
  nonempty_indicies_binary_list = trues(original_size)
  nonempty_indicies_binary_list[empty_indicies] .= false
  nonempty_indicies = findall(nonempty_indicies_binary_list)

  original_solution = zeros(original_size)
  original_solution[nonempty_indicies] = solution[1:length(nonempty_indicies)]

  return original_solution
end

function undo_presolve(
  presolve_info::PresolveInfo,
  primal_solution::Vector{Float64},
  dual_solution::Vector{Float64},
)
  # recover primal solution
  primal_solution = FirstOrderLp.recover_original_solution(
    primal_solution,
    presolve_info.empty_columns,
    presolve_info.original_primal_size,
  )

  projection!(
    primal_solution,
    presolve_info.variable_lower_bound,
    presolve_info.variable_upper_bound,
  )

  # recover dual solution
  dual_solution = FirstOrderLp.recover_original_solution(
    dual_solution,
    presolve_info.empty_rows,
    presolve_info.original_dual_size,
  )
  return primal_solution, dual_solution
end

"""
Rescales a quadratic programming problem by dividing each row and column of the
constraint matrix by the sqrt its respective L2 norm, adjusting the other
problem data accordingly.

# Arguments
- `problem::QuadraticProgrammingProblem`: The input quadratic programming
  problem. This is modified to store the transformed problem.

# Returns

A tuple of vectors `constraint_rescaling`, `variable_rescaling` such that
the original problem is recovered by
`unscale_problem(problem, constraint_rescaling, variable_rescaling)`.
"""
function l2_norm_rescaling(problem::QuadraticProgrammingProblem)
  num_constraints, num_variables = size(problem.constraint_matrix)

  norm_of_rows = vec(l2_norm(problem.constraint_matrix, 2))
  norm_of_columns = vec(l2_norm(problem.constraint_matrix, 1))

  norm_of_rows[iszero.(norm_of_rows)] .= 1.0
  norm_of_columns[iszero.(norm_of_columns)] .= 1.0

  column_rescale_factor = sqrt.(norm_of_columns)
  row_rescale_factor = sqrt.(norm_of_rows)
  scale_problem(problem, row_rescale_factor, column_rescale_factor)

  return row_rescale_factor, column_rescale_factor
end

"""
Uses a modified Ruiz rescaling algorithm to rescale the matrix M=[Q,A';A,0]
where Q is objective_matrix and A is constraint_matrix, and returns the
cumulative scaling vectors. More details of Ruiz rescaling algorithm can be
found at: http://www.numerical.rl.ac.uk/reports/drRAL2001034.pdf.

In the p=Inf case, both matrices approach having all row and column LInf norms
of M equal to 1 as the number of iterations goes to infinity. This convergence
is fast (linear).

In the p=2 case, the goal is all row L2 norms of [Q,A'] equal to 1 and all row
L2 norms of A equal to sqrt(num_variables/(num_constraints+num_variables))
for QP, and all row L2 norms of A equal to
sqrt(num_variables/num_constraints) for LP. Having a different
goal for the row and col norms is required since the sum of squares of the
entries of the A matrix is the same when the sum is grouped by rows or grouped
by columns. In particular, for the LP case, all L2 norms of A must be
sqrt(num_variables/num_constraints) when all row L2 norm of [Q,A'] equal to 1.

The Ruiz rescaling paper (link above) only analyzes convergence in the p < Inf
case when the matrix is square, and it does not preserve the symmetricity of
the matrix, and that is why we need to modify it for p=2 case.

TODO: figure out when this converges.

# Arguments
- `problem::QuadraticProgrammingProblem`: the quadratic programming problem.
  This is modified to store the transformed problem.
- `num_iterations::Int64` the number of iterations to run Ruiz rescaling
  algorithm. Must be positive.
- `p::Float64`: which norm to use. Must be 2 or Inf.

# Returns

A tuple of vectors `constraint_rescaling`, `variable_rescaling` such that
the original problem is recovered by
`unscale_problem(problem, constraint_rescaling, variable_rescaling)`.
"""
function ruiz_rescaling(
  problem::QuadraticProgrammingProblem,
  num_iterations::Int64,
  p::Float64 = Inf,
)
  num_constraints, num_variables = size(problem.constraint_matrix)
  cum_constraint_rescaling = ones(num_constraints)
  cum_variable_rescaling = ones(num_variables)

  for i in 1:num_iterations
    constraint_matrix = problem.constraint_matrix
    objective_matrix = problem.objective_matrix

    if p == Inf
      variable_rescaling = vec(
        sqrt.(
          max.(
            maximum(abs, constraint_matrix, dims = 1),
            maximum(abs, objective_matrix, dims = 1),
          ),
        ),
      )
    else
      @assert p == 2
      variable_rescaling = vec(
        sqrt.(
          sqrt.(
            l2_norm(constraint_matrix, 1) .^ 2 +
            l2_norm(objective_matrix, 1) .^ 2,
          ),
        ),
      )
    end
    variable_rescaling[iszero.(variable_rescaling)] .= 1.0

    if num_constraints == 0
      constraint_rescaling = Float64[]
    else
      if p == Inf
        constraint_rescaling =
          vec(sqrt.(maximum(abs, constraint_matrix, dims = 2)))
      else
        @assert p == 2
        norm_of_rows = vec(l2_norm(problem.constraint_matrix, 2))

        # If the columns all have norm 1 and the row norms are equal they should
        # equal sqrt(num_variables/num_constraints) for LP, and roughly
        # sqrt(num_variables / (num_constraints + num_variables) for QP.
        if iszero(problem.objective_matrix)
          target_row_norm = sqrt(num_variables / num_constraints)
        else
          target_row_norm =
            sqrt(num_variables / (num_constraints + num_variables))
        end
        constraint_rescaling = vec(sqrt.(norm_of_rows / target_row_norm))
      end
      constraint_rescaling[iszero.(constraint_rescaling)] .= 1.0
    end
    scale_problem(problem, constraint_rescaling, variable_rescaling)

    cum_constraint_rescaling .*= constraint_rescaling
    cum_variable_rescaling .*= variable_rescaling
  end

  return cum_constraint_rescaling, cum_variable_rescaling
end

"""
Applies the rescaling proposed by Pock and Cambolle (2011),
"Diagonal preconditioning for first order primal-dual algorithms
in convex optimization"
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.381.6056&rep=rep1&type=pdf

Although presented as a form of diagonal preconditioning, it can be
equivalently implemented by rescaling the problem data.

Each column of the constraint matrix is divided by
sqrt(sum_{elements e in the column} |e|^(2 - alpha))
and each row of the constraint matrix is divided by
sqrt(sum_{elements e in the row} |e|^alpha)

Lemma 2 in Pock and Chambolle demonstrates that this rescaling causes the
operator norm of the rescaled constraint matrix to be less than or equal to
one, which is a desireable property for PDHG.

# Arguments
- `problem::QuadraticProgrammingProblem`: the quadratic programming problem.
  This is modified to store the transformed problem.
- `alpha::Float64`: the exponent parameter. Must be in the interval [0, 2].

# Returns

A tuple of vectors `constraint_rescaling`, `variable_rescaling` such that
the original problem is recovered by
`unscale_problem(problem, constraint_rescaling, variable_rescaling)`.
"""
function pock_chambolle_rescaling(
  problem::QuadraticProgrammingProblem,
  alpha::Float64,
)
  @assert 0 <= alpha <= 2

  constraint_matrix = problem.constraint_matrix

  variable_rescaling = vec(
    sqrt.(
      mapreduce(
        t -> abs(t)^(2 - alpha),
        +,
        constraint_matrix,
        dims = 1,
        init = 0.0,
      ),
    ),
  )
  constraint_rescaling = vec(
    sqrt.(
      mapreduce(t -> abs(t)^alpha, +, constraint_matrix, dims = 2, init = 0.0),
    ),
  )

  variable_rescaling[iszero.(variable_rescaling)] .= 1.0
  constraint_rescaling[iszero.(constraint_rescaling)] .= 1.0

  scale_problem(problem, constraint_rescaling, variable_rescaling)

  return constraint_rescaling, variable_rescaling
end

"""
Rescales `problem` in place. If we let `D = diag(cum_variable_rescaling)` and
`E = diag(cum_constraint_rescaling)`, then `problem` is modified such that:

    objective_matrix = D^-1 objective_matrix D^-1
    objective_vector = D^-1 objective_vector
    objective_constant = objective_constant
    variable_lower_bound = D variable_lower_bound
    variable_upper_bound = D variable_upper_bound
    constraint_matrix = E^-1 constraint_matrix D^-1
    right_hand_side = E^-1 right_hand_side

The scaling vectors must be positive.
"""
function scale_problem(
  problem::QuadraticProgrammingProblem,
  constraint_rescaling::Vector{Float64},
  variable_rescaling::Vector{Float64},
)
  @assert all(t -> t > 0, constraint_rescaling)
  @assert all(t -> t > 0, variable_rescaling)
  problem.objective_vector ./= variable_rescaling
  problem.objective_matrix =
    Diagonal(1 ./ variable_rescaling) *
    problem.objective_matrix *
    Diagonal(1 ./ variable_rescaling)
  problem.variable_upper_bound .*= variable_rescaling
  problem.variable_lower_bound .*= variable_rescaling
  problem.right_hand_side ./= constraint_rescaling
  problem.constraint_matrix =
    Diagonal(1 ./ constraint_rescaling) *
    problem.constraint_matrix *
    Diagonal(1 ./ variable_rescaling)
  return
end

"""
Recovers the original problem from the scaled problem and the scaling vectors
in place. The inverse of `scale_problem`. This function should be only used for
testing.
"""
function unscale_problem(
  problem::QuadraticProgrammingProblem,
  constraint_rescaling::Vector{Float64},
  variable_rescaling::Vector{Float64},
)
  scale_problem(problem, 1 ./ constraint_rescaling, 1 ./ variable_rescaling)
  return
end

"""
    row_permute_in_place(matrix::SparseMatrixCSC{Float64, Int64},
                         old_row_to_new::Vector{Int64})

Permutes the rows of `matrix` in place (without allocating a new matrix)
according to the map `old_row_to_new`. Assumes and does not verify that
`old_row_to_new` is a permutation.
"""
function row_permute_in_place(
  matrix::SparseMatrixCSC{Float64,Int64},
  old_row_to_new::Vector{Int64},
)
  coefficients = SparseArrays.nonzeros(matrix)
  row_indices = SparseArrays.rowvals(matrix)
  row_coef_tuples = Tuple{Int64,Float64}[]
  for col in 1:size(matrix, 2)
    nonzero_range = nzrange(matrix, col)
    empty!(row_coef_tuples)
    sizehint!(row_coef_tuples, length(nonzero_range))
    for index_in_matrix in nonzero_range
      new_row = old_row_to_new[row_indices[index_in_matrix]]
      push!(row_coef_tuples, (new_row, coefficients[index_in_matrix]))
    end
    # SparseMatrixCSC requires row indices to be sorted within a column. So we
    # do this sort and then replace the terms in the column.
    sort!(row_coef_tuples, by = t -> t[1])
    for (index_in_column, index_in_matrix) in enumerate(nonzero_range)
      new_row, new_coefficient = row_coef_tuples[index_in_column]
      row_indices[index_in_matrix] = new_row
      coefficients[index_in_matrix] = new_coefficient
    end
  end
  return
end

"""Preprocesses the original problem, and returns a ScaledQpProblem struct.
Applies L_inf Ruiz rescaling for `l_inf_ruiz_iterations` iterations. If
`l2_norm_rescaling` is true, applies L2 norm rescaling. `problem` is not
modified.
"""
function rescale_problem(
  l_inf_ruiz_iterations::Int,
  l2_norm_rescaling::Bool,
  pock_chambolle_alpha::Union{Float64,Nothing},
  verbosity::Int64,
  original_problem::QuadraticProgrammingProblem,
)
  problem = deepcopy(original_problem)
  if verbosity >= 4
    println("Problem before rescaling:")
    print_problem_details(original_problem)
  end

  num_constraints, num_variables = size(problem.constraint_matrix)
  constraint_rescaling = ones(num_constraints)
  variable_rescaling = ones(num_variables)

  if l_inf_ruiz_iterations > 0
    con_rescale, var_rescale =
      FirstOrderLp.ruiz_rescaling(problem, l_inf_ruiz_iterations, Inf)
    constraint_rescaling .*= con_rescale
    variable_rescaling .*= var_rescale
  end

  if l2_norm_rescaling
    con_rescale, var_rescale = FirstOrderLp.l2_norm_rescaling(problem)
    constraint_rescaling .*= con_rescale
    variable_rescaling .*= var_rescale
  end

  if !isnothing(pock_chambolle_alpha)
    con_rescale, var_rescale =
      FirstOrderLp.pock_chambolle_rescaling(problem, pock_chambolle_alpha)
    constraint_rescaling .*= con_rescale
    variable_rescaling .*= var_rescale
  end

  scaled_problem = ScaledQpProblem(
    original_problem,
    problem,
    constraint_rescaling,
    variable_rescaling,
  )

  if verbosity >= 3
    if l_inf_ruiz_iterations == 0 && !l2_norm_rescaling
      println("No rescaling.")
    else
      print("Problem after rescaling ")
      print("(Ruiz iterations = $l_inf_ruiz_iterations, ")
      println("l2_norm_rescaling = $l2_norm_rescaling):")
      print_problem_details(scaled_problem.scaled_qp)
    end
  end

  return scaled_problem
end


text:
# Example from "Bounds on the reachable sets of nonlinear control systems" Constant Temperature Profile
id = :ScottB2013b

k = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
pL = k
pU = 10.0*k
ScottB2013b_x0(p::Vector{T}) where T = [34.0*one(T); 20.0*one(T)]
function ScottB2013b_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T

    R = 8.134
    A1 = 2400
    A2 = 8800
    E1 = 6.9*10^3
    E2 = 1.69*10^4

    du[1] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[6]*u[6]
    du[2] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[3]*u[3]
    du[3] = p[1]*u[1]*u[2] - p[2]*u[3] - p[3]*u[3]
    du[4] = p[3]*u[3] - p[4]*u[1]*u[6] + p[5]*u[6]
    du[5] = -p[4]*u[4]*u[5] + p[5]*u[6] + p[6]*u[6]
    du[6] = p[4]*u[4]*u[5] - p[5]*u[6] - p[6]*u[6]
    return
end
tspan = (0.0, 0.04)

prob = ODERelaxProb(ScottB2013b_f!, tspan, ScottB2013b_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/pii/S0005109812004839"
source = "ScottB2013"
desc = "3 Series Rxn Network"

ScottB2013b_id = :ScottB2013b
ScottB2013b_desc = "3 Series Rxn Network"

ext: jl


text:
"""
for defining Lattice, Nondispersive, and Dispersive Domains
"""
module Domains

files = (
    "1D/Domains1D.jl",
    "2D/Domains2D.jl",
    # "3D/Domains3D.jl"
)

export AbstractDomain
export LatticeDomain
export NondispersiveDomain
export DispersiveDomain

using ..Boundaries
using ..DielectricFunctions
using ..PumpFunctions
using ..Dispersions
using ..Lattices
using ..Points
using ..Shapes
using RecipesBase

import ..Symmetric, ..Unsymmetric
import LinearAlgebra.norm

"""
    AbstractDomain{N}
"""
abstract type AbstractDomain{N} end

################################################################################
# LATTICE DOMAIN

struct LatticeDomain{N,CLASS,TC,TBND,TLAT,TP} <: AbstractDomain{N}
    boundary::TBND
    lattice::TLAT
    n::ComplexF64
    ε::ComplexF64
    type::Symbol
    name::Symbol
    x::Vector{Point{N,TP}}
    indices::Vector{CartesianIndex{N}}
    interior::BitArray{1}
    bulk::BitArray{1}
    surface::BitArray{1}
    corner::BitArray{1}
    nnm::NTuple{N,Vector{Int}}
    nnp::NTuple{N,Vector{Int}}

    function LatticeDomain{CLASS}(
        boundary::TBND,
        lattice::TLAT,
        n::Number,
        x::Vector{Point{N,TP}},
        indices::Vector{CartesianIndex{N}},
        interior::BitArray{1},
        bulk::BitArray{1},
        surface::BitArray{1},
        corner::BitArray{1},
        nnm::NTuple{N,Vector{Int}},
        nnp::NTuple{N,Vector{Int}};
        type::Symbol = :generic,
        name::Symbol = :anonymous
        ) where {TBND<:Boundary{N},TP,TLAT<:Lattice{N,TC},CLASS} where {N,TC}

        new{N,CLASS,TC,TBND,TLAT,TP}(boundary,lattice,n,n^2,type,name,x,indices,interior,bulk,surface,corner,nnm,nnp)
    end
end

LatticeDomain(lattice::Lattice, boundary::Boundary, n::Number=1; kwargs...) =
    LatticeDomain(boundary, lattice, n; kwargs...)

LatticeDomain(boundary::Boundary, n::Number, lattice::Lattice; kwargs...) =
    LatticeDomain(boundary, lattice, n; kwargs...)

LatticeDomain(lattice::Lattice, n::Number, boundary::Boundary; kwargs...) =
    LatticeDomain(boundary, lattice, n; kwargs...)

LatticeDomain(n::Number, boundary::Boundary, lattice::Lattice; kwargs...) =
    LatticeDomain(boundary, lattice, n; kwargs...)

LatticeDomain(n::Number, lattice::Lattice, boundary::Boundary; kwargs...) =
    LatticeDomain(boundary, lattice, n; kwargs...)

function Base.getproperty(ldom::LatticeDomain, sym::Symbol)
    if sym == :shape
        return getfield(getfield(ldom,:boundary),:shape)
    elseif Base.sym_in(sym, propertynames(getfield(ldom,:lattice)))
        return getproperty(getfield(ldom,:lattice),sym)
    else
        return getfield(ldom,sym)
    end
end

################################################################################
# NONDISPERSIVE DOMAIN

struct NondispersiveDomain{N,TSH,TDF} <: AbstractDomain{N}
    shape::TSH
    dielectric::TDF
    type::Symbol
    name::Symbol

    function NondispersiveDomain(
                shape::TSH,
                dielectric::TDF,
                type::Symbol = :generic,
                name::Symbol = :anonymous
                ) where {TSH<:AbstractShape{N}, TDF<:AbstractDielectricFunction} where N

        new{N,TSH,TDF}(shape,dielectric,type,name)
    end
end

NondispersiveDomain(shape::AbstractShape, n::Number=1, args...) =
    NondispersiveDomain(shape,DielectricFunctions.PiecewiseConstant(n), args...)

NondispersiveDomain(shape::AbstractShape, n1::Number, n2::Number, args...) =
    NondispersiveDomain(shape,DielectricFunctions.PiecewiseConstant(n1,n2), args...)

NondispersiveDomain(boundary::Boundary, args...) = NondispersiveDomain(boundary.shape, args...)

Base.propertynames(ndom::NondispersiveDomain, private=false) = (fieldnames(NondispersiveDomain)..., propertynames(ndom.dielectric,private)...)

function Base.getproperty(ndom::NondispersiveDomain, sym::Symbol)
    dielectric = getfield(ndom,:dielectric)
    if Base.sym_in(sym,propertynames(dielectric,true))
        return getproperty(dielectric,sym)
    else
        return getfield(ndom,sym)
    end
end

function Base.setproperty!(ndom::NondispersiveDomain, sym::Symbol, val)
    if Base.sym_in(sym,propertynames(ndom.dielectric,true))
        return setproperty!(ndom.dielectric,sym,val)
    else
        return setfield!(ndom,sym,val)
    end
end

################################################################################
# DISPERSIVE DOMAIN

struct DispersiveDomain{N,TSH,TCHI,TPF} <: AbstractDomain{N}
    shape::TSH
    χ::TCHI
    pump::TPF
    type::Symbol
    name::Symbol

    function DispersiveDomain(
                shape::TSH,
                χ::TCHI,
                pump::TPF;
                type::Symbol = :generic,
                name::Symbol = :anonymous
                ) where {TSH<:AbstractShape{N}, TPF<:AbstractPumpFunction, TCHI<:AbstractDispersion} where N

        new{N,TSH,TCHI,TPF}(shape,χ,pump,type,name)
    end
end

DispersiveDomain(shape::AbstractShape, χ::AbstractDispersion, F::Number=1; kwargs...) =
    DispersiveDomain(shape, χ, PumpFunctions.PiecewiseConstant(F); kwargs...)

DispersiveDomain(boundary::Boundary, args...; kwargs...) = DispersiveDomain(boundary.shape, args...; kwargs...)

Base.propertynames(ddom::DispersiveDomain, private=false) = (fieldnames(DispersiveDomain)..., propertynames(ddom.pump,private)...)

function Base.getproperty(ddom::DispersiveDomain, sym::Symbol)
    pump = getfield(ddom,:pump)
    if Base.sym_in(sym,propertynames(pump,true))
        return getproperty(pump,sym)
    else
        return getfield(ddom,sym)
    end
end

function Base.setproperty!(ddom::DispersiveDomain, sym::Symbol, val)
    if Base.sym_in(sym,propertynames(ddom.pump,true))
        return setproperty!(ddom.pump,sym,val)
    else
        return setfield!(ddom,sym,val)
    end
end

################################################################################
# load dimensional files (mostly lattice domain)
foreach(include,files)

################################################################################
# Pretty Printing

import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE

function Base.show(io::IO,dom::LatticeDomain{N,CLASS}) where {N,CLASS}
    print(io,"$(N)D ")
    CLASS<:Symmetric ? print(io,"Symmetric ") : nothing
    CLASS<:Unsymmetric ? print(io,"Unsymmetric ") : nothing
    printstyled(io,"LatticeDomain ",color=PRINTED_COLOR_DARK)
    println(io,"(",dom.type,"): ",dom.name)
    print(io,"\tNumber of sites: ")
    printstyled(io,length(dom.x),"\n",color=PRINTED_COLOR_NUMBER)
    print(io,"\tNumber of bulk points: ")
    printstyled(io,sum(dom.bulk),"\n",color=PRINTED_COLOR_NUMBER)
    print(io,"\tNumber of surface points: ")
    printstyled(io,sum(dom.surface),"\n",color=PRINTED_COLOR_NUMBER)
    print(io,"\tNumber of corner points: ")
    printstyled(io,sum(dom.corner),"\n",color=PRINTED_COLOR_NUMBER)
    print(io,"\tBackground Index: ")
    printstyled(io,dom.n;color=PRINTED_COLOR_NUMBER)
    println(io)
    print(io,"\t\t====================")
    println(io)
    println(IOContext(io,:tabbed2=>true),dom.lattice)
    println(io)
    print(IOContext(io,:tabbed2=>true),dom.boundary)
end

function Base.show(io::IO,ndom::NondispersiveDomain)
    printstyled(io,"NondispersiveDomain ",color=PRINTED_COLOR_DARK)
    println(io,"(",ndom.type,"): ",ndom.name)
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    println(io,"\t\t",ndom.shape)
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    show(IOContext(io,:tabbed2=>true), ndom.dielectric)
end

function Base.show(io::IO,ddom::DispersiveDomain)
    printstyled(io,"DispersiveDomain ",color=PRINTED_COLOR_DARK)
    println(io,"(",ddom.type,"): ",ddom.name)
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    println(io,"\t\t",ddom.shape)
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    show(IOContext(io,:tabbed2=>true),ddom.pump)
    println(io)
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    println(io,"\t\t",ddom.χ)
end

end # module


text:
import LazySets: dim, inner
using LazySets: isapproxzero

"""
    AbstractDirections{N, VN}

Abstract type for template direction representations.

### Notes

This type is parameterzed by `N` and `VN`, where:

- `N` stands for the numeric type
- `VN` stands for the vector type with coefficients of type `N`

Each subtype is an iterator over a set of prescribed directions.

All subtypes should implement the standard iterator methods from `Base`, namely
`Base.length` (returns the number of directions in the template), and `Base.iterate`.
Moreover, the following methods should be implemented:

- `dim`    -- return the ambient dimension of the template
- `eltype` -- return the type of each vector in the template

Optionally, subtypes may implement:

- `isbounding`   -- (defaults to `false`) return `true` if an overapproximation with
                    a list of template directions results in a bounded set, given a
                    bounded input set, and `false` otherwise
- `isnormalized` -- (defaults to `false`) returns `true` if each direction in the
                    given template has norm one w.r.t. the usual vector 2-norm
"""
abstract type AbstractDirections{N, VN} end

"""
    dim(ad::AbstractDirections)

Returns the dimension of the generated directions.

### Input

- `ad` -- template directions

### Output

The dimension of the generated directions.
"""
function dim(ad::AbstractDirections) end

"""
    isbounding(ad::Type{<:AbstractDirections})

Checks if an overapproximation with a list of template directions results in a
bounded set, given a bounded input set.

### Input

- `ad` -- template directions

### Output

Given a bounded set ``X``, we can construct an outer approximation of ``X`` by
using the template directions `ad` as normal vectors of the facets.
If this function returns `true`, then the result is again a bounded set (i.e., a
polytope).
Note that the result does not depend on the specific shape of ``X``, as long as
``X`` is bounded.

### Notes

By default, this function returns `false` in order to be conservative.
Custom subtypes of `AbstractDirections` should hence add a method for this
function.
"""
function isbounding(ad::Type{<:AbstractDirections})
    return false
end

isbounding(::AD) where {AD<:AbstractDirections} = isbounding(AD)

"""
    isnormalized(ad::Type{<:AbstractDirections})

Returns whether the given template directions is normalized with respect to the
2-norm.

### Input

- `ad` -- template directions

### Output

`true` if the 2-norm of each element in `ad` is one and `false` otherwise.
"""
function isnormalized(ad::Type{<:AbstractDirections})
    return false
end

isnormalized(::AD) where {AD<:AbstractDirections} = isnormalized(AD)

# ==================================================
# Box directions
# ==================================================

"""
    BoxDirections{N, VN} <: AbstractDirections{N, VN}

Box directions representation.

### Fields

- `n` -- dimension

### Notes

Box directions can be seen as the vectors where only one entry is ±1, and all
other entries are 0. In dimension ``n``, there are ``2n`` such directions.

The default vector representation used in this template is a
`LazySets.Arrays.SingleEntryVector`, although other implementations can be used
such as a regular `Vector` and a sparse vector, `SparseVector`.

### Examples

The template can be constructed by passing the dimension. For example, in dimension
two,

```jldoctest dirs_Box
julia> dirs = BoxDirections(2)
BoxDirections{Float64,LazySets.Arrays.SingleEntryVector{Float64}}(2)

julia> length(dirs)
4
```

By default, each direction is represented in this iterator as a `SingleEntryVector`,
i.e. a vector with only one non-zero element,

```jldoctest dirs_Box
julia> eltype(dirs)
LazySets.Arrays.SingleEntryVector{Float64}
```
In two dimensions, the directions defined by `BoxDirections` are normal to the
facets of a box.

```jldoctest dirs_Box
julia> collect(dirs)
4-element Array{LazySets.Arrays.SingleEntryVector{Float64},1}:
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]
```

The numeric type can be specified as well:

```jldoctest
julia> BoxDirections{Rational{Int}}(10)
BoxDirections{Rational{Int64},LazySets.Arrays.SingleEntryVector{Rational{Int64}}}(10)

julia> length(ans)
20
```
"""
struct BoxDirections{N, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}
    n::Int
end

# convenience constructor for type Float64
BoxDirections(n::Int) = BoxDirections{Float64, SingleEntryVector{Float64}}(n)

# constructor where only N is specified
BoxDirections{N}(n::Int) where {N} = BoxDirections{N, SingleEntryVector{N}}(n)

Base.eltype(::Type{BoxDirections{N, VN}}) where {N, VN} = VN
Base.length(bd::BoxDirections) = 2 * bd.n

# interface functions
dim(bd::BoxDirections) = bd.n
isbounding(::Type{<:BoxDirections}) = true
isnormalized(::Type{<:BoxDirections}) = true

# The idea is that positive states run through vectors with +1 entry,
# and negative states run through vectors with -1 entry
# (1, 0)   state = 1
# (0, 1)   state = 2
# (0, -1)  state = -2
# (-1, 0)  state = -1
function Base.iterate(bd::BoxDirections{N, SingleEntryVector{N}}, state::Int=1) where {N}
    if state == 0
        return nothing
    end
    vec = SingleEntryVector(abs(state), bd.n, convert(N, sign(state)))
    state = (state == bd.n) ? -bd.n : state + 1
    return (vec, state)
end

function Base.iterate(bd::BoxDirections{N, Vector{N}}, state::Int=1) where {N}
    if state == 0
        return nothing
    end
    vec = zeros(N, bd.n)
    vec[abs(state)] = convert(N, sign(state))
    state = (state == bd.n) ? -bd.n : state + 1
    return (vec, state)
end

function Base.iterate(bd::BoxDirections{N, SparseVector{N, Int}}, state::Int=1) where {N}
    if state == 0
        return nothing
    end
    vec = sparsevec([abs(state)], convert(N, sign(state)), bd.n)
    state = (state == bd.n) ? -bd.n : state + 1
    return (vec, state)
end

# ==================================================
# Octagonal directions
# ==================================================

"""
    OctDirections{N, VN} <: AbstractDirections{N, VN}

Octagon directions representation.

### Fields

- `n` -- dimension

### Notes

Octagon directions consist of all vectors that are zero almost everywhere except
in two dimensions ``i``, ``j`` (possibly ``i = j``) where it is ``±1``. In dimension
``n``, there are ``2n^2`` such directions.

### Examples

The template can be constructed by passing the dimension. For example, in dimension
two,

```jldoctest dirs_Oct
julia> dirs = OctDirections(2)
OctDirections{Float64,SparseArrays.SparseVector{Float64,Int64}}(2)

julia> length(dirs) # number of directions
8
```
By default, each direction is represented in this iterator as a sparse vector:

```jldoctest dirs_Oct
julia> eltype(dirs)
SparseArrays.SparseVector{Float64,Int64}
```
In two dimensions, the directions defined by `OctDirections` are normal to
the facets of an octagon.

```jldoctest dirs_Oct
julia> first(dirs)
2-element SparseArrays.SparseVector{Float64,Int64} with 2 stored entries:
  [1]  =  1.0
  [2]  =  1.0

julia> Vector.(collect(dirs))
8-element Array{Array{Float64,1},1}:
 [1.0, 1.0]
 [1.0, -1.0]
 [-1.0, 1.0]
 [-1.0, -1.0]
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]
```

The numeric type can be specified as well:

```jldoctest
julia> OctDirections{Rational{Int}}(10)
OctDirections{Rational{Int64},SparseArrays.SparseVector{Rational{Int64},Int64}}(10)

julia> length(ans)
200
```
"""
struct OctDirections{N, VN} <: AbstractDirections{N, VN}
    n::Int
end

# constructor for type Float64
OctDirections(n::Int) = OctDirections{Float64, SparseVector{Float64, Int}}(n)

# constructor where only N is specified
OctDirections{N}(n::Int) where {N} = OctDirections{N, SparseVector{N, Int}}(n)

Base.eltype(::Type{OctDirections{N, VN}}) where {N, VN} = VN
Base.length(od::OctDirections) = 2 * od.n^2

# interface functions
dim(od::OctDirections) = od.n
isbounding(::Type{<:OctDirections}) = true
isnormalized(::Type{<:OctDirections}) = false

function Base.iterate(od::OctDirections{N, SparseVector{N, Int}}) where {N}
    if od.n == 1
        # fall back to box directions in 1D case
        return iterate(od, 1)
    end
    vec = spzeros(N, od.n)
    vec[1] = one(N)
    vec[2] = one(N)
    return (copy(vec), (vec, 1, 2))
end

# basic idea: modify the vector from the previous iteration
# have two indices i = 1, j = i + 1
# - run j through all indices > i
# - i = i + 1, j = i + 1
# - for any pair (i, j) create four vectors [..., i: ±1, ..., j: ±1, ...]
# in the end continue with box directions
function Base.iterate(od::OctDirections{N, SparseVector{N, Int}}, state::Tuple) where {N}
    # continue with octagon directions
    vec = state[1]
    i = state[2]
    j = state[3]
    if vec[j] > 0
        # change sign at j
        vec[j] = -one(N)
    elseif vec[i] > 0
        # change sign at i and j
        vec[i] = -one(N)
        vec[j] = one(N)
    elseif j < od.n
        # advance j, reset i
        vec[i] = one(N)
        vec[j] = zero(N)
        j = j + 1
        vec[j] = one(N)
    elseif i < od.n - 1
        # advance i
        vec[i] = zero(N)
        vec[j] = zero(N)
        i = i + 1
        j = i + 1
        vec[i] = one(N)
        vec[j] = one(N)
    else
        # continue with box directions
        return iterate(od, 1)
    end
    return (copy(vec), (vec, i, j))
end

function Base.iterate(od::OctDirections{N, SparseVector{N, Int}}, state::Int) where {N}
    # continue with box directions
    return iterate(BoxDirections{N, SparseVector{N, Int}}(od.n), state)
end

# ==================================================
# Box-diagonal directions
# ==================================================

"""
    BoxDiagDirections{N, VN} <: AbstractDirections{N, VN}

Box-diagonal directions representation.

### Fields

- `n` -- dimension

### Notes

Box-diagonal directions can be seen as the union of diagonal directions (all
entries are ±1) and box directions (one entry is ±1, all other entries are 0).
The iterator first enumerates all diagonal directions, and then all box
directions. In dimension ``n``, there are in total ``2^n + 2n`` such directions.

## Examples

The template can be constructed by passing the dimension. For example, in dimension
two,

```jldoctest dirs_BoxDiag
julia> dirs = BoxDiagDirections(2)
BoxDiagDirections{Float64,Array{Float64,1}}(2)

julia> length(dirs) # number of directions
8
```
By default, each direction is represented in this iterator as a regular vector:

```jldoctest dirs_BoxDiag
julia> eltype(dirs)
Array{Float64,1}
```
In two dimensions, the directions defined by `BoxDiagDirections` are normal to
the facets of an octagon.

```jldoctest dirs_BoxDiag
julia> collect(dirs)
8-element Array{Array{Float64,1},1}:
 [1.0, 1.0]
 [-1.0, 1.0]
 [1.0, -1.0]
 [-1.0, -1.0]
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]
```

The numeric type can be specified as well:

```jldoctest
julia> BoxDiagDirections{Rational{Int}}(10)
BoxDiagDirections{Rational{Int64},Array{Rational{Int64},1}}(10)

julia> length(ans)
1044
```
"""
struct BoxDiagDirections{N, VN} <: AbstractDirections{N, VN}
    n::Int
end

# constructor for type Float64
BoxDiagDirections(n::Int) = BoxDiagDirections{Float64, Vector{Float64}}(n)

# constructor where only N is specified
BoxDiagDirections{N}(n::Int) where {N} = BoxDiagDirections{N, Vector{N}}(n)

Base.eltype(::Type{BoxDiagDirections{N, VN}}) where {N, VN} = VN
Base.length(bdd::BoxDiagDirections) = bdd.n == 1 ? 2 : 2^bdd.n + 2 * bdd.n

# interface function
dim(bdd::BoxDiagDirections) = bdd.n
isbounding(::Type{<:BoxDiagDirections}) = true
isnormalized(::Type{<:BoxDiagDirections}) = false

function Base.iterate(bdd::BoxDiagDirections{N, Vector{N}}) where {N}
    return (ones(N, bdd.n), ones(N, bdd.n))
end

function Base.iterate(bdd::BoxDiagDirections{N}, state::Vector{N}) where {N}
    # continue with diagonal directions
    i = 1
    while i <= bdd.n && state[i] < 0
        state[i] = -state[i]
        i = i+1
    end
    if i > bdd.n
        if bdd.n == 1
            # finish here to avoid duplicates
            return nothing
        else
            # continue with box directions
            return iterate(bdd, 1)
        end
    else
        state[i] = -state[i]
        return (copy(state), state)
    end
end

function Base.iterate(bdd::BoxDiagDirections{N, Vector{N}}, state::Int) where {N}
    # continue with box directions
    return iterate(BoxDirections{N, Vector{N}}(bdd.n), state)
end

# ==================================================
# Polar directions
# ==================================================

"""
    PolarDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

Polar directions representation.

### Fields

- `Nφ` -- length of the partition of the polar angle

### Notes

The `PolarDirections` constructor provides a sample of the unit sphere
in ``\\mathbb{R}^2``, which is parameterized by the polar angle
``φ ∈ Dφ := [0, 2π]``; see the wikipedia entry
[Polar coordinate system](https://en.wikipedia.org/wiki/Polar_coordinate_system)
for details.

The integer argument ``Nφ`` defines how many samples of ``Dφ`` are taken. The
Cartesian components of each direction are obtained with

```math
[cos(φᵢ), sin(φᵢ)].
```

### Examples

The integer passed as an argument is used to discretize ``φ``:

```jldoctest; filter = r"2246[0-9]*e-16"
julia> pd = PolarDirections(2);

julia> pd.stack
2-element Array{Array{Float64,1},1}:
 [1.0, 0.0]
 [-1.0, 1.2246467991473532e-16]

julia> length(pd)
2
```
"""
struct PolarDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}
    Nφ::Int
    stack::Vector{VN} # stores the polar directions
end

# convenience constructors
PolarDirections(Nφ::Int) = PolarDirections{Float64, Vector{Float64}}(Nφ)
PolarDirections{N}(Nφ::Int) where {N} = PolarDirections{N, Vector{N}}(Nφ)

function PolarDirections{N, Vector{N}}(Nφ::Int) where {N}
    if Nφ <= 0
        throw(ArgumentError("Nφ = $Nφ is invalid; it shoud be at least 1"))
    end
    stack = Vector{Vector{N}}(undef, Nφ)
    φ = range(N(0), N(2*pi), length=Nφ+1)  # discretization of the polar angle

    @inbounds for i in 1:Nφ  # skip last (repeated) angle
        stack[i] = N[cos(φ[i]), sin(φ[i])]
    end
    return PolarDirections{N, Vector{N}}(Nφ, stack)
end

Base.eltype(::Type{PolarDirections{N, VN}}) where {N, VN} = VN
Base.length(pd::PolarDirections) = pd.Nφ

# interface functions
dim(pd::PolarDirections) = 2
isbounding(pd::PolarDirections) = pd.Nφ > 2
isnormalized(::Type{<:PolarDirections}) = true

function Base.iterate(pd::PolarDirections{N, Vector{N}}, state::Int=1) where {N}
    state == pd.Nφ + 1 && return nothing
    return (pd.stack[state], state + 1)
end

# ==================================================
# Spherical directions
# ==================================================

"""
    SphericalDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

Spherical directions representation.

### Fields

- `Nθ`    -- length of the partition of the azimuthal angle
- `Nφ`    -- length of the partition of the polar angle
- `stack` -- list of computed directions

### Notes

The `SphericalDirections` constructor provides a sample of the unit sphere
in ``\\mathbb{R}^3``, which is parameterized by the azimuthal and polar angles
``θ ∈ Dθ := [0, π]`` and ``φ ∈ Dφ := [0, 2π]`` respectively, see the wikipedia
entry [Spherical coordinate system](https://en.wikipedia.org/wiki/Spherical_coordinate_system)
for details.

The integer arguments ``Nθ`` and ``Nφ`` define how many samples along the domains
``Dθ`` and ``Dφ`` respectively are taken. The Cartesian components of each direction
are obtained with

```math
[sin(θᵢ)*cos(φᵢ), sin(θᵢ)*sin(φᵢ), cos(θᵢ)].
```

The north and south poles are treated separately so that those points
are not considered more than once.

### Examples

A `SphericalDirections` template can be built in different ways. If you pass
only one integer, the same value is used to discretize both ``θ`` and ``φ``:

```jldoctest spherical_directions; filter = r"1232[0-9]*e-17.*2246[0-9]*e-16.*1232[0-9]*e-17"
julia> sd = SphericalDirections(3);

julia> sd.Nθ, sd.Nφ
(3, 3)

julia> length(sd)
4
```

Pass two integers to control the discretization in ``θ`` and in ``φ`` separately:

```jldoctest spherical_directions
julia> sd = SphericalDirections(4, 5);

julia> length(sd)
10

julia> sd = SphericalDirections(4, 8);

julia> length(sd)
16
```
"""
struct SphericalDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}
    Nθ::Int
    Nφ::Int
    stack::Vector{VN} # stores the spherical directions
end

# convenience constructors
SphericalDirections(Nθ::Int) = SphericalDirections(Nθ, Nθ)
SphericalDirections(Nθ::Int, Nφ::Int) = SphericalDirections{Float64, Vector{Float64}}(Nθ::Int, Nφ::Int)
SphericalDirections{N}(Nθ::Int, Nφ::Int) where {N} = SphericalDirections{N, Vector{N}}(Nθ::Int, Nφ::Int)

function SphericalDirections{N, Vector{N}}(Nθ::Int, Nφ::Int) where {N}
    if Nθ <= 1 || Nφ <= 1
        throw(ArgumentError("(Nθ, Nφ) = ($Nθ, $Nφ) is invalid; both shoud be at least 2"))
    end
    stack = Vector{Vector{N}}()
    θ = range(N(0), N(pi), length=Nθ)    # discretization of the azimuthal angle
    φ = range(N(0), N(2*pi), length=Nφ)  # discretization of the polar angle

    # add north pole (θ = 0)
    push!(stack, N[0, 0, 1])

    # add south pole (θ = pi)
    push!(stack, N[0, 0, -1])

    for φᵢ in φ[1:Nφ-1]  # delete repeated angle
        for θⱼ in θ[2:Nθ-1] # delete north and south poles
            d = N[sin(θⱼ)*cos(φᵢ), sin(θⱼ)*sin(φᵢ), cos(θⱼ)]
            push!(stack, d)
        end
    end
    return SphericalDirections{N, Vector{N}}(Nθ, Nφ, stack)
end

Base.eltype(::Type{SphericalDirections{N, VN}}) where {N, VN} = VN
Base.length(sd::SphericalDirections) = length(sd.stack)

# interface functions
dim(::SphericalDirections) = 3
isbounding(sd::SphericalDirections) = sd.Nθ > 2 && sd.Nφ > 2
isnormalized(::Type{<:SphericalDirections}) = true

function Base.iterate(sd::SphericalDirections, state::Int=1)
    state == length(sd.stack) + 1 && return nothing
    return (sd.stack[state], state + 1)
end

# ==================================================
# Custom template directions
# ==================================================

"""
    CustomDirections{N, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

User-defined template directions.

### Fields

- `directions` -- list of template directions
- `n`          -- dimension
- `isbounding` -- boundedness status

### Notes

This struct is a wrapper type for a set of user-defined directions which are
iterated over. It has fields for the list of directions, the set dimension,
and (boolean) cache fields for the boundedness and normalization properties.
The latter are checked by default upon construction.

To check boundedness, we overapproximate the unit ball in the infinity
norm using the given directions and check if the resulting set is bounded.

The dimension will also be determined automatically, unless the empty vector is
passed (in which case the optional argument `n` needs to be specified).

## Examples

Creating a template with box directions in dimension two:

```jldoctest
julia> dirs = CustomDirections([[1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, -1.0]]);

julia> dirs.directions
4-element Array{Array{Float64,1},1}:
 [1.0, 0.0]
 [-1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]

julia> LazySets.Approximations.isbounding(dirs)
true

julia> LazySets.Approximations.isnormalized(dirs)
true
```
"""
struct CustomDirections{N, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}
    directions::Vector{VN}
    n::Int
    bounded::Bool
    normalized::Bool
end

function CustomDirections(directions::Vector{VN},
                          n::Int=_determine_dimension(directions);
                          check_boundedness::Bool=true,
                          check_normalization::Bool=true) where {N, VN<:AbstractVector{N}}
    bounded = check_boundedness ? _isbounding(directions) : false
    normalized = check_normalization ? _isnormalized(directions) : false
    return CustomDirections(directions, n, bounded, normalized)
end

function _determine_dimension(directions)
    isempty(directions) && throw(ArgumentError("empty template directions " *
                                               "need a specified dimension"))
    return length(directions[1])
end

function _isbounding(directions::Vector{VN}) where {N, VN<:AbstractVector{N}}
    isempty(directions) && return false

    # check boundedness of the polyhedron `⋂_a ax <= 1` where `a` is a direction
    P = HPolyhedron([HalfSpace(dir, one(N)) for dir in directions])
    return isbounded(P)
end

function _isnormalized(directions::Vector{VN}) where {N, VN<:AbstractVector{N}}
    return all(x -> _isapprox(x, one(N)), norm.(directions, 2))
end

Base.eltype(::Type{CustomDirections{N, VN}}) where {N, VN} = VN
Base.length(cd::CustomDirections) = length(cd.directions)

# interface functions
dim(cd::CustomDirections) = cd.n
isbounding(cd::Type{<:CustomDirections}) = false # it is not a property of the type
isbounding(cd::CustomDirections) = cd.bounded
isnormalized(cd::Type{<:CustomDirections}) = false # it is not a property of the type
isnormalized(cd::CustomDirections) = cd.normalized

function Base.iterate(cd::CustomDirections{N}, state::Int=1) where {N}
    if state > length(cd.directions)
        return nothing
    end
    vec = cd.directions[state]
    state = state + 1
    return (vec, state)
end


text:
@testset "459.repeated-substring-pattern.jl" begin
    @test repeated_substring_pattern("abab") == true
    @test repeated_substring_pattern("aba") == false
    @test repeated_substring_pattern("abcabcabcabc") == true
end


text:
import Base: sqrt, exp, log, sin, tanh, -, abs, abs2, sign

cuda1 = [
"sqrt",
# "rsqrt",
# "cbrt",
# "rcbrt",
"exp",
# "exp2",
# "exp10",
# "expm1",
"log",
# "log2",
# "log10",
# "log1p",
"sin",
# "cos",
# "tan",
# "sinpi",
# "cospi",
# "asin",
# "acos",
# "atan",
# "sinh",
# "cosh",
"tanh",
# "asinh",
# "acosh",
# "atanh",
# "erf",
# "erfc",
# "erfinv",
# "erfcinv",
# "erfcx",
# "normcdf",
# "normcdfinv",
# "lgamma",
# "tgamma",
# "logb",
# "ilogb",
# "j0",
# "j1",
# "y0",
# "y1",
# "cyl_bessel_i0",
# "cyl_bessel_i1",
# "trunc",
# "round",
# "rint",
# "nearbyint",
# "ceil",
# "floor",
# "lrint",
# "lround",
# "llrint",
# "llround",
("neg", "-", "-xi"),
("invx", "invx", "1/xi"),
("relu", "relu", "(xi>0?xi:0)"),
("sigm", "sigm", "1/(1+exp(-xi))"),
("abs", "abs", "(xi<0?-xi:xi)"),
("abs2", "abs2", "(xi*xi)"),
("sign", "sign", "(xi>0?1:xi<0?-1:0)"),
]

function cuda1def(f, j=f, o...)
    J=Symbol(j)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$S"
        @eval begin
            function $J(x::DexArray{$T})
                y = similar(x)
                ccall(($F,$libcudex),Void,(Cint,Ptr{$T},Ptr{$T}),length(y),x,y)
                return y
            end
        end
    end
end

#if isdefined(:libcudex)
    for f in cuda1
        isa(f,Tuple) || (f=(f,))
        cuda1def(f...)
    end
#end

# Define some common operations as primitives for efficiency:
# 1. Avoid creating intermediate arrays
# 2. Avoid taking derivatives of intermediate operations

for (f,g,y,dx) in ((:invx, :invxback, :(one(T)/x[i]), :(-y[i]*y[i]*dy[i])),
                   (:relu, :reluback, :(max(zero(T),x[i])), :(ifelse(y[i]>0,dy[i],zero(T)))),
                   (:sigm, :sigmback, :(one(T)/(one(T)+exp(-x[i]))), :(dy[i]*y[i]*(one(T)-y[i]))),
                   (:tanx, :tanhback, :(tanh(x[i])), :(dy[i]*(one(T)-y[i]*y[i]))),
                   )
    @eval begin
        function $f{T<:AbstractFloat}(x::Array{T})
            y = similar(x)
            @inbounds for i=1:length(y)
                y[i] = $y
            end
            return y
        end
        function $g{T<:AbstractFloat}(dy::Array{T},y::Array{T})
            dx = similar(dy)
            @inbounds for i=1:length(dx)
                dx[i] = $dx
            end
            return dx
        end
        # @primitive $f(x),dy,y $g(dy,y)
    end
end

# To avoid conflict with AutoGrad:
# @primitive tanh(x::Array),dy,y     tanhback(dy,y)
# @primitive tanh(x::DexArray),dy,y tanhback(dy,y)
# @primitive tanhback(dy,y),ddx  ddx.*(1.-y.*y)  ddx.*(-2.*dy.*y)

# Math for the cross-entropy loss: x is unnormalized input, p is
# target probabilities, q is estimated probabilities. Read left column
# down, right column (loss gradients) back up.

# x			dx = -p + qz/z = -p + exp(logq)
# xmax  = max(x,1)	-sum(db)=0
# logqz = x .- xmax	-p + qz/z
# qz    = exp(logqz)	rep(1/z)
# z     = sum(qz,1)	1/z
# logz  = log(z)	sum(p)=1
# logq  = logqz.-logz	-p
# plogq = p .* logq	-1
# loss  = -sum(plogq)	1

"""

logp(x,[dims]) treats entries in x as as unnormalized logp and returns
normalized logp.  If dims is not specified the normalization is over
the whole x, otherwise the normalization is performed over the given
dimensions.  In particular dims=1 normalizes columns of x and dims=2
normalizes rows of x.

"""
function logp(x,d...)
    x = x .- maximum(x,d...)
    x = x .- log(sum(exp(x),d...))
end

# dy should be -p and y=logq so this should give us -p+q
# @primitive  logp(x,d...),dy,y  (dy - exp(y).*sum(dy,d...))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct StatusMessage <: SwaggerModel
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    error::Any # spec type: Union{ Nothing, ErrorResponse } # spec name: error

    function StatusMessage(;status=nothing, error=nothing)
        o = new()
        validate_property(StatusMessage, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(StatusMessage, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type StatusMessage

const _property_map_StatusMessage = Dict{Symbol,Symbol}(Symbol("status")=>Symbol("status"), Symbol("error")=>Symbol("error"))
const _property_types_StatusMessage = Dict{Symbol,String}(Symbol("status")=>"String", Symbol("error")=>"ErrorResponse")
Base.propertynames(::Type{ StatusMessage }) = collect(keys(_property_map_StatusMessage))
Swagger.property_type(::Type{ StatusMessage }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_StatusMessage[name]))}
Swagger.field_name(::Type{ StatusMessage }, property_name::Symbol) =  _property_map_StatusMessage[property_name]

function check_required(o::StatusMessage)
    true
end

function validate_property(::Type{ StatusMessage }, name::Symbol, val)
end


text:
module TestFill

using Test
using RecordArrays

rawdata = """
fill((a = 1, b = 2))
fill((a = 1, b = 2), 10)
fill((a = 1, b = 2), 2, 3)
fill((a = 1, b = 2), 2, 3, 4)
fill(Some(0))
fill(Some(0), 10)
fill(Some(0), 2, 3)
fill(Some{Union{Int,Missing}}(0))
fill(Some{Union{Int,Missing}}(0), 10)
fill(Some{Union{Int,Missing}}(0), 2, 3)
""" |> x -> split(x, "\n", keepempty = false)

@testset "$code" for code in rawdata
    base = Base.include_string(@__MODULE__, code)
    rec = Base.include_string(@__MODULE__, "RecordArrays.$code")
    @test collect(rec) == base
end

end  # module


text:
abstract type ErrorMeasure end

(f::ErrorMeasure)(sstate::SimulationState) = f(sstate.hstate, sstate.fstate)

function (f::ErrorMeasure)(
    imodel, smodel, omodel, filter;
    num_timesteps, timestep, burnin::Int = 0
)
    sim = Simulation(imodel, smodel, omodel, filter)
    f(sim, num_timesteps, timestep, f, burnin)
end

function (f::ErrorMeasure)(sim::Simulation; num_timesteps, timestep, burnin::Int = 0)
    run(sim, num_timesteps, timestep, f, burnin)
end

struct MSE <: ErrorMeasure end

ComputeMSE = MSE()

function (f::MSE)(hstate::State, fstate::FilterState)
    w = hstate.w
    μ = fstate.μ
    msd(w, μ)
end


text:
"""
Placeholder for a short summary about ShortTermFundingMonitor.
"""
module ShortTermFundingMonitor

using Downloads: request
using JSON3
using DataFrames
using CodecZlib
using Dates

export stfmapicall, stfmmetadataquery, stfmavailablemnemonics, stfmsearchquery
export stfmdatasetdata, stfmseriesdata, stfmsingleseriesdata, stfmseriesspread

const STFM_BASEURL = "https://data.financialresearch.gov/"

const STFM_API_STABLE = "v1"

function stfmapicall(url; timeout=Inf, verbose=false)
    io = IOBuffer()
    res=request(url;output=io, method="GET", timeout=timeout, verbose=verbose)
    gzipcompressed = Pair("content-encoding","gzip") in res.headers
    if gzipcompressed
        js = JSON3.read(String(transcode(GzipDecompressor,take!(io))))
    else
        js = JSON3.read(String(take!(io)))
    end
    close(io)
    return js
end

# function downloadprogress(total,now)
#     iszero(total) && return nothing
#     f = round(Int,100 * now/total)
#     println("$f percent done downloading.")
#     return nothing
# end

function stfmavailablemnemonics(dataset=nothing; parse=true, timeout=Inf, verbose=false)
    point = isnothing(dataset) ? "/metadata/mnemonics?output=by_dataset" : "/metadata/mnemonics?dataset=$dataset"
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,point);timeout,verbose)
    parse == false && return js
    df = isnothing(dataset) ? mapreduce(_pair2table,vcat,pairs(js)) : DataFrame(values(js))
    return df
end

function stfmmetadataquery(mnemonic,fields=""; timeout=Inf, verbose=false)
    querybase="/metadata/query?mnemonic=$mnemonic"
    fieldstring = fields isa String ? string("&fields=",fields) : string("&fields=",join(fields,","))
    query = isempty(fields) ? querybase : string(querybase,fieldstring)
    return stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
end

function stfmsearchquery(s; parse=true, timeout=Inf, verbose=false)
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,"/metadata/search?query=",s);timeout,verbose)
    parse == false && return js
    if isempty(js)
        println("No results returned.")
        return nothing
    else
        df = DataFrame(NamedTuple.(values(js)))
        df[!,"resultnum"] = Int.(keys(js))
        select!(df,"resultnum",:)
        return df
    end
    return js
end

function stfmdatasetdata(dataset::AbstractString; parse=true,timeout=Inf,verbose=false,
    start_date=nothing,end_date=nothing,periodicity=nothing,how=nothing,remove_nulls=nothing,time_format=nothing)
    query = "/series/dataset"
    if isempty(dataset)
        js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
        parse == false && return js
        df = DataFrame(NamedTuple.(values(js)))
        df[!,"dataset"] = string.(keys(js))
        select!(df,"dataset",:)
        return df
    end
    query = string(query,"?dataset=$dataset") * _constructparams(;start_date,end_date,periodicity,how,remove_nulls,time_format)
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
    parse == false && return js

    metadatadict = _constructmetadatadict(js.timeseries)
    df = _constructmultiseriesdata(js.timeseries)

    return (data=df,metadata=metadatadict)
end

function stfmseriesdata(mnemonics;endpoint="multifull",parse=true,timeout=Inf,verbose=false,
    start_date=nothing,end_date=nothing,periodicity=nothing,how=nothing,remove_nulls=nothing,time_format=nothing)
    mstr = mnemonics isa String ? mnemonics : join(mnemonics,",")
    @assert !isempty(mstr) "Requires at least 1 mnemonic as input."
    mflag = isequal(endpoint,"timeseries") ? "mnemonic" : "mnemonics"
    query = string("/series/",endpoint,"?",mflag,"=",mstr) *  _constructparams(;start_date,end_date,periodicity,how,remove_nulls,time_format)
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
    parse == false && return js

    metadatadict = _constructmetadatadict(js)
    df = _constructmultiseriesdata(js)

    return (data=df,metadata=metadatadict)
end

function stfmsingleseriesdata(mnemonic,label=nothing;add_ids=false,parse=true,timeout=Inf,verbose=false,
    start_date=nothing,end_date=nothing,periodicity=nothing,how=nothing,remove_nulls=nothing,time_format=nothing)
    query = string("/series/timeseries?mnemonic=",mnemonic) *  _constructparams(;start_date,end_date,periodicity,how,remove_nulls,time_format)
    if !isnothing(label)
        query *= string("&label=",label)
    end
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
    parse == false && return js

    df = add_ids ? _parsetimeseries(js;mnemonic,label) : _parsetimeseries(js)

    return df
end

function stfmseriesspread(xmnemonic,ymnenomic;add_ids=false,parse=true,timeout=Inf,verbose=false,
    start_date=nothing,end_date=nothing,periodicity=nothing,how=nothing,remove_nulls=nothing,time_format=nothing)
    query = string("/calc/spread?x=",xmnemonic,"&y=",ymnenomic) *  _constructparams(;start_date,end_date,periodicity,how,remove_nulls,time_format)
    js = stfmapicall(string(STFM_BASEURL,STFM_API_STABLE,query);timeout,verbose)
    parse == false && return js

    df = add_ids ? _parsetimeseries(js;mnemonic=xmnemonic,label=ymnenomic) : _parsetimeseries(js)
    add_ids && rename!(df,"mnemonic"=>"xmnemonic","label"=>"ymnenomic")

    return df
end

function _pair2table(p)
    key,values = p
    df = DataFrame(NamedTuple.(values))
    df[!,"dataset"] .= string(key)
    select!(df,"dataset",:)
    return df
end

function _parsetimeseries(x; mnemonic=nothing, label=nothing)
    date = first(first(x)) isa Int ? first.(x) : Date.(first.(x))
    v = replace(last.(x),nothing=>missing)
    v = v * 1.0             # converts to Float64
    if all(!ismissing,v)
        v = disallowmissing(v)
    end
    df = DataFrame((date=date,value=v))
    namevec = names(df)
    if !isnothing(label)
        df[!,"label"] .= label
        pushfirst!(namevec,"label")
    end
    if !isnothing(mnemonic)
        df[!,"mnemonic"] .= mnemonic
        pushfirst!(namevec,"mnemonic")
    end
    select!(df,namevec...)
    return df
end

function _constructparams(;start_date=nothing,end_date=nothing,periodicity=nothing,how=nothing,remove_nulls=nothing,time_format=nothing)
    str = ""
    !isnothing(start_date)      && (str *= string("&start_date=",start_date))
    !isnothing(end_date)        && (str *= string("&end_date=",end_date))
    !isnothing(periodicity)     && (str *= string("&periodicity=",periodicity))
    !isnothing(how)             && (str *= string("&how=",how))
    !isnothing(remove_nulls)    && (str *= string("&remove_nulls=",remove_nulls))
    !isnothing(time_format)     && (str *= string("&time_format=",time_format))
    return str
end

function _constructmultiseriesdata(jst)
    mnemonics = keys(jst)
    df = DataFrame()
    for m in mnemonics
        for k in keys(jst[m])
            if k != :metadata
                for j in keys(jst[m][k])
                    dftmp = _parsetimeseries(jst[m][k][j];mnemonic=string(m),label=string(j))
                    df = vcat(df,dftmp)
                end
            end
        end
    end
    return df
end

function _constructmetadatadict(jst)
    mnemonics = keys(jst)
    metadatadict = Dict{String,Dict{Symbol,Any}}()
    for m in mnemonics
        for k in keys(jst[m])
            if k == :metadata
                push!(metadatadict,string(m) => Dict(jst[m][k]))
            end
        end
    end
    return metadatadict
end

end # module


text:
include("private/core.jl")
include("private/private.jl")

mutable struct Window

    m_width::Int64
    m_height::Int64
    m_title::String
    m_window::_private.J3DWindow

    Window(width::Int64, height::Int64, title::String) = new(width, height, title, _private.J3DWindow())

end

function Init(window::Window)
    ccall(j3dloader.window_create, Cvoid, (Ref{_private.J3DWindow}, Csize_t, Csize_t, Cstring), window.m_window, window.m_width, window.m_height, window.m_title)
    ccall(j3dloader.window_initialize, Cchar, (Ref{_private.J3DWindow}, ), window.m_window)
end

function IsClose(window::Window)::Bool
    ccall(j3dloader.window_is_close, Cchar, (Ref{_private.J3DWindow}, ), window.m_window)
end

function SetTitle(window::Window, title::String)
    window.m_title = title
    ccall(j3dloader.window_set_title, Cvoid, (Ref{_private.J3DWindow}, Cstring), window.m_window, window.m_title)
end

function PollEvent()
    ccall(j3dloader.window_poll_event, Cvoid, ())
end

function Clear(window::Window)
    ccall(j3dloader.window_clear, Cvoid, (Ref{_private.J3DWindow}, ), window.m_window)
end

function Display(window::Window)
    ccall(j3dloader.window_display, Cvoid, (Ref{_private.J3DWindow}, ), window.m_window)
end

function Destroy(window::Window)
    ccall(j3dloader.window_destroy, Cvoid, (Ref{_private.J3DWindow}, ), window.m_window)
end


text:
abstract type AbstractODEProblem end

struct ODEProblem{F,T,U,P} <: AbstractODEProblem
    f::F
    tspan::T
    u0::U
    θ::P
end

abstract type ODESolver end
struct Euler{T} <: ODESolver
    dt::T
end
struct RK2{T} <: ODESolver
    dt::T
end

function f(x,θ)
    α, β, γ, δ = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

function solve(prob::AbstractODEProblem, solver::ODESolver)
    t = prob.tspan[1]; u = prob.u0
    us = [u]; ts = [t]
    while t < prob.tspan[2]
        (u,t) = solver(prob, u, t)
        push!(us,u)
        push!(ts,t)
    end
    ts, reduce(hcat,us)
end

function (solver::Euler)(prob::ODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    (u + dt*f(u,θ), t+dt)
end

function (solver::RK2)(prob::ODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    uh = u + f(u,θ)*dt
    u + dt/2*(f(u,θ) + f(uh,θ)), t+dt
end

θ = [0.1,0.2,0.3,0.2]
u0 = [1.0,1.0]
tspan = (0.,100.)
dt = 0.1
prob = ODEProblem(f,tspan,u0,θ)

t,X=solve(prob, RK2(0.2))

# using Plots
# p1 = plot(t, X[1,:], label="x", lw=3)
# plot!(p1, t, X[2,:], label="y", lw=3)

# display(p1)

#
θ = [0.2,0.2,0.3,0.2]
u0 = [1.0,1.0]
tspan = (0.,100.)
dt = 0.1
prob2 = ODEProblem(f,tspan,u0,θ)

t,X2=solve(prob2, RK2(0.2))

using Optim

function loss(θin,prob::ODEProblem,Y)
    prob.θ.=θin
    t,Xn=solve(prob,RK2(0.2))
    sum((Y.-Xn).^2)
end
θopt = copy(θ)
O=Optim.optimize(θ->loss(θ,prob,X),θopt)
O=Optim.optimize(θ->loss(θ,prob,X),θopt,LBFGS())

using DiffEqFlux
nn=FastDense(2,2)
p = initial_params(nn)
nn([1,2],p)

function fy(x,θ)
    α, β, γ, δ, ω = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂ + ω*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

#
θy = [0.2,0.2,0.3,0.2,0.1]
u0 = [1.0,1.0]
tspan = (0.,100.)
dt = 0.1
proby = ODEProblem(fy,tspan,u0,θy)

t,Xy=solve(proby, RK2(0.2))

# py = plot(t, Xy[1,:], label="x", lw=3)
# plot!(py, t, Xy[2,:], label="y", lw=3)
# savefig("LV_omega.svg")

function fnn(x,θ)
    α, β, γ, δ = θ[1:4]
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]+nn(x,@view θ[5:end])
end

θnn = [0.2,0.2,0.3,0.2,0.01*initial_params(nn)...]
probnn = ODEProblem(fnn,tspan,u0,θnn)

θopt = copy(θnn)
O=Optim.optimize(θ->loss(θ,probnn,Xy),θopt,Optim.Options(iterations=10000))


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

using Markdown
using InteractiveUtils

# ╔═╡ 2d0d4838-eb58-4f2a-9864-40f4336f43f1
#hide
import Pkg

# ╔═╡ e7a181c5-e00e-44cb-8c42-c27983deae25
#hide
Pkg.activate(Base.current_project())

# ╔═╡ 1d0f6811-67cc-4ef1-b2dd-2f2a1e809061
using PlutoTest

# ╔═╡ 64f0c27d-9276-4e9f-b09c-414dd93ffe83
#hide
include("../src/Export.jl")

# ╔═╡ 6f1cd19f-b529-41d7-a9d6-aca3cb2888be
#hide
include("../src/Documenter.jl")

# ╔═╡ 3d2e0115-5643-488a-8932-006912864c5b
md"This module takes care of the creation of the support system document generation such as directory creation, index.yml and mkdocs.yml creation.

Additionally it calls mkdocs commands to build the site"

# ╔═╡ 0c7570c3-8871-494f-91a0-fc04376344c3
md"## makepath"

# ╔═╡ ab5e0c09-223c-46c3-992c-928d42b0cb3b
#export
"""create the given path if not present"""
function makepath(path:: AbstractString)
	if !isdir(path)
		mkdir(path)
	else
		@warn "Path not created becasue it was already present."
	end
end

# ╔═╡ 0fc70e35-557d-4bc1-a0bf-b03b68987d81
Documenter.showdoc(makepath)

# ╔═╡ 23f638ee-4f7d-4480-8c3e-dbd29f5dd39b
#hide
#makepath("../test/docs")

# ╔═╡ 6c17a385-11ba-4989-bedb-d9cb1b94b123
begin
	path = joinpath("..", "test", "dir")
	makepath(joinpath("..", "test", path))
end

# ╔═╡ f7950f63-f57b-427a-a3d7-8d280a21b22a
#noop
@test isdir(path) == true

# ╔═╡ bc3d912d-1d56-4e86-83d7-5578feb21e92
#cleanup
rm(path)

# ╔═╡ fe6028bc-58b8-45d1-bf7b-8f5e31e8d452
md"## mkdocsyml"

# ╔═╡ ad99a75a-2918-4e27-b679-3ab3611cc3fc
#export
"""
> mkdocsyml()--> Creates a mkdocs.yml file
> at the document root and fills it with initial
> content.
"""
function mkdocsyml(path:: AbstractString)
	mkdocyml = joinpath(path, "mkdocs.yml")
	if !isfile(mkdocyml)
	    touch(mkdocyml)
	else
		@warn "File not created because it was already present"
    end

end

# ╔═╡ 18af50bd-bc5f-4e38-ab97-2a8cef56805a
Documenter.showdoc(mkdocsyml)

# ╔═╡ 4b3477a9-cf4a-4ea6-8eac-359cb34d97be
#mkdocs_yml("../docs")

# ╔═╡ df25274d-b3b5-41d4-ae31-8cecc1d8fdea
begin
makepath(path)
mkdocsyml(path)
end

# ╔═╡ 79c67ee4-e136-47ab-ba61-a0fefd402640
#noop
@test isfile(joinpath(path, "mkdocs.yml")) == true

# ╔═╡ 0688f354-b3a4-4c45-a709-17a73861ee44
#cleanup
rm(joinpath(path,"mkdocs.yml"))

# ╔═╡ e13bfb2b-7630-495c-a049-9c7283cbd191
#cleanup
rm(path)

# ╔═╡ eeb046c0-bdcc-11eb-31a1-838f33f977e3
md"## newsitegen"

# ╔═╡ cc200174-9530-438f-8004-3883c372f674
#hide
readdir("../")

# ╔═╡ 4e27289a-5e5c-4a04-820d-e3fb7cc78afb
#export
begin
"""
> newsitegen()--> Create required directory structure for hosting documents and the initial files(index.md and mkdocs.yml) if not already there.
"""
function newsitegen()
	path = joinpath(".","docs")
    makepath(path)
end
end

# ╔═╡ 7f1e88c8-6802-4a1c-be12-115eba510fac
#hide
#newsitegen()

# ╔═╡ c5160350-960f-4f42-a6bd-a0c2d9b03b7f
Documenter.showdoc(newsitegen)

# ╔═╡ 1ed0902f-34ce-4b62-b484-fdf70124394a
#hide
Export.notebook2script(joinpath("..", "nbs"), joinpath("..", "src"))

# ╔═╡ Cell order:
# ╠═2d0d4838-eb58-4f2a-9864-40f4336f43f1
# ╠═e7a181c5-e00e-44cb-8c42-c27983deae25
# ╠═64f0c27d-9276-4e9f-b09c-414dd93ffe83
# ╠═6f1cd19f-b529-41d7-a9d6-aca3cb2888be
# ╠═1d0f6811-67cc-4ef1-b2dd-2f2a1e809061
# ╠═3d2e0115-5643-488a-8932-006912864c5b
# ╠═0c7570c3-8871-494f-91a0-fc04376344c3
# ╠═ab5e0c09-223c-46c3-992c-928d42b0cb3b
# ╠═0fc70e35-557d-4bc1-a0bf-b03b68987d81
# ╠═23f638ee-4f7d-4480-8c3e-dbd29f5dd39b
# ╠═6c17a385-11ba-4989-bedb-d9cb1b94b123
# ╠═f7950f63-f57b-427a-a3d7-8d280a21b22a
# ╠═bc3d912d-1d56-4e86-83d7-5578feb21e92
# ╠═fe6028bc-58b8-45d1-bf7b-8f5e31e8d452
# ╠═ad99a75a-2918-4e27-b679-3ab3611cc3fc
# ╠═18af50bd-bc5f-4e38-ab97-2a8cef56805a
# ╠═4b3477a9-cf4a-4ea6-8eac-359cb34d97be
# ╠═df25274d-b3b5-41d4-ae31-8cecc1d8fdea
# ╠═79c67ee4-e136-47ab-ba61-a0fefd402640
# ╠═0688f354-b3a4-4c45-a709-17a73861ee44
# ╠═e13bfb2b-7630-495c-a049-9c7283cbd191
# ╠═eeb046c0-bdcc-11eb-31a1-838f33f977e3
# ╠═cc200174-9530-438f-8004-3883c372f674
# ╠═4e27289a-5e5c-4a04-820d-e3fb7cc78afb
# ╠═7f1e88c8-6802-4a1c-be12-115eba510fac
# ╠═c5160350-960f-4f42-a6bd-a0c2d9b03b7f
# ╠═1ed0902f-34ce-4b62-b484-fdf70124394a


text:
###############################################################################
#
# SoilPlantAirContinuum system for mono species system
#
###############################################################################
"""
    mutable struct SPACMono{FT}

Struct that mono species SoilPlantAirContinuum system.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct SPACMono{FT<:AbstractFloat}
    "Soil layers bounds `[m]`"
    soil_bounds::Array{FT,1} = FT[0,-0.1,-0.2,-0.3,-0.5,-0.8,-1.2,-2.0]
    "Air layers bounds `[m]`"
    air_bounds ::Array{FT,1} = collect(FT,0:1:20)

    "Root depth `[m]`"
    z_root  ::FT = FT(-1)
    "Canopy maximal height `[m]`"
    z_canopy::FT = FT(10)

    "Plant hydraulic system"
    plant_hs::AbstractPlantOrganism{FT} = create_grass(z_root, z_canopy,
                                                       soil_bounds, air_bounds)
    "Number of canopy layers"
    n_canopy::Int = length(plant_hs.canopy_index_in_air)
    "Number of root layers"
    n_root  ::Int = length(plant_hs.root_index_in_soil)
    "Plant photosynthesis systems"
    plant_ps::Array{CanopyLayer{FT},1} = [CanopyLayer{FT}() for
                                          i in 1:n_canopy]
    "Basal area `[m²]`"
    ba::FT = sum( [plant_hs.roots[i].area for i in 1:n_root] )
    "Ground area `[m²]`"
    ga::FT = ba * 500
    "Leaf area `[m²]`"
    la::FT = sum( plant_hs.leaves[i].area for i in 1:n_canopy )

    # Surrounding AirLayer
    "Air layers"
    envirs::Array{AirLayer{FT},1} = [AirLayer{FT}() for i in 1:n_canopy]

    # Wind related
    "Aerodynamic roughness `[m]`"
    wind_z0::FT = z_canopy * FT(0.07)
    "Zero plane displacement `[m]`"
    wind_d ::FT = z_canopy * 2/3
    "Mean layer height `[m]`"
    wind_zs::Array{FT,1} = [(air_bounds[i]+air_bounds[i+1])/2 for
                            i in 1:length(air_bounds)-1]
    "Wind speed per layer `[m s⁻¹]`"
    winds::Array{FT,1} = [FT(1) for i in 1:length(air_bounds)-1]

    # Soil layers information
    # TODO bridge Soil module later
    "Maximal soil water content"
    mswc  ::Array{FT,1} = [plant_hs.roots[i].sh.Θs for i in 1:n_root]
    "Current soil water content"
    swc   ::Array{FT,1} = [plant_hs.roots[i].sh.Θs for i in 1:n_root]
    "Array of soil matric potential `[MPa]`"
    p_soil::Array{FT,1} = FT[0 for i in 1:n_root]
    "Maximal soil depth `[m]`"
    h_soil::FT = abs(z_root)

    # geography related
    "Latitude `[°]`"
    latitude ::FT = 30
    "Longitude `[°]`"
    longitude::FT = 116
    "Elevation `[m]`"
    elevation::FT = 0

    # photosynthesis mode and stomatal model scheme
    "Photosynthesis parameter set"
    photo_set::AbstractPhotoModelParaSet{FT} = C3CLM(FT)
    "Stomatal behavior scheme"
    stomata_model::AbstractStomatalModel{FT} = ESMBallBerry{FT}()

    # For CanopyLayers module
    "Solar angle container"
    angles::SolarAngles{FT} = SolarAngles{FT}()
    "Canopy4RT container"
    canopy_rt::Canopy4RT{FT} = create_canopy_rt(FT, nLayer=n_canopy, LAI=la/ga)
    "Wave length container"
    wl_set::WaveLengths{FT} = create_wave_length(FT)
    "RT dimensions"
    rt_dim::CanopyLayers.RTDimensions = create_rt_dims(canopy_rt, wl_set);
    "CanopyRads container"
    can_rad::CanopyRads{FT} = create_canopy_rads(FT, rt_dim)
    "CanopyOpticals container"
    can_opt::CanopyOpticals{FT} = create_canopy_opticals(FT, rt_dim)
    "Array of LeafBios container"
    leaves_rt::Array{LeafBios{FT},1} = [create_leaf_bios(FT, rt_dim) for
                                         i in 1:n_canopy]
    "SoilOpticals container"
    soil_opt::SoilOpticals{FT} = create_soil_opticals(wl_set)
    "Incoming radiation container"
    in_rad::IncomingRadiation{FT} = create_incoming_radiation(wl_set)
    "RT container"
    rt_con::RTCache{FT} = create_rt_cache(FT, rt_dim)
    "Container for sunlit leaf area fraction in each layer"
    f_SL::Array{FT,1} = repeat(canopy_rt.lidf, outer=[ canopy_rt.nAzi ]) /
                        canopy_rt.nAzi;

    # local storage for canopy GPP and NPP
    "Canopy GPP per ground area"
    f_gpp::FT = 0
    "Canopy GPP per ground area"
    f_npp::FT = 0
    "Canopy water flux per ground area"
    f_H₂O::FT = 0
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
function resultConverter_v2(x, V, anchoLado, matConexionVertices, vecVertices, vecAlturas)

    template = Int(x[end])
    alt = min(x[1], maximum(vecAlturas))
    theta = x[2]

    # Genera Poligono que corta VolTeor a altura alt
    idAlt = findfirst(y -> y >= alt, vecAlturas) - 1
    vert_poly_0 = Int.(vecVertices[idAlt])
    vert_poly_1 = Int.(matConexionVertices[vert_poly_0,2])
    poly_0 = copy(V[vert_poly_0, 1:2])[:,1,:]
    poly_1 = copy(V[vert_poly_1, 1:2])[:,1,:]
    alfa_ = (vecAlturas[idAlt + 1] - alt) / (vecAlturas[idAlt + 1] - vecAlturas[idAlt])
    polyCorte_alt = alfa_ .* poly_0 .+ (1 - alfa_) .* poly_1

    if template == 1

        alfa = x[3]
        pos_x = x[4]
        pos_y = x[5]
        largo1 = x[6]
        largo2 = x[7]

        R1 = poly2D.rotationMatrix(theta);
        cr1 = [pos_x; pos_y]
        p1_1 = [pos_x; pos_y];
        p2_1 = R1 * ([pos_x + largo1; pos_y] - cr1) + cr1;
        p3_1 = R1 * ([pos_x + largo1; pos_y + anchoLado] - cr1) + cr1;
        p4_1 = R1 * ([pos_x ; pos_y + anchoLado] - cr1) + cr1;
        R2 = poly2D.rotationMatrix(alfa + theta);
        cr2 = [pos_x; pos_y]
        p1_2 = [pos_x; pos_y];
        p2_2 = R2 * ([pos_x + anchoLado; pos_y] - cr2) + cr2;
        p3_2 = R2 * ([pos_x + anchoLado; pos_y + largo2] - cr2) + cr2;
        p4_2 = R2 * ([pos_x; pos_y + largo2] - cr2) + cr2;
        V1 = [p1_1';p2_1';p3_1';p4_1']
        V2 = [p1_2';p2_2';p3_2';p4_2'];

        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        ps_base = polyShape.polyUnion_v2(ps1, ps2)
        areaBasal = polyShape.polyArea_v2(ps_base)
        ps_baseSeparada = PolyShape([V1, V2], 2)

    elseif template == 2
        phi1 = x[3]
        phi2 = x[4]
        pos_x0 = x[5]
        pos_y0 = x[6]
        largo0 = x[7]
        largo1 = x[8]
        largo2 = x[9]

        R_theta = poly2D.rotationMatrix(theta);
        cr_theta  = [pos_x0; pos_y0];

        p1_0  = (R_theta * ([pos_x0; pos_y0] - cr_theta) + cr_theta)'
        p2_0  = (R_theta * ([pos_x0 + largo0; pos_y0] - cr_theta) + cr_theta)'
        p3_0  = (R_theta * ([pos_x0 + largo0; pos_y0 + anchoLado] - cr_theta) + cr_theta)'
        p4_0  = (R_theta * ([pos_x0; pos_y0 + anchoLado] - cr_theta) + cr_theta)'

        R_phi1 = poly2D.rotationMatrix(phi1);
        cr_phi1  = [pos_x0; pos_y0];
        p1_1  = (R_theta * ((R_phi1 * ([pos_x0; pos_y0] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p2_1  = (R_theta * ((R_phi1 * ([pos_x0 + anchoLado; pos_y0] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p3_1  = (R_theta * ((R_phi1 * ([pos_x0 + anchoLado; pos_y0 + largo1] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p4_1  = (R_theta * ((R_phi1 * ([pos_x0; pos_y0 + largo1] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';

        R_phi2 = poly2D.rotationMatrix(-phi2);
        cr_phi2  = [pos_x0 + largo0; pos_y0];
        p1_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0 - anchoLado; pos_y0] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p2_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0; pos_y0] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p3_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0; pos_y0 + largo2] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p4_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0 - anchoLado; pos_y0 + largo2] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';

        V0 = [p1_0;p2_0;p3_0;p4_0]
        V1 = [p1_1;p2_1;p3_1;p4_1]
        V2 = [p1_2;p2_2;p3_2;p4_2]

        ps0 = PolyShape([V0], 1)
        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        ps_base = polyShape.polyUnion_v2(ps0, ps1)
        ps_base = polyShape.polyUnion_v2(ps_base, ps2)
        areaBasal = polyShape.polyArea_v2(ps_base)
        ps_baseSeparada = PolyShape([V0, V1, V2], 3)

    elseif template == 3

    elseif template == 4

        alfa = x[3]
        pos_x = x[4]
        pos_y = x[5]

        largo1 = x[6] #
        largo2 = x[7] #

        R1 = poly2D.rotationMatrix(theta);
        cr1 = [pos_x; pos_y]
        p1_1 = [pos_x; pos_y];
        p2_1 = R1 * ([pos_x + largo1; pos_y] - cr1) + cr1;
        p3_1 = R1 * ([pos_x + largo1; pos_y + anchoLado] - cr1) + cr1;
        p4_1 = R1 * ([pos_x ; pos_y + anchoLado] - cr1) + cr1;
        R2 = poly2D.rotationMatrix(theta - alfa);
        cr2 = [pos_x; pos_y]
        p1_2 = R2 * ([pos_x - anchoLado; pos_y] - cr2) + cr2;
        p2_2 = R2 * ([pos_x; pos_y] - cr2) + cr2;
        p3_2 = R2 * ([pos_x; pos_y + largo2] - cr2) + cr2;
        p4_2 = R2 * ([pos_x - anchoLado; pos_y + largo2] - cr2) + cr2;
        V1 = [p1_1';p2_1';p3_1';p4_1']
        V2 = [p1_2';p2_2';p3_2';p4_2'];

        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        ps_base = polyShape.polyUnion_v2(ps1, ps2)
        areaBasal = polyShape.polyArea_v2(ps_base)
        ps_baseSeparada = PolyShape([V1, V2], 2)

    end

    return alt, areaBasal, ps_base, ps_baseSeparada, psCorte

end


text:
module GrayUtils

using CUDA
using ..Device: array_type, array_device, CUDADevice, CPU
using ..AtmosphericStates
using ..Fluxes

using CLIMAParameters
using CLIMAParameters.Planet: grav, cp_d

export update_profile_lw!, compute_gray_heating_rate!

"""
    update_profile_lw!(
        t_lay,
        t_lev,
        flux_dn,
        flux_net,
        hr_lay,
        flux_grad,
        T_ex_lev,
        Δt,
        nlay,
        nlev,
        ncol,
        ::Type{FT},
    )

Updates t_lay and t_lev based on heating rate.
"""
function update_profile_lw!(
    t_lay,
    t_lev,
    flux_dn,
    flux_net,
    hr_lay,
    flux_grad,
    T_ex_lev,
    Δt,
    nlay,
    nlev,
    ncol,
    ::Type{FT},
) where {FT<:AbstractFloat}
    # updating t_lay and t_lev based on heating rate
    args = (
        t_lay,
        t_lev,
        hr_lay,
        flux_dn,
        flux_net,
        flux_grad,
        T_ex_lev,
        Δt,
        nlay,
        nlev,
    )
    device = array_device(t_lev)
    if device === CUDADevice()
        max_threads = 256
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        @cuda threads = (tx) blocks = (bx) update_profile_lw_CUDA!(
            ncol,
            args...,
        )
    else
        # using julia native multithreading
        Threads.@threads for gcol = 1:ncol
            update_profile_lw_kernel!(args..., gcol)
        end
    end
    #----------------------------------------
end

function update_profile_lw_CUDA!(ncol, args...)
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    if gcol ≤ ncol
        update_profile_lw_kernel!(args..., gcol)
    end
    return nothing
end

function update_profile_lw_kernel!(
    t_lay,
    t_lev,
    hr_lay,
    flux_dn,
    flux_net,
    flux_grad,
    T_ex_lev,
    Δt::FT,
    nlay,
    nlev,
    gcol,
) where {FT<:AbstractFloat}
    # updating t_lay based on heating rate
    for glay = 1:nlay
        @inbounds t_lay[glay, gcol] += Δt * hr_lay[glay, gcol]
    end
    #--------compute t_lev from t_lay--------------------------------------
    for glev = 2:nlay-1
        @inbounds t_lev[glev, gcol] =
            FT(1 / 3) * t_lay[glev-1, gcol] + FT(5 / 6) * t_lay[glev, gcol] -
            FT(1 / 6) * t_lay[glev+1, gcol]
    end

    @inbounds t_lev[nlay, gcol] =
        FT(1 / 3) * t_lay[nlay, gcol] + FT(5 / 6) * t_lay[nlay-1, gcol] -
        FT(1 / 6) * t_lay[nlay-2, gcol]

    @inbounds t_lev[1, gcol] = FT(2) * t_lay[1, gcol] - t_lev[2, gcol]

    @inbounds t_lev[nlay+1, gcol] =
        FT(2) * t_lay[nlay, gcol] - t_lev[nlay, gcol]
    #-----------------------------------------------------------------------
    sbc = FT(Stefan())
    for glev = 1:nlev
        @inbounds T_ex_lev[glev, gcol] =
            (
                (flux_dn[glev, gcol] + (flux_net[glev, gcol] / FT(2))) / sbc
            )^FT(0.25)
    end

    for glev = 2:nlev
        @inbounds flux_grad[glev-1, gcol] =
            abs(flux_net[glev, gcol] - flux_net[glev-1, gcol])
    end
    #-----------------------------------------------------------------------
end

"""
    compute_gray_heating_rate!(
        hr_lay,
        p_lev,
        ncol,
        nlay,
        flux_net,
        param_set,
        ::Type{FT},
    )

Computes the heating rate for the gray radiation simulation.
"""
function compute_gray_heating_rate!(
    hr_lay,
    p_lev,
    ncol,
    nlay,
    flux_net,
    param_set,
    ::Type{FT},
) where {FT}
    cp_d_ = FT(cp_d(param_set))
    grav_ = FT(grav(param_set))
    args = (hr_lay, flux_net, p_lev, grav_, cp_d_, nlay)
    device = array_device(p_lev)
    if device === CUDADevice()
        max_threads = 256
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        @cuda threads = (tx) blocks = (bx) compute_gray_heating_rate_CUDA!(
            ncol,
            args...,
        )
    else
        # Launcing native Julia multithreading Kernel
        Threads.@threads for gcol = 1:ncol
            compute_gray_heating_rate_kernel!(args..., gcol)
        end
    end
    return nothing
end

function compute_gray_heating_rate_CUDA!(ncol, args...)
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    if gcol ≤ ncol
        compute_gray_heating_rate_kernel!(args..., gcol)
    end
    return nothing
end

function compute_gray_heating_rate_kernel!(
    hr_lay,
    flux_net,
    p_lev,
    grav_,
    cp_d_,
    nlay,
    gcol,
)
    for ilay = 1:nlay
        @inbounds hr_lay[ilay, gcol] =
            grav_ * (flux_net[ilay+1, gcol] - flux_net[ilay, gcol]) /
            (p_lev[ilay+1, gcol] - p_lev[ilay, gcol]) / cp_d_
    end
    return nothing
end
end


text:
#--- AbstractMesh

"""
Abastract type AbastractMesh

AstractMesh is a interface defined by overloading the following methods:

- [`get_coord_vector(trian::Triangulation)`](@ref)
- [`get_cell_vector`](trian::Triangulation)`](@ref)
"""

abstract type AbstractMesh <: JalerkinType end

"""
get_coord_vector(Mesh::AbstractMesh) -> AbstractArray{Vector{<:Point}}}
"""
function get_coord_vector(mesh::AbstractMesh)
    @abstractmethod
end

"""
get_cell_vector(Mesh::AbstractMesh) -> AbstractArray{Vector{<:Cell}}}
"""
function get_cell_vector(mesh::AbstractMesh)
    @abstractmethod
end

#--- Mesh
"""Type Point (given dimension, type of real numbers).

Static Vector of Dim real numbers
"""
Point = SVector{Dim, fp_precision()} where {Dim}

"""Type Cell (given number of vértices)

Static vector of indices (identifiers of the cell points)
"""
Cell = SVector{NumVertex, Int64} where {NumVertex}

"""Mesh will contain sets of Cells and Vertices (Points)"""
struct Mesh{CellType, VertexType}  <: AbstractMesh
    "Vector of cells (each cell contains a list of vertex indices)"
    cell_vector :: Vector{CellType}

    "Vector of coordinates of each vertex"
    coord_vector :: Vector{VertexType}
end

"""
This function is part of the interface AbstractMesh
"""
@inline function get_coord_vector(mesh::Mesh{CellType, VertexType}) where {CellType, VertexType}
    mesh.coord_vector
end

"""
This function is part of the interface AbstractMesh
"""
@inline function get_cell_vector(mesh::Mesh{CellType, VertexType}) where {CellType, VertexType}
    mesh.cell_vector
end

"Constructor for a 1D mesh in the interval [x1, x2]"
function Mesh1D(x1, x2, ncells)
    CellType = Cell{2}
    VertexType = Point{1}
    cells = [CellType(i, i+1) for i=1:ncells]
    h = (x2-x1)/ncells
    vertices = [VertexType(x1 + h*(i-1)) for i=1:ncells+1]
    Mesh{Cell{2}, Point{1}}(cells, vertices)
end

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
"""
    ModelEye(assembly::LensAssembly{T}, nsamples::Int = 17; pupil_radius::T = 3.0, detpixels::Int = 1000, transform::RigidBodyTransform{T} = identitytransform(T))

Geometrically accurate model of the human eye focussed at infinity with variable `pupil_radius`.
The eye is added to the provided `assembly` to create a [`CSGOpticalSystem`](@ref) with the retina of the eye as the detector.

The eye can be positioned in the scene using the `transform` argument and the resolution of the detector specified with `detpixels`.
By default the eye is directed along the positive z-axis with the vertex of the cornea at the origin.

`nsamples` determines the resolution at which accelerated surfaces within the eye are triangulated.
"""
function ModelEye(assembly::LensAssembly{T}; nsamples::Int = 17, pupil_radius::T = 3.0, detpixels::Int = 1000, transform::RigidBodyTransform{T} = identitytransform(T), detdatatype::Type{D} = Float32) where {T<:Real,D<:Real}
    cornea_front = AcceleratedParametricSurface(ZernikeSurface(6.9, radius = -7.8, conic = -0.5), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.CORNEA, OpticSim.GlassCat.Air))
    cornea_rear = Plane(SVector{3,T}(0, 0, -1), SVector{3,T}(0, 0, -3.2))
    cornea_csg = csgintersection(cornea_front, cornea_rear)
    cornea = cornea_csg(transform)

    anterior_chamber_front = AcceleratedParametricSurface(ZernikeSurface(6.0, radius = -6.7, conic = -0.3), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.AQUEOUS, OpticSim.GlassCat.EYE.CORNEA))
    anterior_chamber_rear = Plane(SVector{3,T}(0, 0, -1), SVector(0.0, 0.0, -3.2), vishalfsizeu = 6.0, vishalfsizev = 6.0, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.AQUEOUS, OpticSim.GlassCat.EYE.VITREOUS))
    anterior_chamber_csg = csgintersection(anterior_chamber_front, anterior_chamber_rear)
    anterior_chamber = leaf(anterior_chamber_csg, translation(0.0, 0.0, -0.52))(transform)

    pupil = Annulus(pupil_radius, 5.9, rotate(transform, SVector{3,T}(0, 0, 1)), transform * SVector{3,T}(0.0, 0.0, -3.62))

    lens_front = AcceleratedParametricSurface(ZernikeSurface(5.1, radius = -10.0), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_rear = AcceleratedParametricSurface(ZernikeSurface(5.1, radius = 6.0, conic = -3.25), nsamples, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_barrel = Cylinder(5.0, 5.0, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.LENS, OpticSim.GlassCat.EYE.VITREOUS))
    lens_csg = csgintersection(csgdifference(leaf(lens_front), leaf(lens_rear, translation(0.0, 0.0, -3.7))), leaf(lens_barrel, translation(0.0, 0.0, -2.0)))
    lens = leaf(lens_csg, translation(0.0, 0.0, -3.72))(transform)

    vitreous_chamber_csg = csgintersection(csgintersection(leaf(Sphere(11.0, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.VITREOUS, OpticSim.GlassCat.Air, reflectance = 0.0, transmission = 0.0)), translation(0.0, 0.0, -13.138998863513297)), leaf(Plane(0.0, 0.0, 1.0, 0.0, 0.0, -3.72, interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.VITREOUS, OpticSim.GlassCat.EYE.AQUEOUS)))), leaf(Plane(0.0, 0.0, -1.0, 0.0, 0.0, -24 + (11 - sqrt(11^2 - 10^2)), interface = NullInterface(T))))
    vitreous_chamber = vitreous_chamber_csg(transform)

    new_assembly = OpticSim.LensAssembly(assembly, cornea, anterior_chamber, pupil, vitreous_chamber, lens)
    retina = SphericalCap(11.0, asin(10 / 11), rotate(transform, SVector{3,T}(0, 0, 1)), transform * SVector{3,T}(0.0, 0.0, -24.0), interface = FresnelInterface{T}(OpticSim.GlassCat.EYE.VITREOUS, OpticSim.GlassCat.EYE.VITREOUS, reflectance = zero(T), transmission = zero(T)))

    return CSGOpticalSystem(new_assembly, retina, detpixels, detpixels, D)
end
export ModelEye

#! format: off
"""
    ArizonaEye(::Type{T} = Float64; accommodation::T = 0.0)

The popular Arizona eye model taken from [this definition](https://photonengr.com/wp-content/uploads/kbasefiles/ArizonaEyeModel.pdf).
The `accommodation` of the eye can be varied in this model.
Returns a `DataFrame` specifying the prescription of the eye model.
"""
function ArizonaEye(::Type{T} = Float64; accommodation::T = 0.0) where {T<:Real}
    # from https://photonengr.com/wp-content/uploads/kbasefiles/ArizonaEyeModel.pdf also in our documentation/papers directory
    return DataFrame(
        Name = ["Air", "Cornea", "Aqueous", "Lens", "Vitreous", "Retina"],
        Surface = [:Object, 1, 2, 3, 4, :Image],
        Radius = [Inf64, 7.8, 6.5, 12.0 - 0.4*accommodation, -5.22 + 0.2*accommodation, -13.4],
        Conic = [missing, -.25, -.25, -7.52 + 1.29*accommodation, -1.35 - 0.43*accommodation, missing],
        Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.modelglass(1.377, 57.1, 0.0), OpticSim.GlassCat.modelglass(1.337, 61.3, 0.0), OpticSim.GlassCat.modelglass(1.42 + 0.0026 * accommodation - 0.00022 * accommodation^2, 51.9, 0.0), OpticSim.GlassCat.modelglass(1.336, 61.1, 0.0), missing],
        SemiDiameter = [Inf64, 5.3, 5.3, 5.3, 5.3, 5.3],
        Thickness = [Inf64, 0.55, 2.97 - 0.04*accommodation, 3.77 + 0.04accommodation, 16.713, missing],
    )
end
export ArizonaEye
#! format: on


text:
export inverseparameters

"""

Parameters that can be used to invert to `x` from `y`
`parameters(f, types, args...)`

Produces parameters `θ` such that:

`inverseparameters(f, types, f(args...), θ) == args`

These parameters can be used for learning.
"""
function inverseparameters end


text:
module JavisViewer

using Cairo
using Gtk
using GtkReactive
using Javis

end

max_stars_repo_name: gpucce/JavisViewer

max_issues_repo_name: gpucce/JavisViewer

max_forks_repo_name: gpucce/JavisViewer


text:
module Day17
__revise_mode__ = :eval

using Advent: save_input, read_charmatrix

function countoccupied(S, I)
    I1 = oneunit(I)
    count(J -> (I != J) && J in S, I - I1:I + I1)
end

function simulate(S)
    next = deepcopy(S)
    low, high = extrema(S)
    I1 = oneunit(low)
    for I in low - I1:high + I1
        c = countoccupied(S, I)
        if I in S
            c != 2 && c != 3 && pop!(next, I)
        else
            c == 3 && push!(next, I)
        end
    end
    next
end

function read_input(filename, ndims=3)::Set{CartesianIndex}
    layer = read_charmatrix(filename)
    A = reshape(layer, (size(layer)..., ntuple(_ -> 1, ndims - 2)...))
    Set(findall(==('#'), A))
end

function solve(filename, ndims, rounds = 6)
    S = read_input(filename, ndims)
    for _ in 1:rounds
        S = simulate(S)
    end
    length(S)
end

solve_a(filename) = solve(filename, 3)
solve_b(filename) = solve(filename, 4)

solve(filename = save_input(2020, 17)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive