text:
// __BEGIN_LICENSE__
// Copyright (C) 2006-2011 United States Government as represented by
// the Administrator of the National Aeronautics and Space Administration.
// All Rights Reserved.
// __END_LICENSE__

#include <vw/Core/ProgressCallback.h>
#include <vw/Core/Log.h>
#include <vw/FileIO/TemporaryFile.h>
#include <vw/Plate/detail/LocalIndex.h>
#include <vw/Plate/detail/RemoteIndex.h>
#include <vw/Plate/Blob.h>
#include <vw/Plate/BlobManager.h>
#include <vw/Plate/detail/Seed.h>

using namespace vw;
using namespace vw::platefile;
using namespace vw::platefile::detail;

#include <fstream>
#include <boost/foreach.hpp>
#include <boost/filesystem.hpp>
#include <boost/regex.hpp>
#include <algorithm>
namespace fs = boost::filesystem;

LocalIndexPage::LocalIndexPage(std::string filename, uint32 level, uint32 base_col,
                               uint32 base_row, uint32 page_width, uint32 page_height)
  : IndexPage(level, base_col, base_row, page_width, page_height),
    m_filename(filename), m_needs_saving(false)
{
  if (fs::exists(filename)) {
    this->deserialize();
  }
}

LocalIndexPage::~LocalIndexPage() {
  this->sync();
}

// Hijack this method momentartarily to mark the page as "dirty" by
// setting m_needs_saving to true.
void LocalIndexPage::set(TileHeader const& header, IndexRecord const& record) {

  // First call up to the parent class and let the original code run.
  IndexPage::set(header, record);

  // Then, mark this page is 'dirty' so that it gets saved to disk
  // when destroyed.
  m_needs_saving = true;
}

void LocalIndexPage::sync() {
  if (m_needs_saving) {
    this->serialize();
    m_needs_saving = false;
  }
}

void LocalIndexPage::serialize() {

  fs::path pagename(m_filename);
  fs::path parent(pagename.parent_path());

  try {
    fs::create_directories(parent);
  } catch ( const std::exception& e ) {
    vw_throw(IOErr() << "Could not create IndexPage directory: " << e.what());
  }

  std::string tmpname;
  {
    TemporaryFile tmp(parent.file_string(), false);
    IndexPage::serialize(tmp);
    tmpname = tmp.filename();
  }

  if (::rename(tmpname.c_str(), m_filename.c_str()) == -1)
    vw_throw(IOErr() << "LocalIndexPage::serialize(): failed to rename temp page to " << m_filename << ": " << ::strerror(errno));

  m_needs_saving = false;
}

void LocalIndexPage::deserialize() {

  std::ifstream istr(m_filename.c_str());
  if (!istr.good())
    vw_throw(IOErr() << "IndexPage::deserialize() failed.  Could not open " << m_filename << " for reading.");

  // Call up to superclass to finish deserializing.
  try {
    IndexPage::deserialize(istr);
  } catch (const vw::IOErr &e) {
    // Add more useful error reporting.
    vw_throw(IOErr() << "Failed to load page \"" << m_filename << "\": " << e.what());
  }

  m_needs_saving = false;
}

 // ----------------------------------------------------------------------
 //                    LOCAL INDEX PAGE GENERATOR
 // ----------------------------------------------------------------------

 LocalPageGenerator::LocalPageGenerator( std::string filename, uint32 level, uint32 base_col,
                                         uint32 base_row, uint32 page_width, uint32 page_height)
 : m_filename( filename ), m_level(level), m_base_col(base_col), m_base_row(base_row),
   m_page_width(page_width), m_page_height(page_height) {}

boost::shared_ptr<IndexPage>
LocalPageGenerator::generate() const {
 return boost::shared_ptr<IndexPage>(new LocalIndexPage(m_filename, m_level,
       m_base_col, m_base_row,
       m_page_width, m_page_height) );
}

boost::shared_ptr<PageGeneratorBase>
LocalPageGeneratorFactory::create(uint32 level, uint32 base_col, uint32 base_row, uint32 page_width, uint32 page_height)
{
  // Generate a filename
  std::ostringstream filename;
  filename << m_plate_filename
    << "/index"
    << "/" << level
    << "/" << base_row
    << "/" << base_col;

  // Create the proper type of page generator.
  boost::shared_ptr<PageGeneratorBase> page_gen(
      new LocalPageGenerator(filename.str(), level, base_col, base_row,
        page_width, page_height) );

  return page_gen;
}

std::string LocalPageGeneratorFactory::who() const {
  return fs::basename(m_plate_filename);
}

// -------------------------------------------------------------------
//                            LOCAL INDEX
// -------------------------------------------------------------------

// ------------------------ Private Methods --------------------------

std::string LocalIndex::index_filename() const {
  return m_plate_filename + "/plate.index";
}

std::string LocalIndex::log_filename() const {
  return m_plate_filename + "/plate.log";
}

std::vector<std::string> LocalIndex::blob_filenames() const {

  std::vector<std::string> result;

  if ( !fs::exists( m_plate_filename ) )
    vw_throw(IOErr() << "Could not list blob filenames for plate (" << m_plate_filename << "): no such directory");

  boost::regex re("plate_(\\d+)\\.blob");
  typedef fs::directory_iterator iter_t;

  BOOST_FOREACH(const fs::path& elt, std::make_pair(iter_t(m_plate_filename), iter_t())) {
    boost::cmatch matches;
    if (!boost::regex_match(elt.filename().c_str(), matches, re)) {
      vw_out(DebugMessage, "plate") << "Skipping non-blob-filename " << elt.filename() << std::endl;
      continue;
    }

    std::string blob_id_str(matches[1].first, matches[1].second);
    uint32 blob_id = boost::lexical_cast<uint32>(blob_id_str);
    if (result.size() < blob_id+1)
      result.resize(blob_id+1);
    result[blob_id] = elt.string();
  }

  return result;
}

void LocalIndex::save_index_file() const {

  std::ofstream ofstr( this->index_filename().c_str() );
  if (!ofstr.is_open())
    vw_throw(IOErr() << "LocalIndex::save_index_file() -- Could not create index file for writing.");

  // Serialize the header data.
  m_header.SerializeToOstream(&ofstr);
  ofstr.close();

}

 // ------------------------ Public Methods --------------------------

 /// Create a new index.
 LocalIndex::LocalIndex( std::string plate_filename, IndexHeader new_index_info)
  : PagedIndex(boost::shared_ptr<PageGeneratorFactory>( new LocalPageGeneratorFactory(plate_filename) )),
   m_plate_filename(plate_filename), m_blob_manager(boost::shared_ptr<BlobManager>( new BlobManager(m_plate_filename) ))
{
   // Create subdirectory for storing hard copies of index pages.
   std::string base_index_path = plate_filename + "/index";

   if (fs::exists(base_index_path)) {
      open_impl();
#define WARN_IF_DIFFERENT(field) do {if (new_index_info.field() != m_header.field()) vw_out(WarningMessage, "plate") << "Refusing to change default for field " << #field << " from " << m_header.field() << " to " << new_index_info.field() << std::endl;} while(0)
     // Ignore attempts to change the header, but warn if anything is different.
      WARN_IF_DIFFERENT(tile_size);
      WARN_IF_DIFFERENT(tile_filetype);
      WARN_IF_DIFFERENT(pixel_format);
      WARN_IF_DIFFERENT(channel_type);
      WARN_IF_DIFFERENT(type);
      WARN_IF_DIFFERENT(description);
#undef WARN_IF_DIFFERENT
     return;
   }

   fs::create_directories(base_index_path);

   // Start with the new_index_info, which provides the tile size, file
   // type, pixel and channel types, etc.  Then we augment it with a
   // few things to help track the new platefile.
   m_header = new_index_info;

   // Check to make sure we've selected a sane file type
   if (m_header.channel_type() == VW_CHANNEL_FLOAT32) {
     if (m_header.tile_filetype() == "png" || m_header.tile_filetype() == "jpg")
       vw_out(WarningMessage, "plate")
         << "Constructing 32-bit floating point platefile with a non-32-bit file-type ("
         << m_header.tile_filetype() << ").\n";
   }

#define WARN_IF_SET(field) do {if (new_index_info.has_ ## field()) vw_out(ErrorMessage, "plate") << #field << " is a private field. Ignoring your value " << m_header.field() << std::endl;} while(0)
   WARN_IF_SET(platefile_id);
   WARN_IF_SET(version);
   WARN_IF_SET(transaction_read_cursor);
   WARN_IF_SET(transaction_write_cursor);
   WARN_IF_SET(num_levels);
#undef WARN_IF_SET

   // Create a unique (random) platefile_id using the random()
   // function.  It would probably be better to use some sort of fancy
   // MD5 hash of the current time and date, or something like that,
   // but little 'ol random() will probably work just fine for our
   // humble purposes.
   plate_seed_random();
   m_header.set_platefile_id(vw::uint32(random()));

   // Set up the IndexHeader and write it to disk.
   m_header.set_version(VW_PLATE_INDEX_VERSION);
   m_header.set_transaction_read_cursor(0);  // Transaction 0 is the empty mosaic
   m_header.set_transaction_write_cursor(1); // Transaction 1 is the first valid transaction to write into
   m_header.set_num_levels(0);               // Index initially contains zero levels
   this->save_index_file();

   // Create the logging facility
   m_log = boost::shared_ptr<LogInstance>( new LogInstance(this->log_filename()) );
   this->log() << "Created new index: \"" << this->index_filename() << "\n"
               <<  m_header.DebugString() << "\n";
 }

  void LocalIndex::open_impl() {
    // First, check to make sure the platefile directory exists.
    if ( !fs::exists( m_plate_filename ) )
      vw_throw(ArgumentErr() << "Could not open local platefile (" << m_plate_filename << "): no such directory");

    // The load the index file & parse the IndexHeader protobuffer
    std::ifstream ifstr( this->index_filename().c_str() );
    if (!ifstr.is_open())
      vw_throw(IOErr() << "LocalIndex::LocalIndex() -- Could not load index file for reading.");

    bool worked = m_header.ParseFromIstream(&ifstr);
    if (!worked)
      vw_throw(IOErr() << "LocalIndex::LocalIndex() -- Could not parse index information from "
          << this->index_filename());

    ifstr.close();

    // Load Index Levels for PagedIndex
    for (uint32 level = 0; level < this->num_levels(); ++level) {
      boost::shared_ptr<IndexLevel> new_level( new IndexLevel(m_page_gen_factory, level,
            m_page_width, m_page_height,
            m_default_cache_size) );
      m_levels.push_back(new_level);
    }

    m_log = boost::shared_ptr<LogInstance>( new LogInstance(this->log_filename()) );
    this->log() << "Reopened index \"" << this->index_filename() << "\n" << m_header.DebugString() << "\n";
  }

/// Open an existing index from a file on disk.
LocalIndex::LocalIndex(std::string plate_filename) :
  PagedIndex(boost::shared_ptr<PageGeneratorFactory>( new LocalPageGeneratorFactory(plate_filename) ) ),  // superclass constructor
  m_plate_filename(plate_filename), m_blob_manager(boost::shared_ptr<BlobManager>( new BlobManager(m_plate_filename) ))
{
  open_impl();
}

 /// Use this to send data to the index's logfile like this:
 ///
 ///   index_instance.log() << "some text for the log...\n";
 ///
 std::ostream& LocalIndex::log() {
   return (*m_log)(InfoMessage, "console");
 }

