text:
/*
 * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "otbOGRDataSourceWrapper.h"
// standard includes
#include <cassert>
#include <numeric>
#include <algorithm>
#include <clocale> // toupper
#include <boost/bind.hpp>
#include <boost/foreach.hpp>
#include <boost/range/begin.hpp>
#include <boost/range/end.hpp>
// ITK includes
#include "itkMacro.h" // itkExceptionMacro
#include "itkExceptionObject.h"
#include "itksys/SystemTools.hxx"
// OTB includes
#include "otbOGRDriversInit.h"
#include "otbSystem.h"
// OGR includes

/*===========================================================================*/
/*=======================[ construction/destruction ]========================*/
/*===========================================================================*/

bool otb::ogr::DataSource::Clear()
{
  Reset(nullptr);
  return true;
}

void otb::ogr::DataSource::Reset(GDALDataset* source)
{
  if (m_DataSource)
  {
    // OGR makes a pointless check for non-nullity in
    // GDALDataset::DestroyDataSource (pointless because "delete 0" is
    // perfectly valid -> it's a no-op)
    GDALClose(m_DataSource); // void, noexcept
  }
  m_DataSource = source;
}

namespace
{ // Anonymous namespace
  /**\ingroup GeometryInternals
   * \brief Type for associating filename extension with OGR driver names.
   * \since OTB v 3.14.0
   */
struct ExtensionDriverAssociation
{
  char const* extension;
  char const* driverName;
  bool Matches(std::string const& ext) const
  {
    return ext == extension;
  }
};
/**\ingroup GeometryInternals
 * \brief Associative table of filename extension -> OGR driver names.
 * \since OTB v 3.14.0
 */
const ExtensionDriverAssociation k_ExtensionDriverMap[] = {
    {".SHP", "ESRI Shapefile"}, {".TAB", "MapInfo File"}, {".GML", "GML"}, {".GMT", "OGR_GMT"}, {".GPX", "GPX"},
    {".SQLITE", "SQLite"},      {".KML", "KML"},          {".CSV", "CSV"}, {".GPKG", "GPKG"}};
/**\ingroup GeometryInternals
 * \brief Returns the OGR driver name associated to a filename.
 * \since OTB v 3.14.0
 * \note Relies on the driver name associated to the filename extension in \c
 * k_ExtensionDriverMap.
 * \note As a special case, filenames starting with "PG:" are bound to
 * "PostgreSQL".
 */
char const* DeduceDriverName(std::string filename)
{
  std::transform(filename.begin(), filename.end(), filename.begin(), (int (*)(int))toupper);
  if (filename.compare(0, 3, "PG:") == 0)
  {
    return "PostgreSQL";
  }
  const std::string                 extension = itksys::SystemTools::GetFilenameLastExtension(filename);
  ExtensionDriverAssociation const* whichIt =
      std::find_if(boost::begin(k_ExtensionDriverMap), boost::end(k_ExtensionDriverMap), [&](auto const& x) { return x.Matches(extension); });
  if (whichIt == boost::end(k_ExtensionDriverMap))
  {
    return nullptr; // nothing found
  }
  return whichIt->driverName;
}
} // Anonymous namespace

otb::ogr::DataSource::DataSource() : m_DataSource(nullptr), m_LayerOptions(), m_OpenMode(Modes::Update_LayerUpdate), m_FirstModifiableLayerID(0)
{
  Drivers::Init();

  GDALDriver* d = GetGDALDriverManager()->GetDriverByName("Memory");
  assert(d && "OGR Memory driver not found");
  m_DataSource = d->Create("in-memory", 0, 0, 0, GDT_Unknown, 0);

  if (!m_DataSource)
  {
    itkExceptionMacro(<< "Failed to create OGRMemDataSource: " << CPLGetLastErrorMsg());
  }
}

otb::ogr::DataSource::DataSource(GDALDataset* source, Modes::type mode, const std::vector<std::string>& options /*NULL*/)
  : m_DataSource(source), m_LayerOptions(options), m_OpenMode(mode), m_FirstModifiableLayerID(0)
{
  m_FirstModifiableLayerID = GetLayersCount();
}

otb::ogr::DataSource::Pointer otb::ogr::DataSource::OpenDataSource(std::string const& datasourceName, Modes::type mode)
{
  FileNameHelperType::Pointer fileNameHelper = FileNameHelperType::New();
  fileNameHelper->SetExtendedFileName(datasourceName.c_str());
  std::string simpleFileName = fileNameHelper->GetSimpleFileName();

  bool         update = (mode != Modes::Read);
  GDALDataset* source = (GDALDataset*)GDALOpenEx(simpleFileName.c_str(), (update ? GDAL_OF_UPDATE : GDAL_OF_READONLY) | GDAL_OF_VECTOR, NULL,
                                                 otb::ogr::StringListConverter(fileNameHelper->GetGDALOpenOptions()).to_ogr(), NULL);
  if (!source)
  {
    // In read mode, this is a failure
    // In write mode (Overwrite and Update), create the data source transparently
    if (mode == Modes::Read)
    {
      itkGenericExceptionMacro(<< "Failed to open GDALDataset file " << simpleFileName << " : " << CPLGetLastErrorMsg());
    }

    // Hand made factory based on file extension.
    char const* driverName = DeduceDriverName(simpleFileName);
    if (!driverName)
    {
      itkGenericExceptionMacro(<< "No OGR driver known to OTB to create and "
                                  "handle a DataSource named <"
                               << simpleFileName << ">.");
    }

    GDALDriver* d = GetGDALDriverManager()->GetDriverByName(driverName);

    if (!d)
    {
      itkGenericExceptionMacro(<< "Could not create OGR driver " << driverName << ", check your OGR configuration for available drivers.");
    }

    source = d->Create(simpleFileName.c_str(), 0, 0, 0, GDT_Unknown, otb::ogr::StringListConverter(fileNameHelper->GetGDALCreationOptions()).to_ogr());
    if (!source)
    {
      itkGenericExceptionMacro(<< "Failed to create GDALDataset <" << simpleFileName << "> (driver name: <" << driverName << ">: " << CPLGetLastErrorMsg());
    }
  }
  return otb::ogr::DataSource::New(source, mode, fileNameHelper->GetGDALLayerOptions());
}

void DeleteDataSource(std::string const& datasourceName)
{
  otb::OGRExtendedFilenameToOptions::Pointer fileNameHelper = otb::OGRExtendedFilenameToOptions::New();
  fileNameHelper->SetExtendedFileName(datasourceName.c_str());
  std::string simpleFileName = fileNameHelper->GetSimpleFileName();

  // Open dataset
  GDALDataset* poDS = (GDALDataset*)GDALOpenEx(simpleFileName.c_str(), GDAL_OF_UPDATE | GDAL_OF_VECTOR, NULL, NULL, NULL);

  GDALDriver* poDriver = NULL;
  if (poDS)
  {
    poDriver = poDS->GetDriver();
    GDALClose(poDS);
  }
  else
  {
    itkGenericExceptionMacro(<< "Cannot open data source " << simpleFileName << ": " << CPLGetLastErrorMsg());
  }
  if (poDriver)
  {
    OGRErr ret = poDriver->Delete(simpleFileName.c_str());
    if (ret != OGRERR_NONE)
    {
      itkGenericExceptionMacro(<< "Deletion of data source " << simpleFileName << " failed: " << CPLGetLastErrorMsg());
    }
  }
  else
  {
    itkGenericExceptionMacro(<< "Cannot get driver associated with data source " << simpleFileName << ": " << CPLGetLastErrorMsg());
  }
}

otb::ogr::DataSource::Pointer otb::ogr::DataSource::New(std::string const& datasourceName, Modes::type mode)
{
  FileNameHelperType::Pointer fileNameHelper = FileNameHelperType::New();
  fileNameHelper->SetExtendedFileName(datasourceName.c_str());
  std::string simpleFileName = fileNameHelper->GetSimpleFileName();

  if (mode < Modes::Read || mode >= Modes::MAX__)
  {
    itkGenericExceptionMacro(<< "Wrong mode when opening " << simpleFileName);
  }

  Drivers::Init();
  GDALDataset* ds        = (GDALDataset*)GDALOpenEx(simpleFileName.c_str(), GDAL_OF_READONLY | GDAL_OF_VECTOR, NULL, NULL, NULL);
  bool         ds_exists = (ds != nullptr);

  GDALClose(ds);

  if (ds_exists && mode == Modes::Overwrite)
  {
    DeleteDataSource(datasourceName);
  }

  return OpenDataSource(datasourceName, mode);
}

/*static*/
otb::ogr::DataSource::Pointer otb::ogr::DataSource::New(GDALDataset* source, Modes::type mode, const std::vector<std::string>& layerOptions)
{
  Pointer res = new DataSource(source, mode, layerOptions);
  res->UnRegister();
  return res;
}

/*virtual*/ otb::ogr::DataSource::~DataSource()
{
  Clear();
}

/*===========================================================================*/
/*================================[ layers ]=================================*/
/*===========================================================================*/
otb::ogr::DataSource::const_iterator otb::ogr::DataSource::cbegin() const
{
  return const_iterator(*this, 0);
}

otb::ogr::DataSource::const_iterator otb::ogr::DataSource::cend() const
{
  return const_iterator(*this, GetLayersCount());
}

otb::ogr::DataSource::iterator otb::ogr::DataSource::begin()
{
  return iterator(*this, 0);
}

otb::ogr::DataSource::iterator otb::ogr::DataSource::end()
{
  return iterator(*this, GetLayersCount());
}

otb::ogr::Layer otb::ogr::DataSource::CreateLayer(std::string const& name, OGRSpatialReference* poSpatialRef /* = NULL */,
                                                  OGRwkbGeometryType eGType /* = wkbUnknown */, std::vector<std::string> const& papszOptions /* = NULL */)
{
  assert(m_DataSource && "Datasource not initialized");

  // In read mode, try to return an existing layer.
  // and throw if not found
  if (m_OpenMode == Modes::Read)
  {
    otb::ogr::Layer l = GetLayerChecked(name); // will throw if not existing
    itkGenericOutputMacro(<< "Requesting layer creation in read-only "
                             "GDALDataset. Returning the existing layer");
    return l;
  }

  // Other mode : Check if the layer already exists.
  otb::ogr::Layer layer = GetLayer(name); // won't throw on failure

  FileNameHelperType::Pointer layerOptionHelper = FileNameHelperType::GetGDALLayerOptionsHelper(m_LayerOptions);
  layerOptionHelper->AddGDALLayerOptions(papszOptions);
  std::vector<std::string> layerOptions = layerOptionHelper->GetGDALLayerOptions();

  switch (m_OpenMode)
  {
  case Modes::Update_LayerOverwrite:
  {
    // First delete the layer
    if (layer)
    {
      DeleteLayer(GetLayerID(name));
    }

    // Then create it
    OGRLayer* ol = m_DataSource->CreateLayer(name.c_str(), poSpatialRef, eGType, otb::ogr::StringListConverter(layerOptions).to_ogr());

    if (!ol)
    {
      itkGenericExceptionMacro(<< "Failed to create the layer <" << name << "> in the GDALDataset file <" << GetDatasetDescription()
                               << ">: " << CPLGetLastErrorMsg());
    }

    const bool modifiable = true;
    return Layer(ol, modifiable);
  }
  break;
  case Modes::Overwrite:
  case Modes::Update_LayerUpdate:
  {
    if (layer)
    {
      return layer;
    }
    else
    {
      // Then create it
      OGRLayer* ol = m_DataSource->CreateLayer(name.c_str(), poSpatialRef, eGType, otb::ogr::StringListConverter(layerOptions).to_ogr());

      if (!ol)
      {
        itkGenericExceptionMacro(<< "Failed to create the layer <" << name << "> in the GDALDataset file <" << GetDatasetDescription()
                                 << ">: " << CPLGetLastErrorMsg());
      }

      const bool modifiable = true;
      return Layer(ol, modifiable);
    }
  }
  break;
  case Modes::Update_LayerCreateOnly:
  {
    // The layer exists. Returns it
    // It will be non-modifiable if already existing at GDALDataset creation time
    if (layer)
    {
      return layer;
    }

    // Case where the layer does not exists
    OGRLayer* ol = m_DataSource->CreateLayer(name.c_str(), poSpatialRef, eGType, otb::ogr::StringListConverter(layerOptions).to_ogr());

    if (!ol)
    {
      itkGenericExceptionMacro(<< "Failed to create the layer <" << name << "> in the GDALDataset file <" << GetDatasetDescription()
                               << ">: " << CPLGetLastErrorMsg());
    }

    const bool modifiable = true;
    return Layer(ol, modifiable);
  }
  break;
  default:
    assert(false && "Should never be there");
    itkGenericExceptionMacro(<< "GDALDataset opening mode not supported");
    break;
  }

  itkGenericExceptionMacro(<< "GDALDataset opening mode not supported");
  return Layer(nullptr, false); // keep compiler happy
}

otb::ogr::Layer otb::ogr::DataSource::CopyLayer(Layer& srcLayer, std::string const& newName, std::vector<std::string> const& papszOptions /* = NULL */)
{
  assert(m_DataSource && "Datasource not initialized");

  switch (m_OpenMode)
  {
  case Modes::Invalid:
    assert(false && "Invalid GDALDataset opening mode");
    itkGenericExceptionMacro(<< "Invalid GDALDataset opening mode");
    break;
  case Modes::Read:
    itkGenericExceptionMacro(<< "GDALDataset is opened in Read mode : "
                                "cannot create a layer");
    break;
  default:
    break;
  }

  FileNameHelperType::Pointer layerOptionHelper = FileNameHelperType::GetGDALLayerOptionsHelper(m_LayerOptions);
  layerOptionHelper->AddGDALLayerOptions(papszOptions);
  std::vector<std::string> layerOptions = layerOptionHelper->GetGDALLayerOptions();

  OGRLayer* l0 = &srcLayer.ogr();
  OGRLayer* ol = m_DataSource->CopyLayer(l0, newName.c_str(), otb::ogr::StringListConverter(layerOptions).to_ogr());
  if (!ol)
  {
    itkGenericExceptionMacro(<< "Failed to copy the layer <" << srcLayer.GetName() << "> into the new layer <" << newName << "> in the GDALDataset file <"
                             << GetDatasetDescription() << ">: " << CPLGetLastErrorMsg());
  }
  const bool modifiable = true;
  Layer      l(ol, modifiable);
  return l;
}

void otb::ogr::DataSource::DeleteLayer(size_t i)
{
  assert(m_DataSource && "Datasource not initialized");

  switch (m_OpenMode)
  {
  case Modes::Invalid:
    assert(false && "Invalid GDALDataset opening mode");
    itkGenericExceptionMacro(<< "Invalid GDALDataset opening mode");
    break;
  case Modes::Read:
    itkGenericExceptionMacro(<< "GDALDataset is opened in Read mode : "
                                "cannot delete a layer");
    break;
  case Modes::Update_LayerCreateOnly:
    itkGenericExceptionMacro(<< "GDALDataset is opened in "
                                "Update_LayerCreateOnly mode : cannot delete a layer");
    break;
  default:
    break;
  }

  const int nb_layers = GetLayersCount();
  if (int(i) >= nb_layers)
  {
    itkExceptionMacro(<< "Cannot delete " << i << "th layer in the GDALDataset <" << GetDatasetDescription() << "> as it contains only " << nb_layers
                      << "layers.");
  }
  const OGRErr err = m_DataSource->DeleteLayer(int(i));
  if (err != OGRERR_NONE)
  {
    itkExceptionMacro(<< "Cannot delete " << i << "th layer in the GDALDataset <" << GetDatasetDescription() << ">: " << CPLGetLastErrorMsg());
  }
}

bool otb::ogr::DataSource::IsLayerModifiable(size_t i) const
{
  assert(m_DataSource && "Datasource not initialized");
  switch (m_OpenMode)
  {
  case Modes::Read:
    return false;
  case Modes::Update_LayerCreateOnly:
    return int(i) >= m_FirstModifiableLayerID;
  default:
    return true;
  }
}

bool otb::ogr::DataSource::IsLayerModifiable(std::string const& layername) const
{
  assert(m_DataSource && "Datasource not initialized");
  switch (m_OpenMode)
  {
  case Modes::Read:
    return false;
  case Modes::Update_LayerCreateOnly:
  {
    const int id = this->GetLayerIDUnchecked(layername);
    return id >= m_FirstModifiableLayerID;
  }
  default:
    return true;
  }
}

int otb::ogr::DataSource::GetLayerIDUnchecked(std::string const& name) const
{
  assert(m_DataSource && "Datasource not initialized");
  for (int i = 0, N = GetLayersCount(); i < N; i++)
  {
    OGRLayer* raw_layer = GetLayerUnchecked(i);
    // wrapping to give access to GetName() that doesn't exist with all version of gdal/ogr
    Layer layer(raw_layer, false);
    if (layer.GetName() == name)
    {
      return i;
    }
  }
  return -1;
}

size_t otb::ogr::DataSource::GetLayerID(std::string const& name) const
{
  int const id = GetLayerIDUnchecked(name);
  if (id < 0)
  {
    itkExceptionMacro(<< "Cannot fetch any layer named <" << name << "> in the GDALDataset <" << GetDatasetDescription() << ">: " << CPLGetLastErrorMsg());
  }
  return 0; // keep compiler happy
}

otb::ogr::Layer otb::ogr::DataSource::GetLayerChecked(size_t i)
{
  assert(m_DataSource && "Datasource not initialized");
  const int nb_layers = GetLayersCount();
  if (int(i) >= nb_layers)
  {
    itkExceptionMacro(<< "Cannot fetch " << i << "th layer in the GDALDataset <" << GetDatasetDescription() << "> as it contains only " << nb_layers
                      << "layers.");
  }
  OGRLayer* layer_ptr = m_DataSource->GetLayer(int(i));
  if (!layer_ptr)
  {
    itkExceptionMacro(<< "Unexpected error: cannot fetch " << i << "th layer in the GDALDataset <" << GetDatasetDescription() << ">: " << CPLGetLastErrorMsg());
  }
  return otb::ogr::Layer(layer_ptr, IsLayerModifiable(i));
}

OGRLayer* otb::ogr::DataSource::GetLayerUnchecked(size_t i)
{
  assert(m_DataSource && "Datasource not initialized");
  OGRLayer* layer_ptr = m_DataSource->GetLayer(int(i));
  return layer_ptr;
}

otb::ogr::Layer otb::ogr::DataSource::GetLayer(std::string const& name)
{
  assert(m_DataSource && "Datasource not initialized");
  OGRLayer* layer_ptr = m_DataSource->GetLayerByName(name.c_str());
  return otb::ogr::Layer(layer_ptr, layer_ptr ? IsLayerModifiable(name) : false);
}

otb::ogr::Layer otb::ogr::DataSource::GetLayerChecked(std::string const& name)
{
  assert(m_DataSource && "Datasource not initialized");
  OGRLayer* layer_ptr = m_DataSource->GetLayerByName(name.c_str());
  if (!layer_ptr)
  {
    itkExceptionMacro(<< "Cannot fetch any layer named <" << name << "> in the GDALDataset <" << GetDatasetDescription() << ">: " << CPLGetLastErrorMsg());
  }
  return otb::ogr::Layer(layer_ptr, IsLayerModifiable(name));
}

int otb::ogr::DataSource::GetLayersCount() const
{
  assert(m_DataSource && "Datasource not initialized");
  return m_DataSource->GetLayerCount();
}

otb::ogr::Layer otb::ogr::DataSource::ExecuteSQL(std::string const& statement, OGRGeometry* poSpatialFilter, char const* pszDialect)
{
  assert(m_DataSource && "Datasource not initialized");
  const bool modifiable = false;
  OGRLayer*  layer_ptr  = m_DataSource->ExecuteSQL(statement.c_str(), poSpatialFilter, pszDialect);
  if (!layer_ptr)
  {
#if defined(PREFER_EXCEPTION)
    itkExceptionMacro(<< "Unexpected error: cannot execute the SQL request <" << statement << "> in the GDALDataset <" << GetDatasetDescription()
                      << ">: " << CPLGetLastErrorMsg());
#else
    // Cannot use the deleter made for result sets obtained from
    // GDALDataset::ExecuteSQL because it checks for non-nullity....
    // *sigh*
    return otb::ogr::Layer(nullptr, modifiable);
#endif
  }
  return otb::ogr::Layer(layer_ptr, *m_DataSource, modifiable);
}

void otb::ogr::DataSource::SetLayerCreationOptions(const std::vector<std::string>& options)
{
  FileNameHelperType::Pointer helper = FileNameHelperType::New();
  helper->SetGDALLayerOptions(options);
  m_LayerOptions = helper->GetGDALLayerOptions();
  // perf : do we move code from  helper->SetGDALLayerOptions in here?
}

void otb::ogr::DataSource::AddLayerCreationOptions(std::vector<std::string> options)
{
  FileNameHelperType::Pointer helper = FileNameHelperType::New();
  helper->SetGDALLayerOptions(m_LayerOptions);
  helper->AddGDALLayerOptions(options);
  m_LayerOptions = helper->GetGDALLayerOptions();
  // perf : do we move code from helper->AddGDALLayerOptions in here?
}

const std::vector<std::string>& otb::ogr::DataSource::GetLayerCreationOptions() const
{
  return m_LayerOptions;
}
/*===========================================================================*/
/*===============================[ features ]================================*/
/*===========================================================================*/
namespace
{ // Anonymous namespace
  /**\ingroup GeometryInternals
   * \brief %Functor used to accumulate the sizes of the layers in a \c DataSource.
   * \since OTB v 3.14.0
   */
struct AccuLayersSizes
{
  AccuLayersSizes(bool doForceComputation) : m_doForceComputation(doForceComputation)
  {
  }
  int operator()(int accumulated, otb::ogr::Layer const& layer) const
  {
    const int loc_size = layer.GetFeatureCount(m_doForceComputation);
    return loc_size < 0 ? loc_size : loc_size + accumulated;
  }

private:
  bool m_doForceComputation;
};
} // Anonymous namespace

int otb::ogr::DataSource::Size(bool doForceComputation) const
{
  return std::accumulate(begin(), end(), 0, AccuLayersSizes(doForceComputation));
}

/*===========================================================================*/
/*=================================[ Misc ]==================================*/
/*===========================================================================*/

OGREnvelope otb::ogr::DataSource::GetGlobalExtent(bool force /* = false */, std::string* outwkt) const
{
  assert(m_DataSource && "Datasource not initialized");
  const_iterator lit = this->begin();

  if (lit == this->end())
  {
    itkGenericExceptionMacro(<< "Cannot compute global extent because there "
                                "are no layers in the DataSource");
  }

  const OGRSpatialReference* ref_srs = lit->GetSpatialRef();
  OGREnvelope                sExtent = lit->GetExtent(force);

  if (outwkt)
  {
    *outwkt = lit->GetProjectionRef();
  }

  ++lit;

  for (; lit != this->end(); ++lit)
  {
    OGREnvelope cExtent = lit->GetExtent(force); // may throw

    const OGRSpatialReference* current_srs = lit->GetSpatialRef();

    // If both srs are valid and if they are different
    if (ref_srs && current_srs && current_srs->IsSame(ref_srs) == 0)
    {
      // Reproject cExtent in ref_srs
      // OGRCreateCoordinateTransformation is not const-correct
      OGRCoordinateTransformation* coordTransformation =
          OGRCreateCoordinateTransformation(const_cast<OGRSpatialReference*>(current_srs), const_cast<OGRSpatialReference*>(ref_srs));

      coordTransformation->Transform(1, &cExtent.MinX, &cExtent.MinY);
      coordTransformation->Transform(1, &cExtent.MaxX, &cExtent.MaxY);

      const double real_minx = std::min(cExtent.MinX, cExtent.MaxX);
      const double real_miny = std::min(cExtent.MinY, cExtent.MaxY);
      const double real_maxx = std::max(cExtent.MinX, cExtent.MaxX);
      const double real_maxy = std::max(cExtent.MinY, cExtent.MaxY);

      cExtent.MinX = real_minx;
      cExtent.MinY = real_miny;
      cExtent.MaxX = real_maxx;
      cExtent.MaxY = real_maxy;

      OGRCoordinateTransformation::DestroyCT(coordTransformation);
    }
    // else: If srs are invalid, we assume that extent are coherent

    // Merge with previous layers' extent
    sExtent.Merge(cExtent);
  } // for each layer

  return sExtent;
}

std::string otb::ogr::DataSource::GetGlobalExtent(double& ulx, double& uly, double& lrx, double& lry, bool force) const
{
  std::string       outwkt;
  const OGREnvelope sExtent = GetGlobalExtent(force, &outwkt);
  ulx                       = sExtent.MinX;
  uly                       = sExtent.MinY;
  lrx                       = sExtent.MaxX;
  lry                       = sExtent.MaxY;

  return outwkt;
}

/*virtual*/
void otb::ogr::DataSource::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  assert(m_DataSource && "Datasource not initialized");
  BOOST_FOREACH (Layer const& l, *this)
  {
    l.PrintSelf(os, indent);
  }
}

/*virtual*/ void otb::ogr::DataSource::Graft(const itk::DataObject* itkNotUsed(data))
{
  assert(!"Disabled to check if it makes sense...");
}

bool otb::ogr::DataSource::HasCapability(std::string const& capabilityName) const
{
  assert(m_DataSource && "Datasource not initialized");
  return m_DataSource->TestCapability(capabilityName.c_str());
}

void otb::ogr::DataSource::SyncToDisk()
{
  assert(m_DataSource && "Datasource not initialized");
  m_DataSource->FlushCache();
}

std::string otb::ogr::DataSource::GetDatasetDescription() const
{
  std::vector<std::string> files       = otb::ogr::GetFileListAsStringVector(m_DataSource);
  std::string              description = "";
  for (std::vector<std::string>::const_iterator it = files.begin(); it != files.end(); ++it)
    description += (*it) + ", ";

  return description;
}

ext: cxx

max_stars_repo_name: yyxgiser/OTB

max_issues_repo_name: yyxgiser/OTB

max_forks_repo_name: yyxgiser/OTB


text:
#pragma once
#include <boost/program_options.hpp>
#include <string>
#include <iostream>
#include <iomanip>
#include <math.h>
#include <bits/stdc++.h>

/*
README
To avoid confusion in indexing:
  x goes along heigth
  y goes along width
  z goes along depth

All triple nested loops go:
for x in height
  for y in width
    for z in depth

And the 3D dimension is organized as [h, w, d]
*/

namespace utils {

  struct Options {
    // Command line arguments (with default values)
    unsigned num_ipus;
    unsigned num_iterations;
    float alpha;
    std::size_t height;
    std::size_t width;
    std::size_t depth;
    std::string vertex;
    bool cpu;
    // Not command line arguments
    std::size_t side;
    std::size_t tiles_per_ipu = 0;
    std::size_t num_tiles_available = 0;
    std::size_t halo_volume = 0;
    std::vector<std::size_t> splits = {0,0,0};
    std::vector<std::size_t> smallest_slice = {std::numeric_limits<size_t>::max(),1,1};
    std::vector<std::size_t> largest_slice = {0,0,0};
  };

  inline
  Options parseOptions(int argc, char** argv) {
    Options options;
    namespace po = boost::program_options;
    po::options_description desc("Flags");
    desc.add_options()
    ("help", "Show command help.")
    (
      "num-ipus",
      po::value<unsigned>(&options.num_ipus)->default_value(1),
      "Number of IPUs (must be a power of 2)"
    )
    (
      "num-iterations",
      po::value<unsigned>(&options.num_iterations)->default_value(100000),
      "PDE: number of iterations to execute on grid."
    )
    (
      "height",
      po::value<std::size_t>(&options.height)->default_value(0),
      "Heigth of a custom 3D grid"
    )
    (
      "width",
      po::value<std::size_t>(&options.width)->default_value(0),
      "Width of a custom 3D grid"
    )
    (
      "depth",
      po::value<std::size_t>(&options.depth)->default_value(0),
      "Depth of a custom 3D grid"
    )
    (
      "alpha",
      po::value<float>(&options.alpha)->default_value(0.1),
      "PDE: update step size given as a float."
    )
    (
      "vertex",
      po::value<std::string>(&options.vertex)->default_value("HeatEquationOptimized"),
      "Name of vertex (from codelets.cpp) to use for the computation."
    )
    (
      "cpu",
      po::bool_switch(&options.cpu)->default_value(false),
      "Also perform CPU execution to control results from IPU."
    ); // NOTE: remember to remove this semicolon if more options are added in future
    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);
    if (vm.count("help")) {
      std::cout << desc << "\n";
      throw std::runtime_error("Show help");
    }
    po::notify(vm);
    return options;
  }

} // End of namespace Utils

poplar::Device getDevice(utils::Options &options) {
  /* return a Poplar device with the desired number of IPUs */
  std::size_t n = options.num_ipus;
  if (n!=1 && n!=2 && n!=4 && n!=8 && n!=16 && n!=32 && n!=64)
    throw std::runtime_error("Invalid number of IPUs.");

  // Create device manager
  auto manager = poplar::DeviceManager::createDeviceManager();
  auto devices = manager.getDevices(poplar::TargetType::IPU, n);

  // Use the first available device
  for (auto &device : devices)
    if (device.attach())
      return std::move(device);

  throw std::runtime_error("No hardware device available.");
}

std::size_t side_length(std::size_t num_ipus, std::size_t base_length) {
  std::size_t log2_num_ipus = log(num_ipus) / log(2);
  std::size_t side = base_length*pow(1.26, log2_num_ipus);
  return side;
}

inline float randomFloat() {
  return static_cast <float> (rand() / static_cast <float> (RAND_MAX));
}

inline static unsigned index(unsigned x, unsigned y, unsigned z, unsigned width, unsigned depth) {
  return (z) + (y)*(depth) + (x)*(width)*(depth);
}

inline static unsigned block_low(unsigned id, unsigned p, unsigned n) {
  return (id*n)/p;
}

inline static unsigned block_high(unsigned id, unsigned p, unsigned n) {
  return block_low(id+1, p, n);
}

inline static unsigned block_size(unsigned id, unsigned p, unsigned n) {
  return block_high(id, p, n) - block_low(id, p, n);
}

std::size_t volume(std::vector<std::size_t> shape) {
  // return volume of shape vector (3D)
  return shape[0]*shape[1]*shape[2];
}

void workDivision(utils::Options &options) {
  /* Function UPDATES options.splits
   * 1) all tiles will be used, hence options.num_tiles_available must
   *    previously be updated by using the target object
   * 2) the average resulting slice will have shape [height/nh, width/nw, depth/nd'] and
   *    this function chooses nh, nw, nd, so that the surface area is minimized.
   */
  float smallest_surface_area = std::numeric_limits<float>::max();
  std::size_t height = (options.height - 2);
  std::size_t width = (options.width - 2);
  std::size_t depth = (options.depth - 2) / options.num_ipus;
  std::size_t tile_count = options.num_tiles_available / options.num_ipus;
  for (std::size_t i = 1; i <= tile_count; ++i) {
    if (tile_count % i == 0) { // then i is a factor
      // Further, find two other factors, to obtain exactly three factors
      std::size_t other_factor = tile_count/i;
      for (std::size_t j = 1; j <= other_factor; ++j) {
        if (other_factor % j == 0) { // then j is a second factor
          std::size_t k = other_factor/j; // and k is the third factor
          std::vector<std::size_t> splits = {i,j,k};
          if (i*j*k != tile_count) {
            throw std::runtime_error("workDivision(), factorization does not work.");
          }
          for (std::size_t l = 0; l < 3; ++l) {
            for (std::size_t m = 0; m < 3; ++m) {
              for (std::size_t n = 0; n < 3; ++n) {
                if (l != m && l != n && m != n) {
                  float slice_height = float(height)/float(splits[l]);
                  float slice_width = float(width)/float(splits[m]);
                  float slice_depth = float(depth)/float(splits[n]);
                  float surface_area = 2.0*(slice_height*slice_width + slice_depth*slice_width + slice_depth*slice_height);
                  if (surface_area <= smallest_surface_area) {
                    smallest_surface_area = surface_area;
                    options.splits = splits;
                  }
                }
              }
            }
          }
        }
      }
    }
  }
}

std::vector<float> heatEquationCpu(
  const std::vector<float> initial_values,
  const utils::Options &options) {
  /*
   * Heat Equation Vertex, computed on the CPU, unparallelized.
   * The purpose of this function is to serve as the "true"
   * solution.
   * NOTE: can be very slow for large grids/large number of iterations
   */
  const float beta = 1.0 - 6.0*options.alpha;
  unsigned h = options.height;
  unsigned w = options.width;
  unsigned d = options.depth;
  unsigned iter = options.num_iterations;
  std::vector<float> a(initial_values.size());
  std::vector<float> b(initial_values.size());

  // initial copy to include edges
  for (std::size_t i = 0; i < initial_values.size(); ++i)  {
    a[i] = initial_values[i];
    b[i] = initial_values[i];
  }

  // Heat Equation iterations
  for (std::size_t t = 0; t < iter; ++t) {
    for (std::size_t x = 1; x < h - 1; ++x) {
      for (std::size_t y = 1; y < w - 1; ++y) {
        for (std::size_t z = 1; z < d - 1; ++z) {
          a[index(x,y,z,w,d)] = beta*b[index(x,y,z,w,d)] +
            options.alpha*(
              b[index(x+1,y,z,w,d)] +
              b[index(x-1,y,z,w,d)] +
              b[index(x,y+1,z,w,d)] +
              b[index(x,y-1,z,w,d)] +
              b[index(x,y,z+1,w,d)] +
              b[index(x,y,z-1,w,d)]
            );
        }
      }
    }
    for (std::size_t x = 1; x < h - 1; ++x) {
      for (std::size_t y = 1; y < w - 1; ++y) {
        for (std::size_t z = 1; z < d - 1; ++z) {
          b[index(x,y,z,w,d)] = a[index(x,y,z,w,d)];
        }
      }
    }
  }
  return a;
}

void printMeanSquaredError(
  std::vector<float> a,
  std::vector<float> b,
  utils::Options &options) {
  /*
   * Compute the MSE, __only the inner elements__, of two 3D grids
   */
  double squared_error = 0, diff;
  std::size_t h = options.height;
  std::size_t w = options.width;
  std::size_t d = options.depth;
  for (std::size_t x = 1; x < h - 1; ++x) {
    for (std::size_t y = 1; y < w - 1; ++y) {
      for (std::size_t z = 1; z < d - 1; ++z) {
        diff = double(a[index(x,y,z,w,d)] - b[index(x,y,z,w,d)]);
        squared_error += diff*diff;
      }
    }
  }
  double mean_squared_error = squared_error / (double) ((h-2)*(w-2)*(d-2));

  std::cout << "\nMean Squared Error (IPU vs. CPU) = " << mean_squared_error;
  if (mean_squared_error == double(0.0))
    std::cout << " (exactly)";
  std::cout << "\n";
}

void printResults(utils::Options &options, double wall_time) {

  // Calculate metrics
  double inner_volume = (double) options.height * (double) options.width * (double) options.depth;
  double flops_per_element = 8.0;
  double flops = inner_volume * options.num_iterations * flops_per_element / wall_time;
  double internal_communication_ops = 2.0*(double)options.halo_volume*(double)options.num_ipus;
  double external_communication_ops = 4.0*(double)options.height*(double)options.width*(options.num_ipus - 1.0); // 2 load and 2 stores of a slice (partition along depth)
  double bandwidth = (7.0*inner_volume + internal_communication_ops + external_communication_ops)*(double)options.num_iterations*sizeof(float)/wall_time;
  double tflops = flops*1e-12;
  double bandwidth_TB_s = bandwidth*1e-12;

  std::cout << "3D Isotropic Diffusion"
    << "\n----------------------"
    << "\nVertex             = " << options.vertex
    << "\nNo. IPUs           = " << options.num_ipus
    << "\nNo. Tiles          = " << options.num_tiles_available
    << "\nTotal Grid         = " << options.height << "*" << options.width << "*" << options.depth << " = "
                                 << options.height*options.width*options.depth*1e-6 << " million elements"
    << "\nSmallest Sub-grid  = " << options.smallest_slice[0] << "*" << options.smallest_slice[1] << "*" << options.smallest_slice[2]
    << "\nLargest Sub-grid   = " << options.largest_slice[0] << "*" << options.largest_slice[1] << "*" << options.largest_slice[2]
    << "\nalpha              = " << options.alpha
    << "\nNo. Iterations     = " << options.num_iterations
    << "\n"
    << "\nLaTeX Tabular Row"
    << "\n-----------------"
    << "\nNo. IPUs & Grid & No. Iterations & Time [s] & Throughput [TFLOPS] & Minimum Bandwidth [TB/s] \\\\\n"
    << options.num_ipus << " & "
    << "$" << options.height << "\\times " << options.width << "\\times " << options.depth << "$ & "
    << options.num_iterations << " & " << std::fixed
    << std::setprecision(2) << wall_time << " & "
    << std::setprecision(2) << tflops << " & "
    << std::setprecision(2) << bandwidth_TB_s << " \\\\"
    << "\n";
}

void printMultiIpuGridInfo(std::size_t base_length) {
  float base_volume = 0;

  std::cout
    << "\\begin{tabular}{ccc}"
    << "\n\\toprule"
    << "\nNumber of IPUs & Grid Shape & Relative Volume \\\\\\midrule"
    << std::fixed << std::setprecision(2);

  for (std::size_t i = 1; i <= 64; i*=2) {
    std::size_t side = side_length(i, base_length);
    float volume = float(side)*float(side)*float(side);
    if (i == 1) base_volume = volume;
    std::cout << "\n"
      << i << " & $"
      << side << "\\times " << side << "\\times " << side << "$ & "
      << volume/base_volume << "x \\\\";
  }

  std::cout
    << "\n\\bottomrule"
    << "\n\\end{tabular}\n";
}

ext: hpp

max_stars_repo_name: simehaa/IPU

max_issues_repo_name: simehaa/IPU

max_forks_repo_name: simehaa/IPU


text:
//***************************************************************************
//* Copyright (c) 2015 Saint Petersburg State University
//* Copyright (c) 2011-2014 Saint Petersburg Academic University
//* All Rights Reserved
//* See file LICENSE for details.
//***************************************************************************

#ifndef __HAMMER_IT_CONSENSUS_HPP__
#define __HAMMER_IT_CONSENSUS_HPP__

#include "HSeq.hpp"

#include <boost/numeric/ublas/matrix.hpp>
#include <limits>

namespace hammer {
namespace iontorrent {

inline std::pair<hammer::HomopolymerRun, double> consensus(
    const boost::numeric::ublas::matrix<double>& scores) {
  double inf = -std::numeric_limits<double>::infinity();

  double max = inf;
  uint8_t nucl = 0;
  uint8_t len = 1;
  for (uint8_t j = 0; j < 4; ++j)
    for (uint8_t k = 1; k < 64; ++k)
      if (scores(j, k) > max) {
        nucl = j;
        len = k;
        max = scores(j, k);
      }

  return std::make_pair(hammer::HomopolymerRun(nucl, len), max);
}

};  // namespace iontorrent
};  // namespace hammer

#endif  // __HAMMER_IT_CONSENSUS_HPP__

max_stars_repo_name: STRIDES-Codes/Exploring-the-Microbiome-

max_issues_repo_name: STRIDES-Codes/Exploring-the-Microbiome-

max_forks_repo_name: STRIDES-Codes/Exploring-the-Microbiome-


text:
// Copyright (c) 2015-2018, CNRS
// Authors: Justin Carpentier <jcarpent@laas.fr>

#ifndef __multicontact_api_scenario_contact_model_planar_hpp__
#define __multicontact_api_scenario_contact_model_planar_hpp__

#include "multicontact-api/scenario/fwd.hpp"
#include "multicontact-api/serialization/archive.hpp"
#include "multicontact-api/serialization/eigen-matrix.hpp"
#include <iostream>
#include <Eigen/Dense>
#include <pinocchio/spatial/skew.hpp>

namespace multicontact_api {
namespace scenario {

template <typename _Scalar>
struct ContactModelTpl : public serialization::Serializable<ContactModelTpl<_Scalar> > {
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW

  typedef _Scalar Scalar;
  typedef Eigen::Matrix<Scalar, 3, 3> Matrix3;
  typedef Eigen::Matrix<Scalar, 3, Eigen::Dynamic> Matrix3X;
  typedef Eigen::Matrix<Scalar, 6, Eigen::Dynamic> Matrix6X;

  /// \brief Default constructor.
  ContactModelTpl()
      : m_mu(-1.),
        m_contact_type(ContactType::CONTACT_UNDEFINED),
        m_num_contact_points(1),
        m_contact_points_positions(Matrix3X::Zero(3, 1)) {}

  /// \brief Constructor with friction
  ContactModelTpl(const Scalar mu)
      : m_mu(mu),
        m_contact_type(ContactType::CONTACT_UNDEFINED),
        m_num_contact_points(1),
        m_contact_points_positions(Matrix3X::Zero(3, 1)) {}

  /// \brief Full constructor
  ContactModelTpl(const Scalar mu, const ContactType contact_type)
      : m_mu(mu),
        m_contact_type(contact_type),
        m_num_contact_points(1),
        m_contact_points_positions(Matrix3X::Zero(3, 1)) {}

  /// \brief Copy constructor
  template <typename S2>
  explicit ContactModelTpl(const ContactModelTpl<S2>& other)
      : m_mu(other.mu),
        m_contact_type(other.m_contact_type),
        m_contact_points_positions(other.m_contact_points_positions),
        m_num_contact_points(other.m_num_contact_points) {}

  template <typename S2>
  bool operator==(const ContactModelTpl<S2>& other) const {
    return m_mu == other.m_mu && m_contact_type == other.m_contact_type &&
           m_num_contact_points == other.m_num_contact_points &&
           m_contact_points_positions == other.m_contact_points_positions;
  }

  template <typename S2>
  bool operator!=(const ContactModelTpl<S2>& other) const {
    return !(*this == other);
  }

  void disp(std::ostream& os) const {
    os << "ContactType: " << m_contact_type << ", mu: " << m_mu << std::endl
       << "Number of contact points: " << m_num_contact_points << ", positions: " << std::endl
       << m_contact_points_positions << std::endl;
  }

  template <typename S2>
  friend std::ostream& operator<<(std::ostream& os, const ContactModelTpl<S2>& cp) {
    cp.disp(os);
    return os;
  }

  size_t num_contact_points() const { return m_num_contact_points; }

  void num_contact_points(const size_t num) {
    m_num_contact_points = num;
    m_contact_points_positions = Matrix3X::Zero(3, num);
  }

  Matrix3X contact_points_positions() const { return m_contact_points_positions; }

  void contact_points_positions(const Matrix3X& positions) {
    m_contact_points_positions = positions;
    m_num_contact_points = positions.cols();
  }

  /**
   * @brief generatorMatrix Return a 6x(num_contact_points*3) matrix
   * containing the generator used to compute contact forces.
   * @return
   */
  Matrix6X generatorMatrix() const {
    Matrix6X gen = Matrix6X::Zero(6, m_num_contact_points * 3);
    for (size_t i = 0; i < m_num_contact_points; i++) {
      gen.block(0, i * 3, 3, 3) = Matrix3::Identity();
      gen.block(3, i * 3, 3, 3) = pinocchio::skew(m_contact_points_positions.col(i));
    }
    return gen;
  }

  /// \brief Friction coefficient.
  Scalar m_mu;
  /// \brief ZMP radius.
  ContactType m_contact_type;

 private:
  /// \brief The number of contact points used to model this contact
  size_t m_num_contact_points;
  /// \brief 3xnum_contact_points matrix defining the contact points positions in the frame of the contact placement
  Matrix3X m_contact_points_positions;

  // Serialization of the class
  friend class boost::serialization::access;

  template <class Archive>
  void save(Archive& ar, const unsigned int /*version*/) const {
    ar& boost::serialization::make_nvp("mu", m_mu);
    ar& boost::serialization::make_nvp("contact_type", m_contact_type);
    ar& boost::serialization::make_nvp("num_contact_points", m_num_contact_points);
    ar& boost::serialization::make_nvp("contact_points_positions", m_contact_points_positions);
  }

  template <class Archive>
  void load(Archive& ar, const unsigned int /*version*/) {
    ar >> boost::serialization::make_nvp("mu", m_mu);
    ar >> boost::serialization::make_nvp("contact_type", m_contact_type);
    ar >> boost::serialization::make_nvp("num_contact_points", m_num_contact_points);
    ar >> boost::serialization::make_nvp("contact_points_positions", m_contact_points_positions);
  }

  BOOST_SERIALIZATION_SPLIT_MEMBER()
};
}  // namespace scenario
}  // namespace multicontact_api

DEFINE_CLASS_TEMPLATE_VERSION(typename Scalar, multicontact_api::scenario::ContactModelTpl<Scalar>)

#endif  // ifndef __multicontact_api_scenario_contact_model_planar_hpp__

max_stars_repo_name: proyan/multicontact-api

max_issues_repo_name: proyan/multicontact-api

max_forks_repo_name: proyan/multicontact-api


text:
/*    Copyright (c) 2010-2019, Delft University of Technology
 *    All rigths reserved
 *
 *    This file is part of the Tudat. Redistribution and use in source and
 *    binary forms, with or without modification, are permitted exclusively
 *    under the terms of the Modified BSD license. You should have received
 *    a copy of the license with this file. If not, please or visit:
 *    http://tudat.tudelft.nl/LICENSE.
 */

#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MAIN

#include <limits>

#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/random.hpp>

#include "tudat/math/statistics/boostProbabilityDistributions.h"

namespace tudat
{
namespace unit_tests
{

BOOST_AUTO_TEST_SUITE( test_boost_distributions )

//! Test wrappers for continuous probability distributions from boost.
BOOST_AUTO_TEST_CASE( testContinuousBoostDistribution )
{

    // Set independent variables.
    std::vector< double > testVector;
    std::vector< double > positiveTestVector;
    for( unsigned int i = 0; i < 1001; i++ )
    {
        testVector.push_back( -5.0 + static_cast< double >( i ) * 0.01 );
        positiveTestVector.push_back( static_cast< double >( i ) * 0.01 );
    }

    // Set probability vector (used for quantile computation).
    std::vector< double > probabilityVector;
    for( unsigned int i = 1; i < 100; i++ )
    {
        probabilityVector.push_back( static_cast< double >( i ) * 0.01 );
    }

    // Test uniform_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( -1.0 );
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::uniform_boost_distribution, parameters );
        boost::math::uniform_distribution< > manualDistribution( parameters.at( 0 ), parameters.at( 1 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < testVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( testVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, testVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( testVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, testVector.at( i ) ) );
        }

        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }
    }

    // Test normal_boost_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( -1.0 );
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::normal_boost_distribution, parameters );
        boost::math::normal_distribution< > manualDistribution( parameters.at( 0 ), parameters.at( 1 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < testVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( testVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, testVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( testVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, testVector.at( i ) ) );
        }

        // Test inverse cdf (quantile).
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }

    }

    // Test exponential_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::exponential_boost_distribution, parameters );
        boost::math::exponential_distribution< > manualDistribution( parameters.at( 0 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < testVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( positiveTestVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( positiveTestVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
        }

        // Test inverse cdf (quantile).
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }

    }

    // Test gamma_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( 1.0 );
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::gamma_boost_distribution, parameters );
        boost::math::gamma_distribution< > manualDistribution( parameters.at( 0 ), parameters.at( 1 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < testVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( positiveTestVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( positiveTestVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
        }

        // Test inverse cdf (quantile).
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }
    }

    // Test lognormal_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( -1.0 );
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::lognormal_boost_distribution, parameters );
        boost::math::lognormal_distribution< > manualDistribution( parameters.at( 0 ), parameters.at( 1 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < testVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( positiveTestVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( positiveTestVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, positiveTestVector.at( i ) ) );
        }

        // Test inverse cdf (quantile).
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }
    }

    // Test beta_boost_distribution
    {
        std::vector< double > parameters;
        parameters.push_back( 1.0 );
        parameters.push_back( 2.5 );

        std::shared_ptr< statistics::InvertibleContinuousProbabilityDistribution< double > > randomVariable =
                statistics::createBoostRandomVariable(
                    statistics::beta_boost_distribution, parameters );
        boost::math::beta_distribution< > manualDistribution( parameters.at( 0 ), parameters.at( 1 ) );

        // Test pdf and cdf.
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluatePdf( probabilityVector.at( i ) ),
                        boost::math::pdf< double >( manualDistribution, probabilityVector.at( i ) ) );
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateCdf( probabilityVector.at( i ) ),
                        boost::math::cdf< double >( manualDistribution, probabilityVector.at( i ) ) );
        }

        // Test inverse cdf (quantile).
        for( unsigned int i = 0; i < probabilityVector.size( ); i++ )
        {
            BOOST_CHECK_EQUAL(
                        randomVariable->evaluateInverseCdf( probabilityVector.at( i ) ),
                        boost::math::quantile< double >( manualDistribution, probabilityVector.at( i ) ) );
        }
    }
}

BOOST_AUTO_TEST_SUITE_END( )

} // namespace unit_tests
} // namespace tudat

ext: cpp


text:
/**
 * $Id$
 *
 * Software License Agreement (GNU General Public License)
 *
 *  Copyright (C) 2015:
 *
 *    Johann Prankl, prankl@acin.tuwien.ac.at
 *    Aitor Aldoma, aldoma@acin.tuwien.ac.at
 *
 *      Automation and Control Institute
 *      Vienna University of Technology
 *      Gusshausstraße 25-29
 *      1170 Vienn, Austria
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 * @author Johann Prankl, Aitor Aldoma
 *
 */
#ifndef EIGEN_BOOST_SERIALIZATION
#define EIGEN_BOOST_SERIALIZATION
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <boost/serialization/split_free.hpp>
#include <boost/serialization/vector.hpp>

// A macro to get rid of the unused warning
#ifndef UNUSED
#define UNUSED(expr) \
  do {               \
    (void)(expr);    \
  } while (0)
#endif

namespace boost {
namespace serialization {
template <class Archive, typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void save(Archive& ar, const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& m,
          const unsigned int version) {
  UNUSED(version);
  int rows = m.rows(), cols = m.cols();
  ar& rows;
  ar& cols;
  ar& boost::serialization::make_array(m.data(), rows * cols);
}
template <class Archive, typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void load(Archive& ar, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& m,
          const unsigned int version) {
  UNUSED(version);
  int rows, cols;
  ar& rows;
  ar& cols;
  m.resize(rows, cols);
  ar& boost::serialization::make_array(m.data(), rows * cols);
}

template <class Archive, typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
void serialize(Archive& ar, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& m,
               const unsigned int version) {
  split_free(ar, m, version);
}

/*template <class Archive, typename _Scalar>
    void save(Archive & ar, const Eigen::Triplet<_Scalar> & m, const unsigned int version) {
        ar & m.row();
        ar & m.col();
        ar & m.value();
    }
template <class Archive, typename _Scalar>
    void load(Archive & ar, Eigen::Triplet<_Scalar> & m, const unsigned int version) {
        int row,col;
        _Scalar value;
        ar & row;
        ar & col;
        ar & value;
        m = Eigen::Triplet<_Scalar>(row,col,value);
    }

template <class Archive, typename _Scalar>
    void serialize(Archive & ar, Eigen::Triplet<_Scalar> & m, const unsigned int version) {
        split_free(ar,m,version);
    }

template <class Archive, typename _Scalar, int _Options,typename _Index>
    void save(Archive & ar, const Eigen::SparseMatrix<_Scalar,_Options,_Index> & m, const unsigned int version) {
        int innerSize=m.innerSize();
        int outerSize=m.outerSize();
        typedef typename Eigen::Triplet<_Scalar> Triplet;
        std::vector<Triplet> triplets;
        for(int i=0; i < outerSize; ++i) {
            for(typename Eigen::SparseMatrix<_Scalar,_Options,_Index>::InnerIterator it(m,i); it; ++it) {
            triplets.push_back(Triplet(it.row(), it.col(), it.value()));
            }
        }
        ar & innerSize;
        ar & outerSize;
        ar & triplets;
    }
template <class Archive, typename _Scalar, int _Options, typename _Index>
    void load(Archive & ar, Eigen::SparseMatrix<_Scalar,_Options,_Index>  & m, const unsigned int version) {
        int innerSize;
        int outerSize;
        ar & innerSize;
        ar & outerSize;
        int rows = m.IsRowMajor?outerSize:innerSize;
        int cols = m.IsRowMajor?innerSize:outerSize;
        m.resize(rows,cols);
        typedef typename Eigen::Triplet<_Scalar> Triplet;
        std::vector<Triplet> triplets;
        ar & triplets;
        m.setFromTriplets(triplets.begin(), triplets.end());

    }
template <class Archive, typename _Scalar, int _Options, typename _Index>
    void serialize(Archive & ar, Eigen::SparseMatrix<_Scalar,_Options,_Index> & m, const unsigned int version) {
        split_free(ar,m,version);
    }
*/
}  // namespace serialization
}  // namespace boost
#endif


text:
//  Copyright (c) 2020 Andrey Semashev
//
//  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 <boost/atomic/atomic_ref.hpp>
#include <boost/atomic/capabilities.hpp>

#include <boost/config.hpp>
#include <boost/cstdint.hpp>

#include "wait_test_helpers.hpp"

int main(int, char *[])
{
    test_wait_notify_api< atomic_ref_wrapper, boost::uint8_t >(1, 2, 3, BOOST_ATOMIC_HAS_NATIVE_INT8_WAIT_NOTIFY);
    test_wait_notify_api< atomic_ref_wrapper, boost::uint16_t >(1, 2, 3, BOOST_ATOMIC_HAS_NATIVE_INT16_WAIT_NOTIFY);
    test_wait_notify_api< atomic_ref_wrapper, boost::uint32_t >(1, 2, 3, BOOST_ATOMIC_HAS_NATIVE_INT32_WAIT_NOTIFY);
    test_wait_notify_api< atomic_ref_wrapper, boost::uint64_t >(1, 2, 3, BOOST_ATOMIC_HAS_NATIVE_INT64_WAIT_NOTIFY);
#if defined(BOOST_HAS_INT128) && !defined(BOOST_ATOMIC_TESTS_NO_INT128)
    test_wait_notify_api< atomic_ref_wrapper, boost::uint128_type >(1, 2, 3, BOOST_ATOMIC_HAS_NATIVE_INT128_WAIT_NOTIFY);
#endif

    {
        struct_3_bytes s1 = {{ 1 }};
        struct_3_bytes s2 = {{ 2 }};
        struct_3_bytes s3 = {{ 3 }};
        test_wait_notify_api< atomic_ref_wrapper, struct_3_bytes >(s1, s2, s3);
    }
    {
        large_struct s1 = {{ 1 }};
        large_struct s2 = {{ 2 }};
        large_struct s3 = {{ 3 }};
        test_wait_notify_api< atomic_ref_wrapper, large_struct >(s1, s2, s3);
    }

    return boost::report_errors();
}


text:
#include "quote.hpp"

#include <boost/algorithm/string/join.hpp>

namespace configure {

	static char const* unix_special_characters = "`^~<>|;()$*'\\&# \t\n\r";

	static char const* windows_special_characters = "^<>|'&# \t\n\r";

	template<CommandParser target>
	static bool needs_quote(char const c)
	{
		char const* tab = nullptr;
		switch (target)
		{
		case CommandParser::unix_shell:
		case CommandParser::make:
			tab = unix_special_characters;
			break;
		case CommandParser::windows_shell:
		case CommandParser::nmake:
			tab = windows_special_characters;
			break;
		default:
			std::abort();
		}

		while (*tab != '\0' && *tab != c)
			tab++;
		return *tab != '\0';
	}

	template<CommandParser target>
	std::string quote_arg(std::string const& arg)
	{
		bool enclosing_quotes = std::any_of(arg.begin(), arg.end(), &needs_quote<target>);

		std::string res;
		res.reserve(arg.size() + 2);
		if (enclosing_quotes)
			res.push_back('"');

		int backslashes = 0;
		for (auto c: arg)
		{
			if ((target == CommandParser::unix_shell ||
			     target == CommandParser::make) &&
			    (c == '\\' || c == '\'' || c == '`' || c == '$' || c == '"'))
				res.push_back('\\');
			if (target == CommandParser::windows_shell ||
			    target == CommandParser::nmake)
			{
				if (c == '\\') backslashes += 1;
				else if (c == '"')
				{
					while (backslashes > 0)
					{
						backslashes -= 1;
						res.push_back('\\');
					}
					res.push_back('\\');
				}
				else backslashes = 0;
			}

			if (c == '$')
			{
				if (target == CommandParser::make)
					res.push_back('$');
			}

			if (target == CommandParser::nmake ||
			    target == CommandParser::windows_shell)
			{
				if (c == '%')
					res.push_back('%');
			}
			res.push_back(c);
		}

		if (enclosing_quotes)
		{
			while (backslashes > 0)
			{
				backslashes -= 1;
				res.push_back('\\');
			}
			res.push_back('"');
		}
		return res;
	}

	template<CommandParser target>
	std::string quote(std::vector<std::string> const& cmd)
	{
		std::vector<std::string> res;
		for (auto& arg: cmd)
			res.push_back(quote_arg<target>(arg));
		return boost::join(res, " ");
	}

#define INSTANCIATE(e) \
	template \
	std::string quote<CommandParser::e>(std::vector<std::string> const&); \
/**/
	INSTANCIATE(unix_shell);
	INSTANCIATE(windows_shell);
	INSTANCIATE(make);
	INSTANCIATE(nmake);
#undef INSTANCIATE

	std::string quote_arg(CommandParser target, std::string const& arg)
	{
		switch (target)
		{
#define CASE(value) \
		case CommandParser::value: \
			return quote_arg<CommandParser::value>(arg)
		CASE(unix_shell);
		CASE(windows_shell);
		CASE(nmake);
		CASE(make);
#undef CASE
		default:
			std::abort();
		}
	}
}

max_stars_repo_name: hotgloupi/configure

max_issues_repo_name: hotgloupi/configure

max_forks_repo_name: hotgloupi/configure


text:
// g2o - General Graph Optimization
// Copyright (C) 2012 R. Kümmerle
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// * Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
// * 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
// HOLDER 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 <Eigen/Core>
#include <Eigen/Geometry>
#include <iostream>

