text:
module Kernels

using UnstructuredGrids.Helpers

export generate_dual_connections
export generate_cell_to_faces
export generate_face_to_ftype
export generate_face_to_vertices
export find_cell_to_faces
export generate_facet_to_isboundary
export generate_face_to_isboundary
export rewind_ptrs!
export length_to_ptrs!
export generate_data_and_ptrs
export UNSET
export append_ptrs
export refine_grid_connectivity
export generate_tface_to_face
export find_gface_to_face

"""
Given the faces on the boundary of each cell,
find the cells around each face.
"""
function generate_dual_connections(
  cell_to_faces_data::AbstractVector{<:Integer},
  cell_to_faces_ptrs::AbstractVector{<:Integer},
  nfaces = maximum(cell_to_faces_data))
  _face_to_cells(cell_to_faces_data,cell_to_faces_ptrs,nfaces)
end

"""
Given the vertices on the boundary of each cell,
and, for each cell type, the local vertices on
the boundary of each local face, returns the faces
on the boundary of each cell
"""
function generate_cell_to_faces(
  cell_to_vertices_data::AbstractVector{<:Integer},
  cell_to_vertices_ptrs::AbstractVector{<:Integer},
  ctype_to_lface_to_lvertices::Vector{Vector{Vector{Int}}},
  cell_to_ctype::AbstractVector{<:Integer},
  vertex_to_cells_data::AbstractVector{<:Integer},
  vertex_to_cells_ptrs::AbstractVector{<:Integer})
  _cell_to_faces(
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs)
end

function find_cell_to_faces(
  cell_to_vertices_data::AbstractVector{<:Integer},
  cell_to_vertices_ptrs::AbstractVector{<:Integer},
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data::AbstractVector{<:Integer},
  vertex_to_faces_ptrs::AbstractVector{<:Integer})
  _cell_to_faces_from_vertex_to_faces(
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_faces_data,
    vertex_to_faces_ptrs)
end

function generate_facet_to_isboundary(face_to_cells_ptrs::AbstractVector{<:Integer})
  nfaces = length(face_to_cells_ptrs)-1
  face_to_isboundary = fill(false,nfaces)
  _generate_face_to_isboundary_fill!(face_to_isboundary,face_to_cells_ptrs)
  face_to_isboundary
end

function generate_face_to_isboundary(
  facet_to_isboundary::AbstractVector{Bool},
  face_to_facets_data::AbstractVector{<:Integer},
  face_to_facets_ptrs::AbstractVector{<:Integer})
  nobjects = length(face_to_facets_ptrs)-1
  object_to_isboundary = fill(false,nobjects)
  _generate_object_to_isboundary_fill!(
    object_to_isboundary,
    facet_to_isboundary,
    face_to_facets_data,
    face_to_facets_ptrs)
  object_to_isboundary
end

function generate_face_to_ftype(
  cell_to_faces_data::AbstractVector{<:Integer},
  cell_to_faces_ptrs::AbstractVector{<:Integer},
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_lface_to_ftype::AbstractVector{<:AbstractVector{<:Integer}},
  nfaces::Int=maximum(cell_to_faces_data))
  _face_to_ftype(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_ftype,
    nfaces)
end

function generate_face_to_vertices(
  cell_to_vertices_data::AbstractVector{<:Integer},
  cell_to_vertices_ptrs::AbstractVector{<:Integer},
  cell_to_faces_data::AbstractVector{<:Integer},
  cell_to_faces_ptrs::AbstractVector{<:Integer},
  cell_to_ctype::AbstractVector{<:Integer},
  ctype_to_lface_to_lvertices::Vector{Vector{Vector{Int}}},
  nfaces::Int=maximum(cell_to_faces_data))
  _face_to_vertices(
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices,
    nfaces)
end

"""
Rewind the given vector.
"""
function rewind_ptrs!(ptrs::AbstractVector{<:Integer})
  @inbounds for i in (length(ptrs)-1):-1:1
    ptrs[i+1] = ptrs[i]
  end
  ptrs[1] = 1
end

"""
Given a vector of integers, mutate it from length state to pointer state.
"""
function length_to_ptrs!(ptrs::AbstractArray{<:Integer})
  ptrs[1] = 1
  @inbounds for i in 1:(length(ptrs)-1)
    ptrs[i+1] += ptrs[i]
  end
end

"""
Given a vector of vectors compute the corresponding data and and ptrs
"""
function generate_data_and_ptrs(vv::Vector{Vector{T}}) where T
  ptrs = Vector{T}(undef,length(vv)+1)
  _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = Vector{T}(undef,ndata)
  _generate_data_and_ptrs_fill_data!(data,vv)
  (data, ptrs)
end

function append_ptrs(pa::AbstractVector{T},pb::AbstractVector{T}) where T
  na = length(pa)-1
  nb = length(pb)-1
  p = zeros(T,(na+nb+1))
  _append_count!(p,pa,pb)
  length_to_ptrs!(p)
  p
end

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

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

  tcell_to_points_ptrs = zeros(P,ntcells+1)

  _refine_grid_connectivity_count!(
    tcell_to_points_ptrs,
    ncells,
    ltcell_to_lpoints)

  length_to_ptrs!(tcell_to_points_ptrs)

  ndata = tcell_to_points_ptrs[end]-1

  tcell_to_points_data = zeros(T,ndata)

  _refine_grid_connectivity!(
    tcell_to_points_data,
    cell_to_points_data,
    cell_to_points_ptrs,
    ltcell_to_lpoints )

  (tcell_to_points_data, tcell_to_points_ptrs)

end

function generate_tface_to_face(
  cell_to_faces_data::AbstractVector{T},
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lnodes,
  ltface_to_ltnodes,
  lface_to_lnodes,
  ntfaces) where T

  ltcell_to_lfaces =  _generate_local_face_relations(
    ltcell_to_lnodes,
    ltface_to_ltnodes,
    lface_to_lnodes)

  _generate_tface_to_face(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    tcell_to_tfaces_data,
    tcell_to_tfaces_ptrs,
    ltcell_to_lfaces,
    ntfaces)

end

function find_gface_to_face(
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data::AbstractVector{T},
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs) where T

  ngfaces = length(gface_to_nodes_ptrs) - 1
  gface_to_face = zeros(T,ngfaces)
  n = max_cells_arround_vertex(node_to_faces_ptrs)
  faces_around = fill(UNSET,n)
  faces_around_scratch = fill(UNSET,n)

  _fill_gface_to_face!(
    gface_to_face,
    face_to_nodes_data,
    face_to_nodes_ptrs,
    node_to_faces_data,
    node_to_faces_ptrs,
    gface_to_nodes_data,
    gface_to_nodes_ptrs,
    faces_around,
    faces_around_scratch)

  gface_to_face

end

# Helpers

const UNSET = 0

function _append_count!(p,pa,pb)
  na = length(pa)-1
  for ca in 1:na
    p[ca+1] = pa[ca+1] - pa[ca]
  end
  nb = length(pb)-1
  for cb in 1:nb
    p[cb+1+na] = pb[cb+1] - pb[cb]
  end
end

function _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  for (i,v) in enumerate(vv)
    ptrs[i+1] = length(v)
  end
end

function _generate_data_and_ptrs_fill_data!(data,vv)
  k = 1
  for v in vv
    for vi in v
      data[k] = vi
      k += 1
    end
  end
end

function _generate_face_to_isboundary_fill!(
  face_to_isboundary, face_to_cells_ptrs)
  nfaces = length(face_to_isboundary)
  for face in 1:nfaces
    ncells_around = face_to_cells_ptrs[face+1] - face_to_cells_ptrs[face]
    if ncells_around == 1
      face_to_isboundary[face] = true
    elseif ncells_around == 2
      face_to_isboundary[face] = false
    else
      @unreachable
    end
  end
end

function  _generate_object_to_isboundary_fill!(
  object_to_isboundary,
  face_to_isboundary,
  object_to_faces_data,
  object_to_faces_ptrs)

  nobjects = length(object_to_isboundary)
  for object in 1:nobjects
    a = object_to_faces_ptrs[object]-1
    b = object_to_faces_ptrs[object+1]
    nlfaces = b-(a+1)
    for lface in 1:nlfaces
      face = object_to_faces_data[a+lface]
      isboundary = face_to_isboundary[face]
      if isboundary
        object_to_isboundary[object] = true
        break
      end
    end
  end

end

function _face_to_cells(
  cell_to_faces_data::AbstractVector{L},
  cell_to_faces_ptrs::AbstractVector{P},
  nfaces) where {L,P}

  face_to_cells_ptrs = zeros(P,nfaces+1)

  _face_to_cells_count!(
    face_to_cells_ptrs, cell_to_faces_data, cell_to_faces_ptrs)

  length_to_ptrs!(face_to_cells_ptrs)

  ndata = face_to_cells_ptrs[end]-1

  face_to_cells_data = Vector{L}(undef,ndata)

  _face_to_cells_fill!(
    face_to_cells_data, face_to_cells_ptrs,
    cell_to_faces_data, cell_to_faces_ptrs)

  rewind_ptrs!(face_to_cells_ptrs)

  (face_to_cells_data, face_to_cells_ptrs)

end

function _face_to_cells_count!(
    face_to_cells_ptrs, cell_to_faces_data, cell_to_faces_ptrs)

  ncells = length(cell_to_faces_ptrs) - 1
  for cell in 1:ncells
    a = cell_to_faces_ptrs[cell]
    b = cell_to_faces_ptrs[cell+1]-1
    for p in a:b
      face = cell_to_faces_data[p]
      face_to_cells_ptrs[face+1] += 1
    end
  end

end

function _face_to_cells_fill!(
    face_to_cells_data, face_to_cells_ptrs,
    cell_to_faces_data, cell_to_faces_ptrs)

  ncells = length(cell_to_faces_ptrs) - 1
  for cell in 1:ncells
    a = cell_to_faces_ptrs[cell]
    b = cell_to_faces_ptrs[cell+1]-1
    for p in a:b
      face = cell_to_faces_data[p]
      q = face_to_cells_ptrs[face]
      face_to_cells_data[q] = cell
      face_to_cells_ptrs[face] += 1
    end
  end

end