// Clients are expected to make a transaction request whenever
// they start a self-contained chunk of mosaicking work.  .
Transaction LocalIndex::transaction_request(std::string transaction_description,
                                            TransactionOrNeg transaction_id_override_neg) {

   if (transaction_id_override_neg != -1) {
     Transaction transaction_id_override = transaction_id_override_neg.promote();
     // If the user has chosen to override the transaction ID's, then we
     // use the transaction ID they specify.  We also increment the
     // transaction_write_cursor if necessary so that we dole out a
     // reasonable transaction ID if we are ever asked again without an
     // override.
     Transaction max_trans_id = std::max(m_header.transaction_write_cursor(), transaction_id_override+1);
     m_header.set_transaction_write_cursor(max_trans_id);
     this->save_index_file();
     this->log() << "Transaction " << transaction_id_override << " started: " << transaction_description << "\n";
     return transaction_id_override;

   } else {

     // Pick the next transaction ID, increment the cursor, and then save
     // the cursor to a file.  (Saving to a file gurantees that we won't
     // accidentally assign two transactions the same ID if the index
     // server crashes and has to be restarted.
     Transaction transaction_id = m_header.transaction_write_cursor();
     m_header.set_transaction_write_cursor(transaction_id + 1);
     this->save_index_file();
     this->log() << "Transaction " << transaction_id << " started: " << transaction_description << "\n";
     return transaction_id;
   }
 }

 // Once a chunk of work is complete, clients can "commit" their
 // work to the mosaic by issuing a transaction_complete method.
 void LocalIndex::transaction_complete(Transaction transaction_id, bool update_read_cursor) {

   // First we save (sync) the index pages to disk
   //this->sync();

   if ( update_read_cursor ) {
     Transaction max_trans_id = std::max(m_header.transaction_read_cursor(), transaction_id+0);
     m_header.set_transaction_read_cursor(max_trans_id);
     this->save_index_file();
   }

   this->log() << "Transaction " << transaction_id << " complete.  "
               << "[ read_cursor = " << m_header.transaction_read_cursor() << " ]\n";
 }

 // Once a chunk of work is complete, clients can "commit" their
 // work to the mosaic by issuing a transaction_complete method.
 void LocalIndex::transaction_failed(Transaction transaction_id) {

   // Log the failure.
   this->log() << "Transaction " << transaction_id << " FAILED.\n";

 }

 // Return the current location of the transaction cursor.  This
 // will be the last transaction id that refers to a coherent
 // version of the mosaic.
 Transaction LocalIndex::transaction_cursor() {
   return m_header.transaction_read_cursor();
 }

 // Load index entries by iterating through TileHeaders saved in the
 // blob file.  This function essentially rebuilds an index in memory
 // using entries that had been previously saved to disk.
 void LocalIndex::rebuild_index() {

   vw_out(InfoMessage) << "Rebuilding index: " << m_plate_filename <<"\n";

   // index in blob_filenames is the blobfile id
   std::vector<std::string> blob_files = this->blob_filenames();
   for (uint32 blob_id = 0; blob_id < blob_files.size(); ++blob_id) {
     const std::string& name = blob_files[blob_id];
     if (name.empty())
       continue;

     TerminalProgressCallback tpc("plate", "\t --> Rebuild from blob" + vw::stringify(blob_id) + " : ");
     tpc.report_progress(0);

     ReadBlob blob(name);

     IndexRecord rec;
     typedef Blob::iterator iter_t;
     for (iter_t tile = blob.begin(), end = blob.end(); tile != end; ++tile) {
       const BlobTileRecord& hdr = *tile;
       rec.set_blob_id(blob_id);
       rec.set_blob_offset(tile.current_base_offset());
       rec.set_filetype(hdr.hdr.filetype());
       this->write_update(hdr.hdr, rec);
       tpc.report_progress(float(tile.current_base_offset()) / float(blob.size()));
    }
    tpc.report_finished();
  }
}

// -----------------------    I/O      ----------------------

/// Writing, pt. 1: Reserve a blob lock
uint32 LocalIndex::write_request() {
  return m_blob_manager->request_lock();
}

/// Writing, pt. 1: Reserve a blob lock
void LocalIndex::write_update(TileHeader const& header, IndexRecord const& record) {

  // Store the number of tiles that are contained in the mosaic.
  size_t starting_size = m_levels.size();

  // Write the update to the PagedIndex superclass.
  PagedIndex::write_update(header, record);

  // If adding the record resulted in more levels, we save that
  // information to the index header.
  if (m_levels.size() != starting_size) {
    m_header.set_num_levels(boost::numeric_cast<uint32>(m_levels.size()));
    this->save_index_file();
  }
}

/// Writing, pt. 3: Signal the completion
void LocalIndex::write_complete(uint32 blob_id) {
  m_blob_manager->release_lock(blob_id);
}

ext: cc

max_stars_repo_name: digimatronics/ComputerVision

max_issues_repo_name: rkrishnasanka/visionworkbench

max_forks_repo_name: rkrishnasanka/visionworkbench


text:
//
// FILE NAME:       $HeadURL: svn+ssh://svn.cm.aol.com/advertising/adlearn/gen1/trunk/lib/cpp/DataProxy/mock/TestableDataProxyClient.cpp $
//
// REVISION:        $Revision: 220478 $
//
// COPYRIGHT:       (c) 2006 Advertising.com All Rights Reserved.
//
// LAST UPDATED:    $Date: 2011-08-23 14:38:03 -0400 (Tue, 23 Aug 2011) $
// UPDATED BY:      $Author: bhh1988 $

#include "TestableDataProxyClient.hpp"
#include "DatabaseConnectionManager.hpp"
#include "NodeFactory.hpp"
#include <boost/scoped_ptr.hpp>

namespace
{
	boost::scoped_ptr< DatabaseConnectionManager > s_pDatabaseConnectionManager( NULL );
}

TestableDataProxyClient::TestableDataProxyClient()
:	DataProxyClient( true )
{
}

TestableDataProxyClient::~TestableDataProxyClient()
{
}

void TestableDataProxyClient::Initialize( const std::string& i_rConfigFileSpec )
{
	if( s_pDatabaseConnectionManager == NULL )
	{
		s_pDatabaseConnectionManager.reset( new DatabaseConnectionManager( *this ) );
	}

	NodeFactory nodeFactory( *this );
	DataProxyClient::Initialize( i_rConfigFileSpec, nodeFactory, *s_pDatabaseConnectionManager );
}

void TestableDataProxyClient::Initialize( const std::string& i_rConfigFileSpec,
										  INodeFactory& i_rNodeFactory )
{
	if( s_pDatabaseConnectionManager == NULL )
	{
		s_pDatabaseConnectionManager.reset( new DatabaseConnectionManager( *this ) );
	}

	DataProxyClient::Initialize( i_rConfigFileSpec, i_rNodeFactory, *s_pDatabaseConnectionManager );
}

void TestableDataProxyClient::Initialize( const std::string& i_rConfigFileSpec,
										  INodeFactory& i_rNodeFactory,
										  DatabaseConnectionManager& i_rDatabaseConnectionManager )
{
	DataProxyClient::Initialize( i_rConfigFileSpec, i_rNodeFactory, i_rDatabaseConnectionManager );
}

ext: cpp

max_stars_repo_name: aol/DataProxyLibrary

max_issues_repo_name: aol/DataProxyLibrary

max_forks_repo_name: aol/DataProxyLibrary


text:
#ifndef _SRC_UTIL_HPP
#define _SRC_UTIL_HPP
#include <sstream>
#include <iomanip>
#include "FileSystem.hpp"
#include "Color.hpp"

#include <string>
#include <vector>
#include <boost/uuid/sha1.hpp>

namespace Sit {
namespace Util {

/**
 * The struct to represent exception
 */
struct SitException {
	/**
	 * The message intended to be shown to user
	 */
	const std::string message;

	/**
	 * The details for debug
	 */
	const std::string detail;
	SitException(const std::string& Message, const std::string& Detail = "")
		: message(Message), detail(Detail) {}
};

/**
 * Return the SHA1 value of `content`
 */
std::string SHA1sum(const std::string &content);

/**
 * Split the `content` into lines
 */
std::vector<std::string> SplitLine();

/**
 * Return the first line of the string
 */
std::string FirstLine(const std::string& str);

/**
 * Convert file mode string to int
 */
int FileModeToInt(const std::string& mode);

/**
 * Convert file mode to string
 */
std::string FileModeToString(const int mode);

/**
 * Generate author string for committing
 */
std::string AuthorString(const std::string& name, const std::string& email, const std::string& datetime);

/*
 *Complete SHA1 value with exists objects.
 */
std::string SHA1Complete(std::string _id);
}
}
#endif

ext: hpp


text:
//=======================================================================
// Copyright 1997, 1998, 1999, 2000 University of Notre Dame.
// Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek
//
// This file is part of the Boost Graph Library
//
// You should have received a copy of the License Agreement for the
// Boost Graph Library along with the software; see the file LICENSE.
// If not, contact Office of Research, University of Notre Dame, Notre
// Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//=======================================================================
#include <boost/config.hpp>
#include <iostream>
#include <vector>
#include <algorithm>
#include <utility>
#include <boost/graph/adjacency_list.hpp>
#include <boost/pending/disjoint_sets.hpp>
#include <boost/graph/incremental_components.hpp>

/*

  This example shows how to use the disjoint set data structure
  to compute the connected components of an undirected, changing
  graph.

  Sample output:

  An undirected graph:
  0 <--> 1 4
  1 <--> 0 4
  2 <--> 5
  3 <-->
  4 <--> 1 0
  5 <--> 2

  representative[0] = 1
  representative[1] = 1
  representative[2] = 5
  representative[3] = 3
  representative[4] = 1
  representative[5] = 5

  component 0 contains: 4 1 0
  component 1 contains: 3
  component 2 contains: 5 2

 */

using namespace std;

int main(int , char* [])
{
  using namespace boost;
  typedef adjacency_list <vecS, vecS, undirectedS> Graph;
  typedef graph_traits<Graph>::vertex_descriptor Vertex;
  typedef graph_traits<Graph>::vertices_size_type size_type;

  const int N = 6;
  Graph G(N);

  std::vector<Vertex> rank(num_vertices(G));
  std::vector<Vertex> parent(num_vertices(G));
  // VC++ version of std::vector has no ::pointer, so
  // I use ::value_type* instead.
  typedef std::vector<size_type>::value_type* Rank;
  typedef std::vector<Vertex>::value_type* Parent;
  disjoint_sets<Rank, Parent>  ds(&rank[0], &parent[0]);

  initialize_incremental_components(G, ds);
  incremental_components(G, ds);

  graph_traits<Graph>::edge_descriptor e;
  bool flag;
  boost::tie(e,flag) = add_edge(0, 1, G);
  ds.union_set(0,1);

  boost::tie(e,flag) = add_edge(1, 4, G);
  ds.union_set(1,4);

  boost::tie(e,flag) = add_edge(4, 0, G);
  ds.union_set(4,0);

  boost::tie(e,flag) = add_edge(2, 5, G);
  ds.union_set(2,5);

  cout << "An undirected graph:" << endl;
  print_graph(G, get(vertex_index, G));
  cout << endl;

  graph_traits<Graph>::vertex_iterator i,end;
  for (boost::tie(i, end) = vertices(G); i != end; ++i)
    cout << "representative[" << *i << "] = " <<
      ds.find_set(*i) << endl;;
  cout << endl;

  typedef component_index<unsigned int> Components;
  Components components(&parent[0], &parent[0] + parent.size());

  for (Components::size_type c = 0; c < components.size(); ++c) {
    cout << "component " << c << " contains: ";
    Components::value_type::iterator
      j = components[c].begin(),
      jend = components[c].end();
    for ( ; j != jend; ++j)
      cout << *j << " ";
    cout << endl;
  }

  return 0;
}


