text:
const P4EST_2_GRIDAP_FACET_2D  = [ 3, 4, 1, 2 ]
const GRIDAP_2_P4EST_FACET_2D  = [ 3, 4, 1, 2 ]

const P4EST_2_GRIDAP_FACET_3D  = [ 5, 6, 3, 4, 1, 2 ]
const GRIDAP_2_P4EST_FACET_3D  = [ 5, 6, 3, 4, 1, 2 ]

P4est_wrapper.quadrant_data(x::Clong) = reinterpret(P4est_wrapper.quadrant_data, x)

function p4est_get_quadrant_vertex_coordinates(connectivity::Ptr{p4est_connectivity_t},
                                               treeid::p4est_topidx_t,
                                               x::p4est_qcoord_t,
                                               y::p4est_qcoord_t,
                                               level::Int8,
                                               corner::Cint,
                                               vxy::Ptr{Cdouble})

    myself=Ref{p4est_quadrant_t}(
      p4est_quadrant_t(x,y,level,Int8(0),Int16(0),
                       P4est_wrapper.quadrant_data(Clong(0))))
    neighbour=Ref{p4est_quadrant_t}(myself[])
    if corner == 1
       p4est_quadrant_face_neighbor(myself,corner,neighbour)
    elseif corner == 2
       p4est_quadrant_face_neighbor(myself,corner+1,neighbour)
    elseif corner == 3
       p4est_quadrant_corner_neighbor(myself,corner,neighbour)
    end
    # Extract numerical coordinates of lower_left
    # corner of my corner neighbour
    p4est_qcoord_to_vertex(connectivity,
                           treeid,
                           neighbour[].x,
                           neighbour[].y,
                           vxy)
end

function  p8est_get_quadrant_vertex_coordinates(connectivity::Ptr{p8est_connectivity_t},
                                                treeid::p4est_topidx_t,
                                                x::p4est_qcoord_t,
                                                y::p4est_qcoord_t,
                                                z::p4est_qcoord_t,
                                                level::Int8,
                                                corner::Cint,
                                                vxyz::Ptr{Cdouble})

  myself=Ref{p8est_quadrant_t}(
       p8est_quadrant_t(x,y,z,level,Int8(0),Int16(0),
                        P4est_wrapper.quadrant_data(Clong(0))))
  neighbour=Ref{p8est_quadrant_t}(myself[])

  if ( corner == 1 )
    p8est_quadrant_face_neighbor(myself,Cint(1),neighbour)
  elseif ( corner == 2 )
    p8est_quadrant_face_neighbor(myself,Cint(3),neighbour)
  elseif ( corner == 3 )
    p8est_quadrant_edge_neighbor(myself,Cint(11),neighbour)
  elseif ( corner == 4 )
    p8est_quadrant_face_neighbor(myself,Cint(5),neighbour)
  elseif ( corner == 5 )
    p8est_quadrant_edge_neighbor(myself,Cint(7),neighbour)
  elseif ( corner == 6 )
    p8est_quadrant_edge_neighbor(myself,Cint(3),neighbour)
  elseif ( corner == 7 )
    p8est_quadrant_corner_neighbor(myself,Cint(7),neighbour)
  end

  # Extract numerical coordinates of lower_left corner of my corner neighbour
  p8est_qcoord_to_vertex(connectivity,
                         treeid,
                         neighbour[].x,
                         neighbour[].y,
                         neighbour[].z,
                         vxyz)
end

function setup_pXest_connectivity(
  coarse_discrete_model::DiscreteModel{Dc,Dp}) where {Dc,Dp}

  trian=Triangulation(coarse_discrete_model)
  node_coordinates=Gridap.Geometry.get_node_coordinates(trian)
  cell_nodes_ids=Gridap.Geometry.get_cell_node_ids(trian)

  if (Dc==2)
    pconn=p4est_connectivity_new(
        p4est_topidx_t(length(node_coordinates)),         # num_vertices
        p4est_topidx_t(num_cells(coarse_discrete_model)), # num_trees
        p4est_topidx_t(0),
        p4est_topidx_t(0))
  else
    @assert Dc==3
    pconn=p8est_connectivity_new(
        p4est_topidx_t(length(node_coordinates)),         # num_vertices
        p4est_topidx_t(num_cells(coarse_discrete_model)), # num_trees
        p4est_topidx_t(0),
        p4est_topidx_t(0),
        p4est_topidx_t(0),
        p4est_topidx_t(0))
  end

  conn=pconn[]
  vertices=unsafe_wrap(Array, conn.vertices, length(node_coordinates)*3)
  current=1
  for i=1:length(node_coordinates)
    p=node_coordinates[i]
    for j=1:Dp
      vertices[current]=Cdouble(p[j])
      current=current+1
    end
    if (Dp==2)
      vertices[current]=Cdouble(0.0) # Z coordinate always to 0.0 in 2D
      current=current+1
    end
  end
  #print("XXX", vertices, "\n")

  tree_to_vertex=unsafe_wrap(Array, conn.tree_to_vertex, length(cell_nodes_ids)*(2^Dc))
  c=Gridap.Arrays.array_cache(cell_nodes_ids)
  current=1
  for j=1:length(cell_nodes_ids)
     ids=Gridap.Arrays.getindex!(c,cell_nodes_ids,j)
     for id in ids
      tree_to_vertex[current]=p4est_topidx_t(id-1)
      current=current+1
     end
  end

  # /*
  #  * Fill tree_to_tree and tree_to_face to make sure we have a valid
  #  * connectivity.
  #  */
  PXEST_FACES=2*Dc
  tree_to_tree=unsafe_wrap(Array, conn.tree_to_tree, conn.num_trees*PXEST_FACES )
  tree_to_face=unsafe_wrap(Array, conn.tree_to_face, conn.num_trees*PXEST_FACES )
  for tree=1:conn.num_trees
    for face=1:PXEST_FACES
      tree_to_tree[PXEST_FACES * (tree-1) + face] = tree-1
      tree_to_face[PXEST_FACES * (tree-1) + face] = face-1
    end
  end

  if (Dc==2)
    p4est_connectivity_complete(pconn)
    @assert Bool(p4est_connectivity_is_valid(pconn))
  else
    p8est_connectivity_complete(pconn)
    @assert Bool(p8est_connectivity_is_valid(pconn))
  end
  pconn
end

function setup_pXest(::Type{Val{Dc}}, comm, connectivity, num_uniform_refinements) where Dc
   if (Dc==2)
       p4est_new_ext(comm,
                     connectivity,
                     Cint(0), Cint(num_uniform_refinements), Cint(1), Cint(0),
                     C_NULL, C_NULL)
   else
      p8est_new_ext(comm,
                    connectivity,
                    Cint(0), Cint(num_uniform_refinements), Cint(1), Cint(0),
                    C_NULL, C_NULL)
   end
end

function setup_pXest_ghost(::Type{Val{Dc}}, ptr_pXest) where Dc
  if (Dc==2)
    p4est_ghost_new(ptr_pXest,P4est_wrapper.P4EST_CONNECT_FULL)
  else
    p8est_ghost_new(ptr_pXest,P4est_wrapper.P8EST_CONNECT_FULL)
  end
end

function setup_cell_prange(::Type{Val{Dc}},
                           parts::MPIData{<:Integer},
                           ptr_pXest,
                           ptr_pXest_ghost) where Dc
  comm = parts.comm

  pXest_ghost = ptr_pXest_ghost[]
  pXest       = ptr_pXest[]

  # Obtain ghost quadrants
  if (Dc==2)
    ptr_ghost_quadrants = Ptr{p4est_quadrant_t}(pXest_ghost.ghosts.array)
  else
    ptr_ghost_quadrants = Ptr{p8est_quadrant_t}(pXest_ghost.ghosts.array)
  end
  proc_offsets = unsafe_wrap(Array, pXest_ghost.proc_offsets, pXest_ghost.mpisize+1)

  global_first_quadrant = unsafe_wrap(Array,
                                      pXest.global_first_quadrant,
                                      pXest.mpisize+1)

  noids,firstgid,hid_to_gid,hid_to_part=map_parts(parts) do part
    hid_to_gid   = Vector{Int}(undef, pXest_ghost.ghosts.elem_count)
    hid_to_part  = Vector{Int32}(undef, pXest_ghost.ghosts.elem_count)
    k=1
    for i=1:pXest_ghost.mpisize
      for j=proc_offsets[i]:proc_offsets[i+1]-1
        quadrant       = ptr_ghost_quadrants[j+1]
        piggy3         = quadrant.p.piggy3
        hid_to_gid[k]  = global_first_quadrant[i]+piggy3.local_num+1
        hid_to_part[k] = Int32(i)
        k=k+1
      end
    end
    pXest.local_num_quadrants,global_first_quadrant[part]+1,hid_to_gid,hid_to_part
  end
  ngids = pXest.global_num_quadrants
  PRange(
    parts,
    ngids,
    noids,
    firstgid,
    hid_to_gid,
    hid_to_part;
    reuse_parts_rcv=true) # Symmetric communication pattern
end

function setup_pXest_lnodes(::Type{Val{Dc}}, ptr_pXest, ptr_pXest_ghost) where Dc
  if (Dc==2)
    p4est_lnodes_new(ptr_pXest, ptr_pXest_ghost, Cint(1))
  else
    p8est_lnodes_new(ptr_pXest, ptr_pXest_ghost, Cint(1))
  end
end

function generate_cell_vertex_gids(ptr_pXest_lnodes, cell_prange)
  pXest_lnodes=ptr_pXest_lnodes[]

  nvertices = pXest_lnodes.vnodes
  element_nodes = unsafe_wrap(Array,
                              pXest_lnodes.element_nodes,
                              pXest_lnodes.num_local_elements*nvertices)

  nonlocal_nodes = unsafe_wrap(Array,
                               pXest_lnodes.nonlocal_nodes,
                               pXest_lnodes.num_local_nodes-pXest_lnodes.owned_count)

  k = 1
  cell_vertex_gids = map_parts(cell_prange.partition) do indices
    n = length(indices.lid_to_part)
    ptrs = Vector{Int32}(undef,n+1)
    ptrs[1]=1
    for i=1:n
      ptrs[i+1]=ptrs[i]+nvertices
    end
    k=1
    current=1
    data = Vector{Int}(undef, ptrs[n+1]-1)
    for i=1:pXest_lnodes.num_local_elements
      for j=1:nvertices
        l=element_nodes[k+j-1]
        if (l < pXest_lnodes.owned_count)
          data[current]=pXest_lnodes.global_offset+l+1
        else
          data[current]=nonlocal_nodes[l-pXest_lnodes.owned_count+1]+1
        end
        current=current+1
      end
      k=k+nvertices
    end
    PArrays.Table(data,ptrs)
  end
  t=async_exchange!(PArrays._replace,cell_vertex_gids,cell_prange.exchanger)
  map_parts(wait∘schedule,t)
  # map_parts(get_part_ids(cell_vertex_gids)) do part
  #   if (part==2)
  #     println(cell_vertex_gids.part)
  #   end
  # end
  cell_vertex_gids
end

function generate_cell_vertex_lids_nlvertices(cell_vertex_gids)
  map_parts(cell_vertex_gids) do cell_vertex_gids
    g2l=Dict{Int,Int}()
    current=1
    data=Vector{Int}(undef,length(cell_vertex_gids.data))
    for (i,gid) in enumerate(cell_vertex_gids.data)
      if haskey(g2l,gid)
        data[i]=g2l[gid]
      else
        data[i]=current
        g2l[gid]=current
        current=current+1
      end
    end
    (PArrays.Table(data,cell_vertex_gids.ptrs), current-1)
  end
end

function generate_node_coordinates(::Type{Val{Dc}},
                                  cell_vertex_lids,
                                  nlvertices,
                                  ptr_pXest_connectivity,
                                  ptr_pXest,
                                  ptr_pXest_ghost) where Dc

  PXEST_CORNERS=2^Dc
  pXest_ghost = ptr_pXest_ghost[]
  pXest       = ptr_pXest[]

  # Obtain ghost quadrants
  if (Dc==2)
    ptr_ghost_quadrants = Ptr{p4est_quadrant_t}(pXest_ghost.ghosts.array)
  else
    ptr_ghost_quadrants = Ptr{p8est_quadrant_t}(pXest_ghost.ghosts.array)
  end

  tree_offsets = unsafe_wrap(Array, pXest_ghost.tree_offsets, pXest_ghost.num_trees+1)
  dnode_coordinates=map_parts(cell_vertex_lids,nlvertices) do cell_vertex_lids, nl
     node_coordinates=Vector{Point{Dc,Float64}}(undef,nl)
     current=1
     vxy=Vector{Cdouble}(undef,Dc)
     pvxy=pointer(vxy,1)
     cell_lids=cell_vertex_lids.data
     for itree=1:pXest_ghost.num_trees
       if (Dc==2)
         tree = p4est_tree_array_index(pXest.trees, itree-1)[]
       else
         tree = p8est_tree_array_index(pXest.trees, itree-1)[]
       end
       for cell=1:tree.quadrants.elem_count
          if (Dc==2)
            quadrant=p4est_quadrant_array_index(tree.quadrants, cell-1)[]
          else
            quadrant=p8est_quadrant_array_index(tree.quadrants, cell-1)[]
          end
          for vertex=1:PXEST_CORNERS
             if (Dc==2)
               p4est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                     p4est_topidx_t(itree-1),
                                                     quadrant.x,
                                                     quadrant.y,
                                                     quadrant.level,
                                                     Cint(vertex-1),
                                                     pvxy)
             else
               p8est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                     p4est_topidx_t(itree-1),
                                                     quadrant.x,
                                                     quadrant.y,
                                                     quadrant.z,
                                                     quadrant.level,
                                                     Cint(vertex-1),
                                                     pvxy)
             end

            node_coordinates[cell_lids[current]]=Point{Dc,Float64}(vxy...)
            current=current+1
          end
       end
     end

     # Go over ghost cells
     for i=1:pXest_ghost.num_trees
      for j=tree_offsets[i]:tree_offsets[i+1]-1
          quadrant = ptr_ghost_quadrants[j+1]
          for vertex=1:PXEST_CORNERS
            if (Dc==2)
               p4est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                     p4est_topidx_t(i-1),
                                                     quadrant.x,
                                                     quadrant.y,
                                                     quadrant.level,
                                                     Cint(vertex-1),
                                                     pvxy)
            else
              p8est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                     p4est_topidx_t(i-1),
                                                     quadrant.x,
                                                     quadrant.y,
                                                     quadrant.z,
                                                     quadrant.level,
                                                     Cint(vertex-1),
                                                     pvxy)

            end
          #  if (MPI.Comm_rank(comm.comm)==0)
          #     println(vxy)
          #  end
           node_coordinates[cell_lids[current]]=Point{Dc,Float64}(vxy...)
           current=current+1
         end
       end
     end
     node_coordinates
  end
end

function generate_grid_and_topology(::Type{Val{Dc}},
                                    cell_vertex_lids,
                                    nlvertices,
                                    node_coordinates) where {Dc}
  grid,topology=
  map_parts(cell_vertex_lids,nlvertices,node_coordinates) do cell_vertex_lids, nl, node_coordinates
    polytope= Dc==2 ? QUAD : HEX
    scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
    cell_types=collect(Fill(1,length(cell_vertex_lids)))
    cell_reffes=[scalar_reffe]
    cell_vertex_lids_gridap=Gridap.Arrays.Table(cell_vertex_lids.data,cell_vertex_lids.ptrs)
    grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                            cell_vertex_lids_gridap,
                                            cell_reffes,
                                            cell_types,
                                            Gridap.Geometry.NonOriented())

    topology = Gridap.Geometry.UnstructuredGridTopology(node_coordinates,
                                      cell_vertex_lids_gridap,
                                      cell_types,
                                      map(Gridap.ReferenceFEs.get_polytope, cell_reffes),
                                      Gridap.Geometry.NonOriented())
    grid,topology
  end
  grid,topology
end

function generate_face_labeling(parts,
                                cell_prange,
                                coarse_discrete_model::DiscreteModel{Dc,Dp},
                                grid,
                                topology,
                                ptr_pXest,
                                ptr_pXest_ghost) where {Dc,Dp}

  pXest       = ptr_pXest[]
  pXest_ghost = ptr_pXest_ghost[]

  coarse_grid_topology  = Gridap.Geometry.get_grid_topology(coarse_discrete_model)
  coarse_grid_labeling  = Gridap.Geometry.get_face_labeling(coarse_discrete_model)

  coarse_cell_vertices = Gridap.Geometry.get_faces(coarse_grid_topology,Dc,0)
  if (Dc==3)
    coarse_cell_edgets = Gridap.Geometry.get_faces(coarse_grid_topology,Dc,1)
  end
  coarse_cell_facets   = Gridap.Geometry.get_faces(coarse_grid_topology,Dc,Dc-1)

  owned_trees_offset=Vector{Int}(undef,pXest_ghost.num_trees+1)
  owned_trees_offset[1]=0
  for itree=1:pXest_ghost.num_trees
    if Dc==2
      tree = p4est_tree_array_index(pXest.trees, itree-1)[]
    else
      tree = p8est_tree_array_index(pXest.trees, itree-1)[]
    end
    owned_trees_offset[itree+1]=owned_trees_offset[itree]+tree.quadrants.elem_count
  end

 faces_to_entity=map_parts(grid,topology) do grid, topology
     # Iterate over corners
     num_vertices=Gridap.Geometry.num_faces(topology,0)
     vertex_to_entity=zeros(Int,num_vertices)
     cell_vertices=Gridap.Geometry.get_faces(topology,Dc,0)
    #  if part==1
    #    println(cell_vertices)
    #  end

     # Corner iterator callback
     function jcorner_callback(pinfo     :: Ptr{p8est_iter_corner_info_t},
                               user_data :: Ptr{Cvoid})
        info=pinfo[]
        if (Dc==2)
          sides=Ptr{p4est_iter_corner_side_t}(info.sides.array)
        else
          sides=Ptr{p8est_iter_corner_side_t}(info.sides.array)
        end
        nsides=info.sides.elem_count
        tree=sides[1].treeid+1
        # We are on the interior of a tree
        data=sides[1]
        if data.is_ghost==1
           ref_cell=pXest.local_num_quadrants+data.quadid+1
        else
           ref_cell=owned_trees_offset[tree]+data.quadid+1
        end
        corner=sides[1].corner+1
        ref_cornergid=cell_vertices[ref_cell][corner]
        # if (MPI.Comm_rank(comm.comm)==0)
        #   println("XXX ", ref_cell, " ", ref_cornergid, " ", info.tree_boundary, " ", nsides, " ", corner)
        # end
        if (info.tree_boundary!=0 && info.tree_boundary==P4est_wrapper.P4EST_CONNECT_CORNER)
              # The current corner is also a corner of the coarse mesh
              coarse_cornergid=coarse_cell_vertices[tree][corner]
              vertex_to_entity[ref_cornergid]=
                 coarse_grid_labeling.d_to_dface_to_entity[1][coarse_cornergid]
        else
          if vertex_to_entity[ref_cornergid]==0
            # We are on the interior of a tree (if we did not touch it yet)
            vertex_to_entity[ref_cornergid]=coarse_grid_labeling.d_to_dface_to_entity[Dc+1][tree]
          end
        end
        # if (MPI.Comm_rank(comm.comm)==0)
        #   println("YYY ", cell_vertices)
        # end
        nothing
     end

     #  C-callable face callback
     ccorner_callback = @cfunction($jcorner_callback,
                                   Cvoid,
                                   (Ptr{p8est_iter_corner_info_t},Ptr{Cvoid}))

     cell_edgets=Gridap.Geometry.get_faces(topology,Dc,1)
     num_edgets=Gridap.Geometry.num_faces(topology,1)
     edget_to_entity=zeros(Int,num_edgets)
     if (Dc==3)
       # Edge iterator callback
       function jedge_callback(pinfo     :: Ptr{p8est_iter_edge_info_t},
                               user_data :: Ptr{Cvoid})
         info=pinfo[]
         sides=Ptr{p8est_iter_edge_side_t}(info.sides.array)

         nsides=info.sides.elem_count
         tree=sides[1].treeid+1
         # We are on the interior of a tree
         data=sides[1].is.full
         if data.is_ghost==1
           ref_cell=pXest.local_num_quadrants+data.quadid+1
         else
           ref_cell=owned_trees_offset[tree]+data.quadid+1
         end
         edge=sides[1].edge+1

         polytope=HEX
         poly_faces=Gridap.ReferenceFEs.get_faces(polytope)
         poly_edget_range=Gridap.ReferenceFEs.get_dimrange(polytope,1)
         poly_first_edget=first(poly_edget_range)
         poly_facet=poly_first_edget+edge-1

        #  if (MPI.Comm_rank(comm.comm)==0)
        #    coarse_edgetgid=coarse_cell_edgets[tree][edge]
        #    coarse_edgetgid_entity=coarse_grid_labeling.d_to_dface_to_entity[2][coarse_edgetgid]
        #    println("PPP ", ref_cell, " ", edge, " TB ", info.tree_boundary, " ", nsides, " ",coarse_edgetgid, " ",  coarse_edgetgid_entity)
        #  end
         if (info.tree_boundary!=0 && info.tree_boundary==P4est_wrapper.P8EST_CONNECT_EDGE)
          coarse_edgetgid=coarse_cell_edgets[tree][edge]
          coarse_edgetgid_entity=coarse_grid_labeling.d_to_dface_to_entity[2][coarse_edgetgid]
          # We are on the boundary of coarse mesh or inter-octree boundary
          for poly_incident_face in poly_faces[poly_facet]
            if poly_incident_face == poly_facet
              ref_edgetgid=cell_edgets[ref_cell][edge]
              edget_to_entity[ref_edgetgid]=coarse_edgetgid_entity
            else
              ref_cornergid=cell_vertices[ref_cell][poly_incident_face]
              # if (MPI.Comm_rank(comm.comm)==0)
              #    println("CCC ", ref_cell, " ", ref_cornergid, " ", info.tree_boundary, " ", nsides)
              # end
              vertex_to_entity[ref_cornergid]=coarse_edgetgid_entity
            end
          end
         else
          # We are on the interior of the domain if we did not touch the edge yet
          ref_edgetgid=cell_edgets[ref_cell][edge]
          if (edget_to_entity[ref_edgetgid]==0)
            edget_to_entity[ref_edgetgid]=coarse_grid_labeling.d_to_dface_to_entity[Dc+1][tree]
          end
         end
         nothing
       end

       # C-callable edge callback
       cedge_callback = @cfunction($jedge_callback,
                                   Cvoid,
                                   (Ptr{p8est_iter_edge_info_t},Ptr{Cvoid}))
     end

     # Iterate over faces
     num_faces=Gridap.Geometry.num_faces(topology,Dc-1)
     facet_to_entity=zeros(Int,num_faces)
     cell_facets=Gridap.Geometry.get_faces(topology,Dc,Dc-1)

     # Face iterator callback
     function jface_callback(pinfo     :: Ptr{p8est_iter_face_info_t},
                             user_data :: Ptr{Cvoid})
        info=pinfo[]
        if Dc==2
          sides=Ptr{p4est_iter_face_side_t}(info.sides.array)
        else
          sides=Ptr{p8est_iter_face_side_t}(info.sides.array)
        end

        nsides=info.sides.elem_count
        tree=sides[1].treeid+1
        # We are on the interior of a tree
        data=sides[1].is.full
        if data.is_ghost==1
           ref_cell=pXest.local_num_quadrants+data.quadid+1
        else
           ref_cell=owned_trees_offset[tree]+data.quadid+1
        end
        face=sides[1].face+1
        if Dc==2
          gridap_facet=P4EST_2_GRIDAP_FACET_2D[face]
        else
          gridap_facet=P4EST_2_GRIDAP_FACET_3D[face]
        end

        polytope= Dc==2 ? QUAD : HEX
        poly_faces=Gridap.ReferenceFEs.get_faces(polytope)
        poly_facet_range=Gridap.ReferenceFEs.get_dimrange(polytope,Dc-1)
        poly_first_facet=first(poly_facet_range)
        poly_facet=poly_first_facet+gridap_facet-1

        # if (MPI.Comm_rank(comm.comm)==0)
        #   coarse_facetgid=coarse_cell_facets[tree][gridap_facet]
        #   coarse_facetgid_entity=coarse_grid_labeling.d_to_dface_to_entity[Dc][coarse_facetgid]
        #   println("PPP ", ref_cell, " ", gridap_facet, " ", info.tree_boundary, " ", nsides, " ",coarse_facetgid, " ",  coarse_facetgid_entity)
        # end
        if (info.tree_boundary!=0)
          coarse_facetgid=coarse_cell_facets[tree][gridap_facet]
          coarse_facetgid_entity=coarse_grid_labeling.d_to_dface_to_entity[Dc][coarse_facetgid]
          # We are on the boundary of coarse mesh or inter-octree boundary
          for poly_incident_face in poly_faces[poly_facet]
            if poly_incident_face == poly_facet
              ref_facetgid=cell_facets[ref_cell][gridap_facet]
              facet_to_entity[ref_facetgid]=coarse_facetgid_entity
            elseif (Dc==3 && poly_incident_face in Gridap.ReferenceFEs.get_dimrange(polytope,1))
              poly_first_edget=first(Gridap.ReferenceFEs.get_dimrange(polytope,1))
              edget=poly_incident_face-poly_first_edget+1
              ref_edgetgid=cell_edgets[ref_cell][edget]
              edget_to_entity[ref_edgetgid]=coarse_facetgid_entity
            else
              ref_cornergid=cell_vertices[ref_cell][poly_incident_face]
              # if (MPI.Comm_rank(comm.comm)==0)
              #    println("CCC ", ref_cell, " ", ref_cornergid, " ", info.tree_boundary, " ", nsides)
              # end
              vertex_to_entity[ref_cornergid]=coarse_facetgid_entity
            end
          end
        else
          # We are on the interior of the domain
          ref_facegid=cell_facets[ref_cell][gridap_facet]
          facet_to_entity[ref_facegid]=coarse_grid_labeling.d_to_dface_to_entity[Dc+1][tree]
        end
        nothing
     end

    #  C-callable face callback
    cface_callback = @cfunction($jface_callback,
                                 Cvoid,
                                 (Ptr{p8est_iter_face_info_t},Ptr{Cvoid}))

    # Iterate over cells
    num_cells=Gridap.Geometry.num_faces(topology,Dc)
    cell_to_entity=zeros(Int,num_cells)

    # Face iterator callback
    function jcell_callback(pinfo     :: Ptr{p8est_iter_volume_info_t},
                            user_data :: Ptr{Cvoid})
      info=pinfo[]
      tree=info.treeid+1
      cell=owned_trees_offset[tree]+info.quadid+1
      #println("XXX $(tree) $(cell)")
      cell_to_entity[cell]=coarse_grid_labeling.d_to_dface_to_entity[Dc+1][tree]
      nothing
    end
    ccell_callback = @cfunction($jcell_callback,
                                 Cvoid,
                                 (Ptr{p8est_iter_volume_info_t},Ptr{Cvoid}))

    if (Dc==2)
       p4est_iterate(ptr_pXest,ptr_pXest_ghost,C_NULL,C_NULL,cface_callback,C_NULL)
       p4est_iterate(ptr_pXest,ptr_pXest_ghost,C_NULL,ccell_callback,C_NULL,ccorner_callback)
    else
       p8est_iterate(ptr_pXest,ptr_pXest_ghost,C_NULL,C_NULL,cface_callback,C_NULL,C_NULL)
       p8est_iterate(ptr_pXest,ptr_pXest_ghost,C_NULL,C_NULL,C_NULL,cedge_callback,C_NULL)
       p8est_iterate(ptr_pXest,ptr_pXest_ghost,C_NULL,ccell_callback,C_NULL,C_NULL,ccorner_callback)
    end
    if (Dc==2)
      vertex_to_entity, facet_to_entity, cell_to_entity
    else
      vertex_to_entity, edget_to_entity, facet_to_entity, cell_to_entity
    end
 end

 vertex_to_entity  = map_parts(x->x[1]   , faces_to_entity)
 if Dc==3
   edget_to_entity = map_parts(x->x[2]   , faces_to_entity)
 end
 facet_to_entity   = map_parts(x->x[Dc]  , faces_to_entity)
 cell_to_entity    = map_parts(x->x[Dc+1], faces_to_entity)

 function cell_to_faces(grid,topology,cell_dim,face_dim)
   map_parts(grid,topology) do grid,topology
    Gridap.Geometry.get_faces(topology,cell_dim,face_dim)
   end
 end

 polytope = Dc==2 ? QUAD : HEX

 update_face_to_entity_with_ghost_data!(vertex_to_entity,
                                        cell_prange,
                                        num_faces(polytope,0),
                                        cell_to_faces(grid,topology,Dc,0))

 if Dc==3
   update_face_to_entity_with_ghost_data!(edget_to_entity,
                                          cell_prange,
                                          num_faces(polytope,1),
                                          cell_to_faces(grid,topology,Dc,1))
 end

 update_face_to_entity_with_ghost_data!(facet_to_entity,
                                        cell_prange,
                                        num_faces(polytope,Dc-1),
                                        cell_to_faces(grid,topology,Dc,Dc-1))

 update_face_to_entity_with_ghost_data!(cell_to_entity,
                                        cell_prange,
                                        num_faces(polytope,Dc),
                                        cell_to_faces(grid,topology,Dc,Dc))

 if (Dc==2)
  faces_to_entity=[vertex_to_entity,facet_to_entity,cell_to_entity]
 else
  faces_to_entity=[vertex_to_entity,edget_to_entity,facet_to_entity,cell_to_entity]
 end

 face_labeling =
  map_parts(faces_to_entity...) do faces_to_entity...
    # if (part == 1)
    #    println("XXX", faces_to_entity[1])
    #    println("XXX", faces_to_entity[2])
    #    println("XXX", faces_to_entity[3])
    #    #println("XXX", faces_to_entity[4])
    # end
    d_to_dface_to_entity       = Vector{Vector{Int}}(undef,Dc+1)
    d_to_dface_to_entity[1]    = faces_to_entity[1]
    if (Dc==3)
      d_to_dface_to_entity[2]  = faces_to_entity[2]
    end
    d_to_dface_to_entity[Dc]   = faces_to_entity[Dc]
    d_to_dface_to_entity[Dc+1] = faces_to_entity[Dc+1]
    Gridap.Geometry.FaceLabeling(d_to_dface_to_entity,
                                 coarse_grid_labeling.tag_to_entities,
                                 coarse_grid_labeling.tag_to_name)
  end
  face_labeling
end

# function p4est_connectivity_print(pconn::Ptr{p4est_connectivity_t})
#   # struct p4est_connectivity
#   #   num_vertices::p4est_topidx_t
#   #   num_trees::p4est_topidx_t
#   #   num_corners::p4est_topidx_t
#   #   vertices::Ptr{Cdouble}
#   #   tree_to_vertex::Ptr{p4est_topidx_t}
#   #   tree_attr_bytes::Csize_t
#   #   tree_to_attr::Cstring
#   #   tree_to_tree::Ptr{p4est_topidx_t}
#   #   tree_to_face::Ptr{Int8}
#   #   tree_to_corner::Ptr{p4est_topidx_t}
#   #   ctt_offset::Ptr{p4est_topidx_t}
#   #   corner_to_tree::Ptr{p4est_topidx_t}
#   #   corner_to_corner::Ptr{Int8}
#   # end
#   conn = pconn[]
#   println("num_vertices=$(conn.num_vertices)")
#   println("num_trees=$(conn.num_trees)")
#   println("num_corners=$(conn.num_corners)")
#   vertices=unsafe_wrap(Array, conn.vertices, conn.num_vertices*3)
#   println("vertices=$(vertices)")
# end

function _fill_data!(data,entry::Integer,k)
  data[k]=entry
  k=k+1
end

function _fill_data!(data,entry::Vector{<:Integer},k)
  for i=1:length(entry)
    data[k]=entry[i]
    k=k+1
  end
  k
end

function init_cell_to_face_entity(num_faces_x_cell,
                                  cell_to_faces,
                                  face_to_entity)
  ptrs = Vector{Int32}(undef,length(cell_to_faces) + 1)
  ptrs[2:end] .= num_faces_x_cell
  Gridap.Arrays.length_to_ptrs!(ptrs)
  data = Vector{eltype(face_to_entity)}(undef, ptrs[end] - 1)
  cell_to_face_entity = lazy_map(Broadcasting(Reindex(face_to_entity)),cell_to_faces)
  k = 1
  for i = 1:length(cell_to_face_entity)
    for j = 1:length(cell_to_face_entity[i])
      k=_fill_data!(data,cell_to_face_entity[i][j],k)
    end
  end
  return PArrays.Table(data, ptrs)
end

function update_face_to_entity!(face_to_entity, cell_to_faces, cell_to_face_entity)
  for cell in 1:length(cell_to_faces)
      i_to_entity = cell_to_face_entity[cell]
      pini = cell_to_faces.ptrs[cell]
      pend = cell_to_faces.ptrs[cell + 1] - 1
      for (i, p) in enumerate(pini:pend)
        lid = cell_to_faces.data[p]
        face_to_entity[lid] = i_to_entity[i]
      end
  end
end

function update_face_to_entity_with_ghost_data!(
   face_to_entity,cell_prange,num_faces_x_cell,cell_to_faces)

   part_to_cell_to_entity = map_parts(init_cell_to_face_entity,
                                      map_parts(x->num_faces_x_cell,cell_prange.partition),
                                      cell_to_faces,
                                      face_to_entity)

   t=async_exchange!(PArrays._replace,part_to_cell_to_entity,cell_prange.exchanger)
   map_parts(wait∘schedule,t)

   map_parts(update_face_to_entity!,
             face_to_entity,
             cell_to_faces,
             part_to_cell_to_entity)
end

"""
  See P4est_wrapper.jl/src/bindings/sc_common.jl for possible/valid
  argument values for the p4est_verbosity_level parameter
"""
function UniformlyRefinedForestOfOctreesDiscreteModel(
    parts::MPIData{<:Integer},
    coarse_discrete_model::DiscreteModel{Dc,Dp},
    num_uniform_refinements::Int;
    p4est_verbosity_level=P4est_wrapper.SC_LP_DEFAULT) where {Dc,Dp}

  sc_init(parts.comm, Cint(true), Cint(true), C_NULL, p4est_verbosity_level)
  p4est_init(C_NULL, p4est_verbosity_level)

  comm = parts.comm
  ptr_pXest_connectivity=setup_pXest_connectivity(coarse_discrete_model)

  # Create a new forest
  ptr_pXest = setup_pXest(Val{Dc},comm,ptr_pXest_connectivity,num_uniform_refinements)

  # Build the ghost layer
  ptr_pXest_ghost=setup_pXest_ghost(Val{Dc},ptr_pXest)

  cell_prange = setup_cell_prange(Val{Dc},parts,ptr_pXest,ptr_pXest_ghost)

  ptr_pXest_lnodes=setup_pXest_lnodes(Val{Dc}, ptr_pXest, ptr_pXest_ghost)

  cell_vertex_gids=generate_cell_vertex_gids(ptr_pXest_lnodes,cell_prange)

  cell_vertex_lids,nlvertices=generate_cell_vertex_lids_nlvertices(cell_vertex_gids)

  node_coordinates=generate_node_coordinates(Val{Dc},
                                             cell_vertex_lids,
                                             nlvertices,
                                             ptr_pXest_connectivity,
                                             ptr_pXest,
                                             ptr_pXest_ghost)

  grid,topology=generate_grid_and_topology(Val{Dc},
                                           cell_vertex_lids,
                                           nlvertices,
                                           node_coordinates)

  face_labeling=generate_face_labeling(parts,
                       cell_prange,
                       coarse_discrete_model,
                       grid,
                       topology,
                       ptr_pXest,
                       ptr_pXest_ghost)

  discretemodel=
    map_parts(grid,topology,face_labeling) do grid, topology, face_labeling
      Gridap.Geometry.UnstructuredDiscreteModel(grid,topology,face_labeling)
    end

  # Write forest to VTK file
  # #p4est_vtk_write_file(unitsquare_forest, C_NULL, "my_step")

  if (Dc==2)
    # Destroy lnodes
    p4est_lnodes_destroy(ptr_pXest_lnodes)
    # Destroy ghost
    p4est_ghost_destroy(ptr_pXest_ghost)
    # Destroy the forest
    p4est_destroy(ptr_pXest)
    # Destroy the connectivity
    p4est_connectivity_destroy(ptr_pXest_connectivity)
  else
    # Destroy lnodes
    p8est_lnodes_destroy(ptr_pXest_lnodes)
    # Destroy ghost
    p8est_ghost_destroy(ptr_pXest_ghost)
    # Destroy the forest
    p8est_destroy(ptr_pXest)
    # Destroy the connectivity
    p8est_connectivity_destroy(ptr_pXest_connectivity)
  end
  sc_finalize()

  GridapDistributed.DistributedDiscreteModel(discretemodel,cell_prange)
end


text:
using Documenter, DiffusionBayes

makedocs(;
    modules=[DiffusionBayes],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/DiffusionBayes.jl/blob/{commit}{path}#L{line}",
    sitename="DiffusionBayes.jl",
    authors="Marcin Mider",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/DiffusionBayes.jl",
)

ext: jl


text:
# Inspired from JuMP/test/solvers

function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

# Load available solvers
grb = try_import(:Gurobi)
cpx = try_import(:CPLEX)
xpr = try_import(:Xpress)
mos = false && try_import(:Mosek)
cbc = try_import(:Cbc)
if cbc; import Clp; end
import GLPKMathProgInterface
glp = true
ipt = try_import(:Ipopt)
eco = try_import(:ECOS)

# Create LP solver list
lp_solvers = Any[]
grb && push!(lp_solvers, Gurobi.GurobiSolver(OutputFlag=0))
cpx && push!(lp_solvers, CPLEX.CplexSolver(CPX_PARAM_SCRIND=0))
xpr && push!(lp_solvers, Xpress.XpressSolver(OUTPUTLOG=0))
mos && push!(lp_solvers, Mosek.MosekSolver(LOG=0))
cbc && push!(lp_solvers, Clp.ClpSolver())
glp && push!(lp_solvers, GLPKMathProgInterface.GLPKSolverLP())
ipt && push!(lp_solvers, Ipopt.IpoptSolver(print_level=0))
eco && push!(lp_solvers, ECOS.ECOSSolver(verbose=false))


text:
@testset "Bert" begin
  using Transformers.BidirectionalEncoder
  using Transformers.BidirectionalEncoder: bert_cased_tokenizer, bert_uncased_tokenizer, WordPiece
  for f ∈ readdir("./bert/")
    include("./bert/$f")
  end
end


text:
using Agents, Random, DataFrames, LightGraphs
using Distributions: Poisson, DiscreteNonParametric
using DrWatson: @dict
using LinearAlgebra: diagind

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::Int
    days_infected::Int  # number of days since is infected
    status::Symbol  # 1: S, 2: I, 3:R
end

function model_initiation(;
    Ns = [500, 500, 500],
    migration_rate = 0.2,
    death_rate = 0.1,
    β_det = 0.05,
    β_und = 0.1,
    infection_period = 10,
    reinfection_probability = 0.05,
    detection_time = 3,
    Is = ones(Int, length(Ns)),
    seed = nothing,
)

    if !isnothing(seed)
        Random.seed!(seed)
    end

    C = length(Ns)
    β_det = [β_det for i in 1:C]
    β_und = [β_und for i in 1:C]

    migration_rates = reshape([migration_rate for i in 1:(C * C)], C, C)
    # normalize migration_rates
    migration_rates_sum = sum(migration_rates, dims = 2)
    for c in 1:C
        migration_rates[c, :] ./= migration_rates_sum[c]
    end

    properties = @dict(
        Ns,
        β_und,
        β_det,
        migration_rates,
        infection_period,
        reinfection_probability,
        detection_time,
        death_rate,
        Is,
        C,
    )
    space = GraphSpace(complete_digraph(C))
    model = ABM(PoorSoul, space; properties)

    # Add initial individuals
    for city in 1:C, n in 1:Ns[city]
        add_agent!(city, model, 0, :S) # Susceptible
    end
    # add infected individuals
    for city in 1:C
        inds = ids_in_position(city, model)
        for n in 1:Is[city]
            agent = model[inds[n]]
            agent.status = :I # Infected
            agent.days_infected = 1
        end
    end
    return model
end

function agent_step!(agent, model)
    migrate!(agent, model)
    transmit!(agent, model)
    update_status!(agent, model)
    recover_or_die!(agent, model)
end

function migrate!(agent, model)
    pid = agent.pos
    d = DiscreteNonParametric(1:(model.C), model.migration_rates[pid, :])
    m = rand(model.rng, d)
    if m ≠ pid
        move_agent!(agent, m, model)
    end
end

function transmit!(agent, model)
    agent.status == :S && return
    rate = if agent.days_infected < model.detection_time
        model.β_und[agent.pos]
    else
        model.β_det[agent.pos]
    end

    d = Poisson(rate)
    n = rand(model.rng, d)
    n == 0 && return

    for contactID in ids_in_position(agent, model)
        contact = model[contactID]
        if contact.status == :S ||
           (contact.status == :R && rand(model.rng) ≤ model.reinfection_probability)
            contact.status = :I
            n -= 1
            n == 0 && return
        end
    end
end

update_status!(agent, model) = agent.status == :I && (agent.days_infected += 1)

function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand(model.rng) ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end


text:
module Hawkes

export hawkes_simulate, hawkes_loglikelihood

"""
Simulate a univariate Hawkes process with intensity:

`λ(t) = α + sum_{t_k<t}{α*exp(δ*(t-t_k))}`
# Arguments
- `u::Float64`: background intensity
- `α::Float64`: arrival intensity jump-size
- `δ::Float64`: exponential decay rate
- `n::Int`: number of times to simulate
"""
function hawkes_simulate(u::Float64, α::Float64, δ::Float64, n::Int)::Vector{Float64}
    λ = u
    t = 0.0
    ts = Vector{Float64}(undef, n)
    for i in 1:n
        D = 1.0 + δ*log(rand())/(λ - u)
        dt1 = -log(rand())/u
        dt2 = min(dt1, -log(max(0.0,D))/δ)
        dt = ifelse(D < 0.0, dt1, dt2)
        t += dt
        ts[i] = t
        λ = α + (λ - u)*exp(-δ*dt) + u
    end
    ts
end

"""
Simulate a multivariate Hawkes process with intensities:

`λ(t)_i = α_i + sum_{t_k<t}{α_(i,j)*exp(δ_i*(t-t_k))}`
# Arguments
- `u::AbstractVector{Float64}`: background intensities
- `α::AbstractMatrix{Float64}`: arrival intensity jump-sizes
- `δ::AbstractVector{Float64}`: exponential decay rates
- `n::Int`: number of times to simulate
"""
function hawkes_simulate(u::AbstractVector{Float64}, α::AbstractMatrix{Float64}, δ::AbstractVector{Float64}, n::Int)::Vector{Vector{Float64}}
    t = 0.0
    M = length(u)
    λ = copy(u)
    ts = [Vector{Float64}() for _ in 1:M]
    for i in 1:n
        dt_min = Inf
        dt_min_m = 0
        for m in 1:M
            D = 1.0 + δ[m]*log(rand())/(λ[m] - u[m])
            dt1 = -log(rand())/u[m]
            dt2 = min(dt1, -log(max(0.0,D))/δ[m])
            dt = ifelse(D < 0.0, dt1, dt2)
            if dt < dt_min
                dt_min = dt
                dt_min_m = m
            end
        end
        t += dt_min
        push!(ts[dt_min_m], t)
        for m in 1:M
            λ[m] = α[m,dt_min_m] + (λ[m] - u[m])*exp(-δ[m]*dt_min) + u[m]
        end
    end
    ts
end

# Real types are used rather than Float64 so that automatic differentiation methods will work.
"""
Calculate the log-likelihood for a univariate Hawkes process with intensity:

`λ(t) = α + sum_{t_k<t}{α*exp(δ*(t-t_k))}`
# Arguments
- `u::Real`: background intensity
- `α::Real`: arrival intensity jump-size
- `δ::Real`: exponential decay rate
- `ts::AbstractVector{Real}`: arrival times
"""
function hawkes_loglikelihood(u::Real, α::Real, δ::Real, ts::AbstractVector)::Real
    ll1 = 0.0
    ll2 = 0.0
    t = 0.0
    λ = u
    N = length(ts)
    T = ts[N]
    A = 0.0
    for i in 1:N
        dt = ts[i] - t
        t = ts[i]
        ll1 += log(u + α*A)
        ll2 += -1.0 + exp(-δ*(T-t))
        dtn = ts[min(N,i+1)] - ts[i]
        A = exp(-δ*dtn)*(1.0 + A)
    end
    ll1 + α*ll2/δ - u*T
end

# Real types are used rather than Float64 so that automatic differentiation methods will work.
"""
Calculate the log-likelihood for a multivariate Hawkes process with intensities:

`λ(t)_i = α_i + sum_{t_k<t}{α_(i,j)*exp(δ_i*(t-t_k))}`
# Arguments
- `u::AbstractVector{Real}`: background intensities
- `α::AbstractMatrix{Real}`: arrival intensity jump-sizes
- `δ::AbstractVector{Real}`: exponential decay rates
- `ts::AbstractVector{AbstractVector{Real}}`: arrival times
"""
function hawkes_loglikelihood(u::AbstractVector{R},α::AbstractMatrix{R},δ::AbstractVector{R},ts::AbstractVector)::R where R<:Real
    res = 0.0
    M = length(ts)
    T = maximum(last.(ts))
    for m in 1:M
        s = 0.0
        Rdiag = zeros(R, length(ts[m]))
        Rnondiag = zeros(R, length(ts[m]))
        ix = ones(Int, M)
        for i in 2:length(ts[m])
            dt = ts[m][i] - ts[m][i-1]
            Rdiag[i] = (1.0 + Rdiag[i-1])*exp(-δ[m]*dt)
            Rnondiag[i] = Rnondiag[i-1]*exp(-δ[m]*dt)
            for n in 1:M
                m == n && continue
                for j in ix[n]:length(ts[n])
                    if ts[n][j] >= ts[m][i-1]
                        if ts[n][j] < ts[m][i]
                            Rnondiag[i] += exp(-δ[m]*(ts[m][i] - ts[n][j]))
                        else
                            ix[n] = j
                            break
                        end
                    end
                end
            end
        end
        for n in 1:M
            for j in 1:length(ts[n])
                s += α[m,n]*(1.0 - exp(-δ[m]*(T - ts[n][j])))/δ[m]
            end
        end
        res += -u[m]*T - s
        for i in 1:length(ts[m])
            s = u[m]
            for n in 1:M
                s += α[m,n]*ifelse(m==n, Rdiag[i], Rnondiag[i])
            end
            res += log(s)
        end
    end
    res
end

end # module


text:
# core algorithm using IterativeSolvers.jl
function update!(pa, x, attrib, ::UseIterativeSolvers)
	update_prepare!(pa, attrib)

	# prepare b
	for i in eachindex(pa.optm.dvec)
		pa.optm.dvec[i]=pa.optm.obs.d[i]
	end

	# create operator
	A=operator(pa, attrib);

	# put the value from model to x
	model_to_x!(x, pa, attrib)

	IterativeSolvers.lsmr!(x, A, pa.optm.dvec)

	x_to_model!(x, pa, attrib)

	J=func_grad!(nothing, x,  pa, attrib)

	update_finalize!(pa, attrib)

	return J
end

# use Optim.jl (ConjugateGradients)
function update!(pa, x, attrib, ::UseOptim)
	update_prepare!(pa, attrib)

	f =x->func_grad!(nothing, x,  pa, attrib)
	g! =(storage, x)->func_grad!(storage, x,  pa, attrib)

	# put the value from model to x
	model_to_x!(x, pa, attrib)

	res = optimize(f, g!,
		x, ConjugateGradient(alphaguess = LineSearches.InitialQuadratic()),
		Optim.Options(g_tol = 1e-8, iterations = 2000, show_trace = false))
	pa.verbose && println(res)

	x_to_model!(Optim.minimizer(res), pa, attrib)

	update_finalize!(pa, attrib)

	return Optim.minimum(res)
end

#=
# use bounded Optim.jl (ConjugateGradients) (for positive Source Time Functions)
function update!(pa, x, attrib, ::UseOptimSTF)

	if(typeof(attrib) ≠ S)
		error("bounded Optim only for S")
	end

	fill!(pa.sx.lower_x, 0.0)
	fill!(pa.sx.upper_x, Inf)

	update_prepare!(pa, attrib)

	f =x->func_grad!(nothing, x,  pa, attrib)
	g! =(storage, x)->func_grad!(storage, x,  pa, attrib)
	println("FFFFAA", maximum(x))

	# put the value from model to x
	model_to_x!(x, pa, attrib)

	# check if x<0?
	for i in eachindex(x)
		if(x[i] < 0.0)
			x[i] = 0.0
		end
	end

	println("FFFF", maximum(x))

	res = optimize(f, g!, pa.sx.lower_x, pa.sx.upper_x,
		x, Fminbox(ConjugateGradient(alphaguess = LineSearches.InitialQuadratic())),
		Optim.Options(g_tol = 1e-8, iterations = 2000, show_trace = false))
	pa.verbose && println(res)

	x_to_model!(Optim.minimizer(res), pa, attrib)

	update_finalize!(pa, attrib)

	return Optim.minimum(res)
end
=#

"""
* re_init_flag :: re-initialize inversions with random input or not?
"""
function update_all!(pa, io=stdout;
		     max_roundtrips=100,
		     max_reroundtrips=1,
		     ParamAM_func=nothing,
		     roundtrip_tol=1e-6,
		     optim_tols=[1e-6, 1e-6], verbose=true,
		     log_file=false,
		     )

	global to, optG, optS
	reset_timer!(to)

	if(ParamAM_func===nothing)
		ParamAM_func=xx->ParamAM(xx, optim_tols=optim_tols,name="Blind Decon",
				    roundtrip_tol=roundtrip_tol, max_roundtrips=max_roundtrips,
				    max_reroundtrips=max_reroundtrips,
				    min_roundtrips=10,
				    verbose=verbose,
				    reinit_func=xxx->initialize!(pa),
				    #after_roundtrip_func=x->(write(io,string(TimerOutputs.flatten(to),"\n"))),
				    #after_roundtrip_func=x->(err!(pa, nothing)),
				    )
	end

	# create alternating minimization parameters
	f1=x->@timeit to "update S()" update!(pa, pa.sx.x, S(), optS)
	f2=x->@timeit to "update G()" update!(pa, pa.gx.x, G(), optG)
	paam=ParamAM_func([f1, f2])

	if(io===nothing)
		logfilename=joinpath(pwd(),string("XBD",Dates.now(),".log"))
		io=open(logfilename, "a+")
	end

	go(paam, io)  # run alternative minimization

	# save log file
	if(log_file)
		AMlogfile=joinpath(pwd(),string("AM",Dates.now(),".log"))
		CSV.write(AMlogfile, paam.log)
	end

	write(io,string(to))
	write(io, "\n")

	# print errors
	err!(pa, io)
	write(io, "\n")
	if(io === nothing)
		close(io)
	end
end

"""
Given g (sign is ambiguous), update s such that s>0
"""
function update_stf!(pa::BD,)
	#(pa.sxp.n ≠ 2) && error("stf update not possible")
	# update source according to the estimated g from fpr
	update_prepare!(pa, S())
	update!(pa, pa.sx.x, S(), optS)
	update_finalize!(pa, S())
	# save functional and estimated source
	f1 = Misfits.error_squared_euclidean!(nothing, pa.optm.cal.d, pa.optm.obs.d,
				       nothing, norm_flag=false)
	s1 = copy(pa.optm.cal.s)

	# what if the sign of g is otherwise
	rmul!(pa.optm.cal.g, -1.0)
	update_prepare!(pa, S())
	update!(pa, pa.sx.x, S(), optS)
	update_finalize!(pa, S())
	f2 = Misfits.error_squared_euclidean!(nothing, pa.optm.cal.d, pa.optm.obs.d,
			       nothing, norm_flag=false)
	s2 = copy(pa.optm.cal.s)

	# sxp is not used here, because the estimated g after fpr has ambiguous sign
	# as a result, positivity constraint cannot be imposed on s
	if(f1<f2)
		copyto!(pa.optm.cal.s, s1)
		rmul!(pa.optm.cal.g, -1.0) # change sign back
	else
		# no need to change sign of g
		copyto!(pa.optm.cal.s, s2)
	end

	return nothing
end


text:
function isvalid(formula)
    material = falses(8)
    for it in formula
        material[it] = true
    end
    ! ( material[1] & material[2] ) &
    ! ( material[3] & material[4] ) &
    ( material[5] == material[6] ) &
    ( material[7] | material[8] )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@inline function save_noise!(W::NoiseProcess)
  if W.t != W.curt
    push!(W.W,copy(W.curW))
    push!(W.t,copy(W.curt))
    if W.Z != nothing
      push!(W.Z,copy(W.curZ))
    end
  end
end

@inline function accept_step!(W::NoiseProcess,dt,setup_next=true)

  W.curt += W.dt
  W.iter += 1

  if isinplace(W)
    @.. W.curW += W.dW
    if W.Z != nothing
      @.. W.curZ += W.dZ
    end
  else
    W.curW += W.dW
    if W.Z != nothing
      W.curZ += W.dZ
    end
  end

  if W.save_everystep
    push!(W.W,copy(W.curW))
    push!(W.t,copy(W.curt))
    if W.Z != nothing
      push!(W.Z,copy(W.curZ))
    end
  end

  W.dt = dt #dtpropose
  # Setup next step
  if setup_next
    setup_next_step!(W::NoiseProcess)
  end
end

@inline function setup_next_step!(W::NoiseProcess)
  if adaptive_alg(W)==:RSwM3
    ResettableStacks.reset!(W.S₂) #Empty W.S₂
  end
  if adaptive_alg(W)==:RSwM1
    if !isempty(W.S₁)
      W.dt,W.dW,W.dZ = pop!(W.S₁)
    else # Stack is empty
      calculate_step!(W,W.dt)
    end
  elseif adaptive_alg(W)==:RSwM2 || adaptive_alg(W)==:RSwM3
    if !isinplace(W)
      dttmp = zero(W.dt); W.dW = zero(W.dW)
      if W.Z != nothing
        W.dZ = zero(W.dZ)
      end
    else
      dttmp = zero(W.dt); fill!(W.dW,zero(eltype(W.dW)))
      if W.Z != nothing
        fill!(W.dZ,zero(eltype(W.dZ)))
      end
    end
    while !isempty(W.S₁)
      L₁,L₂,L₃ = pop!(W.S₁)
      qtmp = (W.dt-dttmp)/L₁
      if qtmp>1
        dttmp+=L₁
        if isinplace(W)
          @.. W.dW+=L₂
          if W.Z != nothing
            @.. W.dZ+=L₃
          end
        else
          W.dW+=L₂
          if W.Z != nothing
            W.dZ+=L₃
          end
        end
        if adaptive_alg(W)==:RSwM3
          if L₁ > W.rswm.discard_length
            push!(W.S₂,(L₁,L₂,L₃))
          end
        end
      else #Popped too far
        # Generate numbers to bridge and step perfectly
        dttmp += qtmp*L₁
        if isinplace(W)
          W.bridge(W.dWtilde,W,W.curW,L₂,qtmp,L₁,W.rng)
          if W.Z != nothing
            W.bridge(W.dZtilde,W,W.curZ,L₃,qtmp,L₁,W.rng)
          end
        else
          W.dWtilde = W.bridge(W,W.curW,L₂,qtmp,L₁,W.rng)
          if W.Z != nothing
            W.dZtilde = W.bridge(W,W.curZ,L₃,qtmp,L₁,W.rng)
          end
        end
        if isinplace(W)
          @.. W.dW += W.dWtilde
          if W.Z != nothing
            @.. W.dZ += W.dZtilde
          end
        else
          W.dW += W.dWtilde
          if W.Z != nothing
            W.dZ += W.dZtilde
          end
        end
        if (1-qtmp)*L₁ > W.rswm.discard_length
          if isinplace(W)
            @.. L₂ -= W.dWtilde
            if W.Z != nothing
              @.. L₃ -= W.dZtilde
            end
          else
            L₂ -= W.dWtilde
            if W.Z != nothing
              L₃ -= W.dZtilde
            end
          end
          if W.Z == nothing
            push!(W.S₁,((1-qtmp)*L₁,L₂,nothing))
          else
            push!(W.S₁,((1-qtmp)*L₁,L₂,L₃))
          end
          if adaptive_alg(W)==:RSwM3 && qtmp*L₁ > W.rswm.discard_length
            if W.Z == nothing
              ResettableStacks.copyat_or_push!(W.S₂,(qtmp*L₁,W.dWtilde,nothing))
            else
              ResettableStacks.copyat_or_push!(W.S₂,(qtmp*L₁,W.dWtilde,W.dZtilde))
            end
          end
        end
        break
      end
    end #end while empty
    # This is a control variable so do not diff through it
    dtleft = DiffEqBase.ODE_DEFAULT_NORM(W.dt - dttmp,W.curt)
    if dtleft > W.rswm.discard_length #Stack emptied
      if isinplace(W)
        W.dist(W.dWtilde,W,dtleft,W.rng)
        if W.Z != nothing
          W.dist(W.dZtilde,W,dtleft,W.rng)
        end
      else
        W.dWtilde = W.dist(W,dtleft,W.rng)
        if W.Z != nothing
          W.dZtilde = W.dist(W,dtleft,W.rng)
        end
      end
      if isinplace(W)
        @.. W.dW += W.dWtilde
        if W.Z != nothing
          @.. W.dZ += W.dZtilde
        end
      else
        W.dW += W.dWtilde
        if W.Z != nothing
          W.dZ += W.dZtilde
        end
      end
      if adaptive_alg(W)==:RSwM3
        if W.Z == nothing
          ResettableStacks.copyat_or_push!(W.S₂,(dtleft,W.dWtilde,nothing))
        else
          ResettableStacks.copyat_or_push!(W.S₂,(dtleft,W.dWtilde,W.dZtilde))
        end
      end
    end
  end # End RSwM2 and RSwM3
end

@inline function calculate_step!(W::NoiseProcess,dt)
  if isinplace(W)
    W.dist(W.dW,W,dt,W.rng)
    if W.Z != nothing
      W.dist(W.dZ,W,dt,W.rng)
    end
  else
    W.dW = W.dist(W,dt,W.rng)
    if W.Z != nothing
      W.dZ = W.dist(W,dt,W.rng)
    end
  end
  W.dt = dt
end

@inline function reject_step!(W::NoiseProcess,dtnew)
  q = dtnew/W.dt
  if adaptive_alg(W)==:RSwM1 || adaptive_alg(W)==:RSwM2
    if isinplace(W)
      W.bridge(W.dWtilde,W,0,W.dW,q,dtnew,W.rng)
      if W.Z != nothing
        W.bridge(W.dZtilde,W,0,W.dZ,q,dtnew,W.rng)
      end
    else
      W.dWtilde = W.bridge(W,0,W.dW,q,dtnew,W.rng)
      if W.Z != nothing
        W.dZtilde=  W.bridge(W,0,W.dZ,q,dtnew,W.rng)
      end
    end
    cutLength = W.dt-dtnew
    if cutLength > W.rswm.discard_length
      if W.Z == nothing
        push!(W.S₁,(cutLength,W.dW-W.dWtilde,nothing))
      else
        push!(W.S₁,(cutLength,W.dW-W.dWtilde,W.dZ-W.dZtilde))
      end
    end
    if length(W.S₁) > W.maxstacksize
        W.maxstacksize = length(W.S₁)
    end
    if isinplace(W)
      copyto!(W.dW,W.dWtilde)
      if W.Z!=nothing
        copyto!(W.dZ,W.dZtilde)
      end
    else
      W.dW = W.dWtilde
      if W.Z != nothing
        W.dZ = W.dZtilde
      end
    end
    W.dt = dtnew
  else # RSwM3
    if !isinplace(W)
      dttmp = zero(W.dt); W.dWtmp = zero(W.dW)
      if W.Z != nothing
        W.dZtmp = zero(W.dZtmp)
      end
    else
      dttmp = zero(W.dt); fill!(W.dWtmp,zero(eltype(W.dWtmp)))
      if W.Z!= nothing
        fill!(W.dZtmp,zero(eltype(W.dZtmp)))
      end
    end
    if length(W.S₂) > W.maxstacksize2
      W.maxstacksize2= length(W.S₂)
    end
    while !isempty(W.S₂)
      L₁,L₂,L₃ = pop!(W.S₂)
      if dttmp + L₁ < (1-q)*W.dt #while the backwards movement is less than chop off
        dttmp += L₁
        if isinplace(W)
          @.. W.dWtmp += L₂
          if W.Z != nothing
            @.. W.dZtmp += L₃
          end
        else
          W.dWtmp += L₂
          if W.Z != nothing
            W.dZtmp += L₃
          end
        end
        push!(W.S₁,(L₁,L₂,L₃))
      else
        push!(W.S₂,(L₁,L₂,L₃))
        break
      end
    end # end while
    dtK = W.dt - dttmp
    qK = q*W.dt/dtK
    if isinplace(W)
      @.. W.dWtmp = W.dW - W.dWtmp
      if W.Z != nothing
        @.. W.dZtmp = W.dZ - W.dZtmp
      end
    else
      W.dWtmp = W.dW - W.dWtmp
      if W.Z != nothing
        W.dZtmp = W.dZ - W.dZtmp
      end
    end
    if isinplace(W)
      W.bridge(W.dWtilde,W,0,W.dWtmp,qK,dtK,W.rng)
      #W.dWtilde .-= W.curW
      if W.Z != nothing
        W.bridge(W.dZtilde,W,0,W.dZtmp,qK,dtK,W.rng)
        #W.dZtilde .-= W.curZ
      end
    else
      W.dWtilde = W.bridge(W,0,W.dWtmp,qK,dtK,W.rng)# - W.curW
      if W.Z != nothing
        W.dZtilde = W.bridge(W,0,W.dZtmp,qK,dtK,W.rng)# - W.curZ
      end
    end
    # This is a control variable so do not diff through it
    cutLength = DiffEqBase.ODE_DEFAULT_NORM((1-qK)*dtK,W.curt)
    if cutLength > W.rswm.discard_length
      if W.Z == nothing
        push!(W.S₁,(cutLength,W.dWtmp-W.dWtilde,nothing))
      else
        push!(W.S₁,(cutLength,W.dWtmp-W.dWtilde,W.dZtmp-W.dZtilde))
      end
    end
    if length(W.S₁) > W.maxstacksize
        W.maxstacksize = length(W.S₁)
    end
    W.dt = dtnew
    if isinplace(W)
      copyto!(W.dW,W.dWtilde)
      if W.Z != nothing
        copyto!(W.dZ,W.dZtilde)
      end
    else
      W.dW = W.dWtilde
      if W.Z != nothing
        W.dZ = W.dZtilde
      end
    end
  end
end

@inline function interpolate!(W::NoiseProcess,t)
  if t > W.t[end] # Steps past W
    dt = t - W.t[end]
    if isinplace(W)
      W.dist(W.dW,W,dt,W.rng)
      W.curW .+= W.dW
      if W.Z != nothing
        W.dist(W.dZ,W,dt,W.rng)
        W.curZ .+= W.dZ
      end
    else
      W.dW = W.dist(W,dt,W.rng)
      W.curW += W.dW
      if W.Z != nothing
        W.dZ = W.dist(W,dt,W.rng)
        W.curZ += W.dZ
      end
    end
    out1 = copy(W.curW)
    if W.save_everystep
      push!(W.t,t)
      push!(W.W,out1)
    end
    if W.Z != nothing
      out2= copy(W.curZ)
      if W.save_everystep
        push!(W.Z,out2)
      end
    else
      out2 = nothing
    end
    return out1,out2
  else # Bridge
    i = searchsortedfirst(W.t,t)
    if t == W.t[i]
      if isinplace(W)
        W.curW .= W.W[i]
      else
        W.curW = W.W[i]
      end
      if W.Z != nothing
        if isinplace(W)
          W.curZ .= W.Z[i]
        else
          W.curZ = W.Z[i]
        end
        return copy(W.curW),copy(W.curZ)
      else
        return copy(W.curW),nothing
      end
    else
      W0,Wh = W.W[i-1],W.W[i]
      if W.Z != nothing
        Z0,Zh = W.Z[i-1],W.Z[i]
      end
      h = W.t[i]-W.t[i-1]
      q = (t-W.t[i-1])/h
      if isinplace(W)
        new_curW = similar(W.dW)
        W.bridge(new_curW,W,W0,Wh,q,h,W.rng)
        new_curW .+= (1-q)*W0
        if W.Z != nothing
          new_curZ = similar(W.dZ)
          W.bridge(new_curZ,W,Z0,Zh,q,h,W.rng)
          new_curZ .+= (1-q)*Z0
        else
          new_curZ = nothing
        end
      else
        new_curW = W.bridge(W,W0,Wh,q,h,W.rng)
        new_curW += (1-q)*W0
        if W.Z != nothing
          new_curZ = W.bridge(W,Z0,Zh,q,h,W.rng)
          new_curZ += (1-q)*Z0
        else
          new_curZ = nothing
        end
      end
      W.curW = new_curW
      if W.save_everystep
        insert!(W.W,i,new_curW)
        insert!(W.t,i,t)
      end
      if W.Z != nothing
        W.curZ = new_curZ
        if W.save_everystep
          insert!(W.Z,i,new_curZ)
        end
      end
      return new_curW,new_curZ
    end
  end
end

@inline function interpolate!(out1,out2,W::NoiseProcess,t)
  if t > W.t[end] # Steps past W
    dt = t - W.t[end]
    W.dist(W.dW,W,dt,W.rng)
    out1 .+= W.dW
    if W.Z != nothing
      W.dist(W.dZ,W,dt,W.rng)
      out2 .+= W.dZ
    end
    if W.save_everystep
      push!(W.t,t)
      push!(W.W,copy(out1))
      if W.Z != nothing
        push!(W.Z,copy(out2))
      end
    end
  else # Bridge
    i = searchsortedfirst(W.t,t)
    if t == W.t[i]
      out1 .= W.W[i]
      if W.Z != nothing
        out2 .= W.Z[i]
      end
    else
      W0,Wh = W.W[i-1],W.W[i]
      if W.Z != nothing
        Z0,Zh = W.Z[i-1],W.Z[i]
      end
      h = W.t[i]-W.t[i-1]
      q = (t-W.t[i-1])/h
      W.bridge(out1,W,W0,Wh,q,h,W.rng)
      out1 .+= (1-q)*W0
      if W.Z != nothing
        W.bridge(out2,W,Z0,Zh,q,h,W.rng)
        out2 .+= (1-q)*Z0
      end
      W.curW .= out1
      if W.save_everystep
        insert!(W.W,i,copy(out1))
        insert!(W.t,i,t)
      end
      if W.Z != nothing
        W.curZ .= out2
        if W.save_everystep
          insert!(W.Z,i,copy(out2))
        end
      end
    end
  end
end

function resize_stack!(W::NoiseProcess,i)
  for j in eachindex(W.S₂.data)
    resize!(W.S₂.data[j][2],i)
    W.S₂.data[j][3] != nothing && resize!(W.S₂.data[j][3],i)
  end
end

function deleteat_stack!(W::NoiseProcess,i)
  for j in eachindex(W.S₂.data)
    deleteat!(W.S₂.data[j][2],i)
    W.S₂.data[j][3] != nothing && deleteat!(W.S₂.data[j][3],i)
  end
end

#=
function addat_stack!(W::NoiseProcess,i)
  for j in eachindex(W.S₂.data)
    resize!(W.S₂.data[j],i)
    W.S₂.data[j][3] != nothing && resize!(W.S₂.data[j][3],i)
  end
end
=#


text:
#Module for D4M
module D4M

    using LinearAlgebra, SparseArrays, PyPlot, DelimitedFiles

    import SparseArrays: nnz, diag
    import Base: &, ==, >, <, -, *, +, /
    import Base: isless, getindex, isempty, print, size, sum, transpose,
            Array, Matrix, adjoint, broadcast
    import PyPlot: spy
    import JLD: writeas, readas

    export  Assoc,
            StartsWith,
            CatKeyMul, CatValMul,
            CatStr, SplitStr, NumStr,
            col2type, val2col,
            ReadCSV, WriteCSV,
            # writeas, readas,
            print, printFull, printTriple,
            norow, nocol,
            logical, str2num, convertvals,
            sqIn, sqOut,
            putAdj, putRow, putCol, putVal,
            getadj, getrow, getcol, getval, find,
            #saveassoc, loadassoc,
            OutDegree, InDegree, diag,
            bounded, strictbounded, adjbfs

    include("Assoc.jl") # Associative Array

    #Helper functions for parsing
    include("parsinghelpers.jl")
    #Helper functions for working with strings and string arrays
    include("stringarrayhelpers.jl")

    if haskey(ENV,"JAVA_HOME")
        # Database functionality
        include("DB/DBserver.jl")
        include("DB/DBtable.jl")
        # Grapulo Calls
        include("DB/bfs.jl")
        include("DB/tablemult.jl")
        include("DB/jaccard.jl")
        include("DB/nmf.jl")
        include("DB/ktruss.jl")

        using JavaCall
        # be sure to keep adding stuff here, so we don't have to prepend calls w "D4M."
        export dbsetup, ls # DBserver
        export delete, addColCombiner, put, putTriple, getiterator, getsplits, addsplits, nnz # DBtable
        export makedegreetable, adjbfs, edgebfs, singlebfs # BFS
        export tablemult # tablemult
        export jaccard # jaccard
        export nmf # NMF
        export ktrussadj, ktrussedge # ktruss
    else
        println("Not loading database capabilities. If you would like to connect to a database, please set JAVA_HOME.")
        println("After setting JAVA_HOME, run:")
        println("D4Mpkg = Base.PkgId(\"D4M\"); Base.compilecache(D4Mpkg)")
    end
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
module MathTeXEngine
# Adapted from matplotlib mathtext engine
# https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/_mathtext.py

using AbstractTrees
using Automa
using FreeTypeAbstraction
using LaTeXStrings

import Automa.RegExp: @re_str
import DataStructures: Stack
import FreeTypeAbstraction:
    ascender, boundingbox, descender, get_extent, hadvance, inkheight, inkwidth,
    leftinkbound, rightinkbound, topinkbound, bottominkbound
import GeometryBasics: Point2f0
import REPL.REPLCompletions: latex_symbols

const re = Automa.RegExp

export TeXExpr, texparse
export TeXElement, TeXChar, VLine, HLine, generate_tex_elements

include("parser/texexpr.jl")
include("parser/command_data.jl")
include("parser/parser.jl")

include("engine/fonts.jl")
include("engine/texelements.jl")
include("engine/layout.jl")

end # module


text:
"""
Metadata from a thread.
More details [here](https://discord.com/developers/docs/resources/channel#thread-metadata-object).
"""
struct ThreadMetadata
    archived::Bool
    auto_archive_duration::Int
    archive_timestamp::DateTime
    locked::Bool
    invitable::Optional{Bool}
end
@boilerplate ThreadMetadata :constructors :docs :lower :merge :mock


text:
#Figure 2-i
include("fig2b.jl")

#calculating evenness
function pielou(n)
    np = filter(x -> x > eps(), n)
    p = np./sum(np)
    ev = length(p) == 1 ? 0.0 : -sum(p.*log.(p))*(1/log(length(p)))
    return ev
end;
ev = mapslices(pielou, Np[:,2:end]; dims=2);
plot(ev, c=:black, title = "Evenness", xlabel = "Time", ylabel = "Relative abundance (log)", leg = false, ylims = (0,1))

# png("Figure 2/graph_2j.png")


text:
# Simulator
function dynamics(model::Hopper, v1, q1, q2, q3, u, λ, b, w, h)
   (M_func(model, q1) * v1
   - M_func(model, q2) * (SVector{4}(q3) - SVector{4}(q2)) / h
   + h * (w
   + transpose(B_func(model, q3)) * SVector{2}(u)
   + transpose(N_func(model, q3)) * SVector{1}(λ)
   + transpose(P_func(model, q3)) * SVector{2}(b)
   - G_func(model, q2)))
end

function maximum_dissipation(model::Hopper, q2, q3, ψ, η, h)
   ψ_stack = ψ[1] * ones(model.nb)
   return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function no_slip(model::Hopper, q2, q3, λ, h)
   return (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function friction_cone(model::Hopper, λ, b)
   return @SVector [model.μ * λ[1] - sum(b)]
end

mutable struct MOISimulator <: MOI.AbstractNLPEvaluator
   num_var::Int                 # number of decision variables
   num_con::Int                 # number of constraints
   primal_bounds
   constraint_bounds

   model
   slack_penalty
   v1
   q1
   q2
   u
   w
   h
end

function simulator_problem(model, v1, q1, q2, u, w, h; slack_penalty = 1.0e5)
   num_var = model.nq + model.nc + model.nb + model.nc + model.nb + model.ns
   num_con = model.nq + model.nc + model.nc + model.nb +  3

   zl = zeros(num_var)
   zl[1:model.nq] .= -Inf
   zu = Inf * ones(num_var)

   cl = zeros(num_con)
   cu = zeros(num_con)
   cu[model.nq .+ (1:model.nc)] .= Inf
   cu[model.nq + model.nc .+ (1:model.nc)] .= Inf
   cu[model.nq + model.nc + model.nc + model.nb + 1] = Inf
   cu[model.nq + model.nc + model.nc + model.nb + 2] = Inf
   cu[model.nq + model.nc + model.nc + model.nb + 3] = Inf

   MOISimulator(num_var,
	   num_con,
	   (zl, zu),
	   (cl, cu),
	   model,
	   slack_penalty,
	   v1,
	   q1,
	   q2,
	   u,
	   w,
	   h)
end
# prob_sim = simulator_problem(model, ones(nq), ones(nq), ones(nq), zeros(nu), ones(nq), h)
# z0 = rand(prob_sim.num_var)

function MOI.eval_objective(prob::MOISimulator, x)
   return prob.slack_penalty * x[prob.num_var]
end

# MOI.eval_objective(prob_sim, z0)

function MOI.eval_objective_gradient(prob::MOISimulator, grad_f, x)
   grad_f .= 0.0
   grad_f[prob.num_var] = prob.slack_penalty
   return nothing
end

# ∇obj = zeros(prob_sim.num_var)
# MOI.eval_objective_gradient(prob_sim, ∇obj, z0)

function MOI.eval_constraint(prob::MOISimulator, g, x)
   model = prob.model
   q3 = view(x, 1:model.nq)
   λ = view(x, model.nq .+ (1:model.nc))
   b = view(x, model.nq + model.nc .+ (1:model.nb))
   ψ = view(x, model.nq + model.nc + model.nb .+ (1:model.nc))
   η = view(x, model.nq + model.nc + model.nb + model.nc .+ (1:model.nb))
   s = view(x, model.nq + model.nc + model.nb + model.nc + model.nb .+ (1:model.ns))
   # λ = 0
   # b = 0
   g[1:model.nq] = dynamics(model,
	   prob.v1, prob.q1, prob.q2, q3, prob.u, λ, b, prob.w, prob.h)
   g[model.nq .+ (1:model.nc)] = ϕ_func(model, q3)
   g[model.nq + model.nc .+ (1:model.nc)] = friction_cone(model, λ, b)
   g[model.nq + model.nc + model.nc .+ (1:model.nb)] = maximum_dissipation(model,
	   prob.q2, q3, ψ, η, prob.h)
   g[model.nq + model.nc + model.nc + model.nb + 1] = s[1] - (λ' * ϕ_func(model, q3))[1]
   g[model.nq + model.nc + model.nc + model.nb + 2] = s[1] - (ψ' * friction_cone(model, λ, b))[1]
   g[model.nq + model.nc + model.nc + model.nb + 3] = s[1] - (η' * b)[1]

   return nothing
end
# prob_sim.h
# c0 = zeros(prob_sim.num_con)
# MOI.eval_constraint(prob_sim, c0, z0)

function MOI.eval_constraint_jacobian(prob::MOISimulator, jac, x)
   con!(g, z) = MOI.eval_constraint(prob, g, z)
   ForwardDiff.jacobian!(reshape(jac, prob.num_var, prob.num_con), con!, zeros(prob.num_con), x)
   return nothing
end

# ∇c0 = vec(zeros(prob_sim.num_con, prob_sim.num_var))
# MOI.eval_constraint_jacobian(prob_sim, ∇c0, z0)
# #
function sparsity_jacobian(prob::MOISimulator)
   row = []
   col = []

   r = (1:prob.num_con)
   c = (1:prob.num_var)

   row_col!(row, col, r, c)

   return collect(zip(row, col))
end

# sparsity_jacobian(prob_sim)
#
# @time solve(prob_sim, z0)

function step_contact(model, v1, q1, q2, u, w, h)
   prob = simulator_problem(model, v1, q1, q2, u, w, h)
   z0 = [copy(q2); 1.0e-5 * ones(model.nc + model.nb + model.nc + model.nb + model.ns)]
   @time z , info = solve(prob, copy(z0), tol = 1.0e-7, c_tol = 1.0e-7, mapl = 0)

   @assert z[end] < 1.0e-5

   return z[1:model.nq]
end

# h = h̄[1]
# q_sim = [q_proj[1], q_proj[2]]
# v_sim = [(q_proj[2] - q_proj[1]) / h]

# for t = 1:T-1
#    # d x rate
#    _q_sim = [q_sim[end-1], q_sim[end]]
#    _v_sim = [v_sim[end]]
#
#    d = 1
#    for i = 1:d
# 	   _h = h / convert(Float64, d)
# 	   _q = step_contact(model,
# 		   _v_sim[end], _q_sim[end-1], _q_sim[end], u_proj[t][model.idx_u], zeros(model.d), _h)
# 	   _v = (_q - _q_sim[end]) / _h
#
# 	   push!(_q_sim, _q)
# 	   push!(_v_sim, _v)
#    end
#
#    push!(q_sim, _q_sim[end])
#    push!(v_sim, _v_sim[end])
#
#    # same rate
#    # push!(q_sim, step_contact(model,
#    # 	v_sim[end], q_sim[end-1], q_sim[end], ū[t][model.idx_u], zeros(model.d), h))
#    # push!(v_sim, (q_sim[end] - q_sim[end-1]) / h)
# end

# plot(hcat(q_proj...)',
#    color = :red,
#    width = 2.0,
#    label = "")
# plot!(hcat(q_sim...)',
#    color = :black,
#    width = 1.0,
#    label = "")
#
# plot(hcat(v_sim...)',
#    color = :black,
#    width = 1.0,
#    label = "")


text:
################################################################################
# InteratomicPotentials API default implementations for non-trainable potentials
################################################################################

get_parameters(::NonTrainablePotential) = NamedTuple()
set_parameters(p::NonTrainablePotential, ::NamedTuple{()}) = p

serialize_parameters(::NonTrainablePotential) = []
deserialize_parameters(p::NonTrainablePotential, ::AbstractVector) = p

get_hyperparameters(::NonTrainablePotential) = NamedTuple()
set_hyperparameters(p::NonTrainablePotential, ::NamedTuple{()}) = p

serialize_hyperparameters(::NonTrainablePotential) = []
deserialize_hyperparameters(p::NonTrainablePotential, ::AbstractVector) = p


max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.3 General/1.2.3.3 (d+e x^n)^q (a+b x^n+c x^(2 n))^p.jl


text:
using DrWatson
@quickactivate "unfoldjl_dev"

using Random
include("sim_utilities.jl")

##---
paramList = Dict(
    "f" => @formula(dv ~ 1 + condition  + (1+condition|subj)),
            #@formula dv ~ 1 + age * condition  + (1+condition|item) + (1+condition|subj)],
    "θ" => [[1., 1., 1.]],
    "σ" => 1.,
    "β" => [[0., 0.]],
    "blup_method" => ["ranef_scaled","olsranefjf", "olsranef"],
    "residual_method" => [:signflip,:shuffle],
    "nRep" => 1000,
    "nPerm"=> 1000,
)

##--
dl = dict_list(paramList)[1]
simMod = sim_model(dl["f"])
res = run_permutationtest(MersenneTwister(1),simMod,dl["nPerm"],dl["β"],dl["σ"],dl["θ"])

##---
nWorkers=20
for dl = dict_list(paramList)

    res = run_permutationtest_distributed(nWorkers,dl["nRep"],simMod,dl["nPerm"],dl["β"],dl["σ"],dl["θ"],dl["residual_method"],dl["blup_method"])

end
using DataFrames, AlgebraOfGraphics

df = DataFrame(z:=>[],:β=>[],:h1=>[],ranef=>[])
df = DataFrame(:z=>res[1][:],:β=>res[2][:],:h1=>[repeat(["1"],size(res[1],1)); repeat(["0"],size(res[1],1))])

data(df) * mapping(:z,layout_y=:h1) * visual(Hist,bins=0:0.01:1)  |>draw()


text:
# Define the new set structure
struct MyNewSet <: InfiniteOpt.AbstractInfiniteSet
    attr1::Float64 # REPLACE WITH DESIRED TYPE
    attr2::Float64 # REPLACE WITH DESIRED TYPE
    # ADD MORE ATTRIBUTES AS NEEDED
    # constructor
    function MyNewSet(attr1, attr2)
        # PUT APPROPRIATE CHECKS HERE AND OTHER NEEDED OPERATIONS
        return new(convert(Float64, attr1), convert(Float64, attr2))
    end
end

# Extend supports_in_set
function InfiniteOpt.supports_in_set(supports::Union{Number, Vector{<:Number}},
                                     set::MyNewSet)::Bool
    # DETERMINE IF SUPPORTS ARE IN THE DOMAIN OF `set`
    in_set = all(set.attr1 .<= supports .<= set.attr2) # REPLACE WITH ACTUAL CHECK
    return in_set
end

# Extend generate_support_values if possible
function InfiniteOpt.generate_support_values(set::MyNewSet;
                                             num_supports::Int = 10,
                                             sig_fig::Int = 5)::Array
    # REPLACE BELOW WITH METHODS TO GENERATE `num_samples` with `sig_fig`
    supports = collect(range(set.attr1, stop = set.attr2, length = num_supports))
    return round.(supports, sigdigits = sig_fig)
end

# Extend generate_and_add_supports! if supports are not generated for parameters individually
# Here we'll assume generate_support_values is sufficient

# Extend generate_supports_and_coeffs in enable integral() with methods
function InfiniteOpt.MeasureEvalMethods.generate_supports_and_coeffs(
    set::MyNewSet,
    params::Union{InfiniteOpt.ParameterRef, AbstractArray{<:InfiniteOpt.ParameterRef}},
    num_supports::Int,
    lb::Union{Number, JuMPC.SparseAxisArray, Nothing},
    ub::Union{Number, JuMPC.SparseAxisArray, Nothing},
    method::Val; kwargs...
    )::Tuple
    # ADD CHECKS IF NEEDED
    return generate_supports_and_coeffs(IntervalSet(set.attr1, set.attr2), params,
                                        num_supports, lb, ub, method) # REPLACE WITH PROPER DISPATCH
end

# Extend JuMP.has_lower_bound (optional if the answer is always false)
function JuMP.has_lower_bound(set::MyNewSet)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.lower_bound(set::MyNewSet)
    return set.attr1 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_lower_bound if appropriate
function JuMP.set_lower_bound(set::MyNewSet, lower::Number)::MyNewSet
    return MyNewSet(lower, set.attr2) # REPLACE WITH ACTUAL CONSTRUCTOR
end

# Extend JuMP.has_supper_bound (optional if the answer is always false)
function JuMP.has_upper_bound(set::MyNewSet)::Bool
    # INSERT NECESSARY CHECKS IF NEEDED
    has_bound = true # REPLACE WITH ACTUAL RESULT
    return has_bound
end

# Extend JuMP.lower_bound if appropriate
function JuMP.upper_bound(set::MyNewSet)
    return set.attr2 # REPLACE WITH ACTUAL
end

# Extend JuMP.set_upper_bound if appropriate
function JuMP.set_upper_bound(set::MyNewSet, upper::Number)::MyNewSet
    return MyNewSet(set.attr1, upper) # REPLACE WITH ACTUAL CONSTRUCTOR
end


text:
module StateSpace

    using Distributions, ForwardDiff
    using Statistics, LinearAlgebra, Random

    # model types
    export AbstractStateSpaceModel
    export AbstractSSM
    export AbstractGaussianSSM
    export LinearGaussianSSM
    export NonlinearGaussianSSM
    export NonlinearSSM

    # filter types
    export AbstractStateSpaceFilter
    export AbstractKalmanFilter
    export LinearKalmanFilter
    export NonlinearKalmanFilter
    export NonlinearFilter
    export KalmanFilter
    export KF
    export ExtendedKalmanFilter
    export EKF
    export UnscentedKalmanFilter
    export UKF
    export EnsembleKalmanFilter
    export EnKF
    export ParticleFilter
    export FilteredState
    export show
    export process_matrix
    export observation_matrix
    export predict
    export observe
    export update
    export update!
    export filter
    export smooth
    export loglikelihood
    export simulate

    include("matrix_utils.jl")
    include("model_types.jl")
    include("filtered_states.jl")
    include("filter_types.jl")
    include("common.jl")
    include("particle_filter.jl")

end # module


text:
using MultiLine
using Base.Test
using TestSetExtensions

@testset ExtendedTestSet "Running MultiLine tests" begin

    @includetests ARGS

end

max_stars_repo_name: floswald/MultiLine

max_issues_repo_name: floswald/MultiLine

max_forks_repo_name: floswald/MultiLine


text:
using SurrealNumbers

n = 3
x = convert.(SurrealFinite, 0:n)
A = zeros(SurrealFinite, n+1, n+1)
for i=1:n+1
    for j=1:n+1
        A[i,j] = x[i] + x[j]
    end
end
float.(A)
A
generation.(A)
A .== canonicalise.(A)

n = 2
x = convert.(SurrealFinite, (0:2^n)/2^n)
A = zeros(SurrealFinite, 2^n+1, 2^n+1)
for i=1:2^n +1
    for j=1:2^n +1
        A[i,j] = x[i] + x[j]
    end
end
float.(A)
A
generation.(A)
A .== canonicalise.(A)

n = 3
x = convert.(SurrealFinite, 0:n)
A = zeros(SurrealFinite, n+1, n+1)
for i=1:n+1
    for j=1:n+1
        A[i,j] = x[i] * x[j]
    end
end
float.(A)
generation.(A)
A .== canonicalise.(A)

# analyse.(A)

# function analyse(s::SurrealFinite)
#     println("depth = ", generation(s) )

# end


text:
using Merlin
using MLDatasets

include("iddict.jl")

function model(nvocabs::Int)
    T = Float32
    x = Var()
    y = Embedding(T,nvocabs,100)(x)
    y = Linear(T,300,100)(y)
    y = relu(y)
    y = Linear(T,100,nvocabs)(y)
    Graph(y, x)
end

function train()
    train_x, train_y = PTBLM.traindata()
    test_x, test_y = PTBLM.testdata()
    dict = IdDict(["<bos>","<eos>"])
    train_x = map(x -> append!(dict,x), train_x)
    train_y = map(y -> append!(dict,y), train_y)
    test_x = map(x -> append!(dict,x), test_x)
    test_y = map(y -> append!(dict,y), test_y)

    train_x = map(train_x) do x
        x = unshift!(x, 1, 1)
        window(x, (3,), pad=2)
    end
    test_x = map(test_x) do x
        x = unshift!(x, 1, 1)
        window(x, (3,), pad=2)
    end
    train_c = sum(map(length, train_x))

    f = model(length(dict))
    opt = SGD(0.01)
    for epoch = 1:10
        println("epoch: $(epoch)")
        loss = fit(train_x, train_y, f, crossentropy, opt)
        ppl = 2^(loss/train_c)
        println("LOSS: $(loss)")
        println("PPL: $(ppl)")
        println("")
    end
end

train()

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
module operators_lazy

import Base: *, /, +, -
import ..operators

using Base.Cartesian
using ..bases, ..states, ..operators, ..operators_sparse

export LazyOperator, LazyTensor, LazySum, LazyProduct

"""
Abstract base class for lazy operators.

Lazy means in this context that operator-operator operations are not performed
immediately but delayed until the operator is applied to a state, e.g. instead
of (A + B) * x it calculates (A * x) + (B * x).

Lazy operations are implemented for:
    * TensorProduct: LazyTensor
    * Addition: LazySum
    * Multiplication: LazyProduct
"""
abstract LazyOperator <: Operator

"""
Lazy implementation of a tensor product of operators.

The suboperators are stored as values in a dictionary where the key is
the index of the subsystem. Additionally a complex factor is stored in the
"factor" field which allows for fast multiplication with a number.
"""
type LazyTensor <: LazyOperator
    basis_l::CompositeBasis
    basis_r::CompositeBasis
    factor::Complex128
    operators::Dict{Int,Operator}

    function LazyTensor(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int,Operator}, factor::Number=1.)
        N = length(basis_l)
        @assert N==length(basis_r)
        @assert maximum(keys(operators))<=N
        @assert 1<=minimum(keys(operators))
        for (i,op) = operators
            @assert op.basis_l==basis_l.bases[i] && op.basis_r==basis_r.bases[i]
        end
        new(basis_l, basis_r, complex(factor), operators)
    end
end

function LazyTensor{T<:Operator}(basis_l::CompositeBasis, basis_r::CompositeBasis, indices::Vector{Int}, operators::Vector{T}, factor::Number=1.)
    @assert length(indices) == length(Set(indices)) == length(operators)
    LazyTensor(basis_l, basis_r, Dict{Int,Operator}([i=>op for (i,op)=zip(indices, operators)]), factor)
end

LazyTensor(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, operator::Operator, factor::Number=1.) = LazyTensor(basis_l, basis_r, [index], [operator], factor)
LazyTensor(basis::CompositeBasis, indices, operators, factor::Number=1.) = LazyTensor(basis, basis, indices, operators, factor)

"""
Lazy evaluation of sum of operators.

All operators have to be given in respect to the same bases. The field
factors accounts for an additional multiplicative factor for each operator.
"""
type LazySum <: LazyOperator
    basis_l::Basis
    basis_r::Basis
    factors::Vector{Complex128}
    operators::Vector{Operator}

    function LazySum(factors::Vector{Complex128}, operators::Vector{Operator})
        @assert length(operators)>0
        @assert length(operators)==length(factors)
        for i = 2:length(operators)
            @assert operators[1].basis_l == operators[i].basis_l
            @assert operators[1].basis_r == operators[i].basis_r
        end
        new(operators[1].basis_l, operators[1].basis_r, factors, operators)
    end
end
LazySum(operators::Operator...) = LazySum(ones(Complex128, length(operators)), Operator[operators...])

"""
Lazy evaluation of product of operators.

The factors of the product are stored in the "operators" field. Additionally a
complex factor is stored in the "factor" field which allows for fast
multiplication with a number.
"""
type LazyProduct <: LazyOperator
    basis_l::Basis
    basis_r::Basis
    factor::Complex128
    operators::Vector{Operator}

    function LazyProduct(factor::Complex128, operators::Vector{Operator})
        @assert length(operators)>1
        for i = 2:length(operators)
            @assert multiplicable(operators[i-1].basis_r, operators[i].basis_l)
        end
        new(operators[1].basis_l, operators[end].basis_r, factor, operators)
    end
end
LazyProduct(operators::Operator...) = LazyProduct(Complex(1.), Operator[operators...])

function operators.full(x::LazyTensor)
    op_list = DenseOperator[]
    for i=1:length(x.basis_l.bases)
        if i in keys(x.operators)
            push!(op_list, full(x.operators[i]))
        else
            push!(op_list, dense_identity(x.basis_l.bases[i], x.basis_r.bases[i]))
        end
    end
    return x.factor*reduce(tensor, op_list)
end

function operators.full(x::LazySum)
    result = x.factors[1]*full(x.operators[1])
    for i=2:length(x.operators)
        result += x.factors[i]*full(x.operators[i])
    end
    return result
end

function operators.full(x::LazyProduct)
    result = x.factor*full(x.operators[1])
    for i=2:length(x.operators)
        result *= full(x.operators[i])
    end
    return result
end

function operators_sparse.sparse(x::LazyTensor)
    op_list = SparseOperator[]
    for i=1:length(x.basis_l.bases)
        if i in keys(x.operators)
            push!(op_list, sparse(x.operators[i]))
        else
            push!(op_list, identity(x.basis_l.bases[i], x.basis_r.bases[i]))
        end
    end
    return x.factor*reduce(tensor, op_list)
end

function operators_sparse.sparse(x::LazySum)
    result = x.factors[1]*sparse(x.operators[1])
    for i=2:length(x.operators)
        result += x.factors[i]*sparse(x.operators[i])
    end
    return result
end

function operators_sparse.sparse(x::LazyProduct)
    result = x.factor*sparse(x.operators[1])
    for i=2:length(x.operators)
        result *= sparse(x.operators[i])
    end
    return result
end

function operators.dagger(op::LazyTensor)
    D = Dict{Int, Operator}()
    for (i, op_i) in op.operators
        D[i] = dagger(op_i)
    end
    LazyTensor(op.basis_r, op.basis_l, D, conj(op.factor))
end
operators.dagger(op::LazySum) = LazySum(conj(op.factors), Operator[dagger(op_i) for op_i in op.operators])

operators.trace(op::LazyTensor) = op.factor*prod([(haskey(op.operators,i) ? trace(op.operators[i]): prod(op.basis_l.shape)) for i=1:length(op.basis_l)])
operators.trace(op::LazySum) = sum([trace(x) for x=op.operators])

function *(a::LazyTensor, b::LazyTensor)
    check_multiplicable(a.basis_r, b.basis_l)
    a_indices = keys(a.operators)
    b_indices = keys(b.operators)
    D = Dict{Int, Operator}()
    for i=intersect(a_indices, b_indices)
        D[i] = a.operators[i]*b.operators[i]
    end
    for i=setdiff(a_indices, b_indices)
        D[i] = a.operators[i]
    end
    for i=setdiff(b_indices, a_indices)
        D[i] = b.operators[i]
    end
    return LazyTensor(a.basis_l, b.basis_r, D, a.factor*b.factor)
end
*(a::LazyTensor, b::Number) = LazyTensor(a.basis_l, a.basis_r, a.operators, a.factor*b)
*(a::Number, b::LazyTensor) = LazyTensor(b.basis_l, b.basis_r, b.operators, a*b.factor)

# *(a::LazySum, b::LazySum) = LazySum([op1*op2 for op1=a.operators, op2=b.operators]...)
*(a::LazySum, b::LazySum) = LazyProduct(a, b)
*(a::LazySum, b::Number) = LazySum(b*a.factors, a.operators)
*(a::Number, b::LazySum) = LazySum(a*b.factors, b.operators)

*(a::LazyProduct, b::LazyProduct) = LazyProduct(a.factor*b.factor, a.operators..., b.operators...)
*(a::LazyProduct, b::Number) = LazyProduct(a.factor*b, a.operators...)
*(a::Number, b::LazyProduct) = LazyProduct(a*b.factor, b.operators...)

+(a::LazyOperator, b::LazyOperator) = LazySum(a, b)
+(a::LazyOperator, b::Operator) = LazySum(a, b)
+(a::Operator, b::LazyOperator) = LazySum(a, b)
+(a::LazySum, b::LazySum) = LazySum([a.factors, b.factors;], [a.operators, b.operators;])
+(a::LazySum, b::LazyOperator) = LazySum([a.factors, Complex(1.);], [a.operators, b;])
+(a::LazyOperator, b::LazySum) = LazySum([Complex(1.), b.factors;], [a, b.operators;])
+(a::LazySum, b::Operator) = LazySum([a.factors, Complex(1.);], [a.operators, b;])
+(a::Operator, b::LazySum) = LazySum([a.factors, Complex(1.);], [a, b.operators;])

-(a::LazyTensor) = LazyTensor(a.basis_l, a.basis_r, a.operators, -a.factor)
-(a::LazySum) = LazySum(-a.factors, a.operators)
-(a::LazyProduct) = LazyProduct(-a.factor, a.operators)

-(a::LazyOperator, b::LazyOperator) = LazySum([Complex(1.), Complex(-1.)], [a, b])
-(a::LazyOperator, b::Operator) = LazySum([Complex(1.), Complex(-1.)], [a, b])
-(a::Operator, b::LazyOperator) = LazySum([Complex(1.), Complex(-1.)], [a, b])
-(a::LazySum, b::LazySum) = LazySum([a.factors, -b.factors], [a.operators, b.operators])
-(a::LazySum, b::LazyOperator) = LazySum([a.factors, -Complex(1.)], [a.operators, b])
-(a::LazyOperator, b::LazySum) = LazySum([Complex(1.), -a.factors], [a, b.operators])
-(a::LazySum, b::Operator) = LazySum([a.factors, -Complex(1.)], [a.operators, b])
-(a::Operator, b::LazySum) = LazySum([Complex(1.), -a.factors], [a, b.operators])

@generated function _lazytensor_gemv!{RANK, INDEX}(rank::Array{Int, RANK}, index::Array{Int, INDEX},
                                        alpha, op::Operator, b::Ket, beta, result::Ket)
    return quote
        x = Ket(b.basis.bases[$INDEX])
        y = Ket(result.basis.bases[$INDEX])
        indices_others = filter(x->(x!=$INDEX), 1:$RANK)
        shape_others = [b.basis.shape[i] for i=indices_others]
        strides_others = [operators._strides(b.basis.shape)[i] for i=indices_others]
        stride = operators._strides(b.basis.shape)[$INDEX]
        N = b.basis.shape[$INDEX]
        @nexprs 1 d->(I_{$(RANK-1)}=1)
        @nloops $(RANK-1) i d->(1:shape_others[d]) d->(I_{d-1}=I_{d}) d->(I_d+=strides_others[d]) begin
            for j=1:N
                x.data[j] = b.data[I_0+stride*(j-1)]
            end
            operators.gemv!(alpha, op, x, beta, y)
            for j=1:N
                result.data[I_0+stride*(j-1)] = y.data[j]
            end
        end
    end
end

function operators.gemv!(alpha, a::LazyTensor, b::Ket, beta, result::Ket)
    rank = zeros(Int, [0 for i=1:length(a.basis_l.shape)]...)
    bases = [b for b=b.basis.bases]
    for op_index in keys(a.operators)
        index = zeros(Int, [0 for i=1:op_index]...)
        bases[op_index] = a.operators[op_index].basis_l
        tmp = Ket(CompositeBasis(bases...))
        _lazytensor_gemv!(rank, index, Complex(1.), a.operators[op_index], b, Complex(0.), tmp)
        b = tmp
    end
    result.data[:] = beta*result.data[:] + b.data[:]
    nothing
end

@generated function _lazytensor_gemv!{RANK, INDEX}(rank::Array{Int, RANK}, index::Array{Int, INDEX},
                                        alpha, b::Bra, op::Operator, beta, result::Bra)
    return quote
        x = Bra(b.basis.bases[$INDEX])
        y = Bra(result.basis.bases[$INDEX])
        indices_others = filter(x->(x!=$INDEX), 1:$RANK)
        shape_others = [b.basis.shape[i] for i=indices_others]
        strides_others = [operators._strides(b.basis.shape)[i] for i=indices_others]
        stride = operators._strides(b.basis.shape)[$INDEX]
        N = b.basis.shape[$INDEX]
        @nexprs 1 d->(I_{$(RANK-1)}=1)
        @nloops $(RANK-1) i d->(1:shape_others[d]) d->(I_{d-1}=I_{d}) d->(I_d+=strides_others[d]) begin
            for j=1:N
                x.data[j] = b.data[I_0+stride*(j-1)]
            end
            operators.gemv!(alpha, x, op, beta, y)
            for j=1:N
                result.data[I_0+stride*(j-1)] = y.data[j]
            end
        end
    end
end

function operators.gemv!(alpha, a::Bra, b::LazyTensor, beta, result::Bra)
    rank = zeros(Int, [0 for i=1:length(b.basis_r.shape)]...)
    bases = [b for b=a.basis.bases]
    for op_index in keys(b.operators)
        index = zeros(Int, [0 for i=1:op_index]...)
        bases[op_index] = b.operators[op_index].basis_r
        tmp = Bra(CompositeBasis(bases...))
        _lazytensor_gemv!(rank, index, Complex(1.), a, b.operators[op_index], Complex(0.), tmp)
        a = tmp
    end
    result.data[:] = beta*result.data[:] + a.data[:]
    nothing
end

function operators.gemv!(alpha, a::LazySum, b::Ket, beta, result::Ket)
    operators.gemv!(alpha*a.factors[1], a.operators[1], b, beta, result)
    for i=2:length(a.operators)
        operators.gemv!(alpha*a.factors[i], a.operators[i], b, Complex(1.), result)
    end
end

function operators.gemv!(alpha, a::Bra, b::LazySum, beta, result::Bra)
    operators.gemv!(alpha*b.factors[1], a, b.operators[1], beta, result)
    for i=2:length(b.operators)
        operators.gemv!(alpha*b.factors[i], a, b.operators[i], Complex(1.), result)
    end
end

function operators.gemv!(alpha, a::LazyProduct, b::Ket, beta, result::Ket)
    tmp1 = Ket(a.operators[end].basis_l)
    operators.gemv!(Complex(1.)*a.factor, a.operators[end], b, Complex(0.), tmp1)
    for i=length(a.operators)-1:-1:2
        tmp2 = Ket(a.operators[i].basis_l)
        operators.gemv!(Complex(1.), a.operators[i], tmp1, Complex(0.), tmp2)
        tmp1 = tmp2
    end
    operators.gemv!(alpha, a.operators[1], tmp1, beta, result)
end

function operators.gemv!(alpha, a::Bra, b::LazyProduct, beta, result::Bra)
    tmp1 = Bra(b.operators[1].basis_r)
    operators.gemv!(Complex(1.)*b.factor, a, b.operators[1], Complex(0.), tmp1)
    for i=2:length(b.operators)-1
        tmp2 = Bra(b.operators[i].basis_r)
        operators.gemv!(Complex(1.), tmp1, b.operators[i], Complex(0.), tmp2)
        tmp1 = tmp2
    end
    operators.gemv!(alpha, tmp1, b.operators[end], beta, result)
end

end


text:
module SpringCollab2020SpikeJumpThroughController

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/SpikeJumpThroughController" SpikeJumpThroughController(x::Integer, y::Integer, persistent::Bool=false)

const placements = Ahorn.PlacementDict(
    "Spiked Jump Through Controller (Spring Collab 2020)" => Ahorn.EntityPlacement(
        SpikeJumpThroughController,
        "point"
    )
)

sprite = "ahorn/SpringCollab2020/spikeJumpThroughController"

function Ahorn.selection(entity::SpikeJumpThroughController)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SpikeJumpThroughController, room::Maple.Room)
	Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
"""
   HiestFiles is an object that represents a list of file (with paths) that are
   to be processed.
"""
struct HeistFiles
    files::Array{String,1}
    theVector::Cxx.CppValue
end

"""
    HeistFiles(filesArray)

    Creates a `HiestFiles` object representing the list of files to process.

       filesArray is a list of strings that are paths to files. It can also
       contain XrootD URI's.
"""
function HeistFiles(files::Array{String,1})
    fv = icxx"""std::vector<std::string> fv; return fv;"""
    for s ∈ files
        icxx"$fv.push_back($s);"
    end

    HeistFiles(files, fv)
end


text:
import QuantumOptics.timeevolution
import QuantumOptics.stochastic
using QuantumOptics: Basis, Operator, CompositeBasis, Ket, StateVector
using QuantumOptics.timeevolution: DecayRates

#############################################################################################
# Wrap QuantumOptics solvers in order to implement our form of basis check
#############################################################################################

# function timeevolution.master_dynamic(tspan, rho0::T, f::Function;
#             rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing, kwargs...
#         ) where {B<:Basis,T<:Operator{B,B}}
#     check_bases(f(0., 0.).basis_l, rho0.basis_l)
#     timeevolution.master_dynamic(tspan, rho0, f; rates=rates, fout=fout, kwargs...)
# end

# function timeevolution.mcwf_dynamic(tspan, psi0::T, f::Function;
#     seed=rand(UInt), rates::DecayRates=nothing,
#     fout=nothing, display_beforeevent=false, display_afterevent=false,
#     kwargs...) where {T<:Ket}
#     check_bases(f(0., 0.).basis_l, psi0.basis)
#     timeevolution.mcwf_dynamic(tspan, psi0, f; seed=seed, rates=rates, fout=fout,
#             display_beforeevent=display_beforeevent, display_afterevent=display_afterevent,
#             kwargs...
#         )
# end

# function stochastic.schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Function;
#             fout::Union{Function,Nothing}=nothing, noise_processes::Int=0,
#             normalize_state::Bool=false, kwargs...
#         ) where T<:Ket
#     check_bases(fdeterm(0.0, 0.0), fstoch)
#     stochastic.schroedinger_dynamic(tspan, psi0, fdeterm, fstoch;
#             fout=fout, noise_processes=noise_processes, normalize_state=normalize_state,
#             kwargs...
#         )
# end

# function stochastic.master_dynamic(tspan, rho0::T, fdeterm::Function, fstoch::Function;
#             rates::DecayRates=nothing, fout::Union{Function,Nothing}=nothing,
#             noise_processes::Int=0, kwargs...
#         ) where {B<:Basis,T<:Operator{B,B}}
#         check_bases()
#         stochastic.master_dynamic(tspan, rho0, fdeterm, fstoch; rates=rates, fout=fout,
#                 noise_processes=noise_processes, kwargs...
#             )
# end

# function timeevolution.schroedinger_dynamic(tspan, psi0::T, f::Function;
#             fout::Union{Function,Nothing}=nothing, kwargs...
#         ) where T<:StateVector
#     check_bases(f(0., 0.).basis_l, psi0.basis)
#     timeevolution.schroedinger_dynamic(tspan, psi0, f; fout=fout, kwargs...)
# end

# when solving for unitary evolution, it is convenient to switch back and forth between an
# initial state that is pure and one that is mixed
function timeevolution.schroedinger_dynamic(tspan, rho0::T, f::Function;
            fout::Union{Function,Nothing}=nothing, kwargs...
        ) where {B<:Basis,T<:Operator{B,B}}
    check_bases(f(0., 0.).basis_l, rho0.basis_l)
    Jvec = []
    timeevolution.master_dynamic(tspan, rho0, (t,rho) -> (f(t, rho), Jvec, Jvec);
        fout=fout, kwargs...)
end

check_bases(b1::CompositeBasis, b2::CompositeBasis) = b1 == b2 || throw(IncompatibleBases())

mutable struct IncompatibleBases <: Exception end