#include "g2o/core/auto_differentiation.h"
#include "g2o/core/base_unary_edge.h"
#include "g2o/core/base_vertex.h"
#include "g2o/core/optimization_algorithm_factory.h"
#include "g2o/core/sparse_optimizer.h"
#include "g2o/stuff/command_args.h"
#include "g2o/stuff/sampler.h"

using namespace std;

G2O_USE_OPTIMIZATION_LIBRARY(dense);

double errorOfSolution(int numPoints, Eigen::Vector2d* points,
                       const Eigen::Vector3d& circle) {
  Eigen::Vector2d center = circle.head<2>();
  double radius = circle(2);
  double error = 0.;
  for (int i = 0; i < numPoints; ++i) {
    double d = (points[i] - center).norm() - radius;
    error += d * d;
  }
  return error;
}

/**
 * \brief a circle located at x,y with radius r
 */
class VertexCircle : public g2o::BaseVertex<3, Eigen::Vector3d> {
 public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
  VertexCircle() {}

  virtual bool read(std::istream& /*is*/) { return false; }

  virtual bool write(std::ostream& /*os*/) const { return false; }

  virtual void setToOriginImpl() {
    cerr << __PRETTY_FUNCTION__ << " not implemented yet" << endl;
  }

  virtual void oplusImpl(const double* update) {
    Eigen::Vector3d::ConstMapType v(update);
    _estimate += v;
  }
};

/**
 * \brief measurement for a point on the circle
 *
 * Here the measurement is the point which is on the circle.
 * The error function computes the distance of the point to
 * the center minus the radius of the circle.
 */
class EdgePointOnCircle
    : public g2o::BaseUnaryEdge<1, Eigen::Vector2d, VertexCircle> {
 public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  EdgePointOnCircle() {}
  virtual bool read(std::istream& /*is*/) { return false; }
  virtual bool write(std::ostream& /*os*/) const { return false; }

  template <typename T>
  bool operator()(const T* circle, T* error) const {
    typename g2o::VectorN<2, T>::ConstMapType center(circle);
    const T& radius = circle[2];

    error[0] = (measurement().cast<T>() - center).norm() - radius;
    return true;
  }

  G2O_MAKE_AUTO_AD_FUNCTIONS  // use autodiff
};

int main(int argc, char** argv) {
  int numPoints;
  int maxIterations;
  bool verbose;
  g2o::CommandArgs arg;
  arg.param("numPoints", numPoints, 100,
            "number of points sampled from the circle");
  arg.param("i", maxIterations, 10, "perform n iterations");
  arg.param("v", verbose, false, "verbose output of the optimization process");
  arg.parseArgs(argc, argv);

  // generate random data
  Eigen::Vector2d center(4.0, 2.0);
  double radius = 2.0;
  Eigen::Vector2d* points = new Eigen::Vector2d[numPoints];

  g2o::Sampler::seedRand();
  for (int i = 0; i < numPoints; ++i) {
    double r = g2o::Sampler::gaussRand(radius, 0.05);
    double angle = g2o::Sampler::uniformRand(0.0, 2.0 * M_PI);
    points[i].x() = center.x() + r * cos(angle);
    points[i].y() = center.y() + r * sin(angle);
  }

  // setup the solver
  g2o::SparseOptimizer optimizer;
  optimizer.setVerbose(false);

  // allocate the solver
  g2o::OptimizationAlgorithmProperty solverProperty;
  optimizer.setAlgorithm(
      g2o::OptimizationAlgorithmFactory::instance()->construct("lm_dense",
                                                               solverProperty));

  // build the optimization problem given the points
  // 1. add the circle vertex
  VertexCircle* circle = new VertexCircle();
  circle->setId(0);
  circle->setEstimate(
      Eigen::Vector3d(3.0, 3.0, 3.0));  // some initial value for the circle
  optimizer.addVertex(circle);
  // 2. add the points we measured
  for (int i = 0; i < numPoints; ++i) {
    EdgePointOnCircle* e = new EdgePointOnCircle;
    e->setInformation(Eigen::Matrix<double, 1, 1>::Identity());
    e->setVertex(0, circle);
    e->setMeasurement(points[i]);
    optimizer.addEdge(e);
  }

  // perform the optimization
  optimizer.initializeOptimization();
  optimizer.setVerbose(verbose);
  optimizer.optimize(maxIterations);

  if (verbose) cout << endl;

  // print out the result
  cout << "Iterative least squares solution" << endl;
  cout << "center of the circle " << circle->estimate().head<2>().transpose()
       << endl;
  cout << "radius of the cirlce " << circle->estimate()(2) << endl;
  cout << "error " << errorOfSolution(numPoints, points, circle->estimate())
       << endl;
  cout << endl;

  // solve by linear least squares
  // Let (a, b) be the center of the circle and r the radius of the circle.
  // For a point (x, y) on the circle we have:
  // (x - a)^2 + (y - b)^2 = r^2
  // This leads to
  // (-2x -2y 1)^T * (a b c) = -x^2 - y^2   (1)
  // where c = a^2 + b^2 - r^2.
  // Since we have a bunch of points, we accumulate Eqn (1) in a matrix and
  // compute the normal equation to obtain a solution for (a b c).
  // Afterwards the radius r is recovered.
  Eigen::MatrixXd A(numPoints, 3);
  Eigen::VectorXd b(numPoints);
  for (int i = 0; i < numPoints; ++i) {
    A(i, 0) = -2 * points[i].x();
    A(i, 1) = -2 * points[i].y();
    A(i, 2) = 1;
    b(i) = -pow(points[i].x(), 2) - pow(points[i].y(), 2);
  }
  Eigen::Vector3d solution =
      (A.transpose() * A).ldlt().solve(A.transpose() * b);
  // calculate the radius of the circle given the solution so far
  solution(2) = sqrt(pow(solution(0), 2) + pow(solution(1), 2) - solution(2));
  cout << "Linear least squares solution" << endl;
  cout << "center of the circle " << solution.head<2>().transpose() << endl;
  cout << "radius of the cirlce " << solution(2) << endl;
  cout << "error " << errorOfSolution(numPoints, points, solution) << endl;

  // clean up
  delete[] points;

  return 0;
}


text:
//! \file test_expected.cpp

// Copyright Pierre Talbot 2013.
// Copyright Vicente J. Botet Escriba 2013,2014.

// Use, modification and distribution are subject to 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 <experimental/expected.hpp>

#include <boost/detail/lightweight_test.hpp>

namespace stde = std::experimental;
template <class T>
using expected_sc = stde::expected<T, std::error_code>;

struct NoDefaultConstructible
{
  NoDefaultConstructible() = delete;
  NoDefaultConstructible(int)  {}
};

struct NoCopyConstructible
{
  NoCopyConstructible()  {}
  NoCopyConstructible(NoCopyConstructible const&) = delete;
  NoCopyConstructible(NoCopyConstructible &&) = default;
};
struct NoMoveConstructible
{
  NoMoveConstructible()  {}
  NoMoveConstructible(NoMoveConstructible const&) = default;
  NoMoveConstructible(NoMoveConstructible &&) = delete;
  NoMoveConstructible& operator=(NoMoveConstructible const&) = default;
  NoMoveConstructible& operator=(NoMoveConstructible &&) = delete;
};

template<typename T, bool> struct base {
  base() = default;
  base(T v_) : v(v_) {}
  base(base const&) = default;
  base(base&&) = default;
  T v;
};

template<typename T> struct base<T, false> {
   base() = default;
   base(T v_) : v(v_) {}
   base(base const&) = default;
   base(base&&) = delete;
   T v;
};

template<typename T>
class optional : public base<T, std::is_move_constructible<T>::value> {
  using base_type = base<T, std::is_move_constructible<T>::value>;
public:
  optional() = default;
  optional(T v) : base_type(v) {}
  optional(const optional&) = default;
  optional(optional &&) = default;
};

int main()
{
  {
    expected_sc<NoDefaultConstructible> x; // FAILS as expected
  }
#if 0
  {
    NoCopyConstructible ncc;
    expected_sc<NoCopyConstructible> x{ncc}; // FAILS as expected
  }
  {
    NoCopyConstructible ncc;
    expected_sc<NoCopyConstructible> x{std::move(ncc)};
    expected_sc<NoCopyConstructible> y{x}; // FAILS as expected

  }
  {
    NoMoveConstructible nmc;
    NoMoveConstructible nmc2 = std::move(nmc); // FAILS as expected

    // fixme
    expected_sc<NoMoveConstructible> x;
    expected_sc<NoMoveConstructible> y =  std::move(x); // DOESN'T FAIL
  }
#endif
  {
    optional<NoMoveConstructible> x;
    optional<NoMoveConstructible> y {std::move(x)}; //FAILS as expected
    (void)y;
  }
  static_assert(! std::is_default_constructible<NoDefaultConstructible>::value, "");
  static_assert(! std::is_default_constructible<expected_sc<NoDefaultConstructible>>::value, "");

  static_assert(! std::is_copy_constructible<NoCopyConstructible>::value, "");
  static_assert(! std::is_constructible<expected_sc<NoCopyConstructible>, NoCopyConstructible const& >::value, "");
  static_assert(! std::is_constructible<expected_sc<NoCopyConstructible>, expected_sc<NoCopyConstructible> const& >::value, "");
  static_assert(! std::is_copy_constructible<expected_sc<NoCopyConstructible>>::value, "");

  {
    // fixme
    expected_sc<NoMoveConstructible> x;
    expected_sc<NoMoveConstructible> y =  std::move(x); // DOESN'T FAIL as copy is used instead
    (void)y;
  }
  {
    // fixme
    expected_sc<NoMoveConstructible> x;
    expected_sc<NoMoveConstructible> y {expected_sc<NoMoveConstructible>{}}; // DOESN'T FAIL as copy is used instead
    (void)y;
  }
  static_assert(! std::is_move_constructible<NoMoveConstructible>::value, "");
  static_assert( std::is_constructible<expected_sc<NoMoveConstructible>, NoMoveConstructible && >::value, "");
  static_assert( std::is_move_constructible<expected_sc<NoMoveConstructible>>::value, "");

  return ::boost::report_errors();
}

max_stars_repo_name: jwakely/std-make

max_issues_repo_name: jwakely/std-make

max_forks_repo_name: jwakely/std-make


text:
#define CATCH_CONFIG_MAIN
#include "../../vendor/catch/catch.hpp"

#include "boost_defs.hpp"

#include "device_detail.hpp"
#include "thread_utility.hpp"
#include <boost/optional/optional_io.hpp>

TEST_CASE("initialize") {
  krbn::thread_utility::register_main_thread();
}

TEST_CASE("to_json") {
  using namespace std::string_literals;

  {
    krbn::device_detail device_detail(boost::none,
                                      boost::none,
                                      boost::none,
                                      boost::none,
                                      boost::none,
                                      boost::none,
                                      boost::none,
                                      boost::none,
                                      true,
                                      false);
    nlohmann::json json;
    json["is_keyboard"] = true;
    json["is_pointing_device"] = false;
    REQUIRE(device_detail.to_json() == json);
  }
  {
    krbn::device_detail device_detail(krbn::vendor_id(123),
                                      krbn::product_id(234),
                                      krbn::location_id(345),
                                      "m"s,
                                      "p"s,
                                      "s"s,
                                      "t"s,
                                      98765,
                                      false,
                                      true);
    nlohmann::json json;
    json["vendor_id"] = 123;
    json["product_id"] = 234;
    json["location_id"] = 345;
    json["manufacturer"] = "m";
    json["product"] = "p";
    json["serial_number"] = "s";
    json["transport"] = "t";
    json["registry_entry_id"] = 98765;
    json["is_keyboard"] = false;
    json["is_pointing_device"] = true;

    REQUIRE(device_detail.to_json() == json);
  }
}

TEST_CASE("from_json") {
  {
    nlohmann::json json;
    json["vendor_id"] = 123;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["product_id"] = 123;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["location_id"] = 123;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["manufacturer"] = "m";
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["product"] = "p";
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["serial_number"] = "s";
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["transport"] = "t";
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["registry_entry_id"] = 123;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["is_keyboard"] = true;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
  {
    nlohmann::json json;
    json["is_pointing_device"] = true;
    REQUIRE(krbn::device_detail(json).to_json() == json);
  }
}

TEST_CASE("compare") {
  using namespace std::string_literals;

  krbn::device_detail device_detail0(boost::none,
                                     boost::none,
                                     boost::none,
                                     boost::none,
                                     boost::none,
                                     boost::none,
                                     boost::none,
                                     boost::none,
                                     false,
                                     false);

  krbn::device_detail device_detail1(krbn::vendor_id(123),
                                     krbn::product_id(234),
                                     krbn::location_id(345),
                                     "m1"s,
                                     "p1"s,
                                     "s1"s,
                                     "t1"s,
                                     98765,
                                     true,
                                     true);

  krbn::device_detail device_detail2(krbn::vendor_id(123),
                                     krbn::product_id(234),
                                     krbn::location_id(345),
                                     "m2"s,
                                     "p1"s,
                                     "s1"s,
                                     "t1"s,
                                     98765,
                                     true,
                                     false);

  krbn::device_detail device_detail3(krbn::vendor_id(123),
                                     krbn::product_id(234),
                                     krbn::location_id(345),
                                     "m1"s,
                                     "p2"s,
                                     "s1"s,
                                     "t1"s,
                                     98765,
                                     false,
                                     true);

  krbn::device_detail device_detail4(krbn::vendor_id(123),
                                     krbn::product_id(234),
                                     krbn::location_id(345),
                                     "m2"s,
                                     "p2"s,
                                     "s1"s,
                                     "t1"s,
                                     98765,
                                     false,
                                     false);

  REQUIRE(device_detail0.compare(device_detail0) == false);
  REQUIRE(device_detail1.compare(device_detail1) == false);
  // 0,1 (["",""], [p1,m1])
  REQUIRE(device_detail1.compare(device_detail0) == false);
  REQUIRE(device_detail0.compare(device_detail1) == true);
  // 1,2 ([p1,m1], [p1,m2])
  REQUIRE(device_detail1.compare(device_detail2) == true);
  REQUIRE(device_detail2.compare(device_detail1) == false);
  // 1,3 ([p1,m1], [p2,m1])
  REQUIRE(device_detail1.compare(device_detail3) == true);
  REQUIRE(device_detail3.compare(device_detail1) == false);
  // 3,4 ([p2,m1], [p2,m2])
  REQUIRE(device_detail3.compare(device_detail4) == true);
  REQUIRE(device_detail4.compare(device_detail3) == false);
}


text:
/*

Copyright (c) 2005-2015, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * 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.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

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 HOLDER 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.

*/

#ifndef CELLVECDATA_HPP_
#define CELLVECDATA_HPP_

#include <boost/shared_ptr.hpp>
#include <map>
#include <string>
#include <vector>

#include <petscvec.h>

#include "AbstractCellProperty.hpp"
#include "ChasteSerialization.hpp"
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/array.hpp>
#include "Exception.hpp"
#include "ArchiveLocationInfo.hpp"

/**
 * CellVecData class.
 *
 * This cell property allows each cell to store one or more 'named' PETSc Vec variables associated with it,
 * for example solutions to reaction-diffusion PDEs solved within the cell. Other classes may interrogate
 * or modify the values stored in this class.
 */
class CellVecData : public AbstractCellProperty
{
private:

    /**
     * The cell data.
     */
    std::map<std::string, Vec> mCellVecData;

    /**
     * Whether to free in the destructor the Vec stored in the map
     */
    bool mFreeVecs;

    /** Needed for serialization. */
    friend class boost::serialization::access;
    /**
     * Archive the member variables.
     *
     * @param archive the archive
     * @param version the current version of this class
     */
    template<class Archive>
    void serialize(Archive & archive, const unsigned int version)
    {
        archive & boost::serialization::base_object<AbstractCellProperty>(*this);
    }

public:

    /**
     * Default constructor
     */
    CellVecData();

    /**
     * Copy constructor used to perform a deep copy of rAnotherCellVecData
     * @param rAnotherCellVecData  the class to copy
     */
    CellVecData(const CellVecData& rAnotherCellVecData);

    /**
     * Constructor required for serialisation. Avoid using otherwise.
     * @param rCellVecDataMap the internal map which needs to be copied
     */
    CellVecData(const std::map<std::string, Vec>& rCellVecDataMap);

    /**
     * We need the empty virtual destructor in this class to ensure Boost
     * serialization works correctly with static libraries.
     */
    virtual ~CellVecData();

    /**
     * This assigns the cell data.
     *
     * @param rVariableName the name of the data to be set.
     * @param data the value to set it to.
     */
    void SetItem(const std::string& rVariableName, Vec data);

    /**
     * @return data.
     *
     * @param rVariableName the index of the data required.
     * throws if rVariableName has not been stored
     *
     */
    Vec GetItem(const std::string& rVariableName) const;

    /**
     * @return number of data items
     */
    unsigned GetNumItems() const;

    /**
     * @return all keys.
     *
     * According to STL these are sorted in lexicographical/alphabetic order (so that the ordering here is predictable).
     */
    std::vector<std::string> GetKeys() const;
};

#include "SerializationExportWrapper.hpp"
// Declare identifier for the serializer
CHASTE_CLASS_EXPORT(CellVecData)

namespace boost
{
namespace serialization
{
template<class Archive>
inline void save_construct_data(
    Archive & ar, const CellVecData * t, const unsigned int file_version)
{
    unsigned map_size = t->GetNumItems();
    ar << map_size;

    std::vector<std::string> keys = t->GetKeys();

    for (std::vector<std::string>::iterator iter = keys.begin(); iter != keys.end(); ++iter)
    {
        std::string key = *iter;
        ar << key;
        Vec vec_data = t->GetItem(key);

        // Machinery for archiving a sequential PETSc Vec as a boost make_array.
        double *p_vec_data;
        VecGetArray(vec_data, &p_vec_data);
        PetscInt size, local_size;
        VecGetSize(vec_data, &size);
        VecGetLocalSize(vec_data, &local_size);
        // This will fail if we run in parallel and the vector is shared MPI (size is global)
        assert( local_size == size);
        ar << size;
        ar << make_array(p_vec_data, size);
        VecRestoreArray(vec_data, &p_vec_data);

        //std::string archive_filename = ArchiveLocationInfo::GetArchiveDirectory() + key + ".vec";
        //PetscTools::DumpPetscObject(t->GetItem(key), archive_filename);
    }

}

template<class Archive>
inline void load_construct_data(
    Archive & ar, CellVecData * t, const unsigned int file_version)
{
    unsigned map_size;
    ar >> map_size;

    std::map<std::string, Vec> archived_cell_vec_data;
    for (unsigned map_entry = 0; map_entry < map_size; ++map_entry)
    {
        std::string key;
        ar >> key;

        // Un-archive vector data and construct a new PETSc Vec for it
        PetscInt size;
        ar >> size;
        Vec archived_vec = PetscTools::CreateVec(size);
        double *p_archived_vec;
        VecGetArray(archived_vec, &p_archived_vec);
        ar >> make_array<double>(p_archived_vec, size);
        VecRestoreArray(archived_vec, &p_archived_vec);

        //std::string archive_filename = ArchiveLocationInfo::GetArchiveDirectory() + key + ".vec";
        //PetscTools::ReadPetscObject(archived_vec, archive_filename);

        archived_cell_vec_data[key] = archived_vec;
    }

     ::new(t)CellVecData(archived_cell_vec_data);
}
}
}

#endif /* CELLVECDATA_HPP_ */

max_stars_repo_name: ktunya/ChasteMod

max_issues_repo_name: ktunya/ChasteMod

max_forks_repo_name: ktunya/ChasteMod