text:
#include "Common/Common.h"
#include "Demos/Visualization/MiniGL.h"
#include "Demos/Visualization/Selection.h"
#include "GL/glut.h"
#include "Demos/Simulation/TimeManager.h"
#include <Eigen/Dense>
#include "GenericConstraintsModel.h"
#include <iostream>
#include "Demos/Simulation/TimeStepController.h"
#include "Demos/Visualization/Visualization.h"
#include "Demos/Utils/Logger.h"
#include "Demos/Utils/Timing.h"
#include "Demos/Utils/FileSystem.h"

// Enable memory leak detection
#if defined(_DEBUG) && !defined(EIGEN_ALIGN)
	#define new DEBUG_NEW
#endif

INIT_TIMING
INIT_LOGGING

using namespace PBD;
using namespace Eigen;
using namespace std;

void timeStep ();
void buildModel ();
void createMesh();
void render ();
void cleanup();
void reset();
void initShader();
void selection(const Eigen::Vector2i &start, const Eigen::Vector2i &end);
void TW_CALL setTimeStep(const void *value, void *clientData);
void TW_CALL getTimeStep(void *value, void *clientData);
void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData);
void TW_CALL getVelocityUpdateMethod(void *value, void *clientData);
void TW_CALL setStiffness(const void *value, void *clientData);
void TW_CALL getStiffness(void *value, void *clientData);
void TW_CALL setBendingStiffness(const void *value, void *clientData);
void TW_CALL getBendingStiffness(void *value, void *clientData);

GenericConstraintsModel model;
TimeStepController simulation;

const int nRows = 30;
const int nCols = 30;
const Real width = 10.0;
const Real height = 10.0;
bool doPause = true;
std::vector<unsigned int> selectedParticles;
Vector3r oldMousePos;
Shader *shader;
string exePath;
string dataPath;

// main
int main( int argc, char **argv )
{
	REPORT_MEMORY_LEAKS

	std::string logPath = FileSystem::normalizePath(FileSystem::getProgramPath() + "/log");
	FileSystem::makeDirs(logPath);
	logger.addSink(unique_ptr<ConsoleSink>(new ConsoleSink(LogLevel::INFO)));
	logger.addSink(unique_ptr<FileSink>(new FileSink(LogLevel::DEBUG, logPath + "/PBD.log")));

	exePath = FileSystem::getProgramPath();
	dataPath = exePath + "/" + std::string(PBD_DATA_PATH);

	// OpenGL
	MiniGL::init (argc, argv, 1024, 768, 0, 0, "Cloth demo");
	MiniGL::initLights ();
	MiniGL::initTexture();
	MiniGL::setClientIdleFunc (50, timeStep);
	MiniGL::setKeyFunc(0, 'r', reset);
	MiniGL::setSelectionFunc(selection);
	initShader();

	buildModel ();

	MiniGL::setClientSceneFunc(render);
	MiniGL::setViewport (40.0f, 0.1f, 500.0f, Vector3r (5.0, 10.0, 30.0), Vector3r (5.0, 0.0, 0.0));

	TwAddVarRW(MiniGL::getTweakBar(), "Pause", TW_TYPE_BOOLCPP, &doPause, " label='Pause' group=Simulation key=SPACE ");
	TwAddVarCB(MiniGL::getTweakBar(), "TimeStepSize", TW_TYPE_REAL, setTimeStep, getTimeStep, &model, " label='Time step size'  min=0.0 max = 0.1 step=0.001 precision=4 group=Simulation ");
	TwType enumType = TwDefineEnum("VelocityUpdateMethodType", NULL, 0);
	TwAddVarCB(MiniGL::getTweakBar(), "VelocityUpdateMethod", enumType, setVelocityUpdateMethod, getVelocityUpdateMethod, &simulation, " label='Velocity update method' enum='0 {First Order Update}, 1 {Second Order Update}' group=Simulation");
	TwAddVarCB(MiniGL::getTweakBar(), "Stiffness", TW_TYPE_REAL, setStiffness, getStiffness, &model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Distance constraints' ");
	TwAddVarCB(MiniGL::getTweakBar(), "BendingStiffness", TW_TYPE_REAL, setBendingStiffness, getBendingStiffness, &model, " label='Bending stiffness'  min=0.0 step=0.01 precision=4 group=Bending ");

	glutMainLoop ();

	cleanup ();

	Timing::printAverageTimes();

	return 0;
}

void initShader()
{
	std::string vertFile = dataPath + "/shaders/vs_smoothTex.glsl";
	std::string fragFile = dataPath + "/shaders/fs_smoothTex.glsl";
	shader = MiniGL::createShader(vertFile, "", fragFile);

	if (shader == NULL)
		return;

	shader->begin();
	shader->addUniform("modelview_matrix");
	shader->addUniform("projection_matrix");
	shader->addUniform("surface_color");
	shader->addUniform("shininess");
	shader->addUniform("specular_factor");
	shader->end();
}

void cleanup()
{
	delete TimeManager::getCurrent();
	delete shader;
}

void reset()
{
	Timing::printAverageTimes();
	Timing::reset();

	model.reset();
	simulation.reset();
	TimeManager::getCurrent()->setTime(0.0);
}

void mouseMove(int x, int y)
{
	Vector3r mousePos;
	MiniGL::unproject(x, y, mousePos);
	const Vector3r diff = mousePos - oldMousePos;

	TimeManager *tm = TimeManager::getCurrent();
	const Real h = tm->getTimeStepSize();

	ParticleData &pd = model.getParticles();
	for (unsigned int j = 0; j < selectedParticles.size(); j++)
	{
		pd.getVelocity(selectedParticles[j]) += 5.0*diff/h;
	}
	oldMousePos = mousePos;
}

void selection(const Eigen::Vector2i &start, const Eigen::Vector2i &end)
{
	std::vector<unsigned int> hits;
	selectedParticles.clear();
	ParticleData &pd = model.getParticles();
	Selection::selectRect(start, end, &pd.getPosition(0), &pd.getPosition(pd.size() - 1), selectedParticles);
	if (selectedParticles.size() > 0)
		MiniGL::setMouseMoveFunc(GLUT_MIDDLE_BUTTON, mouseMove);
	else
		MiniGL::setMouseMoveFunc(-1, NULL);

	MiniGL::unproject(end[0], end[1], oldMousePos);
}

void timeStep ()
{
	if (doPause)
		return;

	// Simulation code
	for (unsigned int i = 0; i < 8; i++)
		simulation.step(model);

	for (unsigned int i = 0; i < model.getTriangleModels().size(); i++)
		model.getTriangleModels()[i]->updateMeshNormals(model.getParticles());
}

void buildModel ()
{
	TimeManager::getCurrent ()->setTimeStepSize (0.005);

	createMesh();
}

void renderTriangleModels()
{
	// Draw simulation model

	const ParticleData &pd = model.getParticles();
	float surfaceColor[4] = { 0.2f, 0.5f, 1.0f, 1 };

	if (shader)
	{
		shader->begin();
		glUniform3fv(shader->getUniform("surface_color"), 1, surfaceColor);
		glUniform1f(shader->getUniform("shininess"), 5.0f);
		glUniform1f(shader->getUniform("specular_factor"), 0.2f);

		GLfloat matrix[16];
		glGetFloatv(GL_MODELVIEW_MATRIX, matrix);
		glUniformMatrix4fv(shader->getUniform("modelview_matrix"), 1, GL_FALSE, matrix);
		GLfloat pmatrix[16];
		glGetFloatv(GL_PROJECTION_MATRIX, pmatrix);
		glUniformMatrix4fv(shader->getUniform("projection_matrix"), 1, GL_FALSE, pmatrix);
	}

	for (unsigned int i = 0; i < model.getTriangleModels().size(); i++)
	{
		// mesh
		TriangleModel *triModel = model.getTriangleModels()[i];
		const IndexedFaceMesh &mesh = triModel->getParticleMesh();
		Visualization::drawTexturedMesh(pd, mesh, triModel->getIndexOffset(), surfaceColor);
	}
	if (shader)
		shader->end();
}

void render ()
{
	MiniGL::coordinateSystem();

	renderTriangleModels();

	float red[4] = { 0.8f, 0.0f, 0.0f, 1 };
	const ParticleData &pd = model.getParticles();
	for (unsigned int j = 0; j < selectedParticles.size(); j++)
	{
		MiniGL::drawSphere(pd.getPosition(selectedParticles[j]), 0.08f, red);
	}

	MiniGL::drawTime( TimeManager::getCurrent ()->getTime ());
}

/** Create a particle model mesh
*/
void createMesh()
{
	TriangleModel::ParticleMesh::UVs uvs;
	uvs.resize(nRows*nCols);

	const Real dy = width / (Real)(nCols - 1);
	const Real dx = height / (Real)(nRows - 1);

	Vector3r points[nRows*nCols];
	for (int i = 0; i < nRows; i++)
	{
		for (int j = 0; j < nCols; j++)
		{
			const Real y = (Real)dy*j;
			const Real x = (Real)dx*i;
			points[i*nCols + j] = Vector3r(x, 1.0, y);

			uvs[i*nCols + j][0] = x/width;
			uvs[i*nCols + j][1] = y/height;
		}
	}
	const int nIndices = 6 * (nRows - 1)*(nCols - 1);

	TriangleModel::ParticleMesh::UVIndices uvIndices;
	uvIndices.resize(nIndices);

	unsigned int indices[nIndices];
	int index = 0;
	for (int i = 0; i < nRows - 1; i++)
	{
		for (int j = 0; j < nCols - 1; j++)
		{
			int helper = 0;
			if (i % 2 == j % 2)
				helper = 1;

			indices[index] = i*nCols + j;
			indices[index + 1] = i*nCols + j + 1;
			indices[index + 2] = (i + 1)*nCols + j + helper;

			uvIndices[index] = i*nCols + j;
			uvIndices[index + 1] = i*nCols + j + 1;
			uvIndices[index + 2] = (i + 1)*nCols + j + helper;
			index += 3;

			indices[index] = (i + 1)*nCols + j + 1;
			indices[index + 1] = (i + 1)*nCols + j;
			indices[index + 2] = i*nCols + j + 1 - helper;

			uvIndices[index] = (i + 1)*nCols + j + 1;
			uvIndices[index + 1] = (i + 1)*nCols + j;
			uvIndices[index + 2] = i*nCols + j + 1 - helper;
			index += 3;
		}
	}
	model.addTriangleModel(nRows*nCols, nIndices / 3, &points[0], &indices[0], uvIndices, uvs);

	ParticleData &pd = model.getParticles();
	for (unsigned int i = 0; i < pd.getNumberOfParticles(); i++)
	{
		pd.setMass(i, 1.0);
	}

	// Set mass of points to zero => make it static
	pd.setMass(0, 0.0);
	pd.setMass((nRows-1)*nCols, 0.0);

	// init constraints
	for (unsigned int cm = 0; cm < model.getTriangleModels().size(); cm++)
	{
		const unsigned int offset = model.getTriangleModels()[cm]->getIndexOffset();
		IndexedFaceMesh &mesh = model.getTriangleModels()[cm]->getParticleMesh();
		const unsigned int nEdges = mesh.numEdges();
		const IndexedFaceMesh::Edge *edges = mesh.getEdges().data();

		// distance constraints
		for (unsigned int i = 0; i < nEdges; i++)
		{
			const unsigned int v1 = edges[i].m_vert[0] + offset;
			const unsigned int v2 = edges[i].m_vert[1] + offset;

			model.addGenericDistanceConstraint(v1, v2);
		}

		// bending constraints
		const unsigned int *tris = mesh.getFaces().data();
		for (unsigned int i = 0; i < nEdges; i++)
		{
			const int tri1 = edges[i].m_face[0];
			const int tri2 = edges[i].m_face[1];
			if ((tri1 != 0xffffffff) && (tri2 != 0xffffffff))
			{
				// Find the triangle points which do not lie on the axis
				const int axisPoint1 = edges[i].m_vert[0];
				const int axisPoint2 = edges[i].m_vert[1];
				int point1 = -1;
				int point2 = -1;
				for (int j = 0; j < 3; j++)
				{
					if ((tris[3 * tri1 + j] != axisPoint1) && (tris[3 * tri1 + j] != axisPoint2))
					{
						point1 = tris[3 * tri1 + j];
						break;
					}
				}
				for (int j = 0; j < 3; j++)
				{
					if ((tris[3 * tri2 + j] != axisPoint1) && (tris[3 * tri2 + j] != axisPoint2))
					{
						point2 = tris[3 * tri2 + j];
						break;
					}
				}
				if ((point1 != -1) && (point2 != -1))
				{
					const unsigned int vertex1 = point1 + offset;
					const unsigned int vertex2 = point2 + offset;
					const unsigned int vertex3 = edges[i].m_vert[0] + offset;
					const unsigned int vertex4 = edges[i].m_vert[1] + offset;
					model.addGenericIsometricBendingConstraint(vertex1, vertex2, vertex3, vertex4);
				}
			}
		}
	}

	LOG_INFO << "Number of triangles: " << nIndices / 3;
	LOG_INFO << "Number of vertices: " << nRows*nCols;

}

void TW_CALL setTimeStep(const void *value, void *clientData)
{
	const Real val = *(const Real *)(value);
	TimeManager::getCurrent()->setTimeStepSize(val);
}

void TW_CALL getTimeStep(void *value, void *clientData)
{
	*(Real *)(value) = TimeManager::getCurrent()->getTimeStepSize();
}

void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData)
{
	const short val = *(const short *)(value);
	((TimeStepController*)clientData)->setVelocityUpdateMethod((unsigned int)val);
}