function _cell_to_faces(
    cell_to_vertices_data::AbstractVector{L},
    cell_to_vertices_ptrs::AbstractVector{P},
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs) where {L,P}

  cell_to_faces_ptrs = _cell_to_faces_count(
    cell_to_vertices_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(cell_to_faces_ptrs)

  ndata = cell_to_faces_ptrs[end]-1

  cell_to_faces_data = fill(L(UNSET),ndata)

  nvertices = max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  vertices = fill(UNSET,nvertices)
  vertices_scratch = fill(UNSET,nvertices)

  nvertices = max_cells_arround_vertex(vertex_to_cells_ptrs)
  cells_around = fill(UNSET,nvertices)
  cells_around_scratch = fill(UNSET,nvertices)
  lface_of_cells_around = fill(UNSET,nvertices)

  _cell_to_faces_fill!(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs,
    vertices,
    vertices_scratch,
    cells_around,
    cells_around_scratch,
    lface_of_cells_around)

  (cell_to_faces_data, cell_to_faces_ptrs)

end

function  _cell_to_faces_from_vertex_to_faces(
  cell_to_vertices_data::AbstractVector{L},
  cell_to_vertices_ptrs::AbstractVector{P},
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data,
  vertex_to_faces_ptrs) where {L,P}

  cell_to_faces_ptrs = _cell_to_faces_count(
    cell_to_vertices_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(cell_to_faces_ptrs)

  ndata = cell_to_faces_ptrs[end]-1

  cell_to_faces_data = fill(L(UNSET),ndata)

  nvertices = max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  vertices = fill(UNSET,nvertices)

  nvertices = max_cells_arround_vertex(vertex_to_faces_ptrs)
  faces_around = fill(UNSET,nvertices)
  faces_around_scratch = fill(UNSET,nvertices)

  _cell_to_faces_fill!(
    cell_to_faces_ptrs,
    cell_to_faces_data,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_faces_data,
    vertex_to_faces_ptrs,
    vertices,
    faces_around,
    faces_around_scratch)

  (cell_to_faces_data, cell_to_faces_ptrs)
end

function _cell_to_faces_count(
  cell_to_vertices_ptrs::AbstractVector{P},
  cell_to_ctype,
  ctype_to_lface_to_lvertices) where P

  ncells = length(cell_to_vertices_ptrs) - 1

  cell_to_faces_ptrs = zeros(P,ncells+1)

  type_to_nlfaces = _type_to_nlfaces(ctype_to_lface_to_lvertices)

  _cell_to_faces_count!(
    cell_to_faces_ptrs,
    type_to_nlfaces,
    cell_to_ctype)

  cell_to_faces_ptrs

end

function _cell_to_faces_fill!(
  cell_to_faces_ptrs,
  cell_to_faces_data,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data,
  vertex_to_faces_ptrs,
  vertices,
  faces_around,
  faces_around_scratch)

  ncells = length(cell_to_ctype)

  for cell in 1:ncells

    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    nlfaces = length(lface_to_lvertices)
    a = cell_to_faces_ptrs[cell]-1

    for lface in 1:nlfaces

      _fill_vertices_in_lface!(
        vertices,
        lface,
        lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs)

      _find_cells_around_vertices!(
        faces_around,
        faces_around_scratch,
        vertices,
        vertex_to_faces_data,
        vertex_to_faces_ptrs)

      for face in faces_around
        if face != UNSET
          cell_to_faces_data[a+lface] = face
          break
        end
      end

    end

  end

end

function max_nvertices_in_lface(ctype_to_lface_to_lvertices)
  n = 0
  for lface_to_lvertices in ctype_to_lface_to_lvertices
    for lvertices in lface_to_lvertices
      m = length(lvertices)
      n = max(n,m)
    end
  end
  n
end

function max_cells_arround_vertex(vertex_to_cells_ptrs)
  n = 0
  nvertices = length(vertex_to_cells_ptrs)-1
  for vertex in 1:nvertices
    m = vertex_to_cells_ptrs[vertex+1]-vertex_to_cells_ptrs[vertex]
    n = max(n,m)
  end
  n
end

function _type_to_nlfaces(ctype_to_lface_to_lvertices)
  [ length(lface_to_lvertices)
   for lface_to_lvertices in ctype_to_lface_to_lvertices]
end

function _cell_to_faces_count!(
    cell_to_faces_ptrs,
    type_to_nlfaces,
    cell_to_ctype)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    nlfaces = type_to_nlfaces[ctype]
    cell_to_faces_ptrs[cell+1] = nlfaces
  end

end

function  _cell_to_faces_fill!(
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    ctype_to_lface_to_lvertices,
    cell_to_ctype,
    vertex_to_cells_data,
    vertex_to_cells_ptrs,
    vertices,
    vertices_scratch,
    cells_around,
    cells_around_scratch,
    lface_of_cells_around)

  face = 1

  ncells = length(cell_to_ctype)

  for cell in 1:ncells

    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    nlfaces = length(lface_to_lvertices)
    a = cell_to_faces_ptrs[cell]-1

    for lface in 1:nlfaces

      if cell_to_faces_data[a+lface] != UNSET
        continue
      end

      _fill_vertices_in_lface!(
        vertices,
        lface,
        lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs)

      _find_cells_around_vertices!(
        cells_around,
        cells_around_scratch,
        vertices,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)

      _find_lface_of_cells_around_lface!(
        lface_of_cells_around,
        cells_around,
        lface,
        ctype_to_lface_to_lvertices,
        cell,
        cell_to_vertices_data,
        cell_to_vertices_ptrs,
        cell_to_ctype,
        vertices,
        vertices_scratch)

      _fill_face_in_cells_arround!(
        cell_to_faces_data,
        cell_to_faces_ptrs,
        face,
        cells_around,
        lface_of_cells_around)

      face += 1
    end

  end

end

function _fill_face_in_cells_arround!(
  cell_to_faces_data,
  cell_to_faces_ptrs,
  face,
  cells_around,
  lface_of_cells_around)

  for icell_around in 1:length(cells_around)
    cell_around = cells_around[icell_around]
    if cell_around == UNSET
      continue
    end
    lface = lface_of_cells_around[icell_around]
    f = cell_to_faces_ptrs[cell_around]-1
    cell_to_faces_data[f+lface] = face
  end

end

function _find_cells_around_vertices!(
  cells_around,
  cells_around_scratch,
  vertices,
  vertex_to_cells_data,
  vertex_to_cells_ptrs)

  ncells_around = UNSET
  ncells_around_scratch = UNSET

  for ivertex in 1:length(vertices)
    vertex = vertices[ivertex]
    if vertex == UNSET
      continue
    end
    if ivertex == 1
      ncells_around = _fill_cells_around_scratch!(
        cells_around,
        vertex,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)
    else
      ncells_around_scratch = _fill_cells_around_scratch!(
        cells_around_scratch,
        vertex,
        vertex_to_cells_data,
        vertex_to_cells_ptrs)
      _set_intersection!(
        cells_around,cells_around_scratch,
        ncells_around,ncells_around_scratch)
    end
  end

end

function _fill_cells_around_scratch!(
  cells_around_scratch,
  vertex,
  vertex_to_cells_data,
  vertex_to_cells_ptrs)

  cells_around_scratch .= UNSET
  d = vertex_to_cells_ptrs[vertex] - 1
  ncells_around = vertex_to_cells_ptrs[vertex+1] - (d + 1)
  for icell_around in 1:ncells_around
    cell_around = vertex_to_cells_data[d+icell_around]
    cells_around_scratch[icell_around] = cell_around
  end
  ncells_around

end

function _set_intersection!(
  cells_around,cells_around_scratch,
  ncells_around,ncells_around_scratch)
  for i in 1:ncells_around
    if cells_around[i] == UNSET
      continue
    end
    _find_eq!(i,cells_around,cells_around_scratch,ncells_around_scratch)
  end
end

function _find_eq!(i,cells_around,cells_around_scratch,ncells_around_scratch)
  for j in 1:ncells_around_scratch
    if cells_around[i] == cells_around_scratch[j]
      return
    end
  end
  cells_around[i] = UNSET
  return
end

function _find_lface_of_cells_around_lface!(
  lface_of_cells_around,
  cells_around,
  lface,
  ctype_to_lface_to_lvertices,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_ctype,
  vertices,
  vertices_scratch)

  lface_of_cells_around .= UNSET

  for icell_around in 1:length(cells_around)
    cell_around = cells_around[icell_around]
    if cell_around == UNSET
      continue
    end

    lface_of_cell_around = _find_lface_with_same_vertices(
      vertices,
      vertices_scratch,
      cell_around,
      cell_to_vertices_data,
      cell_to_vertices_ptrs,
      cell_to_ctype,
      ctype_to_lface_to_lvertices)

    lface_of_cells_around[icell_around] = lface_of_cell_around

  end

end

function _find_lface_with_same_vertices(
  vertices,
  vertices_scratch,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ctype = cell_to_ctype[cell]
  lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]

  nlfaces = length(lface_to_lvertices)
  c = cell_to_vertices_ptrs[cell]-1

  for lface in 1:nlfaces

    _fill_vertices_in_lface!(
      vertices_scratch,
      lface,
      lface_to_lvertices,
      cell,
      cell_to_vertices_data,
      cell_to_vertices_ptrs)

    if _set_equal(vertices,vertices_scratch)
      return lface
    end

  end

  return UNSET

end

function _fill_vertices_in_lface!(
  vertices,
  lface,
  lface_to_lvertices,
  cell,
  cell_to_vertices_data,
  cell_to_vertices_ptrs)

  vertices .= UNSET
  c = cell_to_vertices_ptrs[cell] - 1
  lvertices = lface_to_lvertices[lface]
  nlfvertex = length(lvertices)
  for lfvertex in 1:nlfvertex
    lvertex = lvertices[lfvertex]
    vertex = cell_to_vertices_data[c+lvertex]
    vertices[lfvertex] = vertex
  end

end

function _set_equal(vertices,vertices_scratch)

  b = _is_subset(vertices,vertices_scratch)
  if b == false; return false; end
  b = _is_subset(vertices_scratch,vertices)
  if b == false; return false; end
  return true

end

function _is_subset(vertices,vertices_scratch)
  for i in 1:length(vertices)
    v = vertices[i]
    if v == UNSET
      continue
    end
    b = _find_eq(v,vertices_scratch)
    if b == false; return false; end
  end
  return true
end

function _find_eq(v,vertices_scratch)
  for vs in vertices_scratch
    if v == vs
      return true
    end
  end
  return false
end

function _face_to_ftype(
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_ftype,
  nfaces)
  face_to_ftype = fill(UNSET,nfaces)
  _face_to_ftype_fill!(
    face_to_ftype,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_ftype)
  face_to_ftype
end

function _face_to_ftype_fill!(
  face_to_ftype,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_ftype)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_ftype = ctype_to_lface_to_ftype[ctype]
    a = cell_to_faces_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      if face_to_ftype[face] != UNSET
        continue
      end
      ftype = lface_to_ftype[lface]
      face_to_ftype[face] = ftype
    end
  end

end

function _face_to_vertices(
  cell_to_vertices_data::AbstractVector{L},
  cell_to_vertices_ptrs::AbstractVector{P},
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices,
  nfaces) where {L,P}

  face_to_vertices_ptrs = fill(P(UNSET),nfaces+1)

  _face_to_vertices_count!(
    face_to_vertices_ptrs,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  length_to_ptrs!(face_to_vertices_ptrs)
  ndata = face_to_vertices_ptrs[end]-1
  face_to_vertices_data = fill(L(UNSET),ndata)

  _face_to_vertices_fill!(
    face_to_vertices_data,
    face_to_vertices_ptrs,
    cell_to_vertices_data,
    cell_to_vertices_ptrs,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    cell_to_ctype,
    ctype_to_lface_to_lvertices)

  (face_to_vertices_data, face_to_vertices_ptrs)

end

function _face_to_vertices_count!(
  face_to_vertices_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    a = cell_to_faces_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      if face_to_vertices_ptrs[face+1] != UNSET
        continue
      end
      nfvertices = length(lface_to_lvertices[lface])
      face_to_vertices_ptrs[face+1] = nfvertices
    end
  end

end

function _face_to_vertices_fill!(
  face_to_vertices_data,
  face_to_vertices_ptrs,
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

  ncells = length(cell_to_ctype)
  for cell in 1:ncells
    ctype = cell_to_ctype[cell]
    lface_to_lvertices = ctype_to_lface_to_lvertices[ctype]
    a = cell_to_faces_ptrs[cell]-1
    c = cell_to_vertices_ptrs[cell]-1
    nlfaces = cell_to_faces_ptrs[cell+1] - (a + 1)
    for lface in 1:nlfaces
      face = cell_to_faces_data[a+lface]
      v = face_to_vertices_ptrs[face]-1
      if face_to_vertices_data[v+1] != UNSET
        continue
      end
      lvertices = lface_to_lvertices[lface]
      nfvertices = length(lvertices)
      for lfvertex in 1:nfvertices
        lvertex = lvertices[lfvertex]
        vertex = cell_to_vertices_data[c+lvertex]
        face_to_vertices_data[v+lfvertex] = vertex
      end
    end
  end

end

function  _refine_grid_connectivity_count!(
    tcell_to_points_ptrs,
    ncells,
    ltcell_to_lpoints)

  tcell = 1

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

end

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

  ncells = length(cell_to_points_ptrs) - 1

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

end

function _generate_local_face_relations(
  ltcell_to_lnodes, ltface_to_ltnodes, lface_to_lnodes)

  nltcells = length(ltcell_to_lnodes)
  nltfaces = length(ltface_to_ltnodes)
  nlfaces = length(lface_to_lnodes)

  ltcell_ltface_to_lface = [ zeros(Int,nltfaces) for ltcell in 1:nltcells  ]

  for ltcell in 1:nltcells

    ltnode_to_lnode = ltcell_to_lnodes[ltcell]

    for ltface in 1:nltfaces
      ltnodes = ltface_to_ltnodes[ltface]
      for lface in 1:nlfaces
        lnodes = lface_to_lnodes[lface]
        allin = true
        for ltnode in ltnodes
          lnode = ltnode_to_lnode[ltnode]
          if !(lnode in lnodes)
            allin = false
            break
          end
        end
        if allin
          ltcell_ltface_to_lface[ltcell][ltface] = lface
          break
        end
      end
    end

  end

  ltcell_ltface_to_lface

end

function _generate_tface_to_face(
  cell_to_faces_data::AbstractVector{T},
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lfaces,
  ntfaces) where T

  tface_to_face = zeros(T,ntfaces)

  _generate_tface_to_face!(
    tface_to_face,
    cell_to_faces_data,
    cell_to_faces_ptrs,
    tcell_to_tfaces_data,
    tcell_to_tfaces_ptrs,
    ltcell_to_lfaces)

  tface_to_face

end

function  _generate_tface_to_face!(
  tface_to_face,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  tcell_to_tfaces_data,
  tcell_to_tfaces_ptrs,
  ltcell_to_lfaces)

  ncells = length(cell_to_faces_ptrs)-1
  nltcells = length(ltcell_to_lfaces)

  tcell = 1
  for cell in 1:ncells
    c = cell_to_faces_ptrs[cell]-1
    for ltcell in 1:nltcells
      a = tcell_to_tfaces_ptrs[tcell]-1
      b = tcell_to_tfaces_ptrs[tcell+1]
      nltfaces = b - (a+1)
      ltface_to_lface = ltcell_to_lfaces[ltcell]
      for ltface in 1:nltfaces
        tface = tcell_to_tfaces_data[a+ltface]
        lface = ltface_to_lface[ltface]
        if lface != 0
          face = cell_to_faces_data[c+lface]
          tface_to_face[tface] = face
        end
      end
      tcell += 1
    end
  end

end

function  _fill_gface_to_face!(
  gface_to_face,
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data,
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs,
  faces_around,
  faces_around_scratch)

  ngfaces = length(gface_to_nodes_ptrs) - 1

  nfaces_around = UNSET
  nfaces_around_scratch = UNSET

  for gface in 1:ngfaces

    a = gface_to_nodes_ptrs[gface]-1
    b = gface_to_nodes_ptrs[gface+1]
    nlnodes = b-(a+1)

    for lnode in 1:nlnodes
      node = gface_to_nodes_data[lnode+a]
      if lnode == 1
        nfaces_around = _fill_cells_around_scratch!(
          faces_around,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
      else
        nfaces_around_scratch = _fill_cells_around_scratch!(
          faces_around_scratch,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
        _set_intersection!(
          faces_around,faces_around_scratch,
          nfaces_around,nfaces_around_scratch)
      end
    end

    for face in faces_around
      if face != UNSET
        gface_to_face[gface] = face
        break
      end
    end

  end

end

end # module Kernels


text:
using PyPlot
using LinearAlgebra
using PyTensorFlow
using PyCall
np = pyimport("numpy")
include("poisson_op.jl")

# Solver parameters
m = 20
n = 20
h = 1.0 # 50 meters
T = 0.01 # 100 days
NT = 1000
Δt = T/NT
x = (1:m)*h|>collect
z = (1:n)*h|>collect
X, Z = np.meshgrid(x, z)

# todo
#=
Krw, Kro -- function mxn tensor to mxn tensor
μw, μo -- known mxn
ρw, ρo -- known mxn
K -- scalar
g -- constant ≈ 9.8?
ϕ -- known mxn
=#
function Krw(Sw)
    return Sw ^ 2
end

function Kro(So)
    return So ^2
end
ρw = constant(1.)
ρo = constant(1.)
μw = constant(1.)
μo = constant(1.)
K = Variable(ones(m,n))
g = constant(1.0)
ϕ = Variable(0.25 .* ones(m,n))

function geto(o::Union{Array,PyObject}, i::Int64, j::Int64)
    if i==-1
        ii = 1:m-2
    elseif i==0
        ii = 2:m-1
    else
        ii = 3:m
    end
    if j==-1
        jj = 1:n-2
    elseif j==0
        jj = 2:n-1
    else
        jj = 3:n
    end
    return o[ii,jj]
end

function G(f, p)
    f1 = (geto(f, 0, 0) + geto(f, 1, 0))/2
    f2 = (geto(f, -1, 0) + geto(f, 0, 0))/2
    f3 = (geto(f,0,1) + geto(f,0,0))/2
    f4 = (geto(f,0,-1) + geto(f,0,0))/2
    rhs = -f1.*(geto(p,1,0)-geto(p,0,0)) +
            f2.*(geto(p,0,0)-geto(p,-1,0)) -
            f3.*(geto(p,0,1)-geto(p,0,0)) +
            f4.*(geto(p,0,0)-geto(p,0,-1))
    local q
    if isa(rhs, Array)
        q = zeros(m, n)
        q[2:m-1, 2:n-1] = rhs/h^2
    else
        q = constant(zeros(m, n))
        q = scatter_add(q, 2:m-1, 2:n-1, rhs/h^2)
        # q[2:m-1, 2:n-1] += rhs/h^2
    end
    q
end

# variables : sw, u, v, p
# (time dependent) parameters: qw, qo, ϕ
function onestep(sw, qw, qo)
    # step 1: update p
    λw = Krw(sw)/μw
    λo = Kro(1-sw)/μo
    λ = λw + λo
    f = λw/λ
    q = qw + qo
    Θ = G((λw*ρw+λo*ρo)*g, Z)
    p = poisson_op(λ.*K, Θ+q, constant(h))

    # step 2: update u, v
    rhs_u = -(geto(K, 0, 0)+geto(K,1,0))/2.0 .* (geto(λ, 0, 0) + get(λ, 1, 0))/2h .* (geto(p, 1, 0) - geto(p, 0, 0))
    rhs_v = -(geto(K, 0, 0)+geto(K,0,1))/2.0 .* (geto(λ, 0, 0) + get(λ, 0, 1))/2h .* (geto(p, 0, 1) - geto(p, 0, 0)) +
            (geto(K, 0, 0)+geto(K,0,1))/2.0 .* (geto(λw*ρw+λo*ρo, 0, 0)+geto(λw*ρw+λo*ρo, 0, 1))/2 * g
    u = constant(zeros(m, n))
    v = constant(zeros(m, n))
    u = scatter_add(u, 2:m-1, 2:n-1, rhs_u)
    v = scatter_add(v, 2:m-1, 2:n-1, rhs_v)

    # step 3: update sw
    rhs = geto(qw, 0, 0) - (geto(f, 1, 0)-geto(f, 0, 0))/h.*geto(u, 0, 0) -
            (geto(f, 0, 1)-geto(f, 0, 0))/h.*geto(v, 0, 0) -
            geto(f, 0, 0) .* (
                (geto(u, 0, 0)-geto(u, -1, 0))/h + (geto(v, 0, 0)-geto(v, 0, -1))/h
            ) -
            geto(G(K.*f.*λo*(ρw-ρo)*g, Z), 0, 0)
    rhs = Δt*rhs/geto(ϕ, 0, 0)
    sw = scatter_add(sw, 2:m-1, 2:n-1, rhs)
    return sw, p
end

"""
solve(qw, qo, sw0, p0)
Solve the two phase flow equation.
`qw` and `qo` -- `NT x m x n` numerical array, `qw[i,:,:]` the corresponding value of qw at i*Δt
`sw0` and `p0` -- initial value for `sw` and `p`. `m x n` numerical array.
"""
function solve(qw, qo, sw0)
    qw_arr = constant(qw) # qw: NT x m x n array
    qo_arr = constant(qo)
    function condition(i, tas...)
        i <= NT
    end
    function body(i, tas...)
        ta_sw, ta_p = tas
        sw, p = onestep(read(ta_sw, i), qw_arr[i], qo_arr[i])
        ta_sw = write(ta_sw, i+1, sw)
        ta_p = write(ta_p, i+1, p)
        i+1, ta_sw, ta_p
    end
    ta_sw, ta_p = TensorArray(NT+1), TensorArray(NT+1)
    ta_sw = write(ta_sw, 1, constant(sw0))
    i = constant(1, dtype=Int32)
    _, ta_sw, ta_p = while_loop(condition, body, [i; ta_sw; ta_p])
    out_sw, out_p = stack(ta_sw), stack(ta_p)
end

function vis(val, args...;kwargs...)
    close("all")
    ns = Int64.(round.(LinRange(1,size(val,1),9)))
    for i = 1:9
        subplot(330+i)
        imshow(val[ns[i],:,:], args...;kwargs...)
        colorbar()
    end
end

qw = zeros(NT, m, n)
qw[:,15,5] .= 1.0
qo = zeros(NT, m, n)
sw0 = zeros(m, n)
out_sw, out_p = solve(qw, qo, sw0)

sess = Session(); init(sess)
S, P = run(sess, [out_sw, out_p])
vis(S)
# figure()
# vis(P)
error("stop")

#=
# Step 1: Assign numerical values to qw, qo, sw0, p0
# qw =
# qo =
# sw0 =
# p0 =
qw = zeros(NT, m, n)
qw[:,15,5] .= 0.0018
qo = zeros(NT, m, n)
sw0 = zeros(m, n)
p0 = 3.0337e+07*ones(m,n)
# # Step 2: Construct Graph
out_sw, out_p = solve(qw, qo, sw0, p0)
# # Step 3: Run
sess = Session()
init(sess)
sw, p = run(sess, [out_sw, out_p])
# # Step 4: Visualize
# vis(sw)
=#

ext: jl


text:
## gaussian_random_fields.jl : Gaussian random field types and generator types

abstract type GaussianRandomFieldGenerator end

struct GaussianRandomField{G<:GaussianRandomFieldGenerator,C<:AbstractCovarianceFunction,
                           P,M,D}
    mean::M
    cov::C
    pts::P
    data::D
end

"""
	GaussianRandomField(mean,cov,method,pts...)
	GaussianRandomField(cov,method,pts...)
	GaussianRandomField(mean,cov,method,nodes,elements)
	GaussianRandomField(cov,method,nodes,elements)

Compute a Gaussian random field with mean `mean` and covariance structure `cov` defined in the points `pts`, and computed using the method `method`.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> mn = ones(length(pts1),length(pts2))

julia> grf = GaussianRandomField(mn,c,Cholesky(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a Cholesky decomposition

```
If no `mean` is specified, a zero-mean Gaussian random field is assumed.
```jldoctest
julia> grf = GaussianRandomField(c,Cholesky(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a Cholesky decomposition

```
The  Gaussian random field sampler `method` can be `Cholesky()`, `Spectral()`, `KarhunenLoeve(n)` (where `n` is the number of terms in the expansion), or `CirculantEmbedding()`. The dimension of the points must match the dimension of the covariance function `cov`. The points can be specified as arguments of type `AbstractVector`, in which case a tensor (Kronecker) product is assumed, or as a Finite Element mesh with node table `nodes` and element table `elements`.
```jldoctest
julia> grf = GaussianRandomField(c,KarhunenLoeve(500),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 500 terms

```
For separable Gaussian random fields with a KL expansion, provide a [`SeparableCovarianceFunction`](@ref). Note that the number of terms `n` in `KarhunenLoeve(n)` refers to the total number of terms.
```jldoctest
julia> e1 = Exponential(0.1); e2 = Exponential(0.01);

julia> scov = SeparableCovarianceFunction(e1,e2)
2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), exponential (λ=0.01, σ=1.0, p=2.0) ]

julia> grf = GaussianRandomField(scov,KarhunenLoeve(500),pts1,pts2)
Gaussian random field with 2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), exponential (λ=0.01, σ=1.0, p=2.0) ] on a 51x51 structured grid, using a KL expansion with 500 terms

julia> plot_eigenfunction(grf,3); show()
[...]

```
Also anisotropic Gaussian random fields can be computed. For example, the anisotropic exponential covariance function needs a positive definite matrix `A`. The size of the off-diagonal elements of `A` determine the degree of anisotropy.
```jldoctest
julia> a = AnisotropicExponential([1 0.8;0.8 1])
anisotropic exponential (A=[1.0 0.8; 0.8 1.0], σ=1.0)

julia> acov = CovarianceFunction(2,a)
2d anisotropic covariance function exponential (A=[1.0 0.8; 0.8 1.0], σ=1.0)

julia> pts = range(0,stop = 10,length = 128)
0.0:0.07874015748031496:10.0

julia> grf = GaussianRandomField(acov,CirculantEmbedding(),pts,pts)
WARNING: negative eigenvalue -1.2245106889248049e-18 detected, Gaussian random field will be approximated (ignoring all negative eigenvalues)
WARNING: increase padding if possible
Gaussian random field with 2d anisotropic covariance function exponential (A=[1.0 0.8; 0.8 1.0], σ=1.0) on a 128x128 structured grid, using a circulant embedding

julia> contourf(grf)
[...]

```
Samples from the random field can be computed using the `sample` function.
```jldoctest
julia> sample(grf)
[...]

```
See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref), [`sample`](@ref)
"""
function GaussianRandomField(mean::Array{<:Real}, cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d
    size(mean) == length.(pts) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))

    method isa CirculantEmbedding && !(pts isa NTuple{d,AbstractRange}) &&
        throw(ArgumentError("can only use circulant embedding on a regular grid, supply ranges for pts"))

    method isa Union{CirculantEmbedding,KarhunenLoeve} && any(pt -> length(pt) < 2, pts) &&
        throw(ArgumentError("must have at least 2 points in each direction to use circulant embedding or KL expansion"))

    _GaussianRandomField(mean, cov, method, pts...; kwargs...)
end

# zero-mean GRF
GaussianRandomField(cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d = GaussianRandomField(zeros(eltype(cov), length.(pts)), cov, method, pts...; kwargs...)

# constant mean GRF
GaussianRandomField(mean::Real, cov::CovarianceFunction{d}, method::GaussianRandomFieldGenerator, pts::Vararg{AbstractVector,d}; kwargs...) where d = GaussianRandomField(fill(convert(eltype(cov), mean), length.(pts)), cov, method, pts...; kwargs...)

# obtain generator used to compute random field
generator(::GaussianRandomField{G}) where G = G()

"""
	sample(grf)
	sample(grf, xi=randn(randdim(grf)))

Take a sample from the Gaussian Random Field `grf` using the (optional) random numbers `xi`. The vector`xi` must have appropriate length. The default value is `randn(randdim(grf))`.

# Examples
```
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

julia> sample(grf)
[...]

julia> sample(grf,xi=2*rand(randdim(grf))-1)
[...]

```
"""
function sample(grf::GaussianRandomField; xi::AbstractArray{<:Real}=randn(randdim(grf)))
    length(xi) == prod(randdim(grf)) || throw(DimensionMismatch("length of random points vector must be equal to $(randdim(grf))"))
    _sample(grf,xi)
end

function Base.show(io::IO, grf::GaussianRandomField)
    print(io, "Gaussian random field with ", grf.cov, " on a ")
    showpoints(io, grf.pts)
    print(io, ", using a ", generator(grf))
end

function showpoints(io::IO, points)
    d = length(points)
    if d == 1
        print(io, length(points[1]), "-point")
    else
        print(io, length(points[1]))
        @inbounds for i in 2:d
            print(io, "×", length(points[i]))
        end
    end
    print(io, " structured grid")
end


text:
using Documenter, Schrute

makedocs(;
    modules=[Schrute],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "../tutorial.html"
    ],
    repo="https://github.com/bradlindblad/Schrute.jl/blob/{commit}{path}#L{line}",
    sitename="Schrute.jl",
    authors="Brad Lindblad",
    assets=String[],
)

deploydocs(;
    repo="github.com/bradlindblad/Schrute.jl",
)


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The diet problem

# Solve the classic "diet problem", also known as the [Stigler diet](https://en.wikipedia.org/wiki/Stigler_diet).
# The code is based on [an example from Gurobi](https://www.gurobi.com/documentation/9.0/examples/diet_cpp_cpp.html).

using JuMP
import GLPK
import Test

function print_solution(is_optimal, foods, buy)
    println("RESULTS:")
    if is_optimal
        for food in foods
            println("  $(food) = $(value(buy[food]))")
        end
    else
        println("The solver did not find an optimal solution.")
    end
end

function example_diet(; verbose = true)
    ## Nutrition guidelines
    categories = ["calories", "protein", "fat", "sodium"]
    category_data = Containers.DenseAxisArray([
        1800 2200;
        91   Inf;
        0    65;
        0    1779
        ], categories, ["min", "max"]
    )
    Test.@test category_data["protein", "min"] == 91.0
    Test.@test category_data["sodium", "max"] == 1779.0
    ## Foods
    foods = [
        "hamburger", "chicken", "hot dog", "fries", "macaroni", "pizza",
        "salad", "milk", "ice cream",
    ]
    cost = Containers.DenseAxisArray(
        [2.49, 2.89, 1.50, 1.89, 2.09, 1.99, 2.49, 0.89, 1.59],
        foods
    )
    food_data = Containers.DenseAxisArray(
        [
            410 24 26 730;
            420 32 10 1190;
            560 20 32 1800;
            380  4 19 270;
            320 12 10 930;
            320 15 12 820;
            320 31 12 1230;
            100  8 2.5 125;
            330  8 10 180
        ], foods, categories
    )
    Test.@test food_data["hamburger", "calories"] == 410.0
    Test.@test food_data["milk", "fat"] == 2.5
    ## Build model
    model = Model(GLPK.Optimizer)
    @variables(model, begin
        ## Variables for nutrition info
        category_data[c, "min"] <= nutrition[c = categories] <= category_data[c, "max"]
        ## Variables for which foods to buy
        buy[foods] >= 0
    end)
    ## Objective - minimize cost
    @objective(model, Min, sum(cost[f] * buy[f] for f in foods))
    ## Nutrition constraints
    @constraint(model, [c in categories],
        sum(food_data[f, c] * buy[f] for f in foods) == nutrition[c]
    )
    ## Solve
    if verbose
        println("Solving original problem...")
    end
    optimize!(model)
    term_status = termination_status(model)
    is_optimal = term_status == MOI.OPTIMAL
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 11.8288 atol = 1e-4
    if verbose
        print_solution(is_optimal, foods, buy)
    end
    ## Limit dairy (note that the problem will become infeasible).
    @constraint(model, buy["milk"] + buy["ice cream"] <= 6)
    if verbose
        println("Solving dairy-limited problem...")
    end
    optimize!(model)
    Test.@test termination_status(model) == MOI.INFEASIBLE
    Test.@test primal_status(model) == MOI.NO_SOLUTION
    if verbose
        print_solution(false, foods, buy)
    end
    return
end

example_diet()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Help Document
# -------------
#
# By convention, variable `i` is used to represent the index (or position) of a
# bit vector and `j` is used to represent the count (or cardinality) of a bit
# vector.

"""
    rank0(rb, i)

Count the number of 0s (`false`s) within `bv[1:i]`.
"""
rank0

"""
    rank1(bv, i)

Count the number of 1s (`true`s) within `bv[1:i]`.
"""
rank1

"""
    rank(x, bv, i)

Count the number of `x`s within `bv[1:i]`.
"""
rank

"""
    select0(bv, j)

Return the position of the `j`-th occurrence of 0 in `bv`.
"""
select0

"""
    select1(bv, j)

Return the position of the `j`-th occurrence of 1 in `bv`.
"""
select1

"""
    select(x, bv, j)

Return the position of the `j`-th occurrence of `x` in `bv`.
"""
select

"""
    search(x, bv, i)

Search the position of the next `x` in `bv` starting from `i`.
"""
search

"""
    search0(bv, i)

Search the position of the next 0 in `bv` starting from `i`.
"""
search0

"""
    search1(bv, i)

Search the position of the next 1 in `bv` starting from `i`.
"""
search1

"""
    rsearch(x, bv, i)

Search the position of the previous `x` in `bv` starting from `i`.
"""
rsearch

"""
    rsearch0(bv, i)

Search the position of the previous 0 in `bv` starting from `i`.
"""
rsearch0

"""
    rsearch1(bv, i)

Search the position of the previous 1 in `bv` starting from `i`.
"""
rsearch1


text:
#Test start to finish

using SimpleNeuralNet
using LinearAlgebra #for norm()
using Printf

#Build the neural net
nNeurons = [4; 5; 2]; #number of neurons by layer
nn = neuralNet(nNeurons);

#Define some important functions
SimpleNeuralNet.nnAct(z) = sigmoid.(z); #activation function (assume the same everywhere)
SimpleNeuralNet.nnDAct(z) = dsigmoid.(z); #derivative of activation function (assume the same everywhere)
SimpleNeuralNet.nnCost(datay,y) = 0.5*norm(y-datay,2)^2;
SimpleNeuralNet.nnDCost(datay,y) = y - datay;

#Initialize the weights and biases to random values
nn.w = rand(length(nn.w));
nn.b = rand(length(nn.b));

#Define the training parameters
lr=1.0;      #learning rate
epochs=1000; #epochs

#Define the dataset
nData = 10;
datax = rand(nn.nNeurons[1],nData);   #draw input dataset at random
datay = rand(nn.nNeurons[end],nData); #draw output dataset at random

#Compute the initial costs
c0 = nnComputeCosts(nn,datax,datay;print=false);

#Train
nnTrain(nn,datax,datay,epochs,lr;verbose=true);

#Compute the final costs
cf = nnComputeCosts(nn,datax,datay;print=false);

#Print a comparison of initial/final costs
for d=1:length(cf)
    @printf("dataset %3d: initial cost=%10.8f final cost=%10.8f\n",d,c0[d],cf[d]);
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
using Devectorize

if !isdefined(:__EXPRESSION_HASHES__)
    __EXPRESSION_HASHES__ = Set{UInt64}()
end

macro runonce(expr)
    h = hash(expr)
    return esc(quote
        if !in($h, __EXPRESSION_HASHES__)
            push!(__EXPRESSION_HASHES__, $h)
            $expr
        end
    end)
end

"""
  Defines the type `Mat{T}` as a convenience wrapper around the
  `AbstractArray{T,2}` type.
"""
typealias Mat{T} AbstractArray{T, 2}

"""
  Defines the type `Vec{T}` as a convenience wrapper around the
  `AbstractArray{T,1}` type.
"""
typealias Vec{T} AbstractArray{T, 1}

"""
  # Boltzmann.normalize_samples (utils.jl)
  ## Function Calls
    `normalize_samples(X::Mat{Float64})`

  ## Description
    Given a matrix, `X`, assume that each column represents a different
    data sample and that each row represents a different feature. In this
    case, `normalize_samples` will normalize each individual sample to
    the range `[0,1]` according to the minimum and maximum features in the
    sample.

  ## Returns
    1. `::Mat{Float64}`, The normalized dataset.

  ### See also...
    `normalize_samples!`
"""
function normalize_samples(x::Mat{Float64})
    samples = size(x,2)
    y = zeros(x)

    for i=1:samples
      s = x[:,i]
      mins = minimum(s)
      maxs = maximum(s)
      rans = maxs-mins
      y[:,i] = (s-mins)/rans
    end

    return y
end

"""
  # Boltzmann.normalize_samples! (utils.jl)
  ## Function Calls
    `normalize_samples!(x::Mat{Float64})`

  ## Description
    Given a matrix, `x`, assume that each column represents a different
    data sample and that each row represents a different feature. In this
    case, `normalize_samples` will normalize each individual sample to
    the range `[0,1]` according to the minimum and maximum features in the
    sample.

  ## Returns
    Nothing, modifies `x` in place.

  ### See also...
    `normalize_samples`
"""
function normalize_samples!(x::Mat{Float64})
    samples = size(x,2)

    for i=1:samples
      s = x[:,i]
      mins = minimum(s)
      maxs = maximum(s)
      rans = maxs-mins
      x[:,i] = (s-mins)/rans
    end
end

"""
  # Boltzmann.normalize (utils.jl)
  ## Function Calls
    `normalize(x::Mat{Float64})`
    `normalize(x::Vec{Float64})`

  ## Description
    Given an array, `x`, normalize the entire array to the range
    `[0,1]` according to the maximum and miniumum values of the array.

  ## Returns
    1. `::Mat{Float64}` *or* `::Vec{Float64}` depending on the input.

  ### See also...
    `normalize!`
"""
function normalize(x::Mat{Float64})
    y = zeros(x)
    minx = minimum(x)
    maxx = maximum(x)
    ranx = maxx-minx

    # We can index a two-dimension array as if it were a vector.
    @simd for i=1:length(x)
      @inbounds y[i] = (x[i]-minx) / ranx
    end

    return y
end
function normalize(x::Vec{Float64})
    y = zeros(x)
    minx = minimum(x)
    maxx = maximum(x)
    ranx = maxx-minx

    @simd for i=1:length(x)
      @inbounds y[i] = (x[i]-minx) / ranx
    end

    return y
end

"""
  # Boltzmann.normalize! (utils.jl)
  ## Function Calls
    `normalize!(x::Mat{Float64})`
    `normalize!(x::Vec{Float64})`

  ## Description
    Given an array, `x`, normalize the entire array to the range
    `[0,1]` according to the maximum and miniumum values of the array.

  ## Returns
    Nothing. Modifies `x` in place.

  ### See also...
    `normalize`
"""
function normalize!(x::Mat{Float64})
    minx = minimum(x)
    maxx = maximum(x)
    ranx = maxx-minx

    @simd for i=1:length(x)
      @inbounds x[i] = (x[i]-minx) / ranx
    end
end
function normalize!(x::Vec{Float64})
    minx = minimum(x)
    maxx = maximum(x)
    ranx = maxx-minx

    @simd for i=1:length(x)
      @inbounds x[i] = (x[i]-minx) / ranx
    end
end

"""
  # Boltzmann.binarize (utils.jl)
  ## Function Calls
    `binarize(x::Mat{Float64}[,threshold=0.0])`
    `binarize(x::Vec{Float64}[,threshold=0.0])`

  ## Description
    Given an array, `x`, assign each element of the array to
    either `0` or `1` depending on specified value of `threshold`.
    This is done according to the following rule...
    ```
        if element <= threshold: element = 0
        if element >  threshold: element = 1
    ```

  ## Returns
    1. `::Mat{Float64}` *or* `::Vec{Float64}`, depending on input.

  ### See also...
    `binarize!`
"""
function binarize(x::Mat{Float64};threshold=0.0)
  s = copy(x)
  @simd for i=1:length(x)
    @inbounds s[i] = x[i] > threshold ? 1.0 : 0.0
  end
  return s
end
function binarize(x::Vec{Float64};threshold=0.0)
  s = copy(x)
  @simd for i=1:length(x)
    @inbounds s[i] = x[i] > threshold ? 1.0 : 0.0
  end
  return s
end

"""
  # Boltzmann.binarize! (utils.jl)
  ## Function Calls
    `binarize!(x::Mat{Float64}[,threshold=0.0])`
    `binarize!(x::Vec{Float64}[,threshold=0.0])`

  ## Description
    Given an array, `x`, assign each element of the array to
    either `0` or `1` depending on specified value of `threshold`.
    This is done according to the following rule...
    ```
        if element <= threshold: element = 0
        if element >  threshold: element = 1
    ```

  ## Returns
    Nothing. Modifies `x` in place.

  ### See also...
    `binarize`
"""
function binarize!(x::Mat{Float64};threshold=0.0)
  @simd for i=1:length(x)
    @inbounds x[i] = x[i] > threshold ? 1.0 : 0.0
  end
end
function binarize!(x::Vec{Float64};threshold=0.0)
  @simd for i=1:length(x)
    @inbounds x[i] = x[i] > threshold ? 1.0 : 0.0
  end
end

"""
  # Boltzmann.logsig (utils.jl)
  ## Function Calls
    `logsig(x::Mat{Float64})`
    `logsig(x::Vec{Float64})`

  ## Description
    Perform the logistic sigmoid element-by-element for the entries
    of `x`. Here, the logistic sigmoid is defined as
                      `y  =  1 / (1+e^-x)`.

  ## Returns
    1. `::Mat{Float64}` *or* `::Vec{Float64}` (depending on input).

  ### See also...
    `logsig!`
"""
function logsig(x::Mat{Float64})
    @devec s = 1 ./ (1 + exp(-x))
    return s
end
function logsig(x::Vec{Float64})
    @devec s = 1 ./ (1 + exp(-x))
    return s
end

"""
  # Boltzmann.logsig! (utils.jl)
  ## Function Calls
    `logsig!(x::Mat{Float64})`
    `logsig!(x::Vec{Float64})`

  ## Description
    Perform the logistic sigmoid element-by-element for the entries
    of `x`. Here, the logistic sigmoid is defined as
                      `y  =  1 / (1+e^-x)`.

  ## Returns
    Nothing. Modifies `x` in place.

  ### See also...
    `logsig`
"""
function logsig!(x::Mat{Float64})
  # We would like to have used `@devec`, here, but this macro by default
  # makes a new instantiation in memory and thus the result is distinct
  # from the original.
  @simd for i=1:length(x)
    @inbounds x[i] = 1 ./ (1 + exp(-x[i]))
  end
end
function logsig!(x::Vec{Float64})
  @simd for i=1:length(x)
    @inbounds x[i] = 1 ./ (1 + exp(-x[i]))
  end
end

"""
  # Boltzmann.random_samplse (utils.jl)
  ## Function Calls
    `random_columns(x::Mat{Float64})`
  ## Description
    Select a random permutation of the columsn of the specified
    matrix, returning the permuted sub-matrix as well as the
    column indices selected from the original matrix.

  ## Returns
    1. `::Mat{Float64}`, the random sub-matrix.
    2. `::Vec{Int}`, the array of selected column indices
"""
function random_columns(x::Mat{Float64},n_select::Int)
    n_features = size(x,1)
    n_columns = size(x,2)
    perm = shuffle!(collect(1:n_columns))[1:n_select]
    y = x[:,perm]

    return y, perm
end


text:
module CellArrayApplyTests

using Test
using CellwiseValues
using ..CellValuesMocks
using TensorValues

l = 10

a = [1,2,3]
v = TestIterCellValue(a,l)
w = apply(-,v,broadcast=true)
o = [ CachedArray(-vi) for vi in v ]
test_iter_cell_array(w,o)

a = [1,2,3]
b = [3,2,1]
u = TestIterCellValue(a,l)
v = TestIterCellValue(b,l)
w = apply(-,u,v,broadcast=true)
o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
test_iter_cell_array(w,o)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)

ax = [rand(1,3,4), 1.0      , [v1,v2,v3], [v1,v2,v3], 1         ]
bx = [rand(2,3,1), rand(2,3), [v2,v3,v1], v1        , [v2,v3,v1]]

for (a,b) in zip(ax,bx)
  u = TestIterCellValue(a,l)
  v = TestIterCellValue(b,l)
  o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
  w = apply(-,u,v,broadcast=true)
  test_iter_cell_array(w,o)
end

a = [1,2,3]
v = TestIndexCellValue(a,l)
w = apply(-,v,broadcast=true)
o = [ CachedArray(-vi) for vi in v ]
test_index_cell_array(w,o)

a = [1,2,3]
b = [3,2,1]
u = TestIndexCellValue(a,l)
v = TestIndexCellValue(b,l)
w = apply(-,u,v,broadcast=true)
o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
test_index_cell_array(w,o)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)

ax = [rand(1,3,4), 1.0      , [v1,v2,v3], [v1,v2,v3], 1         ]
bx = [rand(2,3,1), rand(2,3), [v2,v3,v1], v1        , [v2,v3,v1]]

for (a,b) in zip(ax,bx)
  u = TestIndexCellValue(a,l)
  v = TestIndexCellValue(b,l)
  o = [ CachedArray(ui.-vi) for (ui,vi) in zip(u,v) ]
  w = apply(-,u,v,broadcast=true)
  test_index_cell_array(w,o)
end

end # module


text:
"""
    coast(cmd0::String=""; kwargs...)

Plot continents, shorelines, rivers, and borders on maps.
Plots grayshaded, colored, or textured land-masses [or water-masses] on
maps and [optionally] draws coastlines, rivers, and political
boundaries. A map projection must be supplied.

Full option list at [`coast`]($(GMTdoc)coast.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- **A** | **area** :: [Type => Str or Number]

    Features with an area smaller than min_area in km^2 or of
    hierarchical level that is lower than min_level or higher than
    max_level will not be plotted.
    ($(GMTdoc)coast.html#a)
- $(GMT.opt_B)
- **C** | **river_fill** :: [Type => Str]

    Set the shade, color, or pattern for lakes and river-lakes.
    ($(GMTdoc)coast.html#c)
- **D** | **res** | **resolution** :: [Type => Str]		``Arg = c|l|i|h|f|a``

    Selects the resolution of the data set to use ((f)ull, (h)igh, (i)ntermediate, (l)ow, (c)rude), or (a)uto).
    ($(GMTdoc)coast.html#d)
- **E** | **DCW** :: [Type => Str]

    Select painting or dumping country polygons from the Digital Chart of the World.
    ($(GMTdoc)coast.html#e)
    + Tuple("code", Str); Tuple(code, number); Tuple("code" [,"fill"], (pen)); Tuple((...),(...),...)
    + ex: ("PT",(0.5,"red","--")); (("PT","gblue",(0.5,"red"),("ES",(0.5,"yellow")))
    +     DCW=:PT; DCW=(:PT, 1); DCW=("PT", :red)
- **F** | **box** :: [Type => Str]

    Draws a rectangular border around the map scale or rose.
    ($(GMTdoc)coast.html#f)
- **G** | **land** :: [Type => Str]

    Select filling or clipping of “dry” areas.
    ($(GMTdoc)coast.html#g)
- **I** | **rivers** :: [Type => Str]

    Draw rivers. Specify the type of rivers and [optionally] append pen attributes.
    ($(GMTdoc)coast.html#i)
- **L** | **map_scale** :: [Type => Str]

    Draw a map scale.
    ($(GMTdoc)coast.html#l)
- **M** | **dump** :: [Type => Str]

    Dumps a single multisegment ASCII output. No plotting occurs.
    ($(GMTdoc)coast.html#m)
- **N** | **borders** :: [Type => Str]

    Draw political boundaries. Specify the type of boundary and [optionally] append pen attributes
    ($(GMTdoc)coast.html#n)
- $(GMT.opt_P)
- **clip** :: [Type => Str]		``Arg = land|water|end``

    To clip land do *clip=:land*, *clip=:water* clips water. Use *end* to mark end of existing clip path.
    No projection information is needed.
    ($(GMTdoc)coast.html#q)
- **S** | **water** | **ocean** :: [Type => Str]

    Select filling or clipping of “wet” areas.
    ($(GMTdoc)coast.html#s)
- **Td** | **rose`** :: [Type => Str]

    Draws a map directional rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)coast.html#t)
- **Tm** | **compass** :: [Type => Str]

    Draws a map magnetic rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)coast.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **shore** | **shorelines** | **coast**  :: [Type => Str]
    Draw shorelines [Default is no shorelines]. Append pen attributes.
    ($(GMTdoc)coast.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function coast(cmd0::String=""; clip=nothing, first=true, kwargs...)

	gmt_proggy = (IamModern[1]) ? "coast "  : "pscoast "
	(length(kwargs) == 0 && clip === nothing) && return monolitic(gmt_proggy, cmd0)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd = parse_E_coast(d, [:E, :DCW], "")		# Process first to avoid warning about "guess"
	cmd = add_opt(d, cmd, "M", [:M :dump])
	if (!occursin("-E+l", cmd) && !occursin("-E+L", cmd) && !occursin("-M", cmd))
		cmd, = parse_BJR(d, cmd, "", O, "guess")
	end
	cmd, = parse_common_opts(d, cmd, [:F :JZ :UVXY :bo :c :p :t :params], first)
	cmd  = parse_these_opts(cmd, d, [[:A :area], [:C :river_fill], [:D :res :resolution]])
	cmd  = parse_Td(d, cmd)
	cmd  = parse_Tm(d, cmd)
	cmd  = parse_L(d, cmd)
	cmd  = add_opt_fill(cmd, d, [:G :land], 'G')
	cmd  = add_opt_fill(cmd, d, [:S :water :ocean], 'S')

	if (clip !== nothing)
		clip = string(clip)
		if     (clip == "land")    cmd *= " -Gc"
		elseif (clip == "water" || clip == "ocean") cmd *= " -Sc"
		elseif (clip == "end")     cmd *= " -Q"
		else
			@warn("The 'clip' argument can only be a string with 'land', 'water' or 'end'. Ignoring it.")
		end
	end

	# Parse these three options that can be made to respond to same code
	cmd = parse_INW_coast(d, [[:I :rivers], [:N :borders], [:W :shore :shorelines :coast]], cmd, "INW")
	(show_kwargs[1]) && print_kwarg_opts([:I :rivers],  "NamedTuple | Tuple | Dict | String")
	(show_kwargs[1]) && print_kwarg_opts([:N :borders], "NamedTuple | Tuple | Dict | String")
	(show_kwargs[1]) && print_kwarg_opts([:W :shore :shorelines :coast],   "NamedTuple | Tuple | Dict | String")

	if (!occursin(" -C",cmd) && !occursin(" -E",cmd) && !occursin(" -G",cmd) && !occursin(" -I",cmd) &&
		!occursin(" -M",cmd) && !occursin(" -N",cmd) && !occursin(" -Q",cmd) && !occursin(" -S",cmd) && !occursin(" -W",cmd))
		cmd *= " -W0.5p"
	end
	(!occursin("-D",cmd)) && (cmd *= " -Da")			# Then pick automatic
	finish = !occursin(" -M",cmd) && !occursin("-E+l", cmd) && !occursin("-E+L", cmd) ? true : false	# Otherwise the dump would be redirected to GMTjl_tmp.ps
	get_largest = (!finish && occursin(" -E", cmd) && (find_in_dict(d, [:biggest :largest])[1] !== nothing))
	if (finish)  _cmd, K = finish_PS_nested(d, [gmt_proggy * cmd], K)
	else	     _cmd = [gmt_proggy * cmd]
	end

	R = finish_PS_module(d, _cmd, "", K, O, finish)
	if (get_largest)
		ind = argmax(size.(R))
		R = [R[ind]]		# Keep it a vector to be consistent with the other Dump cases
		R[1].proj4, R[1].geom = prj4WGS84, wkbPolygon
	end
	R
end

# ---------------------------------------------------------------------------------------------------
function parse_INW_coast(d::Dict, symbs::Array{Array{Symbol,2},1}, cmd::String, flags::String)
	# This function is also used by pshistogram (opt -N). Must be length(flags) == length(symbs)
	(length(symbs) != length(flags)) && error("Length of symbs must be equal to number of chars in FLAGS")
	for k = 1:length(symbs)
		if ((val = find_in_dict(d, symbs[k], false)[1]) !== nothing)
			if (isa(val, NamedTuple) || isa(val, Dict) || (isa(val, Tuple) && isa(val[1], NamedTuple)))
				cmd = add_opt(d, cmd, flags[k], symbs[k], (type="/#", level="/#", mode="+p#", pen=("", add_opt_pen)))
			elseif (isa(val, Tuple))  cmd *= " -" * flags[k] * parse_pen(val)
			else                      cmd *= " -" * flags[k] * arg2str(val)	# Includes Str, Number or Symb
			end
			del_from_dict(d, symbs[k])		# Now we can delete the kwarg
		end
	end
	return cmd
end

# ---------------------------------------------------------------------------------------------------
function parse_E_coast(d::Dict, symbs::Vector{Symbol}, cmd::String)
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "NamedTuple | Tuple | Dict | String")
	if ((val = find_in_dict(d, symbs, false)[1]) !== nothing)
		if (isa(val, String) || isa(val, Symbol))	# Simple case, ex E="PT,+gblue" or E=:PT
			t = string(" -E", val)
			!contains(t, "+") && (t *= "+p0.5")		# If only code(s), append pen
			cmd *= t
		elseif (isa(val, NamedTuple) || isa(val, Dict))
			cmd = add_opt(d, cmd, "E", [:DCW :E], (country="", name="", continent="=",
			                                       pen=("+p", add_opt_pen), fill=("+g", add_opt_fill)))
		elseif (isa(val, Tuple))
			cmd = parse_dcw(cmd, val)
		end
		(GMTver >= v"6.1") && (cmd *= " -Vq")		# Suppress annoying warnings regarding filling syntax with +r<dpi>
		del_from_dict(d, symbs)
	end
	return cmd
end

# ---------------------------------------------------------------------------------------------------
function parse_dcw(cmd::String, val::Tuple)::String
	# First test if (code, pen) or (code, fill) are passed in.
	if (isa(val, Tuple) && length(val) == 2)
		(isa(val[2], Real)) && return string(cmd, " -E", val[1], "+p", val[2])
		(isa(val[2], String) || isa(val[2], Symbol)) && return string(cmd, " -E", val[1], "+g", string(val[2]))
	end

	for k = 1:length(val)
		if (isa(val[k], NamedTuple) || isa(val[k], Dict))
			if (isa(val[k], Dict))  val[k] = dict2nt(val[k])  end
			cmd *= add_opt(Dict(:DCW => val[k]), "", "E", [:DCW],
			               (country="", name="", continent="=", pen=("+p", add_opt_pen), fill=("+g", add_opt_fill)))
		elseif (isa(val[k], Tuple))
			cmd *= parse_dcw(val[k])
		else
			cmd *= parse_dcw(val)
			break
		end
	end
	return cmd
end

function parse_dcw(val::Tuple)
    t = string("", " -E", val[1])
	if (length(val) > 1)
		if (isa(val[2], Tuple))  t *= "+p" * parse_pen(val[2])
		else                     t *= string(val[2])
		end
		if (length(val) > 2)
			if (isa(val[3], Tuple))  t *= add_opt_fill("+g", Dict(fill => val[3]), [:fill])
			else                     t *= string(val[3])
			end
		end
	end
	return t
end

# ---------------------------------------------------------------------------------------------------
coast!(cmd0::String=""; clip=nothing, kw...) = coast(cmd0; clip=clip, first=false, kw...)

const pscoast  = coast			# Alias
const pscoast! = coast!


text:
using Printf

function main()
    output_dir_name="post_result"
    input_result_dir="result"

    #vtk_頭の文字
    front="# vtk DataFile Version 2.0"*"\n"
    back="\n"*"ASCII"*"\n"*"DATASET UNSTRUCTURED_GRID"*"\n"
    write_file(output_dir_name,input_result_dir,front,back)
end

function write_file(outdir,inresult_dir,front,back)

    nodes_xmax, nodes_ymax = read_nodenum(1)
    nodes = read_nodes(1)

    elements = read_elements(1)

    inf = readdir(inresult_dir)
    cellnum = (nodes_xmax-3)*(nodes_ymax-3)
    rho = zeros(cellnum)
    u = zeros(cellnum)
    v = zeros(cellnum)
    p = zeros(cellnum)
    T = zeros(cellnum)

    for i in 1:length(inf)
        if occursin(".dat", inf[i]) == true
            dname = replace(inf[i],".dat" => "")
            out_file = replace(inf[i],".dat" => ".vtk")
            print("start writing "*out_file*"\n")

            fff=[]
            open(inresult_dir*"/"*inf[i], "r") do f
                fff=read(f,String)
            end
            fff=split(fff,"\n",keepempty=false)

            for j in 2:length(fff)
                fff[j]=replace(fff[j]," \r" => "")
                temp = split(fff[j]," ")
                rho[j-1] = parse(Float64,temp[1])
                u[j-1] = parse(Float64,temp[2])
                v[j-1] = parse(Float64,temp[3])
                p[j-1] = parse(Float64,temp[4])
                T[j-1] = parse(Float64,temp[5])
            end

            write_points(nodes,out_file,outdir,dname,front,back)
            write_cells(elements,out_file,outdir)
            write_result(rho,out_file,outdir,1)
            write_result(u,out_file,outdir,2)
            write_result(v,out_file,outdir,3)
            write_result(p,out_file,outdir,4)
            write_result(T,out_file,outdir,5)
            write_result_vec(u,v,out_file,outdir)
            print("fin writing "*out_file*"\n")
        end
    end
end

function read_nodenum(skipnum)
    # xmax : 仮想セルも含めたnodeのxの数
    # ymax : 仮想セルも含めたnodeのyの数
    fff=[]
    open("grid/nodesnum", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    temp = split(fff[2]," ")
    xmax = parse(Int64,temp[1])
    ymax = parse(Int64,temp[2])

    return xmax, ymax
end

function read_nodes(skipnum)

    fff=[]
    open("grid/nodes_forvtk", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    nodes=zeros(num_nodes,3)
    for i in 1:num_nodes
        temp=split(fff[i+skipnum]," ")

        # x = parse(Float64,temp[1])
        # y = parse(Float64,temp[2])
        # z = parse(Float64,temp[3])
        x = parse(Float64,temp[2])
        y = parse(Float64,temp[3])
        z = 0.0
        nodes[i,1] = x
        nodes[i,2] = y
        nodes[i,3] = z
    end
    return nodes
end

function read_elements(skipnum)
    fff=[]
    open("grid/element_forvtk", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_elements=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    elements = zeros(num_elements,4)
    for i in 1:num_elements
        temp=split(fff[i+skipnum]," ")

        elements[i,1] = parse(Int64,temp[2])-1
        elements[i,2] = parse(Int64,temp[3])-1
        elements[i,3] = parse(Int64,temp[4])-1
        elements[i,4] = parse(Int64,temp[5])-1
    end
    return elements
end

function  write_points(nodes,out_file,outdir,dname,front,back)
    a = size(nodes)[1]
    a_st = @sprintf("%1.0f", a)

    fff = outdir*"/"*out_file
    open(fff,"w") do f
        write(f,front)
        write(f,dname)
        write(f,back)
        write(f,"POINTS "*a_st*" float\n")
        for i in 1:a
            x = @sprintf("%7.7f", nodes[i,1])
            y = @sprintf("%7.7f", nodes[i,2])
            z = @sprintf("%7.7f", nodes[i,3])
            write(f, x*" "*y*" "*z*"\n")
        end
    end
    println("fin writing points")
end

function write_cells(elements,out_file,outdir)
    a = size(elements)[1]
    b = a*5
    a_st = @sprintf("%1.0f", a)
    b_st = @sprintf("%1.0f", b)

    fff = outdir*"/"*out_file
    open(fff,"a") do f
        write(f, "CELLS "*a_st*" "*b_st*"\n")
        for i in 1:a
            d1 = @sprintf("%1.0f", elements[i,1])
            d2 = @sprintf("%1.0f", elements[i,2])
            d3 = @sprintf("%1.0f", elements[i,3])
            d4 = @sprintf("%1.0f", elements[i,4])
            write(f, "4 "*d1*" "*d2*" "*d3*" "*d4*"\n")
        end
        write(f, "CELL_TYPES "*a_st*"\n")
        for i in 1:a
            write(f, "9\n")          #四角のみ
        end
    end
    println("fin writing cells")
end

function write_result(val,out_file,outdir,k)
    dtype = ["rho","u","v","p","T"]

    a = length(val)
    a_st = @sprintf("%1.0f", a)

    temp1 = "CELL_DATA "*a_st*"\n"
    temp2 = "SCALARS "*dtype[k]*" float\nLOOKUP_TABLE default\n"

    fff = outdir*"/"*out_file
    open(fff,"a") do f
        if k == 1
            write(f, temp1)
        end
        write(f, temp2)
        for i in 1:a
            data = @sprintf("%7.7f", val[i])
            write(f, data)
            write(f, "\n")
        end
    end
end

function write_result_vec(u,v,out_file,outdir)
    a = length(u)
    a_st = @sprintf("%1.0f", a)

    temp1 = "CELL_DATA "*a_st*"\n"
    temp2 = "VECTORS velocity_vec float\n"

    fff = outdir*"/"*out_file
    open(fff,"a") do f
        #write(f, temp1)
        write(f, temp2)
        for i in 1:a
            data = @sprintf("%7.7f", u[i])*" "*@sprintf("%7.7f", v[i])*" 0.0"
            write(f, data)
            write(f, "\n")
        end
    end
end

# main
main()


text:
using SpaceShuttle
using Test

@testset "SpaceShuttle.jl" begin
    @test launch() isa String
end


text:
using Fetch
using Test
using DataDeps

@testset "Fetch.jl" begin
    include("gdrive.jl")
    include("kaggle.jl")
end


text:
#EasyPlotPlots base types & core functions
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#TODO: deprecate:
const scalemap = Dict{Symbol, Symbol}(
	:lin => :identity,
	:log => :log10,
)
#:identity, :ln, :log2, :log10, :asinh, :sqrt

#="linetype" values
:none, :line, :sticks, :bar, :hist,
:scatter, :scatter3d,
:path, :path3d,
:contour, :density, :heatmap, :hexbin,
:surface, :wireframe
:steppost, :steppre,
:hline, :vline
=#

const linestylemap = Dict{Symbol, Symbol}(
	:none    => :auto, #Will modify linetype to cover "none".
	:solid   => :solid,
	:dash    => :dash,
	:dot     => :dot,
	:dashdot => :dashdot,
)
#:auto, :dash, :dashdot, :dot, :solid

const markermap = Dict{Symbol, Symbol}(
	:none    => :none,
	:square    => :rect,
	:diamond   => :diamond,
	:uarrow    => :utriangle, :darrow => :dtriangle,
	:larrow    => :utriangle, :darrow => :dtriangle, #NOT SUPPORTED
	:cross     => :cross, :+ => :cross,
	:diagcross => :xcross, :x => :xcross,
	:circle    => :ellipse, :o => :ellipse,
	:star      => :star5, :* => :star5,
)

#=
:auto, :none,
:cross, :xcross
:ellipse, :rect, :diamond, :dtriangle, :utriangle,
:heptagon, :hexagon, :octagon, :pentagon,
:star4, :star5, :star6, :star7, :star8
=#

#const NOMULTISUPPORT = Set([:gr]) #GR now supported using multiplot model.
const NOMULTISUPPORT = Set([])

#==Base types
===============================================================================#
mutable struct WfrmBuilder <: EasyPlot.AbstractWfrmBuilder
	subplot::Plots.Subplot #Axes reference
	theme::EasyPlot.Theme
	fold::Optional{EasyPlot.FoldedAxis}
end

#Top-level plot (which might contain subplots)
abstract type Figure end

#TODO: Is this workaround still needed??
mutable struct FigureSng <: Figure #For backends that support single plots only (no subplots)
	subplots::Vector{Plots.Plot}
end
FigureSng() = FigureSng(Plots.Plot[])

mutable struct FigureMulti <: Figure #For backends that support subplots
	p::Optional{Plots.Plot} #Optional - So construction does not display anything.
#	subplots::Vector{Plots.Subplot}
end
FigureMulti() = FigureMulti(Plots.plot(overwrite_figure=false))

Figure(backend::Symbol) =
	in(backend, NOMULTISUPPORT) ? FigureSng() : FigureMulti()

#Immutable? would be on stack...
mutable struct WfrmAttributes
	label::Optional{String}

	linetype::Symbol
	linestyle::Symbol
	linewidth::Float64
	linecolor::Colorant
#	linealpha

	markershape::Symbol
	markersize::Float64
	markercolor::Colorant #Fill color
	markeralpha::Float64 #Apperently entire marker
	markerstrokewidth::Float64
	markerstrokecolor::Colorant
	markerstrokealpha::Float64
	markerstrokestyle::Symbol
end

#==Constructor-like functions
===============================================================================#
function addsubplots(fig::FigureSng, ncols::Int, nsubplots::Int)
	for i in 1:nsubplots
		push!(fig.subplots, Plots.plot())
	end
	return fig
end
function addsubplots(fig::FigureMulti, ncols::Int, nsubplots::Int)
	nrows = div(nsubplots-1, ncols) + 1
	fig.p = Plots.plot(layout=(nrows, ncols), overwrite_figure=true)
	return fig
end

#==Helper functions
===============================================================================#
getsubplot(fig::FigureSng, idx::Int) = fig.subplots[idx]
getsubplot(fig::FigureMulti, idx::Int) = fig.p.subplots[idx]

#Linewidth:
maplinewidth(w) = w
maplinewidth(::Nothing) = maplinewidth(1) #default

#Marker size:
mapmarkersize(sz) = 5*sz
mapmarkersize(::Nothing) = mapmarkersize(1)

function maplinetype(v::Symbol) #Only support small subset of linetype values
	if :none == v
		return :none
	else
		return :line
	end
end
maplinetype(::Nothing) = :line #default

function maplinestyle(v::Symbol)
	result = get(linestylemap, v, missing)
	if ismissing(result)
		@info("Line style not supported")
		result = maplinestyle(nothing)
	end
	return result
end
maplinestyle(::Nothing) = :solid #default

function mapmarkershape(v::Symbol)
	result = get(markermap, v, missing)
	if ismissing(result)
		@info("Marker shape not supported")
		result = "o" #Use some supported marker
	end
	return result
end
mapmarkershape(::Nothing) = "" #default (no marker)

function WfrmAttributes(id::String, attr::EasyPlot.WfrmAttributes)
	linewidth = maplinewidth(attr.linewidth)

	#Marker line
	markerstrokecolor = attr.glyphlinecolor
	markerstrokealpha = Colors.alpha(markerstrokecolor)

	#Marker fill:
	markercolor = attr.glyphfillcolor
	markeralpha = Colors.alpha(markercolor)

	#WARNING: In GR backend, markeralpha is applied to marker line!

	return WfrmAttributes(id,
		maplinetype(attr.linestyle),
		maplinestyle(attr.linestyle),
		linewidth,
		attr.linecolor,

		mapmarkershape(attr.glyphshape),
		mapmarkersize(attr.glyphsize),
		markercolor, markeralpha,
		linewidth,
		markerstrokecolor, markerstrokealpha,
		:solid, #markerstrokestyle
	)
end

#==AbstractWfrmBuilder implementation
===============================================================================#
EasyPlot.needsfold(b::WfrmBuilder) = b.fold

#Add DataF1 results:
function _addwfrm(sp::Plots.Subplot, d::DataF1, a::WfrmAttributes)
	kwargs = Any[]
	for attrib in fieldnames(WfrmAttributes)
		v = getfield(a,attrib)

		if v != nothing
			push!(kwargs, tuple(attrib, v))
		end
	end

	wfrm = plot!(sp, d.x, d.y; kwargs...)
end

#Called by EasyPlot, for each individual DataF1 ∈ DataMD.
function EasyPlot.addwfrm(b::WfrmBuilder, d::DataF1, id::String,
	la::EasyPlot.LineAttributes, ga::EasyPlot.GlyphAttributes, strip::Int)
	attr = EasyPlot.WfrmAttributes(b.theme, la, ga) #Apply theme to attributes
	plotsattr = WfrmAttributes(id, attr) #Attributes understood by Plots.jl
	_addwfrm(b.subplot, d, plotsattr)
end

#Last line


text:
using Random

board_size=3

@enum GameState begin
    IN_PROGRESS = 0
    DRAW = 1
    X_WINS = 2
    O_WINS = 3
end

struct Move
    x::Int
    y::Int

    function Move(x::Int, y::Int)
        new(x, y)
    end
end

function transform(value)
    return printable[value]
end

function make_move(board,move,player)
    i = move.x
    j = move.y
    board[i][j] = player
    return board
end

function legal_moves(board)
    moves = []
    for i in range(1,stop=board_size)
        for j in range(1,stop=board_size)
            if board[i][j] == 0
                move = Move(i,j)
                push!(moves, move)
            end
        end
    end
    return shuffle(moves)
end

function player_wins(player)
    if player == 1
        result = X_WINS
    else
        result = O_WINS
    end
    return result
end

function get_game_state(board,player)
    result = IN_PROGRESS

    ## horizonal win
    if (result == IN_PROGRESS)
        for row in board
            if sum(row) == board_size*player
                result = player_wins(player)
            end
        end
    end

    ## vertical win
    if (result == IN_PROGRESS)
        for j in range(1,stop=board_size)
            column = [board[i][j] for i in range(1,stop=board_size)]
            if sum(column) == board_size
                result = player_wins(player)
            end
        end
    end

    ## forward diagonal wins
    if (result == IN_PROGRESS)
        diag = [board[i][i] for i in range(1,stop=board_size)]
        if sum(diag) == board_size
            result = player_wins(player)
        end
    end

    ## backward diagonal wins
    if (result == IN_PROGRESS)
        diag = [board[i][i] for i in board_size:-1:1]
        if sum(diag) == board_size
            result = player_wins(player)
        end
    end

    if (result == IN_PROGRESS)
        ## no moves left
        if length(legal_moves(board)) == 0
            result = DRAW
        end
    end

    return result
end

function print_board(board)
    println()
    for row in board
        println([transform(element) for element in row])
    end
    println()
end

## transform numbers to letters for printed board
printable = Dict{Int64,String}(-1=>"O", 0=>" ",1=>"X")

board=[[0 for j in range(1,stop=board_size)] for i in range(1,stop=board_size)]
##push!(board,row)
##board = [i for i in range(board_size)]
print_board(board)
player = 1

game_state = get_game_state(board,player)
while game_state == IN_PROGRESS
    global game_state
    global board
    global player
    moves = legal_moves(board)
    board = make_move(board,moves[1],player)
    print_board(board)
    game_state = get_game_state(board,player)
    player = -player
end
println(game_state)

max_stars_repo_name: mbreault/julia

max_issues_repo_name: mbreault/julia

max_forks_repo_name: mbreault/julia


text:
module CommonRLInterface

using MacroTools

export
    AbstractEnv,
    reset!,
    actions,
    observe,
    act!,
    terminated

abstract type AbstractEnv end

"""
    reset!(env::AbstractEnv)

Reset `env` to its initial state and return `nothing`.

This is a *required function* that must be provided by every AbstractEnv.
"""
function reset! end

"""
    actions(env::AbstractEnv)

Return a collection of all of the actions available for AbstractEnv `env`. If the environment has multiple agents, this function should return the union of all of their actions.

This is a *required function* that must be provided by every AbstractEnv.

This function is a *static property* of the environment; the value it returns should not change based on the state.

---

    actions(env::AbstractEnv, player_index)

Return a collection of all the actions available to a given player.

This function is a *static property* of the environment; the value it returns should not change based on the state.
"""
function actions end

"""
    observe(env::AbstractEnv)

Return an observation from the environment for the current player.

This is a *required function* that must be provided by every AbstractEnv.
"""
function observe end

"""
    r = act!(env::AbstractEnv, a)

Take action `a` and advance AbstractEnv `env` forward one step, and return rewards for all players.

This is a *required function* that must be provided by every AbstractEnv.

If the environment has a single player, it is acceptable to return a scalar number. If there are multiple players, it should return a container with all rewards indexed by player number.

# Example

## Single Player
```julia
function act!(env::MyMDPEnv, a)
    env.state += a + randn()
    return env.s^2
end
```

## Two Player

```julia
function act!(env::MyMDPEnv, a)
    env.positions[player(env)] += a   # In this game, each player has a position that is updated by his or her action
    rewards = in_goal.(env.positions) # Rewards are +1 for being in a goal region, 0 otherwise
    return rewards                    # returns a vector of rewards for each player
end
```
"""
function act! end

"""
    terminated(env::AbstractEnv)

Determine whether an environment has finished executing.

If `terminated(env)` is true, no further actions should be taken and it is safe to assume that no further rewards will be received.
"""
function terminated end

export
    provided,
    @provide

"""
    provided(f, args...)

Test whether an implementation for `f(args...)` has been provided.

If this returns false, you should assume that the environment does not support the function `f`.

Usage is identical to `Base.applicable`, but unlike `Base.applicable`, `provided` can usually be inferred at compile time. This function is usually automatically implemented by the @provide macro. It should only be implemented by the user in exceptional cases where very fine control is needed.

---

    provided(f, types::Type{<:Tuple})

Alternate version of provided with syntax similar to `Base.hasmethod`.

If this returns true, it means that the function is provided for any set of arguments with the given types. If it returns false, it may be provided for certain objects. For this reason, algorithm writers should call the `provided(f, args...)` version when possible.
"""
function provided end

provided(f::Function, args...) = provided(f, typeof(args))
provided(f::Function, ::Type{<:Tuple}) = false

provided(::typeof(reset!), ::Type{<:Tuple{AbstractEnv}}) = true
provided(::typeof(actions), ::Type{<:Tuple{AbstractEnv}}) = true
provided(::typeof(observe), ::Type{<:Tuple{AbstractEnv}}) = true
provided(::typeof(act!), ::Type{<:Tuple{AbstractEnv, Any}}) = true

"""
    @provide f(x::X) = x^2

Indicate that function `f` has been implemented for arguments of type `X`.

This will automatically implement the appropriate methods of `provided`. Both the long and short function definition forms will work.

# Example
```jldoctest
using CommonRLInterface

struct MyEnv <: AbstractEnv
    s::Int
end

# Initially, clone is not provided
@assert !provided(clone, MyEnv(1))

@provide CommonRLInterface.clone(env::MyEnv) = MyEnv(env.s)

@assert provided(clone, MyEnv(1))

clone(MyEnv(1))

# output

MyEnv(1)

```
"""
macro provide(f)
    def = splitdef(f) # TODO: probably give a better error message that mentions @provide if this fails
    @assert isempty(def[:kwargs]) "@provide does not support keyword args yet."

    func = esc(def[:name])

    argtypes = []
    for arg in def[:args]
        if @capture(arg, name_::T_)
            push!(argtypes, esc(T))
        else
            push!(argtypes, :Any)
        end
    end

    quote
        CommonRLInterface.provided(::typeof($func), ::Type{<:Tuple{$(argtypes...)}}) where {$(map(esc, def[:whereparams])...)} = true

        $(esc(f))
    end
end

export
    clone,
    render,
    state,
    setstate!
include("environment.jl")

export
    observations,
    valid_actions,
    valid_action_mask
include("spaces.jl")

export
    players,
    player,
    all_act!,
    all_observe,
    UtilityStyle,
    ZeroSum,
    ConstantSum,
    GeneralSum,
    IdenticalUtility
include("multiplayer.jl")

export
    Wrappers
include("wrappers.jl")

end


text:
using PkgMPC
using Test

F(3,5)

##For Testing
Data_MPC   = Simulation_Data(60.0, 8.0, 3, [1.0; 1.0], 0.35)
Data_Plant = Model_Par(4.0, 0.12, 0.4545, 0.4, 0.4)
Results = Simulate_Plant(Data_MPC, Data_Plant)
Figure = Plot_Plant(Results, Data_MPC)
Save_Plant(Figure,"//home//amir//Documents//Package_Results//MPC_Results.svg")
Export_Plant(Results,Data_MPC, "//home//amir//Documents//Package_Results//MPC_Results.csv")
##for timing
Data_MPC   = @elapsed Simulation_Data(60.0, 8.0, 3, [1.0; 1.0], 0.35)
Data_Plant = @elapsed Model_Par(4.0, 0.12, 0.4545, 0.4, 0.4)
result = @elapsed Simulate_Plant(Data_MPC, Data_Plant)
Figure = @elapsed Plot_Plant(Results, Data_MPC)
time_save = @elapsed Save_Plant(Figure,"//home//amir//Documents//Package_Results//MPC_Results.svg")
time_export = @elapsed Export_Plant(Results,Data_MPC, "//home//amir//Documents//Package_Results//MPC_Results.csv")

@testset "PkgMPC.jl" begin
    # Write your tests here.
    @test F(2,3) == 13
    @test F(3,0) == 6
end


text:
#=
Computes quantities related to the Gaussian linear state space model

    x_{t+1} = A x_t + C w_{t+1}

        y_t = G x_t

The shocks {w_t} are iid and N(0, I)

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-07-28

References
----------

TODO: Come back and update to match `LinearStateSpace` type from py side
TODO: Add docstrings

http://quant-econ.net/jl/linear_models.html

=#
import Distributions: MultivariateNormal, rand
import Base: ==

#=
    numpy allows its multivariate_normal function to have a matrix of
    zeros for the covariance matrix; Stats.jl doesn't. This type just
    gives a `rand` method when we pass in a matrix of zeros for Sigma_0
    so the rest of the api can work, unaffected

    The behavior of `rand` is to just pass back the mean vector when
    the covariance matrix is zero.
=#
type FakeMVTNorm{T <: Real}
    mu_0::Array{T}
    Sigma_0::Array{T}
end

==(f1::FakeMVTNorm, f2::FakeMVTNorm) =
    (f1.mu_0 == f2.mu_0) && (f1.Sigma_0 == f2.Sigma_0)

Base.rand{T}(d::FakeMVTNorm{T}) = copy(d.mu_0)
"""
A type that describes the Gaussian Linear State Space Model
of the form:

    x_{t+1} = A x_t + C w_{t+1}

        y_t = G x_t

where {w_t} and {v_t} are independent and standard normal with dimensions
k and l respectively.  The initial conditions are mu_0 and Sigma_0 for x_0
~ N(mu_0, Sigma_0).  When Sigma_0=0, the draw of x_0 is exactly mu_0.

#### Fields

- `A::Matrix` Part of the state transition equation.  It should be `n x n`
- `C::Matrix` Part of the state transition equation.  It should be `n x m`
- `G::Matrix` Part of the observation equation.  It should be `k x n`
- `k::Int` Dimension
- `n::Int` Dimension
- `m::Int` Dimension
- `mu_0::Vector` This is the mean of initial draw and is of length `n`
- `Sigma_0::Matrix` This is the variance of the initial draw and is `n x n` and
                    also should be positive definite and symmetric

"""
type LSS
    A::Matrix
    C::Matrix
    G::Matrix
    k::Int
    n::Int
    m::Int
    mu_0::Vector
    Sigma_0::Matrix
    dist::Union{MultivariateNormal, FakeMVTNorm}
end

function LSS(A::ScalarOrArray, C::ScalarOrArray, G::ScalarOrArray,
             mu_0::ScalarOrArray,
             Sigma_0::Matrix=zeros(size(G, 2), size(G, 2)))
    k = size(G, 1)
    n = size(G, 2)
    m = size(C, 2)

    # coerce shapes
    A = reshape(vcat(A), n, n)
    C = reshape(vcat(C), n, m)
    G = reshape(vcat(G), k, n)

    mu_0 = reshape([mu_0;], n)

    # define distribution
    if all(Sigma_0 .== 0.0)   # no variance -- no distribution
        dist = FakeMVTNorm(mu_0, Sigma_0)
    else
        dist = MultivariateNormal(mu_0, Sigma_0)
    end
    LSS(A, C, G, k, n, m, mu_0, Sigma_0, dist)
end

# make kwarg version
function LSS(A::ScalarOrArray, C::ScalarOrArray, G::ScalarOrArray;
             mu_0::Vector=zeros(size(G, 2)),
             Sigma_0::Matrix=zeros(size(G, 2), size(G, 2)))
    return LSS(A, C, G, mu_0, Sigma_0)
end

function simulate(lss::LSS, ts_length=100)
    x = Array(Float64, lss.n, ts_length)
    x[:, 1] = rand(lss.dist)
    w = randn(lss.m, ts_length - 1)
    for t=1:ts_length-1
        x[:, t+1] = lss.A * x[:, t] .+ lss.C * w[:, t]
    end
    y = lss.G * x

    return x, y
end

"""
Simulate num_reps observations of x_T and y_T given x_0 ~ N(mu_0, Sigma_0).

#### Arguments

- `lss::LSS` An instance of the Gaussian linear state space model.
- `t::Int = 10` The period that we want to replicate values for.
- `num_reps::Int = 100` The number of replications we want

#### Returns

- `x::Matrix` An n x num_reps matrix, where the j-th column is the j_th
              observation of x_T
- `y::Matrix` An k x num_reps matrix, where the j-th column is the j_th
              observation of y_T

"""
function replicate(lss::LSS, t::Integer, num_reps::Integer=100)
    x = Array(Float64, lss.n, num_reps)
    for j=1:num_reps
        x_t, _ = simulate(lss, t+1)
        x[:, j] = x_t[:, end]
    end

    y = lss.G * x
    return x, y
end

replicate(lss::LSS; t::Integer=10, num_reps::Integer=100) =
    replicate(lss, t, num_reps)

"""
Create a generator to calculate the population mean and
variance-convariance matrix for both x_t and y_t, starting at
the initial condition (self.mu_0, self.Sigma_0).  Each iteration
produces a 4-tuple of items (mu_x, mu_y, Sigma_x, Sigma_y) for
the next period.

#### Arguments

- `lss::LSS` An instance of the Gaussian linear state space model

"""
function moment_sequence(lss::LSS)
    A, C, G = lss.A, lss.C, lss.G
    mu_x, Sigma_x = copy(lss.mu_0), copy(lss.Sigma_0)
    while true
        mu_y, Sigma_y = G * mu_x, G * Sigma_x * G'
        produce((mu_x, mu_y, Sigma_x, Sigma_y))

        # Update moments of x
        mu_x = A * mu_x
        Sigma_x = A * Sigma_x * A' + C * C'
    end
end

"""
Compute the moments of the stationary distributions of x_t and
y_t if possible.  Computation is by iteration, starting from the
initial conditions lss.mu_0 and lss.Sigma_0

#### Arguments

- `lss::LSS` An instance of the Guassian linear state space model
- `;max_iter::Int = 200` The maximum number of iterations allowed
- `;tol::Float64 = 1e-5` The tolerance level one wishes to achieve

#### Returns

- `mu_x::Vector` Represents the stationary mean of x_t
- `mu_y::Vector`Represents the stationary mean of y_t
- `Sigma_x::Matrix` Represents the var-cov matrix
- `Sigma_y::Matrix` Represents the var-cov matrix

"""
function stationary_distributions(lss::LSS; max_iter=200, tol=1e-5)
    # Initialize iteration
    m = @task moment_sequence(lss)
    mu_x, mu_y, Sigma_x, Sigma_y = consume(m)

    i = 0
    err = tol + 1.

    while err > tol
        if i > max_iter
            error("Convergence failed after $i iterations")
        else
            i += 1
            mu_x1, mu_y, Sigma_x1, Sigma_y = consume(m)
            err_mu = Base.maxabs(mu_x1 - mu_x)
            err_Sigma = Base.maxabs(Sigma_x1 - Sigma_x)
            err = max(err_Sigma, err_mu)
            mu_x, Sigma_x = mu_x1, Sigma_x1
        end
    end

    return mu_x, mu_y, Sigma_x, Sigma_y
end

function geometric_sums(lss::LSS, bet, x_t)
    I = eye(lss.n)
    S_x = (I - bet .* lss.A) \ x_t
    S_y = lss.G * S_x
    return S_x, S_y
end


text:
# ---
# title: 851. Loud and Rich
# id: problem851
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/loud-and-rich/description/>
# hidden: true
# ---
#
# In a group of N people (labelled `0, 1, 2, ..., N-1`), each person has
# different amounts of money, and different levels of quietness.
#
# For convenience, we'll call the person with label `x`, simply "person `x`".
#
# We'll say that `richer[i] = [x, y]` if person `x` definitely has more money
# than person `y`.  Note that `richer` may only be a subset of valid
# observations.
#
# Also, we'll say `quiet[x] = q` if person x has quietness `q`.
#
# Now, return `answer`, where `answer[x] = y` if `y` is the least quiet person
# (that is, the person `y` with the smallest value of `quiet[y]`), among all
# people who definitely have equal to or more money than person `x`.
#
#
#
# **Example 1:**
#
#
#
#     Input: richer = [[1,0],[2,1],[3,1],[3,7],[4,3],[5,3],[6,3]], quiet = [3,2,5,4,6,1,7,0]
#     Output: [5,5,2,5,4,5,6,7]
#     Explanation:
#     answer[0] = 5.
#     Person 5 has more money than 3, which has more money than 1, which has more money than 0.
#     The only person who is quieter (has lower quiet[x]) is person 7, but
#     it isn't clear if they have more money than person 0.
#
#     answer[7] = 7.
#     Among all people that definitely have equal to or more money than person 7
#     (which could be persons 3, 4, 5, 6, or 7), the person who is the quietest (has lower quiet[x])
#     is person 7.
#
#     The other answers can be filled out with similar reasoning.
#
#
# **Note:**
#
#   1. `1 <= quiet.length = N <= 500`
#   2. `0 <= quiet[i] < N`, all `quiet[i]` are different.
#   3. `0 <= richer.length <= N * (N-1) / 2`
#   4. `0 <= richer[i][j] < N`
#   5. `richer[i][0] != richer[i][1]`
#   6. `richer[i]`'s are all different.
#   7. The observations in `richer` are all logically consistent.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## stationary GP options
xall = transD_GP.get_image_prediction_points(img)
nmin, nmax= 2, 100
pnorm = 2.
K = transD_GP.GP.Mat32()
ftrain = img.d[img.select]
fbounds = permutedims([extrema(ftrain)...])
δ = 0.3
λ = [141,141]
demean = true
sdev_prop = 0.05*diff(fbounds, dims=2)[:]
sdev_pos = 0.05*abs.([diff([extrema(xall[1,:])...])[1], diff([extrema(xall[2,:])...])[1]])
updatenonstat = false
needλ²fromlog = false
## Initialize a stationary GP using these options
Random.seed!(12)
opt = transD_GP.OptionsStat(nmin = nmin,
                        nmax = nmax,
                        xbounds = [img.x[1] img.x[end];img.y[1] img.y[end]],
                        fbounds = fbounds,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        demean = demean,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        quasimultid = false,
                        K = K,
                        timesλ = 3,
                        needλ²fromlog = needλ²fromlog,
                        updatenonstat = updatenonstat
                        )

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA