text:
// Copyright (c) 2020 Graphcore Ltd. All rights reserved.
#include "popops/Rearrange.hpp"

#include "poplibs_support/Tracepoint.hpp"
#include <boost/icl/interval_map.hpp>
#include <boost/optional.hpp>
#include <poplibs_support/gcd.hpp>
#include <poplibs_support/logging.hpp>
#include <poputil/DebugInfo.hpp>
#include <poputil/Util.hpp>
#include <poputil/VarStructure.hpp>
#include <poputil/VertexTemplates.hpp>
#include <poputil/exceptions.hpp>

using namespace poplar;
using namespace poplar::program;
using namespace poplibs_support;
using namespace poputil;

static inline Tensor groupTensorAux(const Tensor &t, unsigned rank) {
  return t;
}
static inline Tensor ungroupTensorAux(const Tensor &t, unsigned) { return t; }

template <typename... G>
static inline Tensor groupTensorAux(const Tensor &t, unsigned rank,
                                    const GroupingInfo &g, G &&...gs) {
  return groupTensorAux(t.reshapePartial(g.first, g.first + 1,
                                         {t.dim(g.first) / g.second, g.second})
                            .dimRoll(g.first + 1, rank),
                        rank + 1, std::forward<G>(gs)...);
}

template <typename... G>
static inline Tensor ungroupTensorAux(const Tensor &t, unsigned rank,
                                      const GroupingInfo &g, G &&...gs) {
  return ungroupTensorAux(
      t.dimRoll(rank, g.first + 1).flatten(g.first, g.first + 2), rank,
      std::forward<G>(gs)...);
}

template <typename... G>
static inline Tensor groupTensor(const Tensor &t, G &&...gs) {
  return groupTensorAux(t, t.rank(), std::forward<G>(gs)...);
}

template <typename... G>
static inline Tensor ungroupTensor(const Tensor &t, G &&...gs) {
  return ungroupTensorAux(t, unsigned(t.rank() - sizeof...(gs)),
                          std::forward<G>(gs)...);
}

namespace popops {
namespace rearrange {

static std::vector<Type> getValidTransposeDataTypes() {
  return {HALF, FLOAT, UNSIGNED_SHORT, UNSIGNED_INT, SHORT, INT};
}

bool canUseFastTranspose(const poplar::Target &target, const poplar::Type &type,
                         unsigned numRows, unsigned numColumns,
                         unsigned numTranspositions) {

  bool is2ByteType = (type == HALF || type == UNSIGNED_SHORT || type == SHORT);
  if (!is2ByteType ||
      numTranspositions > std::numeric_limits<unsigned short>::max() ||
      numRows % 4 || numColumns % 4) {
    return false;
  }
  // Check machine limits
  if (numColumns == 4 && numRows == 4) {
    if ((numTranspositions >= 2) &&
        (numTranspositions - 2 > target.getRptCountMax()))
      return false;
  } else if (numColumns == 4) {
    if (((numRows >= 8) && (numRows / 4 - 2 > target.getRptCountMax())) ||
        (numRows / 4u * 3u - 1u > (1u << (target.getNumStrideBits() - 1u)))) {
      return false;
    }
  } else {
    if (((numColumns >= 8) && (numColumns / 4 - 2 > target.getRptCountMax())) ||
        (numColumns / 4u * 3u - 1u >
         (1u << (target.getNumStrideBits() - 1u))) ||
        (numRows / 4u >= (1u << (target.getNumStrideBits() - 1)))) {
      return false;
    }
  }
  return true;
}

static bool switchToWorkerTranspose(const unsigned numTileTranspositions,
                                    const unsigned rows, const unsigned cols) {
  return (numTileTranspositions == 1) ||
         ((rows == 4) && (cols == 4) && (numTileTranspositions <= 4));
}

void addTransposeVertices(
    Graph &graph, const ComputeSet &cs, const Type &dType, unsigned rows,
    unsigned cols, const poplar::Graph::TileToTensorMapping &mapping,
    std::function<std::pair<const poplar::Tensor, const poplar::Tensor>(size_t)>
        getInOut) {
  const auto &validTypes = getValidTransposeDataTypes();
  if (std::find(validTypes.begin(), validTypes.end(), dType) ==
      validTypes.end()) {
    throw poplibs_error("Transposition not supported for data type " +
                        dType.toString());
  }
  if (cols > std::numeric_limits<unsigned short>::max() ||
      rows > std::numeric_limits<unsigned short>::max()) {
    throw poplibs_error("Number of source rows and columns exceed sizes "
                        "supported by Transpose/Transpose2d vertex");
  }
  // Shorthand local function to accumulate total size of a vector of Intervals
  auto accumSize = [](const std::vector<Interval> &vi) {
    return std::accumulate(
        vi.begin(), vi.end(), 0,
        [](size_t acc, const Interval &i) { return acc + i.size(); });
  };
  const auto &target = graph.getTarget();
  for (unsigned tile = 0; tile != mapping.size(); ++tile) {
    // All the transpositions to do on this tile. This is a vector of intervals,
    // each one specifying a set of transpositions.
    const auto &tileTranspositions = mapping[tile];

    // How many transpositions in all for this tile?
    unsigned numTileTranspositions = accumSize(tileTranspositions);
    if (numTileTranspositions > 0) {

      // There are 3 types of vertices that we might use. Default is Supervisor
      enum VertexType { TransposeSupervisor, Transpose, Transpose2d };
      std::map<VertexType, std::string> vertexNames = {
          {TransposeSupervisor, "popops::TransposeSupervisor"},
          {Transpose, "popops::Transpose"},
          {Transpose2d, "popops::Transpose2d"},
      };
      VertexType vertexType = TransposeSupervisor;
      // Will we end up splitting among workers (if not supervisor)?
      bool splitToWorkers = false;
      // Can we really use the Supervisor Vertex to do them all?
      if (canUseFastTranspose(target, dType, rows, cols,
                              numTileTranspositions)) {
        // If we have to do a single matrix (of any size), it's faster to run
        // the 'plain' Transpose instead of TransposeSupervisor.
        // Same is true if we have up to four 4x4 matrix
        if (switchToWorkerTranspose(numTileTranspositions, rows, cols)) {
          vertexType = Transpose;
        }
      } else {
        // Will need to partition to workers. vertexType will be chosen later
        splitToWorkers = true;
      }

      // Local function (as a lambda) to add a vertex for 'tile'.
      //     vType:          What kind of vertex to use
      //     tile:           where the vertex must be mapped
      //     transpositions: all transpositions this vertex has to do
      auto addOneVertex = [&](VertexType vType, unsigned tile,
                              std::vector<poplar::Interval> transpositions) {
        // Build inVec[], outVec[] to contain one element for each transposition
        std::vector<poplar::Tensor> inVec, outVec;
        for (const auto &interval : transpositions) {
          for (auto transposition = interval.begin();
               transposition != interval.end(); ++transposition) {
            poplar::Tensor in, out;
            std::tie(in, out) = getInOut(transposition);
            inVec.push_back(in);
            outVec.push_back(out);
            graph.setTileMapping(out, tile);
          }
        }
        std::string vertexName = vertexNames[vType];
        const auto v = graph.addVertex(cs, templateVertex(vertexName, dType));

        graph.setTileMapping(v, tile);
        if ((vType == Transpose) || (vType == TransposeSupervisor)) {
          graph.connect(v["src"], concat(inVec));
          graph.connect(v["dst"], concat(outVec));
          graph.setInitialValue(v["numSrcColumnsD4"], cols / 4);
          graph.setInitialValue(v["numSrcRowsD4"], rows / 4);
          if (vType == Transpose) {
            graph.setInitialValue(v["numTranspositionsM1"], inVec.size() - 1);
          } else {
            // We will run one supervisor vertex, starting the 6 workers.
            // The first 'workerCount' workers (1<=workerCount<=6) will
            // transpose 'numTranspositions' matrices and (6-workerCount)
            // workers transposing (numTranspositions-1) matrices.
            // Note that (6-workerCount) and/or (numTranspositions-1) might
            // be zero.
            // Note that this is NOT the same split as
            // splitRegionsBetweenWorkers() would do.
            unsigned numWorkerContexts = target.getNumWorkerContexts();
            unsigned workerCount = numWorkerContexts, numTranspositions = 1;
            if (numTileTranspositions <= numWorkerContexts) {
              workerCount = numTileTranspositions;
            } else {
              numTranspositions = numTileTranspositions / workerCount;
              unsigned rem = numTileTranspositions % workerCount;
              if (rem > 0) {
                workerCount = rem;
                numTranspositions += 1;
              }
            }
            graph.setInitialValue(v["numTranspositions"], numTranspositions);
            graph.setInitialValue(v["workerCount"], workerCount);
          }
        } else {
          graph.connect(v["src"], inVec);
          graph.connect(v["dst"], outVec);
          graph.setInitialValue(v["numSrcColumns"], cols);
          graph.setInitialValue(v["numSrcRows"], rows);
        }
      }; // addOneVertex()
      if (!splitToWorkers) {
        // A single vertex will do all the transpositions for this
        // tile
        addOneVertex(vertexType, tile, tileTranspositions);
      } else {
        // Need to split to multiple workers on this tile
        auto perWorkerTranspositions =
            splitRegionsBetweenWorkers(target, tileTranspositions, 1);
        for (const auto &transpositions : perWorkerTranspositions) {
          size_t size = accumSize(transpositions);
          vertexType = canUseFastTranspose(target, dType, rows, cols, size)
                           ? Transpose
                           : Transpose2d;
          addOneVertex(vertexType, tile, transpositions);
        } // for each worker
      }   // cannot use Supervisor variant
    }     // if (numTileTranspositions>0)
  }       // for each tile
}

Tensor partialTranspose(Graph &graph, const Tensor &in, const ComputeSet &cs,
                        const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext, DI_ARGS(in, cs));

  const auto rank = in.rank();
  const auto numSrcRows = in.dim(rank - 2);
  const auto numSrcColumns = in.dim(rank - 1);

  // Get a view on the 'in' tensor that is just a 2D matrix where each
  // row is one of the matrices (the 'rightmost' 2 dimensions) to transpose
  // ('inFlat').
  // I.e. flatten all the leftmost N-2 dimension together and also the last 2
  // dimensions together.
  // Get an equivalent view for the 'out' tensor ('outFlat').
  const auto dType = in.elementType();
  auto outShape = in.shape();
  std::swap(outShape[rank - 2], outShape[rank - 1]);
  auto out = graph.addVariable(dType, outShape, {di, "partialTranspose"});
  auto inFlat = in.reshape({in.numElements() / (numSrcRows * numSrcColumns),
                            numSrcRows * numSrcColumns});
  auto outFlat = out.reshape(inFlat.shape());

  // Get the tile mapping for the first element of each 2D matrix to transpose
  // (i.e. the rows of 'inFlat').
  // This tile is where we will allocate the Transpose vertices and the
  // output transposed matrix.
  const auto transpositionMapping = graph.getTileMapping(inFlat.slice(0, 1, 1));

  addTransposeVertices(graph, cs, dType, numSrcRows, numSrcColumns,
                       transpositionMapping, [&](size_t index) {
                         return std::make_pair(inFlat[index], outFlat[index]);
                       });
  di.addOutput(out);
  return out;
}

unsigned getMinimumRegroupGrainSize(const Type &type) {
  bool is2ByteType = (type == HALF || type == UNSIGNED_SHORT || type == SHORT);
  bool is4ByteType = (type == FLOAT || type == UNSIGNED_INT || type == INT);
  if (is2ByteType) {
    return 4;
  } else if (is4ByteType) {
    return 2;
  }
  return 1;
}

// Returns an updated grouping based on original grouping and tile mapping
static std::tuple<GroupingInfo, GroupingInfo, Graph::TileToTensorMapping>
updateGroupingInternal(const Graph &graph, const Tensor &t,
                       const GroupingInfo &from, const GroupingInfo &to) {
  auto grouped = groupTensor(t, to, from);
  auto groupedFlat = grouped.flatten(0, grouped.rank() - 2).flatten(1, 3);
  const auto tMapping = graph.getTileMapping(groupedFlat);
  const auto numTiles = tMapping.size();
  const auto tilesPerIPU = graph.getTarget().getTilesPerIPU();
  const auto numIPUs = numTiles / tilesPerIPU;
  std::vector<std::size_t> elemsPerIpu(numIPUs);
  for (unsigned tile = 0; tile != numTiles; ++tile) {
    const auto ipu = tile / tilesPerIPU;
    const auto &mapping = tMapping[tile];
    if (!mapping.empty()) {
      for (const auto &r : mapping) {
        elemsPerIpu[ipu] += r.size();
      }
    }
  }

  // Minimum number of elements in a group. Groups are split by a multiple of
  // this
  const unsigned minGroupsSize = getMinimumRegroupGrainSize(t.elementType());

  // find entry with max elements
  auto maxIt = std::max_element(std::begin(elemsPerIpu), std::end(elemsPerIpu));

  // A factor by which the number of transposes can increase by without
  // breaking the constraints set on group size
  auto additionalFactor = groupedFlat.dim(1) / (minGroupsSize * minGroupsSize);

  auto isPrime = [](unsigned num) {
    for (unsigned i = 2; i <= num / 2; ++i) {
      if (num % i == 0) {
        return false;
      }
    }
    return true;
  };

  // This limits the number of transpositions allowed on the IPU
  const auto maxTranspositionsAllowedPerIpu = numTiles;

  // actual transpose factor used. Initialise with 1 which means no additional
  // factor is applied
  unsigned transposeFactor = 1;

  // Estimate the number of transpositions on the IPU which has the maximum
  // number of elements mapped
  auto transpositionsOnIpuEstimate =
      (*maxIt + groupedFlat.dim(1) - 1) / groupedFlat.dim(1);

  bool allowIncrease =
      to.second % minGroupsSize == 0 && from.second % minGroupsSize == 0;
  while (allowIncrease && additionalFactor != 1) {
    unsigned factor = 1;
    // TODO: T12892 This assumes that typical transposes are a multiple of very
    // small primes. Investigate other methods (e.g., dividing into prime
    // factors). A method that should give good results is to find the maximum
    // GCD across different values of transpositions (i.e.
    // maxTranspositionsAllowedPerIpu, maxTranspositionsAllowedPerIpu-1, ...)
    for (unsigned x = 2; x <= additionalFactor; ++x) {
      if (additionalFactor % x == 0 && isPrime(x)) {
        factor = x;
        break;
      }
    }
    if (transpositionsOnIpuEstimate * transposeFactor * factor >
        maxTranspositionsAllowedPerIpu) {
      break;
    }
    if (additionalFactor % factor != 0 || factor == 1) {
      throw poputil::poplibs_error("Invalid factor in regrouping");
    }
    transposeFactor *= factor;
    additionalFactor /= factor;
  }

  auto updatedFrom = from;
  auto updatedTo = to;

  if (transposeFactor != 1) {
    // TODO: T12893 Optimise split once the cost of using a supervisor vertex
    // is known.
    auto factorFrom = gcd(transposeFactor, from.second / minGroupsSize);
    transposeFactor /= factorFrom;
    auto factorTo = gcd(transposeFactor, to.second / minGroupsSize);
    updatedFrom.second /= factorFrom;
    updatedTo.second /= factorTo;
  }
  return std::make_tuple(updatedFrom, updatedTo, std::move(tMapping));
}

std::pair<GroupingInfo, GroupingInfo> updateGrouping(const Graph &graph,
                                                     const Tensor &t,
                                                     const GroupingInfo &from,
                                                     const GroupingInfo &to) {
  const auto result = updateGroupingInternal(graph, t, from, to);
  return std::make_pair(std::get<0>(result), std::get<1>(result));
}

Tensor regroupTensorInternal(Graph &graph, const Tensor &t,
                             std::vector<Copy> &copies,
                             const ComputeSet &transposeCS,
                             const GroupingInfo &from_, const GroupingInfo &to_,
                             const DebugNameAndId &dnai) {
  const std::string fnStr = "internal";
  logging::popops::debug("Regroup: debugstr={}", dnai.getPathName());
  logging::popops::debug("  t      shape={}", t.shape());
  logging::popops::debug("  from   grouping={{{},{}}}", from_.first,
                         from_.second);
  logging::popops::debug("  to     grouping={{{},{}}}", to_.first, to_.second);

  if (t.rank() <= 1) {
    return t;
  }

  const auto &validTypes = getValidTransposeDataTypes();
  if (std::find(validTypes.begin(), validTypes.end(), t.elementType()) ==
      validTypes.end()) {
    throw poplibs_error("Data type " + t.elementType().toString() +
                        " not supported by regroupTensor");
  }
  if (from_.first >= t.rank()) {
    throw poplibs_error("Grouping 'from' specifies dimension " +
                        std::to_string(from_.first) +
                        " which is invalid for tensor to regroup of rank " +
                        std::to_string(t.rank()));
  }
  if (to_.first >= t.rank()) {
    throw poplibs_error("Grouping 'to' specifies dimension " +
                        std::to_string(to_.first) +
                        " which is invalid for tensor to regroup of rank " +
                        std::to_string(t.rank()));
  }
  if (t.dim(from_.first) % from_.second) {
    throw poplibs_error(
        "Size of grouping 'from' " + std::to_string(from_.second) +
        " does not evenly divide dimension " + std::to_string(from_.first) +
        " with size " + std::to_string(t.dim(from_.first)) +
        " of the tensor to regroup");
  }
  if (t.dim(to_.first) % to_.second) {
    throw poplibs_error("Size of grouping 'to' " + std::to_string(to_.second) +
                        " does not evenly divide dimension " +
                        std::to_string(to_.first) + " with size " +
                        std::to_string(t.dim(to_.first)) +
                        " of the tensor to regroup");
  }
  GroupingInfo to, from;
  Graph::TileToTensorMapping tMapping;
  std::tie(from, to, tMapping) = updateGroupingInternal(graph, t, from_, to_);
  auto grouped = groupTensor(t, to, from);
  auto groupedFlat = grouped.flatten(0, grouped.rank() - 2).flatten(1, 3);

  if (!(from == from_ && to == to_)) {
    tMapping = graph.getTileMapping(groupedFlat);
  }

  // Explicitly copy to a single variable in order to force
  // regions to be contiguous. Performing a transpose alone
  // may leave multiple regions per-tile, one for each edge to a
  // transpose vertex.
  auto preRegroup =
      graph.addVariable(t.elementType(), grouped.shape(), {dnai, "preRegroup"});
  auto preRegroupTranspose = preRegroup.flatten(0, preRegroup.rank() - 2);
  auto preRegroupFlat =
      preRegroup.flatten(0, preRegroup.rank() - 2).flatten(1, 3);

  // Build a map giving which intervals are mapped to each
  // IPU. Track which tiles on each IPU have any elements
  // mapped.
  const auto numTiles = tMapping.size();
  const auto tilesPerIPU = graph.getTarget().getTilesPerIPU();
  const auto numIPUs = numTiles / tilesPerIPU;

  using IntervalMap = boost::icl::interval_map<std::size_t, unsigned,
                                               boost::icl::partial_enricher>;
  using Interval = boost::icl::interval<std::size_t>;
  IntervalMap intervalsToIPU;
  for (unsigned tile = 0; tile < numTiles; ++tile) {
    auto ipu = tile / tilesPerIPU;
    const auto &mapping = tMapping[tile];
    if (!mapping.empty()) {
      for (const auto &i : mapping) {
        intervalsToIPU.insert(
            std::make_pair(Interval::right_open(i.begin(), i.end()), ipu));
      }
    }
  }

  // Iterate each transposition, mapping this to an IPU based on the first
  // element in each.
  auto elemsPerTransposition = preRegroupFlat.dim(1);
  std::vector<std::vector<poplar::Interval>> ipuTranspositions(numIPUs);
  for (unsigned t = 0; t < preRegroupFlat.dim(0); ++t) {
    auto it = intervalsToIPU.find(Interval::right_open(
        t * elemsPerTransposition, t * elemsPerTransposition + 1));
    assert(it != intervalsToIPU.end());
    auto ipu = it->second;
    auto &ipuTs = ipuTranspositions[ipu];
    // Try and extend the previous region if possible
    if (!ipuTs.empty() && ipuTs.back().end() == t) {
      ipuTs.back() = poplar::Interval(ipuTs.back().begin(), t + 1);
    } else {
      ipuTs.emplace_back(t, t + 1);
    }
  }

  // Finally map slices of the new tensor to transpose mapped linearly
  // across the tiles on which the original tensor was mapped on the same
  // IPU the elements of the transposition were originally mapped to.
  //
  // FIXME: This currently allows external exchange to be incurred for a
  // given transposition. This should not be allowed as it is not expected
  // but for the timebeing the padding constants added to activations
  // are just mapped to tile 0 which can be a different IPU to the one
  // on which it should actually reside. T6427 is required for this as a
  // primary usage of these regrouping functions.
  for (unsigned ipu = 0; ipu < numIPUs; ++ipu) {
    const auto &transpositions = ipuTranspositions[ipu];
    auto numTranspositions = std::accumulate(
        transpositions.begin(), transpositions.end(), std::size_t(0),
        [](std::size_t t, const poplar::Interval &i) { return t + i.size(); });
    if (!numTranspositions)
      continue;

    // spread across all tiles
    auto numTiles = tilesPerIPU;

    // Map transpositions on this IPU evenly across the tiles on which
    // elements of the source tensor reside.
    auto transpositionsPerTile = (numTranspositions + numTiles - 1) / numTiles;
    auto interval = transpositions.begin();
    unsigned intervalOffset = 0;
    for (unsigned i = 0; i < numTiles; ++i) {
      auto remaining = std::min(transpositionsPerTile, numTranspositions);
      numTranspositions -= remaining;
      while (remaining > 0) {
        auto n = std::min(interval->size() - intervalOffset, remaining);
        auto slice =
            preRegroupFlat.slice(interval->begin() + intervalOffset,
                                 interval->begin() + intervalOffset + n, 0);
        graph.setTileMapping(slice, i);
        remaining -= n;
        intervalOffset += n;
        if (interval->begin() + intervalOffset == interval->end()) {
          ++interval;
          intervalOffset = 0;
        }
      }
    }
  }

  copies.emplace_back(grouped, preRegroup, false, DebugContext(dnai, fnStr));

  // Finally, transpose
  auto partiallyTransposed = popops::rearrange::partialTranspose(
      graph, preRegroup, transposeCS, {dnai, fnStr});

  auto output = ungroupTensor(partiallyTransposed, from, to);
  return output;
}

Tensor regroupTensor(Graph &graph, const Tensor &t, std::vector<Copy> &copies,
                     const ComputeSet &transposeCS, const GroupingInfo &from_,
                     const GroupingInfo &to_,
                     const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext,
                                 DI_ARGS(t, copies, transposeCS, from_, to_));
  auto output =
      regroupTensorInternal(graph, t, copies, transposeCS, from_, to_, {di});

  di.addOutput(output);
  return output;
}

Tensor regroupTensor(Graph &graph, const Tensor &t_,
                     poplar::program::Sequence &prog,
                     const ComputeSet &transposeCS, const GroupingInfo &from_,
                     const GroupingInfo &to_,
                     const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext,
                                 DI_ARGS(t_, transposeCS, from_, to_));
  std::vector<Copy> copies;
  auto t = regroupTensor(graph, t_, copies, transposeCS, from_, to_, {di});

  for (const auto &copy : copies) {
    prog.add(copy);
  }
  di.addOutput(t);
  return t;
}

Tensor regroupIfPossible(Graph &graph, const Tensor &t_,
                         poplar::program::Sequence &prog,
                         const GroupingInfo &to_,
                         const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext, DI_ARGS(t_, to_));

  const auto groupings = detectDimGroupings(graph, t_);
  logging::popops::debug("Regroup if possible {}", debugContext.getPathName());
  logging::popops::debug("  shape {}, groupings {}", t_.shape(), groupings);

  const auto invalid = (t_.dim(to_.first) % to_.second) || (t_.rank() <= 1) ||
                       (to_.first >= t_.rank());
  // If the innermost grouping is already the dimension for which regrouping
  // is requested, nothing needs to be done.
  const auto noRegroupingRequired =
      !groupings.empty() && std::get<0>(groupings.at(0)) == std::get<0>(to_);

  const auto grainSize = getMinimumRegroupGrainSize(t_.elementType());

  // Only check for grouping in the innermost dimension. If the innermost
  // dimension is the `to be regrouped` dimension, noRegroupingRequired = true
  // and no regrouping is performed. If the innermost grouping has the
  // requisite grain size and is not the `to be regrouped` dimension we use it
  // up as the from dimension.
  auto suitableGroupingFound =
      !groupings.empty() && std::get<1>(groupings.at(0)) % grainSize == 0;

  if (invalid || noRegroupingRequired || !suitableGroupingFound) {
    logging::popops::debug("  Regrouping not possible: invalid ? {}, "
                           "noRegroupingRequired {}, suitableGroupingFound {}",
                           invalid, noRegroupingRequired,
                           suitableGroupingFound);
    di.addOutput(t_);
    return t_;
  }

  std::vector<Copy> copies;
  auto transposeCS =
      graph.addComputeSet(debugContext.getPathName() + "/regroupTensor");
  auto t = regroupTensorInternal(graph, t_, copies, transposeCS,
                                 groupings.at(0), to_, {di});
  for (const auto &copy : copies) {
    prog.add(copy);
  }
  prog.add(Execute(transposeCS));
  di.addOutput(t);
  return t;
}

Tensor regroupIfBeneficial(Graph &graph, const Tensor &in_, const Tensor &ref,
                           std::vector<Copy> &preTranspose,
                           ComputeSet transposeCS,
                           const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext,
                                 DI_ARGS(in_, ref, preTranspose, transposeCS));

  logging::popops::debug("Regroup if beneficial: debugstr={}",
                         debugContext.getPathName());
  logging::popops::debug("  input      shape={}", in_.shape());
  logging::popops::debug("  reference  shape={}", ref.shape());
  auto in = in_;

  if (in.rank() <= 1 || ref.rank() <= 1) {
    return in;
  }

  // If we can't transpose this data type, it's not beneficial to try
  const auto &transposableTypes = getValidTransposeDataTypes();
  if (std::find(transposableTypes.begin(), transposableTypes.end(),
                in.elementType()) == transposableTypes.end()) {
    return in;
  }

  if (in.shape() != ref.shape()) {
    throw poplibs_error("Input and reference tensors should be of "
                        "the same shape");
  }

  const auto inGrouping = detectDimGroupings(graph, in);
  const auto refGrouping = detectDimGroupings(graph, ref);
  logging::popops::debug("  input      groupings={}", inGrouping);
  logging::popops::debug("  reference  groupings={}", refGrouping);

  // TODO: T10360 Consider avoiding regrouping float inputs.
  auto grainSize = getMinimumRegroupGrainSize(in.elementType());
  if (!inGrouping.empty() && !refGrouping.empty() &&
      inGrouping[0].first != refGrouping[0].first &&
      (inGrouping[0].second % grainSize) == 0 &&
      (refGrouping[0].second % grainSize) == 0) {
    logging::popops::debug("  regrouped");
    in = regroupTensorInternal(graph, in, preTranspose, transposeCS,
                               inGrouping[0], refGrouping[0], {di});
  }
  di.addOutput(in);
  return in;
}

Tensor regroupIfBeneficial(Graph &graph, const Tensor &in_, const Tensor &ref,
                           Sequence &prog,
                           const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext, DI_ARGS(in_, ref));

  std::vector<Copy> preTranspose;
  ComputeSet transposeCS = graph.addComputeSet({di, "Transpose"});

  auto in =
      regroupIfBeneficial(graph, in_, ref, preTranspose, transposeCS, {di});

  for (const auto &copy : preTranspose) {
    prog.add(copy);
  }
  prog.add(Execute(transposeCS, {di}));
  di.addOutput(in);
  return in;
}

Tensor regroupIfBeneficial(Graph &graph, const Tensor &in_,
                           std::size_t preferredGrouping_, Sequence &prog,
                           const poplar::DebugContext &debugContext) {
  POPOPS_TRACEPOINT();
  poputil::PoplibsOpDebugInfo di(debugContext,
                                 DI_ARGS(in_, preferredGrouping_));
  logging::popops::debug("Regroup if beneficial (preferred): debugstr={}",
                         debugContext.getPathName());
  logging::popops::debug("  input        shape={}", in_.shape());
  logging::popops::debug("  preferred grouping={}", preferredGrouping_);
  auto in = in_;

  // If we can't transpose this data type, it's not beneficial to try
  const auto &transposableTypes = getValidTransposeDataTypes();
  if (std::find(transposableTypes.begin(), transposableTypes.end(),
                in.elementType()) == transposableTypes.end()) {
    return in;
  }

  if (in.dim(in.rank() - 1) % preferredGrouping_ != 0) {
    throw poplibs_error("Input tensor's innermost dimension is not "
                        "divisible by the given preferred grouping (" +
                        std::to_string(preferredGrouping_) + ")");
  }

  const auto inGrouping = detectDimGroupings(graph, in);
  logging::popops::debug("  input     grouping={}", inGrouping);
  const auto preferredGrouping =
      GroupingInfo{in.rank() - 1, preferredGrouping_};

  // TODO: T10360 Consider avoiding regrouping float inputs.
  auto grainSize = getMinimumRegroupGrainSize(in.elementType());
  logging::popops::trace(
      "beneficialRegroup decision: !empty {}, groupingChange "
      "{}, in matches grainsize {}, out matches grainsize {} "
      "(grainSize {})",
      !inGrouping.empty(),
      inGrouping.empty() ? 0 : inGrouping[0].first != preferredGrouping.first,
      inGrouping.empty() ? 0 : inGrouping[0].second % grainSize == 0,
      preferredGrouping.second % grainSize == 0, grainSize);
  if (!inGrouping.empty() && inGrouping[0].first != preferredGrouping.first &&
      inGrouping[0].second % grainSize == 0 &&
      preferredGrouping.second % grainSize == 0) {
    logging::popops::debug("  regrouped");
    ComputeSet transposeCS = graph.addComputeSet({di, "Transpose"});
    in = regroupTensor(graph, in, prog, transposeCS, inGrouping[0],
                       preferredGrouping, {di});
    prog.add(Execute(transposeCS, {di}));
  }

  return in;
}

} // end namespace rearrange
} // end namespace popops

max_stars_repo_name: graphcore/poplibs

max_issues_repo_name: graphcore/poplibs

max_forks_repo_name: graphcore/poplibs


text:
#include <actions/actions.h>
#include <dependency_graph/graph.h>
#include <dependency_graph/rtti.h>
#include <possumwood_sdk/app.h>

#include <boost/test/unit_test.hpp>
#include <dependency_graph/node_base.inl>
#include <dependency_graph/nodes_iterator.inl>
#include <dependency_graph/port.inl>
#include <vector>

#include "common.h"

using namespace dependency_graph;
using nlohmann::json;

namespace {

dependency_graph::Nodes::const_iterator findNode(const dependency_graph::Network& net, const std::string& name) {
	dependency_graph::Nodes::const_iterator result = net.nodes().end();

	for(auto it = net.nodes().begin(); it != net.nodes().end(); ++it)
		if(it->name() == name) {
			result = it;
			break;
		}

	return result;
}

}  // namespace

BOOST_AUTO_TEST_CASE(clipboard) {
	std::vector<::nlohmann::json> data;
	std::vector<std::function<void(const dependency_graph::Network&)>> tests;

	//////////////
	// a simple network of nodes
	data.push_back("{\"connections\":[],\"nodes\":{}, \"name\": \"network\", \"type\": \"network\"}"_json);

	data.back()["nodes"]["addition_0"] = nlohmann::json{
	    {"name", "add"},
	    {"type", "addition"},
	    {"ports", {{"input_1", 2.0}, {"input_2", 4.0}}},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}}};

	data.back()["nodes"]["multiplication_0"] = nlohmann::json{
	    {"name", "mult_0"},
	    {"type", "multiplication"},
	    {"ports", {{"input_2", 5.0}}},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}}};

	data.back()["nodes"]["multiplication_1"] = nlohmann::json{
	    {"name", "mult_1"},
	    {"type", "multiplication"},
	    {"ports", {{"input_1", 0.0}, {"input_2", 0.0}}},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}}};

	data.back()["connections"] = nlohmann::json{{
	    {"out_node", "addition_0"},
	    {"out_port", "output"},
	    {"in_node", "multiplication_0"},
	    {"in_port", "input_1"},
	}};

	tests.push_back([&](const dependency_graph::Network& net) {
		auto m0 = findNode(net, "mult_0");
		BOOST_REQUIRE(m0 != net.nodes().end());

		auto m1 = findNode(net, "mult_1");
		BOOST_REQUIRE(m1 != net.nodes().end());

		auto a0 = findNode(net, "add");
		BOOST_REQUIRE(a0 != net.nodes().end());

		BOOST_CHECK_EQUAL(a0->port(0).get<float>(), 2.0f);
		BOOST_CHECK_EQUAL(a0->port(1).get<float>(), 4.0f);
		BOOST_CHECK_EQUAL(a0->port(2).get<float>(), 6.0f);

		BOOST_CHECK_EQUAL(m0->port(0).get<float>(), 6.0f);
		BOOST_CHECK_EQUAL(m0->port(1).get<float>(), 5.0f);
		BOOST_CHECK_EQUAL(m0->port(2).get<float>(), 30.0f);
	});

	///////////////
	// previous network as a subnetwork (with no inputs or outputs)
	data.push_back(nlohmann::json());
	data.back()["nodes"]["network_0"] =
	    json{{"name", "test_network"},
	         {"type", "network"},
	         // {"ports", {}},
	         {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}},
	         {"nodes", data[0]["nodes"]},
	         {"connections", data[0]["connections"]}};
	data.back()["connections"] = "[]"_json;
	data.back()["name"] = "network";
	data.back()["type"] = "network";

	  tests.push_back([&](const dependency_graph::Network& net) {
		  auto s = findNode(net, "test_network");
		  BOOST_REQUIRE(s != net.nodes().end());

		  BOOST_REQUIRE(s->is<dependency_graph::Network>());
		  auto& subnet = s->as<dependency_graph::Network>();

		  auto m0 = findNode(subnet, "mult_0");
		  BOOST_REQUIRE(m0 != subnet.nodes().end());

		  auto m1 = findNode(subnet, "mult_1");
		  BOOST_REQUIRE(m1 != subnet.nodes().end());

		  auto a0 = findNode(subnet, "add");
		  BOOST_REQUIRE(a0 != subnet.nodes().end());

		  BOOST_CHECK_EQUAL(a0->port(0).get<float>(), 2.0f);
		  BOOST_CHECK_EQUAL(a0->port(1).get<float>(), 4.0f);
		  BOOST_CHECK_EQUAL(a0->port(2).get<float>(), 6.0f);

		  BOOST_CHECK_EQUAL(m0->port(0).get<float>(), 6.0f);
		  BOOST_CHECK_EQUAL(m0->port(1).get<float>(), 5.0f);
		  BOOST_CHECK_EQUAL(m0->port(2).get<float>(), 30.0f);
	  });

	///////////////
	// previous network as a subnetwork (with inputs and outputs, but not connected)
	//   -> contains unconnected void ports
	data.push_back(data.back());
	data.back()["nodes"]["network_0"]["nodes"]["input_0"] = json{
	    {"name", "input_1"},
	    {"type", "input"},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}},
	};
	data.back()["nodes"]["network_0"]["nodes"]["input_1"] = json{
	    {"name", "input_2"},
	    {"type", "input"},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}},
	};
	data.back()["nodes"]["network_0"]["nodes"]["output_0"] = json{
	    {"name", "output_1"},
	    {"type", "output"},
	    {"blind_data", {{"type", unmangledTypeId<possumwood::NodeData>()}, {"value", "test blind data"}}},
	};

	tests.push_back([&](const dependency_graph::Network&) {

	});

	////////////
	// previous network, but with connected inputs and output
	data.push_back(data.back());

	// add network port connections
	data.back()["nodes"]["network_0"]["connections"] = json{{
	                                                            {"out_node", "addition_0"},
	                                                            {"out_port", "output"},
	                                                            {"in_node", "output_0"},
	                                                            {"in_port", "data"},
	                                                        },
	                                                        {
	                                                            {"out_node", "addition_0"},
	                                                            {"out_port", "output"},
	                                                            {"in_node", "multiplication_0"},
	                                                            {"in_port", "input_1"},
	                                                        },
	                                                        {
	                                                            {"out_node", "input_0"},
	                                                            {"out_port", "data"},
	                                                            {"in_node", "addition_0"},
	                                                            {"in_port", "input_1"},
	                                                        },
	                                                        {
	                                                            {"out_node", "input_1"},
	                                                            {"out_port", "data"},
	                                                            {"in_node", "addition_0"},
	                                                            {"in_port", "input_2"},
	                                                        }};

	// connected ports will no longer be explicitly listed in addition_0 ports
	data.back()["nodes"]["network_0"]["nodes"]["addition_0"].erase(
	    data.back()["nodes"]["network_0"]["nodes"]["addition_0"].find("ports"));
	// but they will be exposed out on the network
	data.back()["nodes"]["network_0"]["ports"] = json{
	    {"input_1", 7},
	    {"input_2", 8},
	};

	tests.push_back([&](const dependency_graph::Network& net) {
		auto s = findNode(net, "test_network");
		BOOST_REQUIRE(s != net.nodes().end());

		BOOST_REQUIRE(s->is<dependency_graph::Network>());
		auto& subnet = s->as<dependency_graph::Network>();

		auto m0 = findNode(subnet, "mult_0");
		BOOST_REQUIRE(m0 != subnet.nodes().end());

		auto m1 = findNode(subnet, "mult_1");
		BOOST_REQUIRE(m1 != subnet.nodes().end());

		auto a0 = findNode(subnet, "add");
		BOOST_REQUIRE(a0 != subnet.nodes().end());

		BOOST_CHECK_EQUAL(a0->port(0).get<float>(), 7.0f);
		BOOST_CHECK_EQUAL(a0->port(1).get<float>(), 8.0f);
		BOOST_CHECK_EQUAL(a0->port(2).get<float>(), 15.0f);

		BOOST_CHECK_EQUAL(m0->port(0).get<float>(), 15.0f);
		BOOST_CHECK_EQUAL(m0->port(1).get<float>(), 5.0f);
		BOOST_CHECK_EQUAL(m0->port(2).get<float>(), 75.0f);

		BOOST_CHECK_EQUAL(s->port(0).get<float>(), 7.0f);
		BOOST_CHECK_EQUAL(s->port(1).get<float>(), 8.0f);
		BOOST_CHECK_EQUAL(s->port(2).get<float>(), 15.0f);

		// and test pull from outside
		BOOST_REQUIRE_NO_THROW(s->port(0).set(9.0f));

		BOOST_CHECK_EQUAL(s->port(0).get<float>(), 9.0f);
		BOOST_CHECK_EQUAL(s->port(2).get<float>(), 17.0f);
	});

	/////////////////////////////////////////////////////////////////////////////////////////

	// make sure the static handles are initialised
	additionNode();
	multiplicationNode();

	auto test_it = tests.begin();

	for(const auto& pasted_data : data) {
		// make the app "singleton"
		possumwood::AppCore app;

		// test that the undo stack is empty
		BOOST_REQUIRE(app.undoStack().empty());

		// paste the data
		dependency_graph::Selection selection;
		BOOST_REQUIRE_NO_THROW(possumwood::actions::fromJson(app.graph(), selection, pasted_data));

		// check that theres one undo step now
		BOOST_CHECK_EQUAL(app.undoStack().undoActionCount(), 1u);
		BOOST_CHECK_EQUAL(app.undoStack().redoActionCount(), 0u);

		// check that the pasted data are correct
		{
			::json json;
			BOOST_REQUIRE_NO_THROW(json = possumwood::actions::toJson());
			BOOST_CHECK_EQUAL(json, pasted_data);
		}

		// run the appropriate test
		(*test_it)(app.graph());

		// // perform undo
		// BOOST_REQUIRE_NO_THROW(app.undoStack().undo());

		// // graph should be empty now
		// BOOST_REQUIRE(app.graph().nodes().empty());

		// // check that theres one redo step now
		// BOOST_CHECK_EQUAL(app.undoStack().undoActionCount(), 0u);
		// BOOST_CHECK_EQUAL(app.undoStack().redoActionCount(), 1u);

		// // perform redo
		// BOOST_REQUIRE_NO_THROW(app.undoStack().redo());

		// // and check that the result is right
		// {
		// 	::json json;
		// 	BOOST_REQUIRE_NO_THROW(json = possumwood::actions::toJson());
		// 	BOOST_CHECK_EQUAL(json, pasted_data);
		// }

		++test_it;
	}
}

ext: cpp

max_stars_repo_name: martin-pr/possumwood

max_issues_repo_name: martin-pr/possumwood

max_forks_repo_name: martin-pr/possumwood


text:
/*
 *  Copyright (C) 2018 Naomasa Matsubayashi
 *  Licensed under MIT license, see file LICENSE in this source tree.
 */
#ifndef UWG_PEER_IDENTIFIER_HPP
#define UWG_PEER_IDENTIFIER_HPP

#include <cstdint>
#include <string>
#include <boost/asio.hpp>

namespace uwg {
  using kxid_t = uint32_t;
  struct peer_identifier {
    peer_identifier() : kxid( 0 ) {}
    peer_identifier( const boost::asio::ip::udp::endpoint &e, kxid_t i ) : endpoint( e ), kxid( i ) {}
    boost::asio::ip::udp::endpoint endpoint;
    kxid_t kxid;
  };

  bool operator==( const peer_identifier &l, const peer_identifier &r ) {
    return l.endpoint == r.endpoint && l.kxid == r.kxid;
  }
  bool operator!=( const peer_identifier &l, const peer_identifier &r ) {
    return l.endpoint != r.endpoint || l.kxid == r.kxid;
  }
  bool operator<( const peer_identifier &l, const peer_identifier &r ) {
    if( l.endpoint == r.endpoint ) return l.kxid < r.kxid;
    else return l.endpoint < r.endpoint;
  }
  bool operator>( const peer_identifier &l, const peer_identifier &r ) {
    if( l.endpoint == r.endpoint ) return l.kxid > r.kxid;
    else return l.endpoint > r.endpoint;
  }
  bool operator<=( const peer_identifier &l, const peer_identifier &r ) {
    if( l.endpoint == r.endpoint ) return l.kxid <= r.kxid;
    else return l.endpoint < r.endpoint;
  }
  bool operator>=( const peer_identifier &l, const peer_identifier &r ) {
    if( l.endpoint == r.endpoint ) return l.kxid >= r.kxid;
    else return l.endpoint < r.endpoint;
  }
}

#endif

ext: hpp


text:
/*
 * This file is part of hipSYCL, a SYCL implementation based on CUDA/HIP
 *
 * Copyright (c) 2018 Aksel Alpay
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "clang/Basic/SourceManager.h"
#include "clang/Lex/HeaderSearch.h"

#include <boost/filesystem.hpp>

#include "CL/sycl/detail/debug.hpp"

#include <iostream>
namespace hipsycl {
namespace transform {

class RewriteSelectorWhitelist
{
public:
  static RewriteSelectorWhitelist&
  get()
  {
    static RewriteSelectorWhitelist list;
    return list;
  }

  void addToWhitelist(const std::string& path)
  {
    whiteList.push_back(path);
  }

  bool isPathAccepted(const std::string& path) const
  {
    if(!isActivated)
      return true;

    for(const auto& dir : whiteList)
    {
      if(pathContainsFile(dir, path))
        return true;
    }

    return false;
  }

  void enable(bool enabled = true)
  {
    this->isActivated = enabled;
  }

  bool isEnabled() const
  {
    return isActivated;
  }
private:
  bool pathContainsFile(const std::string& dir,
                        const std::string& filename) const
  {
    auto checkIfFilepathStartsWithDir =
      [](std::string p, const std::string& f) -> bool{

      if(p.empty())
        return true;

      if(f.empty())
        return true;

      if(p.back() != '/')
        p += '/';

      return f.find(p) == 0;
    };

    try
    {
      auto canonicalDir = boost::filesystem::canonical(dir);
      auto canonicalFile = boost::filesystem::canonical(filename);

      return checkIfFilepathStartsWithDir(canonicalDir.string(),
                                          canonicalFile.string());
    }
    catch(...)
    {
      // Happens if dir or filename do not exist -- in this case
      // we just check if the filename starts with dir
      return checkIfFilepathStartsWithDir(dir, filename);
    }
  }

  RewriteSelectorWhitelist()
  : isActivated{false}
  {}

  bool isActivated;
  std::vector<std::string> whiteList;
};

class RewriteSelector
{
public:

  bool shouldIncludeBeRewritten(
      const clang::SourceManager& mgr,
      clang::FileID includedId,
      clang::SrcMgr::CharacteristicKind FileType,
      const clang::DirectoryLookup *DirLookup) const
  {
    // Never rewrite the sycl.hpp include file
    if(mgr.getFileEntryForID(includedId)->getName().endswith("CL/sycl.hpp"))
      return false;
#if defined(WIN32) || defined(_WIN32) || defined(__WIN32) && !defined(__CYGWIN__)
    if(mgr.getFileEntryForID(includedId)->tryGetRealPathName().endswith("CL\\sycl.hpp"))
      return false;
#endif

    auto headerName = mgr.getFileEntryForID(includedId)->getName();
    if(boost::filesystem::path{headerName}.is_relative())
    {
      HIPSYCL_DEBUG_WARNING
          << "hipsycl_rewrite_includes: Encountered "
          << "relative include path " << std::string{headerName}
          << ", this may lead to incorrect inlining of includes "
          << "if --rewrite-include-paths is used." << std::endl;
    }

    // Don't rewrite the standard library headers
    if(FileType == clang::SrcMgr::CharacteristicKind::C_System ||
       FileType == clang::SrcMgr::CharacteristicKind::C_ExternCSystem)
    {
      if(endsWithStdHeader(headerName))
        return false;
    }

    if(RewriteSelectorWhitelist::get().isPathAccepted(headerName))
      return true;

    return false;
  }

private:
  bool endsWithStdHeader(const llvm::StringRef& path) const
  {
    for(const std::string& header : getStdHeaders())
    {
      if(path.endswith(header))
        return true;
#if defined(WIN32) || defined(_WIN32) || defined(__WIN32) && !defined(__CYGWIN__)
      std::string winHeader = header;
      for(std::size_t i = 0; i winHeader.size(); ++i)
        if(winHeader[i] == '/')
          winHeader[i] = '\\';
      if(path.endswith(winHeader))
        return true;
#endif

    }
    return false;
  }

  static
  const std::vector<std::string>& getStdHeaders()
  {
    static std::vector<std::string> headers{
      "assert.h",
      "complex.h",
      "ctype.h",
      "errno.h",
      "fenv.h",
      "float.h",
      "inttypes.h",
      "iso646.h",
      "limits.h",
      "locale.h",
      "math.h",
      "setjmp.h",
      "signal.h",
      "stdalign.h",
      "stdarg.h",
      "stdatomic.h",
      "stdbool.h",
      "stddef.h",
      "stdint.h",
      "stdio.h",
      "stdlib.h",
      "stdnoreturn.h",
      "string.h",
      "tgmath.h",
      "threads.h",
      "time.h",
      "uchar.h",
      "wchar.h",
      "wctype.h",
      "iostream",
      "unordered_set",
      "type_traits",
      "cstddef",
      "cerrno",
      "queue",
      "stack",
      "streambuf",
      "fenv.h",
      "cstdint",
      "future",

      "backward/hash_set",
      "backward/auto_ptr.h",
      "backward/backward_warning.h",
      "backward/strstream",
      "backward/hashtable.h",
      "backward/hash_fun.h",
      "backward/binders.h",
      "backward/hash_map",
      "iomanip",
      "memory",

      "tr1/unordered_set",
      "tr1/type_traits",
      "tr1/fenv.h",
      "tr1/cstdint",
      "tr1/memory",
      "tr1/cstdarg",
      "tr1/utility",
      "tr1/shared_ptr.h",
      "tr1/beta_function.tcc",
      "tr1/tuple",
      "tr1/cwchar",
      "tr1/stdint.h",
      "tr1/ell_integral.tcc",
      "tr1/random",
      "tr1/cfloat",
      "tr1/hashtable_policy.h",
      "tr1/unordered_set.h",
      "tr1/wctype.h",
      "tr1/cinttypes",
      "tr1/hashtable.h",
      "tr1/exp_integral.tcc",
      "tr1/complex",
      "tr1/complex.h",
      "tr1/stdlib.h",
      "tr1/stdbool.h",
      "tr1/unordered_map.h",
      "tr1/cmath",
      "tr1/cwctype",
      "tr1/ctype.h",
      "tr1/cfenv",
      "tr1/poly_hermite.tcc",
      "tr1/array",
      "tr1/random.h",
      "tr1/float.h",
      "tr1/bessel_function.tcc",
      "tr1/legendre_function.tcc",
      "tr1/stdarg.h",
      "tr1/ctime",
      "tr1/random.tcc",
      "tr1/ccomplex",
      "tr1/functional_hash.h",
      "tr1/tgmath.h",
      "tr1/regex",
      "tr1/functional",
      "tr1/riemann_zeta.tcc",
      "tr1/inttypes.h",
      "tr1/cctype",
      "tr1/cstdbool",
      "tr1/limits.h",
      "tr1/cstdio",
      "tr1/special_function_util.h",
      "tr1/ctgmath",
      "tr1/climits",
      "tr1/hypergeometric.tcc",
      "tr1/unordered_map",
      "tr1/gamma.tcc",
      "tr1/modified_bessel_func.tcc",
      "tr1/cstdlib",
      "tr1/poly_laguerre.tcc",
      "tr1/math.h",
      "tr1/stdio.h",
      "tr1/wchar.h",
      "deque",
      "cstdarg",
      "list",
      "utility",
      "tuple",
      "scoped_allocator",
      "charconv",
      "cwchar",
      "iterator",
      "random",
      "tr2",
      "tr2/type_traits",
      "tr2/bool_set",
      "tr2/dynamic_bitset",
      "tr2/bool_set.tcc",
      "tr2/ratio",
      "tr2/dynamic_bitset.tcc",
      "cstdalign",
      "cfloat",
      "codecvt",
      "decimal",
      "decimal/decimal",
      "decimal/decimal.h",
      "thread",
      "new",
      "fstream",
      "cinttypes",
      "string_view",
      "bitset",
      "ostream",
      "locale",
      "stdexcept",
      "atomic",
      "istream",
      "chrono",
      "complex",
      "complex.h",
      "stdlib.h",
      "ciso646",
      "filesystem",
      "ios",
      "cxxabi.h",
      "cmath",
      "ratio",
      "cwctype",
      "algorithm",
      "cfenv",
      "variant",
      "any",
      "valarray",
      "cassert",
      "exception",
      "set",
      "string",
      "limits",

      "debug/unordered_set",
      "debug/safe_local_iterator.tcc",
      "debug/stl_iterator.h",
      "debug/deque",
      "debug/list",
      "debug/macros.h",
      "debug/safe_sequence.tcc",
      "debug/safe_iterator.h",
      "debug/formatter.h",
      "debug/safe_unordered_container.tcc",
      "debug/bitset",
      "debug/multimap.h",
      "debug/safe_local_iterator.h",
      "debug/multiset.h",
      "debug/debug.h",
      "debug/functions.h",
      "debug/safe_iterator.tcc",
      "debug/map.h",
      "debug/set",
      "debug/string",
      "debug/array",
      "debug/map",
      "debug/assertions.h",
      "debug/set.h",
      "debug/safe_sequence.h",
      "debug/safe_base.h",
      "debug/safe_container.h",
      "debug/safe_unordered_container.h",
      "debug/vector",
      "debug/helper_functions.h",
      "debug/unordered_map",
      "debug/forward_list",
      "debug/safe_unordered_base.h",
      "numeric",
      "array",
      "cuchar",

      "bits/stl_multimap.h",
      "bits/locale_facets.tcc",
      "bits/regex_automaton.tcc",
      "bits/stl_construct.h",
      "bits/valarray_array.tcc",
      "bits/fs_path.h",
      "bits/move.h",
      "bits/parse_numbers.h",
      "bits/stl_iterator.h",
      "bits/quoted_string.h",
      "bits/stl_vector.h",
      "bits/valarray_array.h",
      "bits/regex_automaton.h",
      "bits/node_handle.h",
      "bits/shared_ptr.h",
      "bits/cpp_type_traits.h",
      "bits/stl_tempbuf.h",
      "bits/deque.tcc",
      "bits/algorithmfwd.h",
      "bits/fs_fwd.h",
      "bits/regex_scanner.tcc",
      "bits/std_function.h",
      "bits/slice_array.h",
      "bits/hashtable_policy.h",
      "bits/stl_algo.h",
      "bits/fs_ops.h",
      "bits/stl_queue.h",
      "bits/regex_constants.h",
      "bits/unordered_set.h",
      "bits/stl_set.h",
      "bits/invoke.h",
      "bits/stl_multiset.h",
      "bits/basic_ios.h",
      "bits/fs_dir.h",
      "bits/regex_executor.h",
      "bits/postypes.h",
      "bits/atomic_base.h",
      "bits/shared_ptr_atomic.h",
      "bits/ptr_traits.h",
      "bits/stringfwd.h",
      "bits/stl_algobase.h",
      "bits/istream.tcc",
      "bits/hashtable.h",
      "bits/atomic_lockfree_defines.h",
      "bits/exception.h",
      "bits/basic_string.tcc",
      "bits/unique_ptr.h",
      "bits/list.tcc",
      "bits/ios_base.h",
      "bits/regex_compiler.h",
      "bits/stl_list.h",
      "bits/std_mutex.h",
      "bits/stl_deque.h",
      "bits/cxxabi_forced.h",
      "bits/nested_exception.h",
      "bits/unordered_map.h",
      "bits/stl_heap.h",
      "bits/locale_classes.h",
      "bits/allocated_ptr.h",
      "bits/localefwd.h",
      "bits/regex_scanner.h",
      "bits/indirect_array.h",
      "bits/specfun.h",
      "bits/fstream.tcc",
      "bits/streambuf.tcc",
      "bits/predefined_ops.h",
      "bits/locale_facets_nonio.tcc",
      "bits/stl_iterator_base_funcs.h",
      "bits/atomic_futex.h",
      "bits/valarray_before.h",
      "bits/regex_executor.tcc",
      "bits/regex_error.h",
      "bits/random.h",
      "bits/stl_tree.h",
      "bits/allocator.h",
      "bits/stl_stack.h",
      "bits/concept_check.h",
      "bits/regex.h",
      "bits/stl_iterator_base_types.h",
      "bits/mask_array.h",
      "bits/locale_facets_nonio.h",
      "bits/vector.tcc",
      "bits/alloc_traits.h",
      "bits/stl_map.h",
      "bits/random.tcc",
      "bits/regex.tcc",
      "bits/range_access.h",
      "bits/functional_hash.h",
      "bits/locale_classes.tcc",
      "bits/locale_conv.h",
      "bits/enable_special_members.h",
      "bits/uniform_int_dist.h",
      "bits/basic_string.h",
      "bits/hash_bytes.h",
      "bits/exception_ptr.h",
      "bits/streambuf_iterator.h",
      "bits/stl_raw_storage_iter.h",
      "bits/std_abs.h",
      "bits/stl_numeric.h",
      "bits/stl_uninitialized.h",
      "bits/stl_relops.h",
      "bits/regex_compiler.tcc",
      "bits/functexcept.h",
      "bits/locale_facets.h",
      "bits/c++0x_warning.h",
      "bits/ostream.tcc",
      "bits/stream_iterator.h",
      "bits/boost_concept_check.h",
      "bits/forward_list.h",
      "bits/string_view.tcc",
      "bits/forward_list.tcc",
      "bits/stl_pair.h",
      "bits/memoryfwd.h",
      "bits/exception_defines.h",
      "bits/basic_ios.tcc",
      "bits/ostream_insert.h",
      "bits/stl_bvector.h",
      "bits/gslice.h",
      "bits/refwrap.h",
      "bits/shared_ptr_base.h",
      "bits/sstream.tcc",
      "bits/uses_allocator.h",
      "bits/codecvt.h",
      "bits/stl_function.h",
      "bits/cxxabi_init_exception.h",
      "bits/gslice_array.h",
      "bits/valarray_after.h",
      "bits/char_traits.h",
      "iosfwd",
      "map",
      "typeindex",
      "ctime",
      "optional",
      "ccomplex",
      "mutex",
      "initializer_list",
      "vector",
      "tgmath.h",
      "regex",
      "functional",
      "clocale",
      "system_error",
      "cctype",
      "cstdbool",
      "csignal",
      "parallel",
      "parallel/compatibility.h",
      "parallel/random_shuffle.h",
      "parallel/par_loop.h",
      "parallel/tags.h",
      "parallel/algorithmfwd.h",
      "parallel/base.h",
      "parallel/merge.h",
      "parallel/unique_copy.h",
      "parallel/multiway_mergesort.h",
      "parallel/find_selectors.h",
      "parallel/losertree.h",
      "parallel/omp_loop.h",
      "parallel/search.h",
      "parallel/equally_split.h",
      "parallel/types.h",
      "parallel/balanced_quicksort.h",
      "parallel/multiway_merge.h",
      "parallel/algorithm",
      "parallel/settings.h",
      "parallel/numeric",
      "parallel/set_operations.h",
      "parallel/workstealing.h",
      "parallel/algobase.h",
      "parallel/multiseq_selection.h",
      "parallel/for_each.h",
      "parallel/partial_sum.h",
      "parallel/quicksort.h",
      "parallel/omp_loop_static.h",
      "parallel/sort.h",
      "parallel/features.h",
      "parallel/iterator.h",
      "parallel/partition.h",
      "parallel/random_number.h",
      "parallel/basic_iterator.h",
      "parallel/for_each_selectors.h",
      "parallel/queue.h",
      "parallel/find.h",
      "parallel/numericfwd.h",
      "parallel/algo.h",
      "parallel/checkers.h",
      "parallel/list_partition.h",
      "parallel/compiletime_settings.h",
      "parallel/parallel.h",
      "cstdio",
      "ctgmath",
      "climits",
      "unordered_map",
      "forward_list",
      "sstream",
      "cstring",
      "cstdlib",
      "math.h",
      "experimental",
      "experimental/unordered_set",
      "experimental/type_traits",
      "experimental/source_location",
      "experimental/memory",
      "experimental/deque",
      "experimental/list",
      "experimental/utility",
      "experimental/tuple",
      "experimental/iterator",
      "experimental/random",
      "experimental/string_view",
      "experimental/chrono",
      "experimental/memory_resource",
      "experimental/filesystem",
      "experimental/ratio",
      "experimental/algorithm",
      "experimental/propagate_const",
      "experimental/any",
      "experimental/set",
      "experimental/string",
      "experimental/numeric",
      "experimental/array",
      "experimental/bits",
      "experimental/bits/fs_path.h",
      "experimental/bits/shared_ptr.h",
      "experimental/bits/fs_fwd.h",
      "experimental/bits/fs_ops.h",
      "experimental/bits/fs_dir.h",
      "experimental/bits/erase_if.h",
      "experimental/bits/lfts_config.h",
      "experimental/bits/string_view.tcc",
      "experimental/map",
      "experimental/optional",
      "experimental/vector",
      "experimental/regex",
      "experimental/functional",
      "experimental/system_error",
      "experimental/unordered_map",
      "experimental/forward_list",
      "shared_mutex",
      "typeinfo",
      "csetjmp",
      "condition_variable"
    };

    return headers;
  }
};

}
}

max_stars_repo_name: nazavode/hipSYCL

max_issues_repo_name: nazavode/hipSYCL

max_forks_repo_name: nazavode/hipSYCL


text:
// Copyright Louis Dionne 2013-2016
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)

#include <boost/hana/assert.hpp>
#include <boost/hana/equal.hpp>
#include <boost/hana/insert.hpp>
#include <boost/hana/integral_constant.hpp>
#include <boost/hana/tuple.hpp>

#include <string>
namespace hana = boost::hana;
using namespace hana::literals;
using namespace std::literals;

int main() {
    auto xs = hana::make_tuple("Hello"s, "world!"s);

    BOOST_HANA_RUNTIME_CHECK(
        hana::insert(xs, 1_c, " "s) == hana::make_tuple("Hello"s, " "s, "world!"s)
    );
}

max_stars_repo_name: nxplatform/nx-mobile

max_issues_repo_name: nxplatform/nx-mobile

max_forks_repo_name: nxplatform/nx-mobile


text:
#include <functional>
#include <iostream>
#include <sstream>
#include <stdlib.h>
#include <string.h>

#include <boost/random.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>

#include "biggles/detail/random.hpp"
#include "biggles/sampling/metropolis_hastings.hpp"

extern "C" {
#include <ccan/tap/tap.h>
}

typedef boost::tuples::tuple<float, float> point;

inline float rosenbrock(const point& p)
{
    float x, y;
    boost::tuples::tie(x,y) = p;
    return (1.f-x)*(1.f-x) + 100.f*(y-x*x)*(y-x*x);
}

// a cooked-up PDF which has a maximum where the Rosenbrock function has a minimum.
inline float rosenbrock_log_pdf(const point& p)
{
    return -rosenbrock(p);
}

typedef boost::tuples::tuple<point, float> proposal_result_t;

inline proposal_result_t propose_gaussian(const point& p)
{
    static boost::normal_distribution<float> norm(0.f, 0.2f);
    static boost::variate_generator<boost::mt19937&, boost::normal_distribution<float> > variate(
        biggles::detail::random_generator, norm);

    point new_p(p.get<0>() + variate(), p.get<1>() + variate());

    return proposal_result_t(new_p, 0.f);
}

using namespace biggles::sampling::metropolis_hastings;

int main(int argc, char** argv)
{
    biggles::detail::seed_prng(0xdeadbeef);

    bool verbose = (argc > 1) && (0 == strcmp(argv[1], "-v"));

    plan_tests(4);

    typedef std::pointer_to_unary_function<const point&, float> pdf_t;
    typedef std::pointer_to_unary_function<const point&, proposal_result_t> proposal_func_t;

    // MH sampler for Rosenbrock function. An example of using a function pointer directly.
    sampler<pdf_t, proposal_func_t> s(std::ptr_fun(rosenbrock_log_pdf), std::ptr_fun(propose_gaussian));

    point optimum(s.last_sample());
    float optimum_log_density(s.current_sample_log_density());

    for(int i=0; i<4096; ++i)
    {
        s.draw();

        if(s.current_sample_log_density() > optimum_log_density)
        {
            optimum_log_density = s.current_sample_log_density();
            optimum = s.last_sample();
        }

        if(verbose)
        {
            std::stringstream ss;
            ss << "sample: " << s.last_sample();
            diag("%s", ss.str().c_str());
        }
    }

    {
        std::stringstream ss;
        ss << "optimum: " << optimum << ", log-pdf: " << optimum_log_density
        << ", alpha: " << s.acceptance_rate() << ", rosenbrock: " << rosenbrock(optimum);
        diag("%s", ss.str().c_str());
    }

    diag("acceptance rate = %.2f%%", s.acceptance_rate()*100.f);
    //ok(fabs(s.acceptance_rate() - 0.25f) < 0.2f, "acceptance rate within 0.2 of 0.25");
    ok(fabs(optimum.get<0>() - 1.f) < 0.1f, "x-co-ordinate within 10%% of optimum");
    ok(fabs(optimum.get<1>() - 1.f) < 0.1f, "y-co-ordinate within 10%% of optimum");

    // the same thing using the optimise function
    point optimum2 = optimise(std::ptr_fun(rosenbrock_log_pdf), std::ptr_fun(propose_gaussian));

    {
        std::stringstream ss;
        ss << "optimum via optimise(): " << optimum2 << ", log-pdf: " << rosenbrock_log_pdf(optimum2)
        << ", rosenbrock: " << rosenbrock(optimum2);
        diag("%s", ss.str().c_str());
    }

    ok(fabs(optimum2.get<0>() - 1.f) < 0.1f, "x-co-ordinate within 10%% of optimum");
    ok(fabs(optimum2.get<1>() - 1.f) < 0.1f, "y-co-ordinate within 10%% of optimum");

    return exit_status();
}

max_stars_repo_name: fbi-octopus/biggles

max_issues_repo_name: fbi-octopus/biggles

max_forks_repo_name: fbi-octopus/biggles


text:
/*=============================================================================
    Copyright (c) 2001-2011 Joel de Guzman

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/
#if !defined(FUSION_REVERSE_07212005_1230)
#define FUSION_REVERSE_07212005_1230

#include <boost/fusion/support/config.hpp>
#include <boost/fusion/support/is_sequence.hpp>
#include <boost/fusion/view/reverse_view/reverse_view.hpp>
#include <boost/utility/enable_if.hpp>

namespace boost {
namespace fusion {
namespace result_of {
template <typename Sequence> struct reverse {
  typedef reverse_view<Sequence> type;
};
} // namespace result_of

template <typename Sequence>
BOOST_CONSTEXPR BOOST_FUSION_GPU_ENABLED inline
    typename enable_if<traits::is_sequence<Sequence>,
                       reverse_view<Sequence const>>::type
    reverse(Sequence const &view) {
  return reverse_view<Sequence const>(view);
}
} // namespace fusion
} // namespace boost

#endif

max_stars_repo_name: henrywarhurst/matrix

max_issues_repo_name: henrywarhurst/matrix

max_forks_repo_name: henrywarhurst/matrix


text:
/*****************************
 * VIO used in mapping and locate.
 * Stereo-VIO is more complicated and unstable than mono-VIO,but it's effect is better.
 * Use stereo's depth information to improve accuracy of estimatation.
 * TODO:to compare pnp and icp. 2018-4-12.
*****************************/
#include<iostream>
#include<algorithm>
#include<fstream>
#include<chrono>

#include <ros/ros.h>
#include <cv_bridge/cv_bridge.h>

#include<opencv2/core/core.hpp>

#include"../../include/System.h"

#include "MsgSync/MsgSynchronizer.h"

#include "../../src/IMU/imudata.h"
#include "../../src/IMU/configparam.h"
#include <rosbag/bag.h>
#include <rosbag/view.h>

#include <boost/foreach.hpp>

using namespace std;

static volatile int keepRunning = 1;
void sig_handler( int sig )
{
    if ( sig == SIGINT )
    {
        keepRunning = 0;
    }
}

#define DEF_320x240 0

class ImageGrabber
{
public:
    ImageGrabber(ORB_SLAM2::System* pSLAM):mpSLAM(pSLAM){}

    void GrabStereo(const sensor_msgs::ImageConstPtr& msgLeft,const sensor_msgs::ImageConstPtr& msgRight);

    ORB_SLAM2::System* mpSLAM;
    bool do_rectify;
    cv::Mat M1l,M2l,M1r,M2r;
};

int main(int argc,char **argv)
{
    ros::init(argc,argv,"Stereo_mapping");
    ros::start();

    signal( SIGINT, sig_handler );

    if ( 2 != argc )
    {
        std::cerr << std::endl << "Usage: ./robot path_to_settings" << std::endl;
        return -1;
    }
    cv::FileStorage top_setting( argv[1], cv::FileStorage::READ );
        if ( ! top_setting.isOpened( ) )
        {
            std::cerr << "Failed to open settings file: " << argv[1] << std::endl;
            return -1;
        }

    #if DEF_320x240
        // ORB_SLAM2 settings
        std::string str_orb_setting_320( top_setting["orb_setting_320"] );
        cv::FileStorage orb_setting_320( str_orb_setting_320, cv::FileStorage::READ );
        if ( ! orb_setting_320.isOpened( ) )
        {
            std::cerr << "Failed to open settings file: " << str_orb_setting_320 << std::endl;
            return -1;
        }

        // ORB Vocabulary Path
        std::string str_orb_vocabulary( orb_setting_320["orb_vocabulary"] );
    #else
        // ORB_SLAM2 settings
        std::string str_orb_setting_640( top_setting["orb_setting_640"] );
        cv::FileStorage orb_setting_640( str_orb_setting_640, cv::FileStorage::READ );
        if ( ! orb_setting_640.isOpened( ) )
        {
            std::cerr << "Failed to open settings file: " << str_orb_setting_640 << std::endl;
            return -1;
        }
        // ORB Vocabulary Path
        std::string str_orb_vocabulary( orb_setting_640["orb_vocabulary"] );
    #endif

    // Sensor settings
    std::string str_sensor_setting( top_setting["sensor_setting"] );
    cv::FileStorage sensor_setting( str_sensor_setting, cv::FileStorage::READ );
    if ( ! sensor_setting.isOpened( ) )
    {
        std::cerr << "Failed to open settings file: " << str_sensor_setting << std::endl;
        return -1;
    }

    ORB_SLAM2::System SLAM(argv[1],argv[2],ORB_SLAM2::System::STEREO,true);
    ORB_SLAM2::ConfigParam config(argv[2]);

    double imageMsgDelaySec = config.GetImageDelayToIMU();
    ORBVIO::MsgSynchronizer msgsync(imageMsgDelaySec);
    ros::NodeHandle nh;
    ros::Subscriber imagesub0;
    ros::Subscriber imagesub1;
    ros::Subscriber imusub;
    if(ORB_SLAM2::ConfigParam::GetRealTimeFlag())
    {
         imagesub0 = nh.subscribe(config._imageLeftTopic, /*200*/ 2, &ORBVIO::MsgSynchronizer::imageleftCallback, &msgsync);
         imagesub1 = nh.subscribe(config._imageRightTopic, /*200*/ 2, &ORBVIO::MsgSynchronizer::imagerightCallback, &msgsync);
         imusub = nh.subscribe(config._imuTopic, 200, &ORBVIO::MsgSynchronizer::imuCallback, &msgsync);
    }
    sensor_msgs::ImageConstPtr imageLeftMsg;
    sensor_msgs::ImageConstPtr imageRightMsg;
    std::vector<sensor_msgs::ImuConstPtr> vimuMsg;

    // 3dm imu output per g. 1g=9.80665 according to datasheet
    const double g3dm = 9.80665;
    const bool bAccMultiply98 = config.GetAccMultiply9p8();

    ros::Rate r(1000);

    if(!ORB_SLAM2::ConfigParam::GetRealTimeFlag())
    {
        ROS_WARN("Run not-realtime");

        std::string bagfile = config._bagfile;
        rosbag::Bag bag;
        bag.open(bagfile,rosbag::bagmode::Read);

        std::vector<std::string> topics;
        std::string imutopic = config._imuTopic;
        std::string imagelefttopic = config._imageLeftTopic;
        std::string imagerighttopic = config._imageRightTopic;
        topics.push_back(imagelefttopic);
        topics.push_back(imagerighttopic);
        topics.push_back(imutopic);

        rosbag::View view(bag, rosbag::TopicQuery(topics));
        //while(ros::ok())
        BOOST_FOREACH(rosbag::MessageInstance const m, view)
        {
            sensor_msgs::ImuConstPtr simu = m.instantiate<sensor_msgs::Imu>();
            if(simu!=NULL)
                msgsync.imuCallback(simu);
            sensor_msgs::ImageConstPtr simage = m.instantiate<sensor_msgs::Image>();
            if(simage!=NULL)
                msgsync.imageleftCallback(simage);
            if(simage!=NULL)
                msgsync.imagerightCallback(simage);
            bool bdata = msgsync.getRecentMsgs(imageLeftMsg,vimuMsg);

            if(bdata)
            {
                std::vector<ORB_SLAM2::IMUData> vimuData;
                //ROS_INFO("image time: %.3f",imageMsg->header.stamp.toSec());
                for(unsigned int i=0;i<vimuMsg.size();i++)
                {
                    sensor_msgs::ImuConstPtr imuMsg = vimuMsg[i];
                    double ax = imuMsg->linear_acceleration.x;
                    double ay = imuMsg->linear_acceleration.y;
                    double az = imuMsg->linear_acceleration.z;
                    if(bAccMultiply98)
                    {
                        ax *= g3dm;
                        ay *= g3dm;
                        az *= g3dm;
                    }
                    ORB_SLAM2::IMUData imudata(imuMsg->angular_velocity.x,imuMsg->angular_velocity.y,imuMsg->angular_velocity.z,
                                    ax,ay,az,imuMsg->header.stamp.toSec());
                    vimuData.push_back(imudata);
                    //ROS_INFO("imu time: %.3f",vimuMsg[i]->header.stamp.toSec());
                }

                // Copy the ros image message to cv::Mat.
                cv_bridge::CvImageConstPtr cv_ptr0;
                cv_bridge::CvImageConstPtr cv_ptr1;
                try
                {
                    cv_ptr0 = cv_bridge::toCvShare(imageLeftMsg);
                    cv_ptr1 = cv_bridge::toCvShare(imageRightMsg);
                }
                catch (cv_bridge::Exception& e)
                {
                    ROS_ERROR("cv_bridge exception: %s", e.what());
                    return -1;
                }

                // Consider delay of image message
                //SLAM.TrackMonocular(cv_ptr->image, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                cv::Mat iml = cv_ptr0->image.clone();
                cv::Mat imr = cv_ptr1->image.clone();
                {
                    // To test relocalization
                    static double startT=-1;
                    if(startT<0)
                        startT = imageMsg->header.stamp.toSec();
                    // Below to test relocalizaiton
                    //if(imageMsg->header.stamp.toSec() > startT+25 && imageMsg->header.stamp.toSec() < startT+25.3)
                    if(imageMsg->header.stamp.toSec() < startT+config._testDiscardTime)
                        im = cv::Mat::zeros(im.rows,im.cols,im.type());
                }
                SLAM.TrackMonoVI(im, vimuData, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                //SLAM.TrackMonoVI(cv_ptr->image, vimuData, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                //cv::imshow("image",cv_ptr->image);

                // Wait local mapping end.
                bool bstop = false;
                while(!SLAM.bLocalMapAcceptKF())
                {
                    if(!ros::ok())
                    {
                        bstop=true;
                    }
                };
                if(bstop)
                    break;

            }

            //cv::waitKey(1);

            ros::spinOnce();
            r.sleep();
            if(!ros::ok())
                break;
        }

    }
    else
    {
        ROS_WARN("Run realtime");
        while( keepRunning )
        {
            bool bdata = msgsync.getRecentMsgs(imageMsg,vimuMsg);

            if(bdata)
            {
                std::vector<ORB_SLAM2::IMUData> vimuData;
                //ROS_INFO("image time: %.3f",imageMsg->header.stamp.toSec());
                for(unsigned int i=0;i<vimuMsg.size();i++)
                {
                    sensor_msgs::ImuConstPtr imuMsg = vimuMsg[i];
                    double ax = imuMsg->linear_acceleration.x;
                    double ay = imuMsg->linear_acceleration.y;
                    double az = imuMsg->linear_acceleration.z;
                    if(bAccMultiply98)
                    {
                        ax *= g3dm;
                        ay *= g3dm;
                        az *= g3dm;
                    }
                    ORB_SLAM2::IMUData imudata(imuMsg->angular_velocity.x,imuMsg->angular_velocity.y,imuMsg->angular_velocity.z,
                                    ax,ay,az,imuMsg->header.stamp.toSec());
                    vimuData.push_back(imudata);
                    //ROS_INFO("imu time: %.3f",vimuMsg[i]->header.stamp.toSec());
                }

                // Copy the ros image message to cv::Mat.
                cv_bridge::CvImageConstPtr cv_ptr;
                try
                {
                    cv_ptr = cv_bridge::toCvShare(imageMsg);
                }
                catch (cv_bridge::Exception& e)
                {
                    ROS_ERROR("cv_bridge exception: %s", e.what());
                    return -1;
                }

                // Consider delay of image message
                //SLAM.TrackMonocular(cv_ptr->image, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                cv::Mat im = cv_ptr->image.clone();
                {
                    // To test relocalization
                    static double startT=-1;
                    if(startT<0)
                        startT = imageMsg->header.stamp.toSec();
                    // Below to test relocalizaiton
                    //if(imageMsg->header.stamp.toSec() > startT+25 && imageMsg->header.stamp.toSec() < startT+25.3)
                    if(imageMsg->header.stamp.toSec() < startT+config._testDiscardTime)
                        im = cv::Mat::zeros(im.rows,im.cols,im.type());
                }
                SLAM.TrackMonoVI(im, vimuData, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                //SLAM.TrackMonoVI(cv_ptr->image, vimuData, imageMsg->header.stamp.toSec() - imageMsgDelaySec);
                //cv::imshow("image",cv_ptr->image);

            }

            //cv::waitKey(1);

            ros::spinOnce();
            r.sleep();
            if(!ros::ok())
                break;
        }
    }

    SLAM.Shutdown();

    std::string map_file_updated_name = orb_settig_640["map_file_updated_name"];
    SLAM.SaveMap( map_file_updated_name );

    std::string path_to_save = orb_setting_640["path_to_save"];
    SLAM.SaveKeyFrameTrajectoryTUM( path_to_save );

    ros::shutdown();
    return 0;
}

void ImageGrabber::GrabStereo(const sensor_msgs::ImageConstPtr& msgLeft,const sensor_msgs::ImageConstPtr& msgRight)
{
    // Copy the ros image message to cv::Mat.
    cv_bridge::CvImageConstPtr cv_ptrLeft;
    try
    {
        cv_ptrLeft = cv_bridge::toCvShare(msgLeft);
    }
    catch (cv_bridge::Exception& e)
    {
        ROS_ERROR("cv_bridge exception: %s", e.what());
        return;
    }

    cv_bridge::CvImageConstPtr cv_ptrRight;
    try
    {
        cv_ptrRight = cv_bridge::toCvShare(msgRight);
    }
    catch (cv_bridge::Exception& e)
    {
        ROS_ERROR("cv_bridge exception: %s", e.what());
        return;
    }

    if(do_rectify)
    {
        cv::Mat imLeft, imRight;
        cv::remap(cv_ptrLeft->image,imLeft,M1l,M2l,cv::INTER_LINEAR);
        cv::remap(cv_ptrRight->image,imRight,M1r,M2r,cv::INTER_LINEAR);
        mpSLAM->TrackStereo(imLeft,imRight,cv_ptrLeft->header.stamp.toSec());
    }
    else
    {
        mpSLAM->TrackStereo(cv_ptrLeft->image,cv_ptrRight->image,cv_ptrLeft->header.stamp.toSec());
    }

}

max_stars_repo_name: michaelczhou/vinsrobot

max_issues_repo_name: michaelczhou/vinsrobot

max_forks_repo_name: michaelczhou/vinsrobot


text:
/*
    Copyright 2013 Adobe
    Distributed under the Boost Software License, Version 1.0.
    (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
/**************************************************************************************************/

#include <adobe/eve_evaluate.hpp>

#include <mutex>

#include <boost/bind.hpp>

#include <adobe/algorithm/sort.hpp>
#include <adobe/array.hpp>
#include <adobe/dictionary.hpp>
#include <adobe/name.hpp>
#include <adobe/static_table.hpp>
#include <adobe/string.hpp>
#include <adobe/virtual_machine.hpp>
#include <adobe/once.hpp>

/**************************************************************************************************/

using namespace std;

/**************************************************************************************************/

namespace {

/**************************************************************************************************/

typedef std::pair<adobe::name_t*, adobe::name_t*> reflected_table_range_t;
typedef adobe::static_table<adobe::name_t, adobe::layout_attributes_t::alignment_t, 7>
alignment_table_t;
typedef adobe::static_table<adobe::name_t, adobe::layout_attributes_placement_t::placement_t, 3>
placement_table_t;

/**************************************************************************************************/

using namespace adobe::literals;

adobe::static_name_t key_spacing = "spacing"_name;
adobe::static_name_t key_indent = "indent"_name;
adobe::static_name_t key_margin = "margin"_name;

adobe::static_name_t key_placement = "placement"_name;

adobe::static_name_t key_horizontal = "horizontal"_name;
adobe::static_name_t key_vertical = "vertical"_name;

adobe::static_name_t key_child_horizontal = "child_horizontal"_name;
adobe::static_name_t key_child_vertical = "child_vertical"_name;

adobe::static_name_t key_align_left = "align_left"_name;
adobe::static_name_t key_align_right = "align_right"_name;
adobe::static_name_t key_align_top = "align_top"_name;
adobe::static_name_t key_align_bottom = "align_bottom"_name;
adobe::static_name_t key_align_center = "align_center"_name;
adobe::static_name_t key_align_proportional = "align_proportional"_name;
adobe::static_name_t key_align_fill = "align_fill"_name;

adobe::static_name_t key_place_row = "place_row"_name;
adobe::static_name_t key_place_column = "place_column"_name;
adobe::static_name_t key_place_overlay = "place_overlay"_name;

adobe::static_name_t key_guide_mask = "guide_mask"_name;
adobe::static_name_t key_guide_balance = "guide_balance"_name;

adobe::static_name_t key_guide_baseline = "guide_baseline"_name;
adobe::static_name_t key_guide_label = "guide_label"_name;

// blank(guide_attach: {left: @guide_baseline});

alignment_table_t* alignment_table_g;
placement_table_t* placement_table_g;
reflected_table_range_t* reflected_range_g;

/**************************************************************************************************/

void init_once_() {
    static alignment_table_t alignment_table_s = {
        {alignment_table_t::entry_type(key_align_left, adobe::layout_attributes_t::align_left),
         alignment_table_t::entry_type(key_align_right, adobe::layout_attributes_t::align_right),
         alignment_table_t::entry_type(key_align_top, adobe::layout_attributes_t::align_top),
         alignment_table_t::entry_type(key_align_bottom, adobe::layout_attributes_t::align_bottom),
         alignment_table_t::entry_type(key_align_center, adobe::layout_attributes_t::align_center),
         alignment_table_t::entry_type(key_align_proportional,
                                       adobe::layout_attributes_t::align_proportional),
         alignment_table_t::entry_type(key_align_fill, adobe::layout_attributes_t::align_fill)}};

    static placement_table_t placement_table_s = {
        {placement_table_t::entry_type(key_place_row,
                                       adobe::layout_attributes_placement_t::place_row),
         placement_table_t::entry_type(key_place_column,
                                       adobe::layout_attributes_placement_t::place_column),
         placement_table_t::entry_type(key_place_overlay,
                                       adobe::layout_attributes_placement_t::place_overlay)}};

    alignment_table_s.sort();
    placement_table_s.sort();

    alignment_table_g = &alignment_table_s;
    placement_table_g = &placement_table_s;

    static adobe::name_t reflected[] = {key_align_left,   key_align_right,  key_align_top,
                                        key_align_bottom, key_align_center, key_align_proportional,
                                        key_align_fill,   key_place_row,    key_place_column,
                                        key_place_overlay};
    static reflected_table_range_t reflected_table_range_s;

    adobe::sort(reflected);

    reflected_table_range_s.first = boost::begin(reflected);
    reflected_table_range_s.second = boost::end(reflected);

    reflected_range_g = &reflected_table_range_s;
}

/**************************************************************************************************/

once_flag flag;
void init_once() { call_once(flag, &init_once_); }

/**************************************************************************************************/

adobe::dictionary_t evaluate_named_arguments(adobe::virtual_machine_t& evaluator,
                                             const adobe::array_t& arguments) {
    evaluator.evaluate(arguments);

    adobe::dictionary_t result(std::move(evaluator.back().cast<adobe::dictionary_t>()));

    evaluator.pop_back();
    return result;
}

/**************************************************************************************************/

void add_cell(adobe::sheet_t& sheet, adobe::eve_callback_suite_t::cell_type_t type,
              adobe::name_t name, const adobe::line_position_t& position,
              const adobe::array_t& init_or_expr) {
    switch (type) {
    case adobe::eve_callback_suite_t::constant_k:
        sheet.add_constant(name, position, init_or_expr);
        break;
    case adobe::eve_callback_suite_t::logic_k:
        sheet.add_logic(name, position, init_or_expr);
        break;
    default:
        assert(false && "cell type not supported.");
    }
}

/**************************************************************************************************/

// REVISIT (sparent) : Copy & Paste from adam_evaluate.cpp

void add_relation(adobe::sheet_t& sheet, const adobe::line_position_t& position,
                  const adobe::array_t& conditional,
                  const adobe::eve_callback_suite_t::relation_t* first,
                  const adobe::eve_callback_suite_t::relation_t* last) {
    typedef std::vector<adobe::sheet_t::relation_t> relation_buffer_t;

    relation_buffer_t relations;

    relations.reserve(relation_buffer_t::size_type(last - first));

    /*
        REVISIT (sparent) : It would be nice to find a simple way to handle a transformed copy
        in a generic fashion when multiple members are needed.
    */

    while (first != last) // copy
    {
        relations.push_back(
            adobe::sheet_t::relation_t(first->name_set_m, first->position_m, first->expression_m));
        ++first;
    }

    sheet.add_relation(position, conditional, &relations[0], &relations[0] + relations.size());
}

/**************************************************************************************************/

} // namespace

/**************************************************************************************************/

namespace adobe {

/**************************************************************************************************/

eve_callback_suite_t bind_layout(const bind_layout_proc_t& proc, sheet_t& sheet,
                                 virtual_machine_t& evaluator) {
    init_once();

    eve_callback_suite_t suite;

    suite.add_view_proc_m = boost::bind(
        proc, _1, _3, boost::bind(&evaluate_named_arguments, boost::ref(evaluator), _4));
    suite.add_cell_proc_m = boost::bind(&add_cell, boost::ref(sheet), _1, _2, _3, _4);
    suite.add_relation_proc_m = boost::bind(&add_relation, boost::ref(sheet), _1, _2, _3, _4);
    suite.add_interface_proc_m =
        boost::bind(&adobe::sheet_t::add_interface, boost::ref(sheet), _1, _2, _3, _4, _5, _6);
    suite.finalize_sheet_proc_m = boost::bind(&sheet_t::update, boost::ref(sheet));

    return suite;
}

/**************************************************************************************************/

void apply_layout_parameters(layout_attributes_t& data, const dictionary_t& parameters) {
    init_once();

    get_value(parameters, key_indent, data.indent_m);

    if (parameters.count(key_horizontal)) {
        (*alignment_table_g)(get_value(parameters, key_horizontal).cast<name_t>(),
                             data.slice_m[extents_slices_t::horizontal].alignment_m);
    }

    if (parameters.count(key_vertical)) {
        (*alignment_table_g)(get_value(parameters, key_vertical).cast<name_t>(),
                             data.slice_m[extents_slices_t::vertical].alignment_m);
    }

    // REVISIT (sparent) If we had named guides then named guide suppression would go here.

    if (parameters.count(key_guide_mask)) // an empty array would allow baselines out.
    {
        // Turn on all guides - then selectively suppress
        data.slice_m[extents_slices_t::vertical].suppress_m = false;
        data.slice_m[extents_slices_t::horizontal].suppress_m = false;

        array_t guide_mask(get_value(parameters, key_guide_mask).cast<array_t>());

        for (array_t::const_iterator iter(guide_mask.begin()), last(guide_mask.end()); iter != last;
             ++iter) {
            if (iter->cast<name_t>() == key_guide_baseline)
                data.slice_m[extents_slices_t::vertical].suppress_m = true;
            if (iter->cast<name_t>() == key_guide_label)
                data.slice_m[extents_slices_t::horizontal].suppress_m = true;
        }
    }

    // REVISIT (sparent) : If we had named guides then named guide balancing would go here...

    /*
        Balanced guides must be supressed to avoid having them attach to outside guides which could
        overconstrain the system.
    */

    if (parameters.count(key_guide_balance)) {
        data.slice_m[extents_slices_t::vertical].balance_m = false;
        data.slice_m[extents_slices_t::horizontal].balance_m = false;

        array_t guide_balance(get_value(parameters, key_guide_balance).cast<array_t>());

        for (array_t::const_iterator iter(guide_balance.begin()), last(guide_balance.end());
             iter != last; ++iter) {
            if (iter->cast<name_t>() == key_guide_baseline) {
                data.slice_m[extents_slices_t::vertical].balance_m = true;
                data.slice_m[extents_slices_t::vertical].suppress_m = true;
            }
            if (iter->cast<name_t>() == key_guide_label) {
                data.slice_m[extents_slices_t::horizontal].balance_m = true;
                data.slice_m[extents_slices_t::horizontal].suppress_m = true;
            }
        }
    }

    // REVISIT (sparent) : I'm seeing a pattern here - with three cases this could be factored

    {
        dictionary_t::const_iterator iter(parameters.find(key_placement));

        if (iter != parameters.end()) {
            (*placement_table_g)(iter->second.cast<name_t>(), data.placement_m);
        }

        // Adjust defaults

        // Specifying a row from the parameters implies enabling baselines unless otherwise
        // specified.

        if (iter != parameters.end() &&
            data.placement_m == layout_attributes_placement_t::place_row &&
            !parameters.count(key_guide_mask)) {
            data.slice_m[extents_slices_t::vertical].suppress_m = false;
        }
    }

    {
        dictionary_t::const_iterator iter(parameters.find(key_child_horizontal));

        if (iter != parameters.end()) {
            (*alignment_table_g)(iter->second.cast<name_t>(),
                                 data.slice_m[extents_slices_t::horizontal].child_alignment_m);
        }
    }

    {
        dictionary_t::const_iterator iter(parameters.find(key_child_vertical));

        if (iter != parameters.end()) {
            (*alignment_table_g)(iter->second.cast<name_t>(),
                                 data.slice_m[extents_slices_t::vertical].child_alignment_m);
        }
    }

    // spacing
    {
        dictionary_t::const_iterator iter(parameters.find(key_spacing));
        if (iter != parameters.end()) {
            if (iter->second.type_info() == typeid(array_t)) {
                const array_t& spacing_array = iter->second.cast<array_t>();
                data.spacing_m.resize(spacing_array.size() + 1);

                layout_attributes_t::spacing_t::iterator dest_iter(data.spacing_m.begin() + 1);

                for (array_t::const_iterator iter(spacing_array.begin());
                     iter != spacing_array.end(); ++iter) {
                    *dest_iter = iter->cast<long>();
                    ++dest_iter;
                }
            } else {
                double tmp(data.spacing_m[1]);
                iter->second.cast(tmp); // Try getting as number
                data.spacing_m[1] = long(tmp);
            }
        }
    }

    // margin
    {
        dictionary_t::const_iterator iter(parameters.find(key_margin));
        if (iter != parameters.end()) {
            if (iter->second.type_info() == typeid(array_t)) {
                const array_t& margin_array = iter->second.cast<array_t>();

                data.vertical().margin_m.first = margin_array[0].cast<long>();
                data.horizontal().margin_m.first = margin_array[1].cast<long>();
                data.vertical().margin_m.second = margin_array[2].cast<long>();
                data.horizontal().margin_m.second = margin_array[3].cast<long>();
            } else {
                long margin = iter->second.cast<long>();

                data.vertical().margin_m.first = margin;
                data.horizontal().margin_m.first = margin;
                data.vertical().margin_m.second = margin;
                data.horizontal().margin_m.second = margin;
            }
        }
    }
}

/**************************************************************************************************/

adobe::any_regular_t layout_variables(adobe::sheet_t& layout_sheet, adobe::name_t name) {
    adobe::name_t* found(
        std::lower_bound(reflected_range_g->first, reflected_range_g->second, name));

    if (found != reflected_range_g->second && *found == name)
        return adobe::any_regular_t(name);

    return layout_sheet.get(name);
}

/**************************************************************************************************/

} // namespace adobe

/**************************************************************************************************/


text:
/*
Part of the Fluid Corpus Manipulation Project (http://www.flucoma.org/)
Copyright 2017-2019 University of Huddersfield.
Licensed under the BSD-3 License.
See license.md file in the project root for full license information.
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 725899).
*/

// Jonathan Driedger, Thomas Prätzlich, and Meinard Müller
// Let It Bee — Towards NMF-Inspired Audio Mosaicing
// Proceedings of ISMIR 2015

#pragma once

#include "STFT.hpp"
#include "../util/FluidEigenMappings.hpp"
#include "../../data/FluidIndex.hpp"
#include "../../data/TensorTypes.hpp"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <iostream>
#include <vector>

namespace fluid {
namespace algorithm {

using _impl::asEigen;
using _impl::asFluid;
using Eigen::Array;
using Eigen::ArrayXd;
using Eigen::ArrayXXd;
using Eigen::Matrix;
using Eigen::MatrixXd;
using Eigen::VectorXd;

class NMFCross
{

public:
  // pass iteration number; returns true if able to continue (i.e. not
  // cancelled)
  using ProgressCallback = std::function<bool(index)>;

  NMFCross(index nIterations) : mIterations(nIterations) {}

  static void synthesize(const RealMatrixView h, const ComplexMatrixView w,
                         ComplexMatrixView out)
  {
    using namespace Eigen;
    using namespace _impl;
    MatrixXd  H = asEigen<Matrix>(h);
    MatrixXcd W = asEigen<Matrix>(w);
    MatrixXcd V = H * W;
    out <<= asFluid(V);
  }

  void process(const RealMatrixView X, RealMatrixView H1, RealMatrixView W0,
               index r, index p, index c) const
  {
    index nFrames = X.extent(0);
    index nBins = X.extent(1);
    index rank = W0.extent(0);
    nBins = W0.extent(1);
    MatrixXd W = asEigen<Matrix>(W0).transpose();
    MatrixXd H;
    H = MatrixXd::Random(rank, nFrames) * 0.5 +
        MatrixXd::Constant(rank, nFrames, 0.5);
    MatrixXd V = asEigen<Matrix>(X).transpose();
    multiplicativeUpdates(V, W, H, r, p, c);
    MatrixXd HT = H.transpose();
    H1 <<= asFluid(HT);
  }

  void addProgressCallback(ProgressCallback&& callback)
  {
    mCallbacks.emplace_back(std::move(callback));
  }

private:
  index                         mIterations;
  std::vector<ProgressCallback> mCallbacks;

  std::vector<index> topC(Eigen::VectorXd vec, index c) const
  {
    using namespace std;
    vector<double> stdVec(vec.data(), vec.data() + vec.size());
    sort(stdVec.begin(), stdVec.end());
    vector<index> idx(asUnsigned(vec.size()));
    iota(idx.begin(), idx.end(), 0);
    sort(idx.begin(), idx.end(),
         [&vec](index i1, index i2) { return vec[i1] > vec[i2]; });
    auto result = std::vector<index>(idx.begin(), idx.begin() + c);
    return result;
  }

  Eigen::MatrixXd promoteContinuity(MatrixXd& H, index size) const
  {
    index    halfSize = (size - 1) / 2;
    MatrixXd kernel = MatrixXd::Identity(size, size);
    MatrixXd padded = MatrixXd::Zero(H.rows() + size, H.cols() + size);
    MatrixXd output = MatrixXd::Zero(H.rows(), H.cols());
    padded.block(halfSize, halfSize, H.rows(), H.cols()) = H;
    for (index i = 0; i < H.rows(); i++)
    {
      for (index j = 0; j < H.cols(); j++)
      {
        output(i, j) =
            padded.block(i, j, size, size).cwiseProduct(kernel).sum();
      }
    }
    return output;
  }

  Eigen::MatrixXd enforceTemporalSparseness(MatrixXd& H, index size,
                                            index iteration) const
  {
    index    halfSize = (size - 1) / 2;
    MatrixXd padded = MatrixXd::Zero(H.rows(), H.cols() + size);
    MatrixXd output = MatrixXd::Zero(H.rows(), H.cols());
    padded.block(0, halfSize, H.rows(), H.cols()) = H;
    for (index i = 0; i < H.rows(); i++)
    {
      for (index j = 0; j < H.cols(); j++)
      {
        VectorXd        neighborhood = padded.row(i).segment(j, size);
        VectorXd::Index maxIndex{0};
        neighborhood.maxCoeff(&maxIndex);
        if (int(maxIndex) != halfSize)
        { output(i, j) = H(i, j) * (1 - ((iteration + 1) / mIterations)); }
        else
        {
          output(i, j) = H(i, j);
        }
      }
    }
    return output;
  }

  Eigen::MatrixXd restrictPolyphony(MatrixXd& H, ArrayXd& energyInW, index size,
                                    index iteration) const
  {
    MatrixXd output = MatrixXd::Zero(H.rows(), H.cols());
    for (index k = 0; k < H.cols(); k++)
    {
      ArrayXd wCol = H.col(k).array() * energyInW.array();
      output.col(k) = H.col(k) * (1 - ((iteration + 1) / mIterations));
      auto top = topC(wCol, size);
      for (auto t : top) { output(t, k) = H(t, k); }
    }
    return output;
  }
  void multiplicativeUpdates(MatrixXd& V, MatrixXd& W, MatrixXd& H, index r,
                             index p, index c) const
  {
    using namespace std;
    using namespace Eigen;
    double const epsilon = std::numeric_limits<double>::epsilon();
    MatrixXd     ones = MatrixXd::Ones(V.rows(), V.cols());
    W = W.array().max(epsilon).matrix();
    // ArrayXd wNorm = W.colwise().sum();
    // W.array().rowwise() /= wNorm.transpose());
    ArrayXd energyInW = W.array().square().colwise().sum();
    for (index i = 0; i < mIterations; i++)
    {
      if ((i % 1) == 0)
      { // TODO: original version seems to work better with one in 5 iterations
        H = enforceTemporalSparseness(H, r, i);
        H = restrictPolyphony(H, energyInW, p, i);
        H = promoteContinuity(H, c);
      }
      ArrayXXd V2 = (W * H).array().max(epsilon);
      ArrayXXd hnum = (W.transpose() * (V.array() / V2).matrix()).array();
      ArrayXXd hden = (W.transpose() * ones).array();
      H = (H.array() * hnum / hden.max(epsilon)).matrix();
      // MatrixXd R = W * H;
      // R = R.cwiseMax(epsilon);
      // double divergence = (V.cwiseProduct(V.cwiseQuotient(R)) - V + R).sum();
      for (auto& cb : mCallbacks)
        if (!cb(i + 1)) return;
    }
    V = W * H;
  }
};
} // namespace algorithm
} // namespace fluid