void TW_CALL getVelocityUpdateMethod(void *value, void *clientData)
{
	*(short *)(value) = (short)((TimeStepController*)clientData)->getVelocityUpdateMethod();
}

void TW_CALL setStiffness(const void *value, void *clientData)
{
	const Real val = *(const Real *)(value);
	((GenericConstraintsModel*)clientData)->setClothStiffness(val);
}

void TW_CALL getStiffness(void *value, void *clientData)
{
	*(Real *)(value) = ((GenericConstraintsModel*)clientData)->getClothStiffness();
}

void TW_CALL setBendingStiffness(const void *value, void *clientData)
{
	const Real val = *(const Real *)(value);
	((GenericConstraintsModel*)clientData)->setClothBendingStiffness(val);
}

void TW_CALL getBendingStiffness(void *value, void *clientData)
{
	*(Real *)(value) = ((GenericConstraintsModel*)clientData)->getClothBendingStiffness();
}


text:
//  (C) Copyright Edward Diener 2011-2015
//  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 <boost/vmd/elem.hpp>
#include <boost/vmd/is_empty.hpp>
#include <boost/detail/lightweight_test.hpp>
#include <boost/preprocessor/list/at.hpp>
#include <boost/preprocessor/seq/elem.hpp>
#include <boost/preprocessor/tuple/elem.hpp>

int main()
  {

#if BOOST_PP_VARIADICS

  #define BOOST_VMD_REGISTER_ggh (ggh)

  #define ANIDENTIFIER ggh
  #define ANUMBER 249
  #define ASEQ (25)(26)(27)
  #define ATUPLE (0,1,2,3,((a,b))((c,d))((e))((f,g,h)))
  #define ATUPLE2 (10,11,12,13)
  #define ATUPLE3 (100,101,102,103)
  #define ATUPLE5 (200,201,202,203)
  #define ALIST (0,(1,(2,(3,BOOST_PP_NIL))))
  #define ANARRAY (3,(a,b,38))
  #define ANARRAY2 (5,(c,d,133,22,15))
  #define ASEQUENCE ANUMBER ALIST ATUPLE ANIDENTIFIER ANARRAY ASEQ
  #define ASEQUENCE2 ANIDENTIFIER ANARRAY2 ALIST ASEQ ATUPLE2
  #define ASEQUENCE3 ASEQ ANIDENTIFIER ATUPLE3 ALIST
  #define ASEQUENCE4
  #define ASEQUENCE5 ALIST ASEQ ATUPLE5 ANIDENTIFIER

  #define A_TUPLE_PLUS (mmf,34,^^,!) 456
  #define PLUS_ATUPLE yyt (j,ii%)
  #define JDATA ggh
  #define KDATA (a,b) name
  #define A_SEQ (25)(26)(27)

#if !BOOST_VMD_MSVC_V8

  #define AN_EMPTY_TUPLE_PLUS () 83

#endif

  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(2,ASEQUENCE,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(4,ASEQUENCE2,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST_EQ(BOOST_PP_LIST_AT(BOOST_VMD_ELEM(2,ASEQUENCE3,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE),3),3);
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,ASEQUENCE4,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(2,ASEQUENCE5,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));

  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,anything,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST_EQ(BOOST_VMD_ELEM(0,A_TUPLE_PLUS,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE),456);
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,PLUS_ATUPLE,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,JDATA,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST(!BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,KDATA,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));
  BOOST_TEST(BOOST_VMD_IS_EMPTY(BOOST_VMD_ELEM(0,A_SEQ,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE)));

#if !BOOST_VMD_MSVC_V8

  BOOST_TEST_EQ(BOOST_VMD_ELEM(0,AN_EMPTY_TUPLE_PLUS,BOOST_VMD_RETURN_ONLY_AFTER,BOOST_VMD_TYPE_TUPLE),83);

#endif

#else

BOOST_ERROR("No variadic macro support");

#endif

  return boost::report_errors();

  }

ext: cxx


text:
#ifndef SKYLARK_HASH_TRANSFORM_COMBBLAS_HPP
#define SKYLARK_HASH_TRANSFORM_COMBBLAS_HPP

#include <map>
#include <boost/serialization/map.hpp>

#include "../utility/external/FullyDistMultiVec.hpp"

namespace skylark { namespace sketch {

/* Specialization: FullyDistMultiVec for input, output */
template <typename IndexType,
          typename ValueType,
          template <typename> class IdxDistributionType,
          template <typename> class ValueDistribution>
struct hash_transform_t <FullyDistMultiVec<IndexType, ValueType>,
                         FullyDistMultiVec<IndexType, ValueType>,
                         IdxDistributionType,
                         ValueDistribution > :
        public hash_transform_data_t<IdxDistributionType,
                                     ValueDistribution> {
    typedef IndexType index_type;
    typedef ValueType value_type;
    typedef FullyDistMultiVec<IndexType, ValueType> matrix_type;
    typedef FullyDistMultiVec<IndexType, ValueType> output_matrix_type;
    typedef FullyDistVec<IndexType, ValueType> mpi_vector_t;
    typedef FullyDistMultiVec<IndexType, ValueType> mpi_multi_vector_t;
    typedef hash_transform_data_t<IdxDistributionType,
                                  ValueDistribution> data_type;

    /**
     * Regular constructor
     */
    hash_transform_t (int N, int S, base::context_t& context) :
        data_type (N, S, context) {}

    /**
     * Copy constructor
     */
    template <typename InputMatrixType,
              typename OutputMatrixType>
    hash_transform_t (hash_transform_t<InputMatrixType,
                                       OutputMatrixType,
                                       IdxDistributionType,
                                       ValueDistribution>& other) :
        data_type(other) {}

    /**
     * Constructor from data
     */
    hash_transform_t (hash_transform_data_t<IdxDistributionType,
                                            ValueDistribution>& other_data) :
        data_type(other_data) {}

    template <typename Dimension>
    void apply (const mpi_multi_vector_t &A,
        mpi_multi_vector_t &sketch_of_A,
        Dimension dimension) const {
        try {
            apply_impl (A, sketch_of_A, dimension);
        } catch(boost::mpi::exception e) {
            SKYLARK_THROW_EXCEPTION (
                base::mpi_exception()
                    << base::error_msg(e.what()) );
        } catch (std::string e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e) );
        } catch (std::logic_error e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e.what()) );
        }
    }

    int get_N() const { return this->_N; } /**< Get input dimension. */
    int get_S() const { return this->_S; } /**< Get output dimension. */

    const sketch_transform_data_t* get_data() const { return this; }

private:
    void apply_impl_single (const mpi_vector_t& a_,
                            mpi_vector_t& sketch_of_a,
                            columnwise_tag) const {
        std::vector<value_type> sketch_term(data_type::_S,0);

        // We are essentially doing a 'const' access to a, but the neccessary,
        // 'const' option is missing from the interface.
        mpi_vector_t &a = const_cast<mpi_vector_t&>(a_);

        /** Accumulate the local sketch vector */
        /** FIXME: Lot's of random access --- not good for performance */
        DenseVectorLocalIterator<index_type, value_type> local_iter(a);
        while(local_iter.HasNext()) {
            index_type idx = local_iter.GetLocIndex();
            index_type global_idx = local_iter.LocalToGlobal(idx);
            index_type global_sketch_idx = data_type::row_idx[global_idx];
            sketch_term[global_sketch_idx] +=
                (local_iter.GetValue()*data_type::row_value[global_idx]);
            local_iter.Next();
        }

        /** Accumulate the global sketch vector */
        /** FIXME: Only need to scatter ... don't need everything everywhere */
        MPI_Allreduce(MPI_IN_PLACE,
            &(sketch_term[0]),
            data_type::_S,
            MPIType<value_type>(),
            MPI_SUM,
            a.commGrid->GetWorld());

        /** Fill in .. SetElement() is dummy for non-local sets, so it's ok */
        for (index_type i=0; i<data_type::_S; ++i) {
            sketch_of_a.SetElement(i,sketch_term[i]);
        }
    }

    void apply_impl (const mpi_multi_vector_t& A,
        mpi_multi_vector_t& sketch_of_A,
        columnwise_tag) const {
        const index_type num_rhs = A.size;
        if (sketch_of_A.size != num_rhs) { /** error */; return; }
        if (A.dim != data_type::_N) { /** error */; return; }
        if (sketch_of_A.dim != data_type::_S) { /** error */; return; }

        /** FIXME: Can sketch all the vectors in one shot */
        for (index_type i=0; i<num_rhs; ++i) {
            apply_impl_single (A[i], sketch_of_A[i], columnwise_tag());
        }
    }
};

/* Specialization: SpParMat for input, output */
template <typename IndexType,
          typename ValueType,
          template <typename> class IdxDistributionType,
          template <typename> class ValueDistribution>
struct hash_transform_t <
    SpParMat<IndexType, ValueType, SpDCCols<IndexType, ValueType> >,
    SpParMat<IndexType, ValueType, SpDCCols<IndexType, ValueType> >,
    IdxDistributionType,
    ValueDistribution > :
        public hash_transform_data_t<IdxDistributionType,
                                     ValueDistribution> {
    typedef IndexType index_type;
    typedef ValueType value_type;
    typedef SpDCCols< index_type, value_type > col_t;
    typedef FullyDistVec< index_type, value_type> mpi_vector_t;
    typedef SpParMat< index_type, value_type, col_t > matrix_type;
    typedef SpParMat< index_type, value_type, col_t > output_matrix_type;
    typedef hash_transform_data_t<IdxDistributionType,
                                  ValueDistribution> data_type;

    /**
     * Regular constructor
     */
    hash_transform_t (int N, int S, base::context_t& context) :
        data_type(N, S, context) {

    }

    /**
     * Copy constructor
     */
    template <typename InputMatrixType,
              typename OutputMatrixType>
    hash_transform_t (hash_transform_t<InputMatrixType,
                                       OutputMatrixType,
                                       IdxDistributionType,
                                       ValueDistribution>& other) :
        data_type(other) {}

    /**
     * Constructor from data
     */
    hash_transform_t (hash_transform_data_t<IdxDistributionType,
                                            ValueDistribution>& other_data) :
        data_type(other_data) {}

    template <typename Dimension>
    void apply (const matrix_type &A,
        output_matrix_type &sketch_of_A,
        Dimension dimension) const {
        try {
            apply_impl (A, sketch_of_A, dimension);
        } catch(boost::mpi::exception e) {
            SKYLARK_THROW_EXCEPTION (
                base::mpi_exception()
                    << base::error_msg(e.what()) );
        } catch (std::string e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e) );
        } catch (std::logic_error e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e.what()) );
        } catch (std::bad_alloc e) {
            SKYLARK_THROW_EXCEPTION (
                base::skylark_exception()
                    << base::error_msg("bad_alloc: out of memory") );
        }
    }

    int get_N() const { return this->_N; } /**< Get input dimension. */
    int get_S() const { return this->_S; } /**< Get output dimension. */

    const sketch_transform_data_t* get_data() const { return this; }

private:
    /**
     * Apply the sketching transform that is described in by the sketch_of_A.
     */
    template <typename Dimension>
    void apply_impl (const matrix_type &A_,
        output_matrix_type &sketch_of_A,
        Dimension dist) const {

        // We are essentially doing a 'const' access to A, but the necessary,
        // 'const' option is missing from the interface
        matrix_type &A = const_cast<matrix_type&>(A_);

        const size_t rank = A.getcommgrid()->GetRank();

        // extract columns of matrix
        col_t &data = A.seq();

        const size_t ncols = sketch_of_A.getncol();

        const size_t my_row_offset = utility::cb_my_row_offset(A);
        const size_t my_col_offset = utility::cb_my_col_offset(A);

        size_t comm_size = A.getcommgrid()->GetSize();
        std::vector< std::set<size_t> > proc_set(comm_size);

        // pre-compute processor targets of local sketch application
        for(typename col_t::SpColIter col = data.begcol();
            col != data.endcol(); col++) {
            for(typename col_t::SpColIter::NzIter nz = data.begnz(col);
                nz != data.endnz(col); nz++) {

                // compute global row and column id, and compress in one
                // target position index
                const index_type rowid = nz.rowid()  + my_row_offset;
                const index_type colid = col.colid() + my_col_offset;
                const size_t pos       = getPos(rowid, colid, ncols, dist);

                // compute target processor for this target index
                const size_t target =
                    utility::owner(sketch_of_A, pos / ncols, pos % ncols);

                if(proc_set[target].count(pos) == 0) {
                    assert(target < comm_size);
                    proc_set[target].insert(pos);
                }
            }
        }

        // constructing arrays for one-sided access
        std::vector<index_type> proc_start_idx(comm_size + 1, 0);
        for(size_t i = 1; i < comm_size + 1; ++i)
            proc_start_idx[i] = proc_start_idx[i-1] + proc_set[i-1].size();

        std::vector<index_type> indicies(proc_start_idx[comm_size], 0);
        std::vector<value_type> values(proc_start_idx[comm_size], 0);

        // Apply sketch for all local values. Note that some of the resulting
        // values might end up on a different processor. The data structure
        // fills values (sorted by processor id) in one continuous array.
        // Subsequently, one-sided operations can be used to access values for
        // each processor.
        for(typename col_t::SpColIter col = data.begcol();
            col != data.endcol(); col++) {
            for(typename col_t::SpColIter::NzIter nz = data.begnz(col);
                nz != data.endnz(col); nz++) {

                // compute global row and column id, and compress in one
                // target position index
                const index_type rowid = nz.rowid()  + my_row_offset;
                const index_type colid = col.colid() + my_col_offset;
                const size_t pos       = getPos(rowid, colid, ncols, dist);

                // compute target processor for this target index
                const size_t proc =
                    utility::owner(sketch_of_A, pos / ncols, pos % ncols);

                // get offset in array for current element
                const size_t ar_idx = proc_start_idx[proc] +
                    std::distance(proc_set[proc].begin(), proc_set[proc].find(pos));

                indicies[ar_idx] = pos;
                values[ar_idx]  += nz.value() *
                                   data_type::getValue(rowid, colid, dist);
            }
        }

        // Creating windows for all relevant arrays
        boost::mpi::communicator comm = utility::get_communicator(A);

        // tell MPI that we will not use locks
        MPI_Info info;
        MPI_Info_create(&info);
        MPI_Info_set(info, "no_locks", "true");

        MPI_Win start_offset_win, idx_win, val_win;

        MPI_Win_create(&proc_start_idx[0], sizeof(size_t) * (comm_size + 1),
                       sizeof(size_t), info, comm, &start_offset_win);

        MPI_Win_create(&indicies[0], sizeof(index_type) * indicies.size(),
                       sizeof(index_type), info, comm, &idx_win);

        MPI_Win_create(&values[0], sizeof(value_type) * values.size(),
                       sizeof(value_type), info, comm, &val_win);

        MPI_Info_free(&info);

        // Synchronize epoch, no subsequent put operations (read only) and no
        // preceding fence calls.
        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOPRECEDE, start_offset_win);
        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOPRECEDE, idx_win);
        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOPRECEDE, val_win);

        // creating a local structure to hold sparse data triplets
        std::vector< tuple< index_type, index_type, value_type > > sketch_data;

        std::map<index_type, size_t> tuple_idx;
        typedef typename std::map<index_type, size_t>::iterator itr_t;

        // in order to compute local indices in sketch_of_A we need to know
        // row and column offset on each rank.
        const size_t res_row_offset = utility::cb_my_row_offset(sketch_of_A);
        const size_t res_col_offset = utility::cb_my_col_offset(sketch_of_A);

        for(size_t p = 0; p < comm_size; ++p) {

            // get the start/end offset
            std::vector<size_t> offset(2);
            MPI_Get(&(offset[0]), 2, boost::mpi::get_mpi_datatype<size_t>(),
                    p, rank, 2, boost::mpi::get_mpi_datatype<size_t>(),
                    start_offset_win);

            MPI_Win_fence(MPI_MODE_NOPUT, start_offset_win);
            size_t num_values = offset[1] - offset[0];

            // and fill indices/values.
            std::vector<index_type> add_idx(num_values);
            std::vector<value_type> add_val(num_values);
            MPI_Get(&(add_idx[0]), num_values,
                    boost::mpi::get_mpi_datatype<index_type>(), p, offset[0],
                    num_values, boost::mpi::get_mpi_datatype<index_type>(),
                    idx_win);

            MPI_Get(&(add_val[0]), num_values,
                    boost::mpi::get_mpi_datatype<value_type>(), p, offset[0],
                    num_values, boost::mpi::get_mpi_datatype<value_type>(),
                    val_win);

            MPI_Win_fence(MPI_MODE_NOPUT, idx_win);
            MPI_Win_fence(MPI_MODE_NOPUT, val_win);

            // finally, add data to local CombBLAS buffer (if we have any).
            //XXX: We cannot have duplicated (row/col) pairs in sketch_data.
            //     Check again with CombBLAS if there is a more suitable way
            //     to create sparse dist matrices.
            for(size_t i = 0; i < num_values; ++i) {

                index_type lrow = add_idx[i] / ncols - res_row_offset;
                index_type lcol = add_idx[i] % ncols - res_col_offset;

                itr_t itr = tuple_idx.find(add_idx[i]);
                if(itr == tuple_idx.end()) {
                    tuple_idx.insert(
                        std::make_pair(add_idx[i], sketch_data.size()));
                    sketch_data.push_back(make_tuple(lrow, lcol, add_val[i]));
                } else {
                    get<2>(sketch_data[itr->second]) += add_val[i];
                }
            }
        }

        // this pointer will be freed in the destructor of col_t (see below).
        //FIXME: verify with Valgrind
        SpTuples<index_type, value_type> *tmp_tpl =
            new SpTuples<index_type, value_type> (
                sketch_data.size(), sketch_of_A.getlocalrows(),
                sketch_of_A.getlocalcols(), &sketch_data[0]);

        // create temporary matrix (no further communication should be
        // required because all the data are local) and assign (deep copy).
        col_t *sp_data = new col_t(*tmp_tpl, false);

        //FIXME: is there a direct way to set the "data buffer" for the output
        //       matrix? Currently this method creates a temporary matrix and
        //       then assigns to the output matrix (deep copy).
        sketch_of_A = output_matrix_type(sp_data, sketch_of_A.getcommgrid());

        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOSUCCEED, start_offset_win);
        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOSUCCEED, idx_win);
        MPI_Win_fence(MPI_MODE_NOPUT | MPI_MODE_NOSUCCEED, val_win);

        MPI_Win_free(&start_offset_win);
        MPI_Win_free(&idx_win);
        MPI_Win_free(&val_win);
    }

    inline index_type getPos(index_type rowid, index_type colid, size_t ncols,
        columnwise_tag) const {
        return colid + ncols * data_type::row_idx[rowid];
    }

    inline index_type getPos(index_type rowid, index_type colid, size_t ncols,
        rowwise_tag) const {
        return rowid * ncols + data_type::row_idx[colid];
    }
};

/* Specialization: SpParMat for input, Local SpMat output */
template <typename IndexType,
          typename ValueType,
          template <typename> class IdxDistributionType,
          template <typename> class ValueDistribution>
struct hash_transform_t <
    SpParMat<IndexType, ValueType, SpDCCols<IndexType, ValueType> >,
    base::sparse_matrix_t<ValueType>,
    IdxDistributionType,
    ValueDistribution > :
        public hash_transform_data_t<IdxDistributionType,
                                     ValueDistribution> {
    typedef IndexType index_type;
    typedef ValueType value_type;
    typedef SpDCCols< index_type, value_type > col_t;
    typedef FullyDistVec< index_type, value_type> mpi_vector_t;
    typedef SpParMat< index_type, value_type, col_t > matrix_type;
    typedef base::sparse_matrix_t< value_type > output_matrix_type;
    typedef hash_transform_data_t<IdxDistributionType,
                                  ValueDistribution> data_type;

    /**
     * Regular constructor
     */
    hash_transform_t (int N, int S, base::context_t& context) :
        data_type(N, S, context) {

    }

    /**
     * Copy constructor
     */
    template <typename InputMatrixType,
              typename OutputMatrixType>
    hash_transform_t (hash_transform_t<InputMatrixType,
                                       OutputMatrixType,
                                       IdxDistributionType,
                                       ValueDistribution>& other) :
        data_type(other) {}

    /**
     * Constructor from data
     */
    hash_transform_t (hash_transform_data_t<IdxDistributionType,
                                            ValueDistribution>& other_data) :
        data_type(other_data) {}

    template <typename Dimension>
    void apply (const matrix_type &A,
        output_matrix_type &sketch_of_A,
        Dimension dimension) const {
        try {
            apply_impl (A, sketch_of_A, dimension);
        } catch(boost::mpi::exception e) {
            SKYLARK_THROW_EXCEPTION (
                base::mpi_exception()
                    << base::error_msg(e.what()) );
        } catch (std::string e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e) );
        } catch (std::logic_error e) {
            SKYLARK_THROW_EXCEPTION (
                base::combblas_exception()
                    << base::error_msg(e.what()) );
        } catch (std::bad_alloc e) {
            SKYLARK_THROW_EXCEPTION (
                base::skylark_exception()
                    << base::error_msg("bad_alloc: out of memory") );
        }
    }

    int get_N() const { return this->_N; } /**< Get input dimension. */
    int get_S() const { return this->_S; } /**< Get output dimension. */

    const sketch_transform_data_t* get_data() const { return this; }

private:
    /**
     * Apply the sketching transform that is described in by the sketch_of_A.
     */
    template <typename Dimension>
    void apply_impl (const matrix_type &A_,
        output_matrix_type &sketch_of_A,
        Dimension dist) const {

        // We are essentially doing a 'const' access to A, but the necessary,
        // 'const' option is missing from the interface
        matrix_type &A = const_cast<matrix_type&>(A_);

        const size_t rank = A.getcommgrid()->GetRank();

        // extract columns of matrix
        col_t &data = A.seq();

        const size_t my_row_offset = utility::cb_my_row_offset(A);
        const size_t my_col_offset = utility::cb_my_col_offset(A);

        int n_res_cols = A.getncol();
        int n_res_rows = A.getnrow();
        data_type::get_res_size(n_res_rows, n_res_cols, dist);

        // Apply sketch for all local values. Subsequently, all values are
        // gathered on processor 0 and the local matrix is populated.
        // XXX: For now sort on each processor to ease the creation of the
        // gathered local matrix. Most likely possible can be avoided by
        // traversing the matrix in the correct order.
        typedef std::map<index_type, value_type> col_values_t;
        col_values_t col_values;
        for(typename col_t::SpColIter col = data.begcol();
            col != data.endcol(); col++) {
            for(typename col_t::SpColIter::NzIter nz = data.begnz(col);
                nz != data.endnz(col); nz++) {

                index_type rowid = nz.rowid()  + my_row_offset;
                index_type colid = col.colid() + my_col_offset;

                const value_type value =
                    nz.value() * data_type::getValue(rowid, colid, dist);
                data_type::finalPos(rowid, colid, dist);
                col_values[colid * n_res_rows + rowid] += value;
            }
        }

        boost::mpi::communicator world(A.getcommgrid()->GetWorld(),
                                       boost::mpi::comm_duplicate);

        std::vector< std::map<index_type, value_type > > result;
        boost::mpi::gather(world, col_values, result, 0);

        // unpack into temp structure
        //FIXME:
        int size_estimate = A.getnnz();
        if(rank == 0) {
            create_local_sp_mat(result, n_res_rows, n_res_cols, size_estimate,
                                sketch_of_A);
        }
    }

    void create_local_sp_mat(
            std::vector< std::map<index_type, value_type > > &result,
            int n_res_rows, int n_res_cols, int size_estimate,
            output_matrix_type &sketch_of_A) const {

        int nnz = 0;
        int *indptr_new = new int[n_res_cols + 1];
        std::vector<int> final_rows(size_estimate);
        std::vector<value_type> final_vals(size_estimate);

        indptr_new[0] = 0;

        typedef typename std::map<index_type, value_type>::iterator itr_t;
        std::vector<itr_t> proc_iters(result.size());
        for(size_t i = 0; i < result.size(); ++i)
            proc_iters[i] = result[i].begin();

        std::vector<index_type> idx_map(n_res_rows, -1);

        for(int col = 0; col < n_res_cols; ++col) {

            // collect all values for column 'col' of all procs
            for(size_t i = 0; i < result.size(); ++i) {

                itr_t cur_itr = proc_iters[i];

                for(; cur_itr != result[i].end() &&
                        static_cast<int>(cur_itr->first / n_res_rows) == col;
                        cur_itr++) {

                    int row    = cur_itr->first % n_res_rows;
                    double val = cur_itr->second;

                    if(idx_map[row] == -1) {
                        idx_map[row] = nnz;
                        final_rows[nnz] = row;
                        final_vals[nnz] = val;
                        nnz++;
                    } else {
                        final_vals[idx_map[row]] += val;
                    }
                }

                proc_iters[i] = cur_itr;
            }

            indptr_new[col + 1] = nnz;

            // reset idx_map
            for(int i = indptr_new[col]; i < nnz; ++i)
                idx_map[final_rows[i]] = -1;
        }

        int *indices_new = new int[nnz];
        std::copy(final_rows.begin(), final_rows.begin() + nnz, indices_new);

        double *values_new = new double[nnz];
        std::copy(final_vals.begin(), final_vals.begin() + nnz, values_new);

        sketch_of_A.attach(indptr_new, indices_new, values_new,
                nnz, n_res_rows, n_res_cols, true);
    }
};

} } /** namespace skylark::sketch */

#endif // SKYLARK_HASH_TRANSFORM_COMBBLAS_HPP

max_issues_repo_name: cjiyer/libskylark

max_forks_repo_name: cjiyer/libskylark


text: #include <boost/fusion/functional/invocation/invoke_procedure.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
// uses boost::synchronized_value to link variable and mutex
#define BOOST_THREAD_VERSION 4
#include <boost/thread.hpp>
#include <boost/thread/synchronized_value.hpp> // include required
#include <boost/asio.hpp>
#include <string>
#include <sstream>
#include <iostream>
#include <functional>

std::string get_robotstxt(const std::string &host)
{
	using namespace boost::asio;
	io_service ioservice;
	ip::tcp::resolver resolver(ioservice);
	ip::tcp::resolver::query query(host, "http");
	auto it = resolver.resolve(query);
	ip::tcp::socket socket(ioservice);
	socket.connect(*it);
	std::string request = "GET /robots.txt HTTP/1.1\r\nHost: " + host + "\r\n\r\n";
	write(socket, buffer(request));
	streambuf response;
	boost::system::error_code ec;
	read(socket, response, transfer_all(), ec);
	if (ec == error::eof)
	{
		std::ostringstream os;
		os << &response;
		std::string s = os.str();
		std::size_t idx = s.find("\r\n\r\n");
		if (idx != std::string::npos)
			s.erase(0, idx + 4);
		return s;
	}
	return "";
}

// lines is a boost::synchronized_value with a vector of strings and a mutex
boost::synchronized_value<std::vector<std::string>, boost::mutex> lines;

void store_lines_from_robotstxt(const std::string &host)
{
	std::string robotstxt = get_robotstxt(host);
	std::istringstream is(robotstxt);
	std::string s;
	while (std::getline(is, s))
	{
		// synchronize() returns a locked pointer of type boost::strict_lock_ptr
		// which provides synchronized access to the vector of strings; the destructor
		// of boost::strict_lock_ptr unlocks the mutex automatically
		boost::strict_lock_ptr<std::vector<std::string>, boost::mutex> p = lines.synchronize();
		p->emplace_back(s);
	}
}

int main()
{
	boost::thread_group group;
	group.create_thread(std::bind(store_lines_from_robotstxt, "theboostcpplibraries.com"));
	group.create_thread(std::bind(store_lines_from_robotstxt, "dieboostcppbibliotheken.de"));
	group.join_all();
	// boost::synchronized_value provides a member function value() to access
	// the vector of strings without having to use the mutex (which is not required
	// here as all threads have been terminated)
	for (auto &line : lines.value())
		std::cout << line << '\n';
}


text:
/*
 * Copyright 2015 David A. Boyuka II
 *
 * 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.
 */

/*
 * timing.h
 *
 *  Created on: Jun 20, 2014
 *      Author: David A. Boyuka II
 */
#ifndef TIMING_H_
#define TIMING_H_

#include <utility> // std::move
#include <boost/timer/timer.hpp>

class TimeBlock {
private:
	static constexpr double SEC_PER_NANO = 1e-9;
public:
	TimeBlock(double *add_wallclock, double *add_cpu = nullptr) : add_wallclock(add_wallclock), add_cpu(add_cpu), timer() {}
	~TimeBlock() {
		timer.stop();
		const boost::timer::cpu_times times = timer.elapsed();
		if (add_wallclock)
			*add_wallclock	+= times.wall * SEC_PER_NANO;
		if (add_cpu)
			*add_cpu 		+= times.user * SEC_PER_NANO;
	}
	TimeBlock(const TimeBlock &) = delete;
	TimeBlock(TimeBlock &&other) : add_wallclock(other.add_wallclock), add_cpu(other.add_cpu), timer(std::move(other.timer)) {
		other.add_wallclock = nullptr;
		other.add_cpu = nullptr;
	}

	TimeBlock & operator=(const TimeBlock &) = delete;
	TimeBlock & operator=(TimeBlock &&other) {
		this->timer = std::move(other.timer);
		this->add_wallclock = other.add_wallclock;
		this->add_cpu = other.add_cpu;
		other.add_wallclock = nullptr;
		other.add_cpu = nullptr;
		return *this;
	}
private:
	double *add_wallclock, *add_cpu;
	boost::timer::cpu_timer timer;
};

#define TIME_BLOCK_BEGIN(wc, cpu) { TimeBlock tb = TimeBlock(wc, cpu);
#define TIME_BLOCK_END() }

#endif /* TIMING_H_ */

max_stars_repo_name: daboyuka/PIQUE

max_issues_repo_name: daboyuka/PIQUE

max_forks_repo_name: daboyuka/PIQUE


text:
// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "Domain/CoordinateMaps/BulgedCube.hpp"

#include <boost/none.hpp>
#include <boost/optional.hpp>
#include <cmath>
#include <exception>
#include <functional>  // for std::reference_wrapper
#include <limits>
#include <pup.h>

#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "ErrorHandling/Assert.hpp"
#include "NumericalAlgorithms/RootFinding/TOMS748.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/DereferenceWrapper.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/Gsl.hpp"
#include "Utilities/MakeWithValue.hpp"

namespace {
class RootFunction {
 public:
  RootFunction(const double radius, const double sphericity,
               const double physical_r_squared, const double x_sq,
               const double y_sq, const double z_sq) noexcept
      : radius_(radius),
        sphericity_(sphericity),
        physical_r_squared_(physical_r_squared),
        x_sq_(x_sq),
        y_sq_(y_sq),
        z_sq_(z_sq) {
    ASSERT(physical_r_squared_ > 0.0,
           "The RootFunction for the inverse map of BulgedCube is only valid "
           "for a radius greater than zero. The zero-radius case should be "
           "handled explicitly without call the root finder.");
  }

  double operator()(const double rho) const noexcept {
    const double x_sq_over_r_sq = x_sq_ / physical_r_squared_;
    const double y_sq_over_r_sq = y_sq_ / physical_r_squared_;
    const double z_sq_over_r_sq = z_sq_ / physical_r_squared_;
    return sqrt(physical_r_squared_) -
           radius_ * rho *
               (1.0 / sqrt(3.0) +
                sphericity_ *
                    (1.0 / sqrt(1.0 + square(rho) *
                                          (x_sq_over_r_sq + y_sq_over_r_sq)) +
                     1.0 / sqrt(1.0 + square(rho) *
                                          (x_sq_over_r_sq + z_sq_over_r_sq)) +
                     1.0 / sqrt(1.0 + square(rho) *
                                          (y_sq_over_r_sq + z_sq_over_r_sq)) -
                     1.0 / sqrt(2.0 + square(rho) * x_sq_over_r_sq) -
                     1.0 / sqrt(2.0 + square(rho) * y_sq_over_r_sq) -
                     1.0 / sqrt(2.0 + square(rho) * z_sq_over_r_sq)));
  }
  double get_r_sq() const noexcept { return physical_r_squared_; }

 private:
  const double radius_;
  const double sphericity_;
  const double physical_r_squared_;
  const double x_sq_;
  const double y_sq_;
  const double z_sq_;
};

boost::optional<double> scaling_factor(RootFunction&& rootfunction) noexcept {
  const double physical_r_squared = rootfunction.get_r_sq();
  try {
    constexpr double tol = 10.0 * std::numeric_limits<double>::epsilon();
    // Use a small nonzero number since the inverse map is singular at r==0 and
    // that case is handled separately.
    constexpr double lower_bound = std::numeric_limits<double>::min();
    // upper_bound = sqrt(3) + tol
    constexpr double upper_bound = 1.7320508075688772 + tol;
    double rho =
        // NOLINTNEXTLINE(clang-analyzer-core)
        RootFinder::toms748(rootfunction, lower_bound, upper_bound, tol, tol);
    rho /= sqrt(physical_r_squared);
    return rho;
  } catch (std::exception& exception) {
    return boost::none;
  }
}
}  // namespace

namespace domain {
namespace CoordinateMaps {
BulgedCube::BulgedCube(const double radius, const double sphericity,
                       const bool use_equiangular_map) noexcept
    : radius_(radius),
      sphericity_(sphericity),
      use_equiangular_map_(use_equiangular_map),
      is_identity_(radius_ == sqrt(3.0) and sphericity_ == 0.0 and
                   not use_equiangular_map_) {
  ASSERT(radius > 0.0, "The radius of the cube must be greater than zero");
  ASSERT(sphericity >= 0.0 and sphericity < 1.0,
         "The sphericity must be strictly less than one.");
}

template <typename T>
std::array<tt::remove_cvref_wrap_t<T>, 3> BulgedCube::operator()(
    const std::array<T, 3>& source_coords) const noexcept {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  const auto physical_coordinates = [this](
      const ReturnType& cap_xi, const ReturnType& cap_eta,
      const ReturnType& cap_zeta) noexcept {
    const auto one_over_rho_xi = 1.0 / sqrt(2.0 + square(cap_xi));
    const auto one_over_rho_eta = 1.0 / sqrt(2.0 + square(cap_eta));
    const auto one_over_rho_zeta = 1.0 / sqrt(2.0 + square(cap_zeta));
    const auto one_over_rho_xi_eta =
        1.0 / sqrt(1.0 + square(cap_xi) + square(cap_eta));
    const auto one_over_rho_xi_zeta =
        1.0 / sqrt(1.0 + square(cap_xi) + square(cap_zeta));
    const auto one_over_rho_eta_zeta =
        1.0 / sqrt(1.0 + square(cap_eta) + square(cap_zeta));
    const ReturnType radial_scaling_factor =
        radius_ * (1.0 / sqrt(3.0) +
                   sphericity_ * (one_over_rho_eta_zeta + one_over_rho_xi_zeta +
                                  one_over_rho_xi_eta - one_over_rho_xi -
                                  one_over_rho_eta - one_over_rho_zeta));

    ReturnType physical_x = radial_scaling_factor * cap_xi;
    ReturnType physical_y = radial_scaling_factor * cap_eta;
    ReturnType physical_z = radial_scaling_factor * cap_zeta;
    return std::array<ReturnType, 3>{
        {std::move(physical_x), std::move(physical_y), std::move(physical_z)}};
  };

  if (use_equiangular_map_) {
    return physical_coordinates(
        tan(M_PI_4 * dereference_wrapper(source_coords[0])),
        tan(M_PI_4 * dereference_wrapper(source_coords[1])),
        tan(M_PI_4 * dereference_wrapper(source_coords[2])));
  }
  return physical_coordinates(dereference_wrapper(source_coords[0]),
                              dereference_wrapper(source_coords[1]),
                              dereference_wrapper(source_coords[2]));
}

boost::optional<std::array<double, 3>> BulgedCube::inverse(
    const std::array<double, 3>& target_coords) const noexcept {
  const double& physical_x = target_coords[0];
  const double& physical_y = target_coords[1];
  const double& physical_z = target_coords[2];
  const double x_sq = square(physical_x);
  const double y_sq = square(physical_y);
  const double z_sq = square(physical_z);
  const double physical_r_squared = x_sq + y_sq + z_sq;
  // For physical_r_squared==0 we know that we are at the origin x=y=z=0, and we
  // know analytically that this map maps the origin to itself. Handling the
  // case r==0 case separately simplifies the root finding procedure for the
  // scaling_factor.
  if (UNLIKELY(physical_r_squared == 0.0)) {
    return {{{0.0, 0.0, 0.0}}};
  }

  // We are not at the origin, find the scaling factor (does a root-find)
  const auto scaling_factor =
      // NOLINTNEXTLINE(clang-analyzer-core)
      ::scaling_factor(
      RootFunction{radius_, sphericity_, physical_r_squared, x_sq, y_sq, z_sq});
  if (not scaling_factor) {
    return boost::none;
  }
  if (use_equiangular_map_) {
    return {{{2.0 * M_2_PI * atan(physical_x * scaling_factor.get()),
              2.0 * M_2_PI * atan(physical_y * scaling_factor.get()),
              2.0 * M_2_PI * atan(physical_z * scaling_factor.get())}}};
  }
  return {
      {{physical_x * scaling_factor.get(), physical_y * scaling_factor.get(),
        physical_z * scaling_factor.get()}}};
}

template <typename T>
std::array<tt::remove_cvref_wrap_t<T>, 3> BulgedCube::xi_derivative(
    const std::array<T, 3>& source_coords) const noexcept {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  const auto derivative_lambda = [this](
      const ReturnType& cap_xi, const ReturnType& cap_eta,
      const ReturnType& cap_zeta, const auto& cap_xi_deriv) {
    const ReturnType one_over_rho_xi_cubed =
        pow<3>(1.0 / sqrt(2.0 + square(cap_xi)));
    const ReturnType one_over_rho_eta = 1.0 / sqrt(2.0 + square(cap_eta));
    const ReturnType one_over_rho_zeta = 1.0 / sqrt(2.0 + square(cap_zeta));
    const ReturnType one_over_rho_xi_eta_cubed =
        pow<3>(1.0 / sqrt(1.0 + square(cap_xi) + square(cap_eta)));
    const ReturnType one_over_rho_xi_zeta_cubed =
        pow<3>(1.0 / sqrt(1.0 + square(cap_xi) + square(cap_zeta)));
    const ReturnType one_over_rho_eta_zeta =
        1.0 / sqrt(1.0 + square(cap_eta) + square(cap_zeta));
    const ReturnType common_factor =
        sphericity_ * radius_ * cap_xi * cap_xi_deriv *
        (one_over_rho_xi_cubed - one_over_rho_xi_eta_cubed -
         one_over_rho_xi_zeta_cubed);

    const ReturnType physical_x =
        radius_ * cap_xi_deriv *
        (1.0 / sqrt(3.0) +
         sphericity_ *
             (((1.0 + square(cap_eta)) * one_over_rho_xi_eta_cubed +
               (1.0 + square(cap_zeta)) * one_over_rho_xi_zeta_cubed -
               2.0 * one_over_rho_xi_cubed) +
              one_over_rho_eta_zeta - one_over_rho_eta - one_over_rho_zeta));
    const ReturnType physical_y = cap_eta * common_factor;
    const ReturnType physical_z = cap_zeta * common_factor;

    return std::array<ReturnType, 3>{{physical_x, physical_y, physical_z}};
  };
  if (use_equiangular_map_) {
    return derivative_lambda(
        tan(M_PI_4 * dereference_wrapper(source_coords[0])),
        tan(M_PI_4 * dereference_wrapper(source_coords[1])),
        tan(M_PI_4 * dereference_wrapper(source_coords[2])),
        ReturnType{M_PI_4 *
                   (1.0 + square(tan(M_PI_4 *
                                     dereference_wrapper(source_coords[0]))))});
  }
  return derivative_lambda(dereference_wrapper(source_coords[0]),
                           dereference_wrapper(source_coords[1]),
                           dereference_wrapper(source_coords[2]), 1.0);
}

template <typename T>
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> BulgedCube::jacobian(
    const std::array<T, 3>& source_coords) const noexcept {
  const auto dX_dxi = xi_derivative(source_coords);
  const auto dX_deta = xi_derivative(
      std::array<std::reference_wrapper<const tt::remove_cvref_wrap_t<T>>, 3>{
          {std::cref(dereference_wrapper(source_coords[1])),
           std::cref(dereference_wrapper(source_coords[0])),
           std::cref(dereference_wrapper(source_coords[2]))}});
  const auto dX_dzeta = xi_derivative(
      std::array<std::reference_wrapper<const tt::remove_cvref_wrap_t<T>>, 3>{
          {std::cref(dereference_wrapper(source_coords[2])),
           std::cref(dereference_wrapper(source_coords[1])),
           std::cref(dereference_wrapper(source_coords[0]))}});
  auto jacobian_matrix =
      make_with_value<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>>(
          dereference_wrapper(source_coords[0]), 0.0);

  get<0, 0>(jacobian_matrix) = dX_dxi[0];
  get<0, 1>(jacobian_matrix) = dX_deta[1];
  get<0, 2>(jacobian_matrix) = dX_dzeta[2];
  get<1, 0>(jacobian_matrix) = dX_dxi[1];
  get<1, 1>(jacobian_matrix) = dX_deta[0];
  get<1, 2>(jacobian_matrix) = dX_dzeta[1];
  get<2, 0>(jacobian_matrix) = dX_dxi[2];
  get<2, 1>(jacobian_matrix) = dX_deta[2];
  get<2, 2>(jacobian_matrix) = dX_dzeta[0];
  return jacobian_matrix;
}

template <typename T>
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>
BulgedCube::inv_jacobian(const std::array<T, 3>& source_coords) const noexcept {
  const auto jac = jacobian(source_coords);
  return determinant_and_inverse(jac).second;
}

void BulgedCube::pup(PUP::er& p) noexcept {
  p | radius_;
  p | sphericity_;
  p | use_equiangular_map_;
  p | is_identity_;
}

bool operator==(const BulgedCube& lhs, const BulgedCube& rhs) noexcept {
  return lhs.radius_ == rhs.radius_ and lhs.sphericity_ == rhs.sphericity_ and
         lhs.use_equiangular_map_ == rhs.use_equiangular_map_ and
         lhs.is_identity_ == rhs.is_identity_;
}

bool operator!=(const BulgedCube& lhs, const BulgedCube& rhs) noexcept {
  return not(lhs == rhs);
}

// Explicit instantiations
/// \cond
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)

#define INSTANTIATE(_, data)                                                  \
  template std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3> BulgedCube::   \
  operator()(const std::array<DTYPE(data), 3>& source_coords) const noexcept; \
  template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>  \
  BulgedCube::jacobian(const std::array<DTYPE(data), 3>& source_coords)       \
      const noexcept;                                                         \
  template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>  \
  BulgedCube::inv_jacobian(const std::array<DTYPE(data), 3>& source_coords)   \
      const noexcept;

GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector,
                                      std::reference_wrapper<const double>,
                                      std::reference_wrapper<const DataVector>))

#undef DTYPE
#undef INSTANTIATE
/// \endcond
}  // namespace CoordinateMaps
}  // namespace domain

max_stars_repo_name: tomwlodarczyk/spectre

max_issues_repo_name: tomwlodarczyk/spectre

max_forks_repo_name: tomwlodarczyk/spectre


text:
/*
 * libasiotap - A portable TAP adapter extension for Boost::ASIO.
 * Copyright (C) 2010-2011 Julien KAUFFMANN <julien.kauffmann@freelan.org>
 *
 * This file is part of libasiotap.
 *
 * libasiotap 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.
 *
 * libasiotap 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/>.
 *
 * In addition, as a special exception, the copyright holders give
 * permission to link the code of portions of this program with the
 * OpenSSL library under certain conditions as described in each
 * individual source file, and distribute linked combinations
 * including the two.
 * You must obey the GNU General Public License in all respects
 * for all of the code used other than OpenSSL.  If you modify
 * file(s) with this exception, you may extend this exception to your
 * version of the file(s), but you are not obligated to do so.  If you
 * do not wish to do so, delete this exception statement from your
 * version.  If you delete this exception statement from all source
 * files in the program, then also delete it here.
 *
 * If you intend to use libasiotap in a commercial software, please
 * contact me : we may arrange this for a small fee or no fee at all,
 * depending on the nature of your project.
 */

/**
 * \file posix_route_manager.cpp
 * \author Julien KAUFFMANN <julien.kauffmann@freelan.org>
 * \brief The POSIX route manager class.
 */

#include "posix/posix_route_manager.hpp"

#include "error.hpp"

#include <sstream>

#include <boost/lexical_cast.hpp>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/algorithm/string.hpp>

#ifdef LINUX
#include <netlinkplus/manager.hpp>
#endif

#include <executeplus/posix_system.hpp>

namespace asiotap
{
	posix_route_manager::route_type posix_route_manager::get_route_for(const boost::asio::ip::address& host)
	{
#ifdef MACINTOSH
		const std::vector<std::string> real_args { "/sbin/route", "-n", "get", boost::lexical_cast<std::string>(host) };

		std::stringstream ss;
		executeplus::checked_execute(real_args, executeplus::get_current_environment(), &ss);

		//The output is like:
		/*
			   route to: 8.8.8.8
			destination: default
				   mask: default
				gateway: 10.7.0.254
			  interface: en0
				  flags: <UP,GATEWAY,DONE,STATIC,PRCLONING>
			 recvpipe  sendpipe  ssthresh  rtt,msec    rttvar  hopcount      mtu     expire
				   0         0         0         0         0         0      1500         0
		*/
		std::map<std::string, std::string> values;
		std::string line;

		while (std::getline(ss, line))
		{
			const auto colon_pos = line.find_first_of(':');

			if (colon_pos == std::string::npos)
			{
				continue;
			}

			std::string key = line.substr(0, colon_pos);
			std::string value = line.substr(colon_pos + 1, std::string::npos);

			boost::algorithm::trim(key);
			boost::algorithm::trim(value);

			values[key] = value;
		}

		if (values.find("interface") == values.end())
		{
			throw boost::system::system_error(make_error_code(asiotap_error::external_process_output_parsing_error));
		}

		const std::string interface = values["interface"];
		boost::optional<boost::asio::ip::address> gw;

		if (values.find("gateway") != values.end())
		{
			gw = boost::asio::ip::address::from_string(values["gateway"]);
		}
#else
#ifdef FREELAN_DISABLE_NETLINK
		const std::vector<std::string> real_args { "/bin/ip", "route", "get", boost::lexical_cast<std::string>(host) };

		std::stringstream ss;
		executeplus::checked_execute(real_args, executeplus::get_current_environment(), &ss);

		// The output is like:
		/*
			8.8.8.8 via 37.59.15.254 dev eth0  src 46.105.57.112
				cache
		*/

		std::string host_addr;

		if (!(ss >> host_addr))
		{
			throw boost::system::system_error(make_error_code(asiotap_error::external_process_output_parsing_error));
		}

		std::map<std::string, std::string> values;
		std::string key, value;

		while (ss >> key >> value)
		{
			boost::algorithm::trim(key);
			boost::algorithm::trim(value);

			values[key] = value;
		}

		if (values.find("dev") == values.end())
		{
			throw boost::system::system_error(make_error_code(asiotap_error::external_process_output_parsing_error));
		}

		const std::string interface = values["dev"];
		boost::optional<boost::asio::ip::address> gw;

		if (values.find("via") != values.end())
		{
			gw = boost::asio::ip::address::from_string(values["via"]);
		}
#else
		const auto route_info = m_netlink_manager.get_route_for(host);
		const auto interface = route_info.output_interface.name();
		const auto gw = route_info.gateway;
#endif
#endif

		const auto route = to_ip_route(to_network_address(host), gw);
		const posix_route_manager::route_type route_entry = { interface, route, 0 };

		return route_entry;
	}

	void posix_route_manager::ifconfig(const std::string& interface, const ip_network_address& address)
	{
#if defined(MACINTOSH) || defined(FREELAN_DISABLE_NETLINK)
		const std::vector<std::string> real_args { "/sbin/ifconfig", interface, boost::lexical_cast<std::string>(address) };

		executeplus::checked_execute(real_args, executeplus::get_current_environment());
#else
		m_netlink_manager.add_interface_address(netlinkplus::interface_entry(interface), to_ip_address(address), to_prefix_length(address));
#endif
	}

	void posix_route_manager::ifconfig(const std::string& interface, const ip_network_address& address, const boost::asio::ip::address& remote_address)
	{
#if defined(MACINTOSH) || defined(FREELAN_DISABLE_NETLINK)
		const std::vector<std::string> real_args { "/sbin/ifconfig", interface, boost::lexical_cast<std::string>(address), boost::lexical_cast<std::string>(remote_address) };

		executeplus::checked_execute(real_args, executeplus::get_current_environment());
#else
		m_netlink_manager.add_interface_address(netlinkplus::interface_entry(interface), remote_address, to_prefix_length(address), to_ip_address(address));
#endif
	}

	void posix_route_manager::set_route(route_action action, const std::string& interface, const ip_network_address& dest)
	{
#if defined(MACINTOSH) || defined(FREELAN_DISABLE_NETLINK)
		const std::string net_host = is_unicast(dest) ? "-host" : "-net";
#ifdef MACINTOSH
		const std::string command = action == route_action::add ? "add" : "delete";
		const std::vector<std::string> real_args { "/sbin/route", "-n", command, net_host, boost::lexical_cast<std::string>(dest), "-interface", interface };
#else
		const std::string command = action == route_action::add ? "add" : "del";
		const std::vector<std::string> real_args { "/sbin/route", "-n", command, net_host, boost::lexical_cast<std::string>(dest), "dev", interface };
#endif

		executeplus::checked_execute(real_args, executeplus::get_current_environment());
#else
		if (action == route_action::add)
		{
			m_netlink_manager.add_route(netlinkplus::interface_entry(interface), to_ip_address(dest), to_prefix_length(dest));
		}
		else
		{
			m_netlink_manager.remove_route(netlinkplus::interface_entry(interface), to_ip_address(dest), to_prefix_length(dest));
		}
#endif
	}

	void posix_route_manager::set_route(route_action action, const std::string& interface, const ip_network_address& dest, const boost::asio::ip::address& gateway)
	{
#if defined(MACINTOSH) || defined(FREELAN_DISABLE_NETLINK)
		const std::string net_host = is_unicast(dest) ? "-host" : "-net";
		const std::string command = action == route_action::add ? "add" : "del";
#ifdef MACINTOSH
		static_cast<void>(interface);
		const std::vector<std::string> real_args { "/sbin/route", "-n", command, net_host, boost::lexical_cast<std::string>(dest), boost::lexical_cast<std::string>(gateway) };
#else
		const std::vector<std::string> real_args { "/sbin/route", "-n", command, net_host, boost::lexical_cast<std::string>(dest), "gw", boost::lexical_cast<std::string>(gateway), "dev", interface };
#endif

		executeplus::checked_execute(real_args, executeplus::get_current_environment());
#else
		if (action == route_action::add)
		{
			m_netlink_manager.add_route(netlinkplus::interface_entry(interface), to_ip_address(dest), to_prefix_length(dest), gateway);
		}
		else
		{
			m_netlink_manager.remove_route(netlinkplus::interface_entry(interface), to_ip_address(dest), to_prefix_length(dest), gateway);
		}
#endif
	}

	void posix_route_manager::register_route(const route_type& route_entry)
	{
		const auto _gateway = gateway(route_entry.route);
		const auto ina = network_address(route_entry.route);

		if (_gateway)
		{
				set_route(route_action::add, route_entry.interface, ina, *_gateway);
		}
		else
		{
				set_route(route_action::add, route_entry.interface, ina);
		}
	}

	void posix_route_manager::unregister_route(const route_type& route_entry)
	{
		const auto _gateway = gateway(route_entry.route);
		const auto ina = network_address(route_entry.route);

		if (_gateway)
		{
				set_route(route_action::remove, route_entry.interface, ina, *_gateway);
		}
		else
		{
				set_route(route_action::remove, route_entry.interface, ina);
		}
	}
}

max_stars_repo_name: krattai/AEBL

max_issues_repo_name: krattai/AEBL

max_forks_repo_name: krattai/AEBL