text:
// kate: replace-tabs off; indent-width 4; indent-mode normal
// vim: ts=4:sw=4:noexpandtab

#define EIGEN_DONT_VECTORIZE
#define EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
#include <Eigen/Eigen>
#include <Eigen/LU>
#include <iostream>
#include <fstream>
#include <cmath>
#include <limits>
#include <vector>
#include <map>

using namespace std;

static double inlierRatio(0.8);
static int restartCount(1);
static int generationCount(64);

double uniformRand()
{
	return double(rand())/RAND_MAX;
}

double gaussianRand(double mean, double sigm)
{
	// Generation using the Polar (Box-Mueller) method.
	// Code inspired by GSL, which is a really great math lib.
	// http://sources.redhat.com/gsl/
	// C++ wrapper available.
	// http://gslwrap.sourceforge.net/
	double r, x, y;

	// Generate random number in unity circle.
	do
	{
		x = uniformRand()*2 - 1;
		y = uniformRand()*2 - 1;
		r = x*x + y*y;
	}
	while (r > 1.0 || r == 0);

	// Box-Muller transform.
	return sigm * y * sqrt (-2.0 * log(r) / r) + mean;
}

struct TrainingEntry
{
	double timeStamp;
	Eigen::Vector3d odom_tr;
	Eigen::Quaterniond odom_rot;
	Eigen::Vector3d icp_tr;
	Eigen::Quaterniond icp_rot;

	TrainingEntry(){}
	TrainingEntry(std::istream& is)
	{
		is >> timeStamp;
		double t_x = 0, t_y = 0, t_z = 0, q_x = 0, q_y = 0, q_z = 0, q_w = 1;
		is >> t_x;
		is >> t_y;
		is >> t_z;
		icp_tr = Eigen::Vector3d(t_x, t_y, t_z);
		is >> q_x;
		is >> q_y;
		is >> q_z;
		is >> q_w;
		icp_rot = Eigen::Quaterniond(q_w, q_x, q_y, q_z).normalized();
		is >> t_x;
		is >> t_y;
		is >> t_z;
		odom_tr = Eigen::Vector3d(t_x, t_y, t_z);
		is >> q_x;
		is >> q_y;
		is >> q_z;
		is >> q_w;
		odom_rot = Eigen::Quaterniond(q_w, q_x, q_y, q_z).normalized();
		//odom_tr = odom_rot*odom_tr;
		//cerr << icp_rot.x() << " " << icp_rot.y() << " " << icp_rot.z() << " " << icp_rot.w() <<endl;
		//cerr << icp_tr.x() << " " << icp_tr.y() << " " << icp_tr.z() << endl;
		// FIXME: bug in Eigen ?
		//icp_tr = icp_rot*icp_tr;
		//cerr << icp_tr.x() << " " << icp_tr.y() << " " << icp_tr.z() << endl;
	}

	void dump(ostream& stream) const
	{
		stream <<
			timeStamp << " : " <<
			icp_tr.x() << " " << icp_tr.y() << " " << icp_tr.z() << " " <<
			icp_rot.x() << " " << icp_rot.y() << " " << icp_rot.z() << " " << icp_rot.w() << " " <<
			odom_tr.x() << " " << odom_tr.y() << " " << odom_tr.z() << " " <<
			odom_rot.x() << " " << odom_rot.y() << " " << odom_rot.z() << " " << odom_rot.w();
	}
};

struct TrainingSet: public std::vector<TrainingEntry>
{
	void dump()
	{
		for (TrainingSet::const_iterator it(begin()); it != end(); ++it)
		{
			const TrainingEntry& entry(*it);
			entry.dump(cout);
			cout << "\n";
		}
	}
};

TrainingSet trainingSet;

struct Params
{
	Eigen::Vector3d tr;
	Eigen::Quaterniond rot;

	Params():tr(0,0,0),rot(1,0,0,0) {}
	Params(const Eigen::Vector3d& _tr, const Eigen::Quaterniond& _rot):
		tr(_tr),
		rot(_rot)
		{}
	Params(const double transVariance):
		tr(
			gaussianRand(0, transVariance),
			gaussianRand(0, transVariance),
			gaussianRand(0, transVariance)
		),
		rot(
			Eigen::Quaterniond(1,0,0,0) *
			Eigen::Quaterniond(Eigen::AngleAxisd(uniformRand() * M_PI*2, Eigen::Vector3d::UnitX())) *
			Eigen::Quaterniond(Eigen::AngleAxisd(uniformRand() * M_PI*2, Eigen::Vector3d::UnitY())) *
			Eigen::Quaterniond(Eigen::AngleAxisd(uniformRand() * M_PI*2, Eigen::Vector3d::UnitZ()))
		)
		{}

	// add random noise
	const Params& mutate(const double amount = 1.0)
	{
		const double count = fabs(gaussianRand(0, 1));
		for (int i = 0; i < int(count + 1); i++)
		{
			const int toMutate = rand() % 6;
			switch (toMutate)
			{
				case 0: tr.x() += gaussianRand(0, 0.1) * amount; break;
				case 1: tr.y() += gaussianRand(0, 0.1) * amount; break;
				case 2: tr.z() += gaussianRand(0, 0.1) * amount; break;
				case 3: rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitX())); break;
				case 4: rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitY())); break;
				case 5: rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitZ())); break;
				default: break;
			};
		}
		/*
		tr.x() += gaussianRand(0, 0.1) * amount;
		tr.y() += gaussianRand(0, 0.1) * amount;
		tr.z() += gaussianRand(0, 0.1) * amount;
		rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitX()));
		rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitY()));
		rot *= Eigen::Quaterniond(Eigen::AngleAxisd(gaussianRand(0, M_PI / 8) * amount, Eigen::Vector3d::UnitZ()));
		*/
		return *this;
	}

	// renormalize quaternion
	void normalize()
	{
		rot.normalize();
		if (rot.x() < 0)
		{
			rot.x() = -rot.x();
			rot.y() = -rot.y();
			rot.z() = -rot.z();
			rot.w() = -rot.w();
		}
	}

	// dump content of this transform
	void dumpTransform(ostream& stream) const
	{
		stream << tr.x() << " " << tr.y() << " " << tr.z() << " " <<
				rot.x() << " " << rot.y() << " " << rot.z() << " " << rot.w();
	}

	// dump what to copy/paste
	void dumpCopyPaste(ostream& stream) const
	{
		/*stream <<
			"translation: x=" << tr.x() << " y=" << tr.y() << " z=" << tr.z() << " " <<
			"quaternion: x=" << rot.x() << " y=" << rot.y() << " z=" << rot.z() << " w=" << rot.w();
		*/
		stream << "<node pkg=\"tf\" type=\"static_transform_publisher\" name=\"base_link_to_kinect\" args=\"";
		dumpTransform(stream);
		stream << " /base_link /kinect 100\"/>"
		" \n\n OR \n\n"
		"rosrun tf static_transform_publisher ";
		dumpTransform(stream);
		stream << " /base_link /kinect 100\n";
	}
};

// we have two typedefs for the same physical object because they are conceptually different
typedef vector<Params> Genome;
typedef vector<Params> ParamsVector;

void normalizeParams(ParamsVector& params)
{
	for (size_t i = 0; i < params.size(); ++i)
		params[i].normalize();
}

void dumpParamsStats(ostream& stream, const ParamsVector& params)
{
	Eigen::Vector3d trMean(0,0,0);
	Eigen::Quaterniond::Coefficients rotMean(0,0,0,0);
	for (size_t i = 0; i < params.size(); ++i)
	{
		trMean += params[i].tr;
		rotMean += params[i].rot.coeffs();
	}
	trMean /= double(params.size());
	rotMean /= double(params.size());

	Params mean(trMean,Eigen::Quaterniond(rotMean));
	mean.normalize();
	mean.dumpCopyPaste(stream); stream << "\n";

	Eigen::Vector3d trVar(0,0,0);
	Eigen::Quaterniond::Coefficients rotVar(0,0,0,0);
	for (size_t i = 0; i < params.size(); ++i)
	{
		trVar += (params[i].tr - trMean).cwiseProduct(params[i].tr - trMean);
		rotVar += (params[i].rot.coeffs() - rotMean).cwiseProduct(params[i].rot.coeffs() - rotMean);
	}
	trVar /= double(params.size());
	rotVar /= double(params.size());

	stream << "\nVariance:\n";
	Params var(trVar,Eigen::Quaterniond(rotVar));
	var.dumpTransform(stream); stream << "\n";

	stream << "\nValues:\n";
	for (size_t i = 0; i < params.size(); ++i)
	{
		params[i].dumpTransform(stream); stream << "\n";
	}
	stream << endl;
}

double computeError(const Params& p, const TrainingEntry& e)
{
	/*
	// version with Eigen::Matrix4d
	Eigen::Matrix4d blk(Eigen::Matrix4d::Identity());
	blk.corner(Eigen::TopLeft,3,3) = p.rot.toRotationMatrix();
	blk.corner(Eigen::TopRight,3,1) = p.tr;

	Eigen::Matrix4d odom(Eigen::Matrix4d::Identity());
	odom.corner(Eigen::TopLeft,3,3) = e.odom_rot.toRotationMatrix();
	odom.corner(Eigen::TopRight,3,1) = e.odom_tr;

	Eigen::Matrix4d blk_i(blk.inverse());

	//const Eigen::Matrix4d pred_icp = blk * odom * blk_i;
	const Eigen::Matrix4d pred_icp = blk_i * odom * blk;

	const Eigen::Matrix3d pred_icp_rot_m = pred_icp.corner(Eigen::TopLeft,3,3);
	const Eigen::Quaterniond pred_icp_rot = Eigen::Quaterniond(pred_icp_rot_m);
	const Eigen::Vector3d pred_icp_tr = pred_icp.corner(Eigen::TopRight,3,1);
	*/

	// version with Eigen::Transform3d
	typedef Eigen::Transform<double, 3, Eigen::Affine> Transform3d;
	typedef Eigen::Translation<double, 3> Translation3d;

	const Transform3d blk = Translation3d(p.tr) * p.rot;
	const Transform3d blk_i = Transform3d(blk.inverse(Eigen::Isometry));
	const Transform3d odom = Translation3d(e.odom_tr) * e.odom_rot;
	//const Eigen::Transform3d pred_icp = blk * odom * blk_i;
	const Transform3d pred_icp = blk_i * odom * blk;

	const Eigen::Matrix3d pred_icp_rot_m = pred_icp.matrix().topLeftCorner(3,3);
	const Eigen::Quaterniond pred_icp_rot = Eigen::Quaterniond(pred_icp_rot_m);
	const Eigen::Vector3d pred_icp_tr = pred_icp.translation();

	// identity checked ok
	//cerr << "mat:\n" << blk.matrix() * blk_i.matrix() << endl;

	//cout << "dist: " << (e.icp_tr - pred_icp.translation()).norm() << endl;
	//cout << "ang dist: " << e.icp_rot.angularDistance(pred_icp_rot) << endl;
	//cout << "tr pred icp:\n" << pred_icp.translation() << "\nicp:\n" << e.icp_tr << "\n" << endl;
	//cout << "rot pred icp:\n" << pred_icp_rot << "\nicp:\n" << e.icp_rot << "\n" << endl;
	// FIXME: tune coefficient for rot vs trans

	const double e_tr((e.icp_tr - pred_icp_tr).norm());
	const double e_rot(e.icp_rot.angularDistance(pred_icp_rot));
	/*if (e_tr < 0)
		abort();
	if (e_rot < 0)
		abort();*/
	//cerr << e_tr << " " << e_rot << endl;
	return e_tr + e_rot;
}

// compute errors with outlier rejection
double computeError(const Params& p)
{
	vector<double> errors;
	errors.reserve(trainingSet.size());
	for (TrainingSet::const_iterator it(trainingSet.begin()); it != trainingSet.end(); ++it)
	{
		const TrainingEntry& entry(*it);
		errors.push_back(computeError(p, entry));
	}
	sort(errors.begin(), errors.end());
	const size_t inlierCount(errors.size() * inlierRatio);
	double error = 0;
	for (size_t i = 0; i < inlierCount; ++i)
		error += errors[i];
	return error;
}

double evolveOneGen(Genome& genome, double annealing = 1.0, Params* bestParams = 0)
{
	typedef multimap<double, Params> EvaluationMap;
	typedef EvaluationMap::iterator EvaluationMapIterator;
	EvaluationMap evalutationMap;

	double totalError = 0;
	double bestError = numeric_limits<double>::max();
	int bestInd = 0;
	for (size_t ind = 0; ind < genome.size(); ind++)
	{
		const double error = computeError(genome[ind]);
		if (error < bestError)
		{
			bestError = error;
			bestInd = ind;
		}

		totalError += error;
		evalutationMap.insert(make_pair(error, genome[ind]));

		/*cout << "E " << ind << " : ";
		genome[ind].dump(cout);
		cout << " = " << error << "\n";*/
	}

	if (bestParams)
		*bestParams = genome[bestInd];

	assert((genome.size() / 4) * 4 == genome.size());

	size_t ind = 0;
	for (EvaluationMapIterator it = evalutationMap.begin(); ind < genome.size() / 4; ++it, ++ind )
	{
		//cout << "S " << it->first << "\n";
		genome[ind * 4] = it->second;
		genome[ind * 4 + 1] = it->second.mutate(annealing);
		genome[ind * 4 + 2] = it->second.mutate(annealing);
		genome[ind * 4 + 3] = it->second.mutate(annealing);
	}

	return bestError;
}

int main(int argc, char** argv)
{
	srand(time(0));

	if (argc < 2)
	{
		cerr << "Usage " << argv[0] << " LOG_FILE_NAME [INLIER_RATIO] [RESTART_COUNT] [GEN_COUNT]\n";
		cerr << "  LOG_FILE_NAME   name of file to load containing delta tf\n";
		cerr << "  INLIER_RATIO    ratio of inlier to use for error computation (range: ]0:1], default: 0.8)\n";
		cerr << "  RESTART_COUNT   number of random restart (default: 1)\n";
		cerr << "  GEN_COUNT       number of generation for the ES (default: 64)\n";
		cerr << endl;
		return 1;
	}
	if (argc >= 3)
	{
		inlierRatio = atof(argv[2]);
		if (inlierRatio <= 0 || inlierRatio > 1)
		{
			cerr << "Invalid inline ratio: " << inlierRatio << ", must be within range ]0:1]" << endl;
			return 2;
		}
	}
	cout << "Inlier ratio: " << inlierRatio << endl;
	if (argc >= 4)
	{
		restartCount = atoi(argv[3]);
		if (restartCount < 1)
		{
			cerr << "Invalid restart count: " << restartCount << ", must be greater than 0" << endl;
			return 3;
		}
	}
	cout << "Restart count: " << restartCount << endl;
	if (argc >= 5)
	{
		generationCount = atoi(argv[4]);
		if (generationCount < 1)
		{
			cerr << "Invalid generation count: " << generationCount << ", must be greater than 0" << endl;
			return 4;
		}
	}
	cout << "Generation count: " << generationCount << endl;

	ifstream ifs(argv[1]);
	while (ifs.good())
	{
		TrainingEntry e(ifs);
		if (ifs.good())
		{
			trainingSet.push_back(e);
			//e.dump(cerr); cerr << endl;
		}
	}
	cout << "Loaded " << trainingSet.size() << " training entries" << endl;
	//trainingSet.dump();

	ParamsVector bests;
	for (int i = 0; i < restartCount; ++i)
	{
		cout << "Starting " << i << " restart on " << restartCount << endl;
		Genome genome(1024);
		for (size_t g = 0; g < genome.size(); ++g)
			genome[g] = Params(0.5);
		for (int g = 0; g < generationCount; ++g)
		{
			cout << "\r" << g << " best has error " << evolveOneGen(genome, 2. * (double)(generationCount - g) / (double)(generationCount)) << "            ";
			cout.flush();
		}
		Params best;
		cout << "\r  Best error of restart " << i << ": " << evolveOneGen(genome, 1.0, &best) << endl;
		bests.push_back(best);
	}

	cout << "\nOptimization completed, code to COPY-PASTE in to use the transformation:\n\n";
	if (restartCount == 1)
	{
		bests[0].dumpCopyPaste(cout);
	}
	else
	{
		normalizeParams(bests);
		dumpParamsStats(cout, bests);
	}

	return 0;
}

max_stars_repo_name: KRudraksh/Relocalisation

max_issues_repo_name: KRudraksh/Relocalisation

max_forks_repo_name: KRudraksh/Relocalisation


text:
#include <Eigen/Core>
#include <Eigen/Dense>

namespace tyco
{
namespace P3
{

template<typename T, typename Cs>
struct row_vector
{
    Eigen::Matrix<T, 1, 4> raw_;
};

template<typename T, typename Cs>
struct column_vector
{
    Eigen::Matrix<T, 4, 1> raw_;
};

template<typename T, typename CsLeft, typename CsRight>
struct homography
{
    Eigen::Matrix<T, 4, 4> raw_;
};

// M^-1
template<typename T, typename CsLeftIn, typename CsRightIn>
homography<T, CsRightIn, CsLeftIn> inverse(const homography<T, CsLeftIn, CsRightIn>& M)
{
    return {M.raw_.inverse()};
}

// M * p
template<typename T, typename CsLeftIn, typename CsRightIn>
column_vector<T, CsLeftIn> operator*(
    const homography<T, CsLeftIn, CsRightIn>& M,
    const column_vector<T, CsRightIn>& p)
{
    return {M.raw_ * p.raw_};
}

// p^T * M
template<typename T, typename CsLeftIn, typename CsRightIn>
row_vector<T, CsRightIn> operator*(
    const row_vector<T, CsLeftIn>& p,
    const homography<T, CsLeftIn, CsRightIn>& M)
{
    return {p.raw_ * M.raw_};
}

// M1 * M2
template<typename T, typename CsLeft, typename CsMiddle, typename CsRight>
const homography<T, CsLeft, CsRight> operator*(
    const homography<T, CsLeft, CsMiddle>& M1,
    const homography<T, CsMiddle, CsRight>& M2)
{
    return {M1.raw_ * M2.raw_};
}

} // namespace P3
} // namespace tyco

ext: hpp

max_stars_repo_name: mabur/tyco

max_issues_repo_name: mabur/tyco

max_forks_repo_name: mabur/tyco


text:
#pragma once

#include <vector>
#include <boost/graph/graph_concepts.hpp>
#include <caffe/ultinous/ImageClassificationModel.h>
#include <caffe/ultinous/FeatureMap.hpp>
#include <caffe/ultinous/AbstractTripletGenerator.hpp>
#include <caffe/ultinous/FeatureCollectorTripletGenerator.hpp>
#include "caffe/proto/caffe.pb.h"

namespace caffe {
namespace ultinous {

template <typename Dtype>
class HardTripletGenerator : public AbstractTripletGenerator
{
public:
  HardTripletGenerator(HardTripletParameter const &htp, BasicModel const &basicModel)
    : m_sampler(basicModel)
    , m_classesInSample(htp.sampledclasses())
    , m_imagesInSampleClass(htp.sampledpictures())
    , m_margin(htp.margin())
    , m_indexInSample(m_classesInSample*m_imagesInSampleClass)
    , m_featureMap(FeatureMapContainer<Dtype>::instance(htp.featuremapid()))
    , m_tooHardTriplets(htp.toohardtriplets())
    , m_hardestPositive(htp.hardestpositive())
    , m_hardestNegative(htp.hardestnegative())
    , m_isLastTripletHard(false)
  {
    CHECK_GT( m_classesInSample, 0 );
    CHECK_GT( m_imagesInSampleClass, 0 );
    ImageSampler::initSample( m_classesInSample, m_imagesInSampleClass, m_sample );

    for( size_t i = 0; i < m_classesInSample*m_imagesInSampleClass; ++i)
      m_shuffle.push_back(i);

    m_numImagesInModel = 0;
    for( int i = 0; i < basicModel.size(); ++i )
      m_numImagesInModel += basicModel[i].images.size();

    m_featureMap.resize( m_numImagesInModel );
    FeatureCollectorTripletGenerator<Dtype>::init( basicModel );
  }
private:
  typedef size_t ImageIndex;
  typedef size_t ClassIndex;
  typedef size_t SampleIndex;
  typedef ImageClassificationModel::ImageIndexes ImageIndexes;
public:

  Triplet nextTriplet()
  {
    static bool featuresCollected = false;

    if( !featuresCollected )
    {
      if( m_featureMap.numFeatures() != m_numImagesInModel )
      {
        m_isLastTripletHard = true;
        return FeatureCollectorTripletGenerator<Dtype>::getInstance().nextTriplet();
      }

      LOG(INFO) << "All features are collected!";
      featuresCollected = true;
    }

    if(classIndex(m_indexInSample) >= m_classesInSample)
      resample();

    Triplet t;
    m_isLastTripletHard = false;

    SampleIndex anchorIndex = m_shuffle[m_indexInSample];

    t.push_back(image(anchorIndex)); // anchor

    const Vec& dvec = m_distances[anchorIndex];
    SampleIndex posSampleBegin = classIndex(anchorIndex)*m_imagesInSampleClass;
    SampleIndex posSampleEnd = posSampleBegin + m_imagesInSampleClass;
    SampleIndex posIndex = 0;

    if( m_hardestPositive )
    {
      Dtype maxPosDistance = 0;
      for(size_t posSample = posSampleBegin; posSample<posSampleEnd; ++posSample)
      {
        if(posSample==anchorIndex)
          continue;
        if(dvec[posSample] > maxPosDistance)
        {
          maxPosDistance = dvec[posSample];
          posIndex = posSample;
        }
      }
    }
    else
    {
      posIndex = 1+anchorIndex;
      if( (posIndex % m_imagesInSampleClass) == 0 )
        posIndex = posSampleBegin;
    }

    t.push_back(image(posIndex)); // hard positive

    Dtype posDistance = dvec[posIndex];
    Dtype posDistanceWithMargin = posDistance+m_margin;

    size_t nSample = m_classesInSample*m_imagesInSampleClass;
    Dtype closeNegDistance = std::numeric_limits<Dtype>::max();
    size_t negIndex = std::numeric_limits<size_t>::max();

    std::vector<size_t> hardNegIndices;

    for(size_t negSample = 0; negSample<nSample; ++negSample)
    {
      if( negSample == posSampleBegin )
      {
        negSample = posSampleEnd-1;
        continue;
      }

      if(negIndex == std::numeric_limits<size_t>::max())
        negIndex = negSample; // the first negative selected as "random"

      if( dvec[negSample] >= posDistance && dvec[negSample] < posDistanceWithMargin )
      {
          m_isLastTripletHard = true;
          if( m_hardestNegative )
          {
            if( dvec[negSample] < closeNegDistance )
            {
              closeNegDistance = dvec[negSample];
              negIndex = negSample;
            }
          }
          else
          {
            hardNegIndices.push_back( negSample );
          }
      }
    }

    if( m_tooHardTriplets && !m_isLastTripletHard )
    {
      closeNegDistance = 0;
      for(size_t negSample = 0; negSample<nSample; ++negSample)
      {
        if( negSample == posSampleBegin )
        {
          negSample = posSampleEnd-1;
          continue;
        }

        if( dvec[negSample] < posDistance )
        {
          m_isLastTripletHard = true;
          if( m_hardestNegative )
          {
            if( dvec[negSample] > closeNegDistance )
            {
              closeNegDistance = dvec[negSample];
              negIndex = negSample;
            }
          }
          else
          {
            hardNegIndices.push_back( negSample );
          }
        }
      }
    }

    if( !m_hardestNegative )
    {
      if( hardNegIndices.size() > 0 )
      {
        negIndex = hardNegIndices[rand()%hardNegIndices.size()];
        m_isLastTripletHard = true;
      }
      else
      {
        negIndex = rand()%(nSample-m_imagesInSampleClass);
        if( negIndex >= posSampleBegin )
          negIndex += m_imagesInSampleClass;
        m_isLastTripletHard = false;
      }
    }

    t.push_back(image(negIndex)); // hard negative

    ++m_indexInSample;

    return t;
  }

  bool isLastTripletHard( ) const
  {
    return m_isLastTripletHard;
  }

  const FeatureMap<Dtype>& getFeatureMap( ) const
  {
    return m_featureMap;
  }

  Dtype getMargin( ) const
  {
    return m_margin;
  }
private:
  ClassIndex classIndex(SampleIndex idx) const { return idx/m_imagesInSampleClass; }
  ImageIndex imageIndex(SampleIndex idx) const { return idx%m_imagesInSampleClass; }
  const ImageIndexes& images(SampleIndex idx) const { return m_sample[classIndex(idx)].images; }
  ImageIndex image(SampleIndex idx) const { return images(idx)[imageIndex(idx)]; }
private:
  void resample()
  {
    m_sampler.sample(m_sample);
    recalcDistances();
    shuffle( m_shuffle.begin(), m_shuffle.end() );
    m_indexInSample = 0;
  }

  void recalcDistancesGPU(); // src/caffe/ultinous/HardTripletGenerator.cu

  void recalcDistancesCPU()
  {
    size_t const nSample = m_classesInSample * m_imagesInSampleClass;

    typename FeatureMap<Dtype>::FeatureVec sqr;

    for( size_t i = 0; i < nSample-1; i++ )
    {
      const typename FeatureMap<Dtype>::FeatureVec& feat1 = m_featureMap.getFeatureVec( image(i) );
      CHECK_GT( feat1.size(), 0);

      sqr.resize( feat1.size() );

      for( size_t j = i+1; j < nSample; j++ )
      {
        const typename FeatureMap<Dtype>::FeatureVec& feat2 = m_featureMap.getFeatureVec( image(j) );

        CHECK_GT( feat2.size(), 0 );
        CHECK_EQ( feat1.size(), feat2.size() );

        caffe_sub( feat1.size(), &(feat1[0]), &(feat2[0]), &(sqr[0]) );
        Dtype dist = caffe_cpu_dot( sqr.size(), &(sqr[0]), &(sqr[0]) );

        m_distances[i][j] = dist;
        m_distances[j][i] = dist;
      }
    }
  }
  void recalcDistances()
  {
    size_t const nSample = m_classesInSample * m_imagesInSampleClass;

    CHECK_GT( nSample, 0 );

    if( m_distances.size() != nSample )
      m_distances = Mat( nSample, Vec(nSample, 0) );

    CHECK_EQ( m_distances.size(), nSample );
    for( size_t i = 0; i < nSample; i++ )
      CHECK_EQ( m_distances[i].size(), nSample );

    #ifdef CPU_ONLY
    recalcDistancesCPU();
    #else
    recalcDistancesGPU();
    #endif
  }

private:
  typedef std::vector<Dtype> Vec;
  typedef std::vector<Vec> Mat;
  typedef ImageSampler::Sample Sample;
private:
  ImageSampler m_sampler;
  size_t m_classesInSample;
  size_t m_imagesInSampleClass;
  Mat m_distances;
  Dtype m_margin;
  Sample m_sample;
  SampleIndex m_indexInSample;
  FeatureMap<Dtype>& m_featureMap;
  bool m_tooHardTriplets;
  bool m_hardestPositive;
  bool m_hardestNegative;
  bool m_isLastTripletHard;

  std::vector<SampleIndex> m_shuffle;

  shared_ptr<SyncedMemory> m_syncedFeatures;
  shared_ptr<SyncedMemory> m_syncedDistances;

  size_t m_numImagesInModel;
};

} // namespace ultinous
} // namespace caffe

max_stars_repo_name: Ultinous/caffe

max_issues_repo_name: Ultinous/caffe

max_forks_repo_name: Ultinous/caffe


text:
extern "C" {
#include <ruby.h>
#include <ruby/intern.h>
}
#include <ruby/encoding.h>
#include <stdio.h>
#include <thrift/protocol/TBinaryProtocol.h>
#include <thrift/protocol/TCompactProtocol.h>
#include <boost/make_shared.hpp>
#include <functional>
#include <map>
#include <vector>
#include "serializer.h"

using namespace std;
using ::apache::thrift::protocol::TType;
using ::apache::thrift::transport::TMemoryBuffer;
using namespace ::apache::thrift;

VALUE sym_for_class;
VALUE sym_for_binary;
VALUE sym_for_type;
VALUE sym_for_element;
VALUE sym_for_key;
VALUE sym_for_value;
VALUE sym_for_optional;
VALUE sym_for_name;
VALUE intern_for_FIELDS;
VALUE intern_for_new;
VALUE intern_for_keys;
VALUE intern_for_values;
ID intern_for_setfield_ivar;
VALUE intern_for_to_a;
VALUE klass_for_union;
VALUE klass_for_integer;
VALUE klass_for_float;
VALUE klass_for_string;
VALUE klass_for_hash;
VALUE klass_for_set;
VALUE klass_for_array;

extern VALUE Sparsam;
extern VALUE SparsamNativeError;
VALUE SparsamTypeMismatchError;
VALUE SparsamMissingMandatory;
VALUE SparsamUnionException;
VALUE SparsamUnknownTypeException;

VALUE SetClass;

KlassFieldsCache klassCache;  // consider the memory leaked.
std::unordered_set<VALUE> unions;

void *serializer_create() { return (void *)(new ThriftSerializer()); }

void serializer_free(void *data) {
  ThriftSerializer *ts = (ThriftSerializer *)(data);
  delete ts;
}

void initialize_constants() {
  sym_for_class = ID2SYM(rb_intern("class"));
  sym_for_binary = ID2SYM(rb_intern("binary"));
  sym_for_type = ID2SYM(rb_intern("type"));
  sym_for_element = ID2SYM(rb_intern("element"));
  sym_for_key = ID2SYM(rb_intern("key"));
  sym_for_value = ID2SYM(rb_intern("value"));
  sym_for_optional = ID2SYM(rb_intern("optional"));
  sym_for_name = ID2SYM(rb_intern("name"));
  intern_for_FIELDS = rb_intern("FIELDS");
  intern_for_new = rb_intern("new");
  intern_for_keys = rb_intern("keys");
  intern_for_values = rb_intern("values");
  intern_for_to_a = rb_intern("to_a");
  intern_for_setfield_ivar = rb_intern("@setfield");
  klass_for_set = rb_const_get_at(rb_cObject, rb_intern("Set"));
  klass_for_integer = rb_const_get_at(rb_cObject, rb_intern("Integer"));
  klass_for_float = rb_const_get_at(rb_cObject, rb_intern("Float"));
  klass_for_string = rb_const_get_at(rb_cObject, rb_intern("String"));
  klass_for_hash = rb_const_get_at(rb_cObject, rb_intern("Hash"));
  klass_for_array = rb_const_get_at(rb_cObject, rb_intern("Array"));
}

void initialize_runtime_constants() {
  klass_for_union = rb_const_get_at(Sparsam, rb_intern("Union"));
  SparsamMissingMandatory =
      rb_const_get_at(Sparsam, rb_intern("MissingMandatory"));
  SparsamTypeMismatchError =
      rb_const_get_at(Sparsam, rb_intern("TypeMismatch"));
  SparsamUnionException = rb_const_get_at(Sparsam, rb_intern("UnionException"));
  SparsamUnknownTypeException =
      rb_const_get_at(Sparsam, rb_intern("UnknownTypeException"));
  SetClass = rb_const_get_at(rb_cObject, rb_intern("Set"));
}

void serializer_init(void *serializer, int protocol, void *str_arg1,
                     uint32_t len) {
  using ::apache::thrift::protocol::TBinaryProtocol;
  using ::apache::thrift::protocol::TCompactProtocol;
  using ::apache::thrift::protocol::TProtocol;
  using ::boost::shared_ptr;

  ThriftSerializer *ts = (ThriftSerializer *)(serializer);
  shared_ptr<TMemoryBuffer> strBuffer;
  if (str_arg1 != NULL) {
    strBuffer = boost::make_shared<TMemoryBuffer>(
        (uint8_t *)str_arg1, len, TMemoryBuffer::TAKE_OWNERSHIP);
  } else {
    strBuffer = boost::make_shared<TMemoryBuffer>();
  }
  Proto proto = static_cast<Proto>(protocol);
  if (proto == compact) {
    ts->tprot = shared_ptr<TProtocol>(new TCompactProtocol(strBuffer));
  } else if (proto == binary) {
    ts->tprot = shared_ptr<TProtocol>(new TBinaryProtocol(strBuffer));
  } else {
    rb_raise(SparsamNativeError, "Unknown protocol %d", proto);
  }
  ts->tmb = strBuffer;
}

#define get_ts()                          \
  void *self_data = NULL;                 \
  Data_Get_Struct(self, void, self_data); \
  ThriftSerializer *ts = (ThriftSerializer *)(self_data);

#define watch_for_texcept() try {
#define catch_thrift_and_reraise()                \
  }                                               \
  catch (::apache::thrift::TException e) {        \
    rb_raise(SparsamNativeError, "%s", e.what()); \
    return Qnil;                                  \
  }

static inline VALUE make_ruby_string(const string &val) {
  return rb_enc_str_new(val.c_str(), val.size(), rb_utf8_encoding());
}

static inline VALUE make_ruby_binary(const string &val) {
  return rb_str_new(val.c_str(), val.size());
}

static void raise_exc_with_struct_and_field_names(VALUE exc_class,
                                                  VALUE msg_prefix,
                                                  VALUE outer_struct_class,
                                                  VALUE field_sym) {
  VALUE struct_name = rb_class_name(outer_struct_class);
  VALUE field_name = rb_sym_to_s(field_sym);

  VALUE msg = rb_sprintf("%s (in %s#%s)", RSTRING_PTR(msg_prefix),
                         RSTRING_PTR(struct_name), RSTRING_PTR(field_name));
  VALUE args[3] = {msg, struct_name, field_name};
  VALUE e = rb_class_new_instance(3, args, exc_class);
  rb_exc_raise(e);
}

static void raise_type_mismatch(VALUE outer_struct, VALUE field_sym,
                                const char *expected, VALUE actual) {
  VALUE actual_name = rb_class_name(CLASS_OF(actual));
  VALUE msg = rb_sprintf(
      "Mismatched type (expected to be compatible with: %s, found: %s)",
      expected, RSTRING_PTR(actual_name));

  raise_exc_with_struct_and_field_names(SparsamTypeMismatchError, msg,
                                        CLASS_OF(outer_struct), field_sym);
}

static void raise_type_mismatch(VALUE outer_struct, VALUE field_sym, int ttype,
                                VALUE actual) {
  raise_type_mismatch(outer_struct, field_sym, TTypeName((size_t)ttype).c_str(),
                      actual);
}

static void raise_type_mismatch(VALUE outer_struct, VALUE field_sym,
                                VALUE klass, VALUE actual) {
  VALUE expected_name = rb_class_name(klass);
  raise_type_mismatch(outer_struct, field_sym, RSTRING_PTR(expected_name),
                      actual);
}

static inline long raise_type_mismatch_as_value(VALUE outer_struct,
                                                VALUE field_sym, int ttype,
                                                VALUE actual) {
  raise_type_mismatch(outer_struct, field_sym, ttype, actual);
  return 0;
}

static inline void Sparsam_Check_Type(VALUE x, int t, VALUE outer_struct,
                                      VALUE field_sym) {
  if (!(RB_TYPE_P(x, t))) {
    raise_type_mismatch(outer_struct, field_sym, t, x);
  }
}

static inline VALUE make_ruby_bool(bool val) { return val ? Qtrue : Qfalse; }

void ThriftSerializer::skip_n_type(uint32_t n, TType ttype) {
  for (uint32_t i = 0; i < n; i++) {
    this->tprot->skip(ttype);
  }
}

void ThriftSerializer::skip_n_pair(uint32_t n, TType type_a, TType type_b) {
  for (uint32_t i = 0; i < n; i++) {
    this->tprot->skip(type_a);
    this->tprot->skip(type_b);
  }
}

// Blatantly copied protobuf's design
// https://git.io/vHuUn
// CONVERT is new here, because we're targeting ruby
#define HANDLE_TYPE(TYPE, CPPTYPE, READ_METHOD, CONVERT) \
  case protocol::T_##TYPE: {                             \
    CPPTYPE value;                                       \
    this->tprot->read##READ_METHOD(value);               \
    ret = CONVERT(value);                                \
    break;                                               \
  }

VALUE ThriftSerializer::readAny(TType ttype, FieldInfo *field_info) {
  VALUE ret = Qnil;
  switch (ttype) {
    // Handle all the non-container types by marco
    HANDLE_TYPE(I16, int16_t, I16, INT2FIX)
    HANDLE_TYPE(I32, int32_t, I32, INT2FIX)
    HANDLE_TYPE(I64, int64_t, I64, LL2NUM)
    HANDLE_TYPE(BOOL, bool, Bool, make_ruby_bool)
    HANDLE_TYPE(DOUBLE, double, Double, DBL2NUM)
    HANDLE_TYPE(BYTE, int8_t, Byte, INT2FIX)

    case protocol::T_STRING: {
      string value;
      if (field_info->isBinaryString) {  // if (field_info[:binary])
        this->tprot->readBinary(value);
        ret = make_ruby_binary(value);
      } else {
        this->tprot->readString(value);
        ret = make_ruby_string(value);
      }
      break;
    }

    case protocol::T_LIST: {
      TType element_type;
      uint32_t size;

      this->tprot->readListBegin(element_type, size);
      if (field_info->elementType == NULL ||
          element_type != field_info->elementType->ftype) {
        this->skip_n_type(size, element_type);
        break;
      }
      ret = rb_ary_new2(size);

      for (uint32_t i = 0; i < size; i++) {
        rb_ary_store(ret, i,
                     this->readAny(element_type, field_info->elementType));
      }
      this->tprot->readListEnd();

      break;
    }

    case protocol::T_SET: {
      TType element_type;
      uint32_t size;

      this->tprot->readSetBegin(element_type, size);
      if (field_info->elementType == NULL ||
          element_type != field_info->elementType->ftype) {
        this->skip_n_type(size, element_type);
        break;
      }
      VALUE ary = rb_ary_new2(size);

      for (uint32_t i = 0; i < size; i++) {
        rb_ary_store(ary, i,
                     this->readAny(element_type, field_info->elementType));
      }
      ret = rb_class_new_instance(1, &ary, klass_for_set);
      this->tprot->readSetEnd();
      break;
    }

    case protocol::T_STRUCT: {
      string cname;
      this->tprot->readStructBegin(cname);
      if (unions.count(field_info->klass) == 1) {
        ret = this->readUnion(field_info->klass);
      } else {
        ret = this->readStruct(field_info->klass);
      }
      this->tprot->readStructEnd();
      break;
    }

    case protocol::T_MAP: {
      TType key_type, value_type;
      uint32_t size;
      VALUE k, v;

      this->tprot->readMapBegin(key_type, value_type, size);

      if (field_info->keyType == NULL ||
          field_info->elementType == NULL) {  // no type check to be consistent
        skip_n_pair(size, key_type, value_type);
        break;
      }

      ret = rb_hash_new();
      for (uint32_t i = 0; i < size; i++) {
        k = this->readAny(key_type, field_info->keyType);
        v = this->readAny(value_type, field_info->elementType);
        rb_hash_aset(ret, k, v);
      }
      this->tprot->readMapEnd();
      break;
    }

    default:
      this->tprot->skip(ttype);
      rb_raise(SparsamUnknownTypeException, "Received unknown type with id: %d",
               ttype);
      break;
  }

  return ret;
}

#undef HANDLE_TYPE

VALUE ThriftSerializer::readStruct(VALUE klass) {
  string cname;
  FieldBegin fieldBegin;
  TType typeId;
  FieldInfo *fieldInfo;
  VALUE ret = rb_class_new_instance(0, NULL, klass);  // ret = &klass.new
  auto fields = FindOrCreateFieldInfoMap(klass);

  while (true) {
    this->tprot->readFieldBegin(cname, fieldBegin.ftype, fieldBegin.fid);
    if (fieldBegin.ftype == protocol::T_STOP) {
      break;
    }
    auto iter = fields->find(fieldBegin.fid);

    if (iter == fields->end()) {
      this->tprot->skip(fieldBegin.ftype);
      this->tprot->readFieldEnd();
      continue;
    }

    fieldInfo = iter->second;

    typeId = fieldInfo->ftype;

    if (typeId != fieldBegin.ftype) {
      raise_exc_with_struct_and_field_names(
          SparsamTypeMismatchError,
          rb_sprintf("Mismatched type (definition: %s, found: %s)",
                     TTypeName(fieldBegin.ftype).c_str(),
                     TTypeName(typeId).c_str()),
          klass, fieldInfo->symName);
    }

    VALUE rb_value = this->readAny(fieldBegin.ftype, iter->second);
    if (!NIL_P(rb_value)) {
      rb_ivar_set(ret, fieldInfo->ivarName, rb_value);
    }

    this->tprot->readFieldEnd();
  }
  return ret;
}

VALUE ThriftSerializer::readUnion(VALUE klass) {
  string cname;
  FieldBegin fieldBegin;

  VALUE ret = rb_class_new_instance(0, NULL, klass);  // ret = &klass.new
  auto fields = FindOrCreateFieldInfoMap(klass);

  VALUE key, rb_value;

  this->tprot->readFieldBegin(cname, fieldBegin.ftype, fieldBegin.fid);
  if (fieldBegin.ftype == protocol::T_STOP) {
    return ret;
  }

  auto iter = fields->find(fieldBegin.fid);

  if (iter == fields->end()) {
    this->tprot->skip(fieldBegin.ftype);
    this->tprot->readFieldEnd();
    return ret;
  }

  rb_value = this->readAny(fieldBegin.ftype, iter->second);
  if (!NIL_P(rb_value)) {
    rb_ivar_set(ret, intern_for_setfield_ivar, iter->second->symName);
    rb_ivar_set(ret, iter->second->ivarName, rb_value);
  }

  this->tprot->readFieldEnd();
  this->tprot->readFieldBegin(cname, fieldBegin.ftype, fieldBegin.fid);

  if (fieldBegin.ftype != protocol::T_STOP) {
    rb_raise(SparsamUnionException, "More than one element in union.");
    return Qnil;
  }

  return ret;
}

// for the unary `+` before lambda:
// https://stackoverflow.com/a/18889029/4944625
// explicit cast to work with signature: (int (*)(...))
#define HASH_FOREACH_BEGIN(hash, ...) \
  void *_args[] = {__VA_ARGS__};      \
  rb_hash_foreach(hash, (int (*)(ANYARGS))(+[](VALUE k, VALUE v, VALUE args) { \
    void **argv = (void **) args;
#define HASH_FOREACH_RET() return (int)ST_CONTINUE;

#define HASH_FOREACH_ABORT() return (int)ST_STOP;

#define HASH_FOREACH_END() \
  HASH_FOREACH_RET()       \
  }), (VALUE) _args);

#define HANDLE_TYPE(TYPE, WRITE_METHOD, CONVERT) \
  case protocol::T_##TYPE: {                     \
    this->tprot->write##WRITE_METHOD(CONVERT);   \
    break;                                       \
  }

static inline long raise_bignum_range_error_as_value() {
  rb_raise(rb_eRangeError, "bignum too big to convert");
  return 0;
}

#define CONVERT_FIXNUM(CONVERT, INTENDED)                             \
  ((FIXNUM_P(actual))                                                 \
       ? CONVERT(actual)                                              \
       : ((RB_TYPE_P(actual, T_BIGNUM))                               \
              ? raise_bignum_range_error_as_value()                   \
              : raise_type_mismatch_as_value(outer_struct, field_sym, \
                                             INTENDED, actual)))

#define CONVERT_I64                                                   \
  ((FIXNUM_P(actual))                                                 \
       ? (LONG_LONG)FIX2LONG(actual)                                  \
       : ((RB_TYPE_P(actual, T_BIGNUM))                               \
              ? rb_big2ll(actual)                                     \
              : raise_type_mismatch_as_value(outer_struct, field_sym, \
                                             protocol::T_I64, actual)))

#ifdef RB_FLOAT_TYPE_P
#define FLOAT_TYPE_P(x) RB_FLOAT_TYPE_P(x)
#else
#define FLOAT_TYPE_P(x) RB_TYPE_P(x, T_FLOAT)
#endif

#define CONVERT_FLOAT(CONVERT)                                 \
  ((FLOAT_TYPE_P(actual))                                      \
       ? CONVERT(actual)                                       \
       : raise_type_mismatch_as_value(outer_struct, field_sym, \
                                      protocol::T_DOUBLE, actual))

static inline bool convertBool(VALUE actual, VALUE outer_struct,
                               VALUE field_sym) {
  switch (actual) {
    case Qtrue:
      return true;
    case Qfalse:
      return false;
    default:
      raise_type_mismatch(outer_struct, field_sym, protocol::T_BOOL, actual);
  }

  /* unreachable */
  return false;
}

#ifdef FIX2SHORT
#define SHORT_CONVERT(x) FIX2SHORT(x)
#else
#define SHORT_CONVERT(x) ((short)FIX2INT(x))
#endif

static inline char byte_convert(VALUE x) {
  short s = SHORT_CONVERT(x);

  if (s <= 127 && s >= -128) {
    return (char)s;
  } else {
    rb_raise(rb_eRangeError, "integer %d out of range for char", s);
  }

  /* unreachable */
  return 0;
}

void ThriftSerializer::writeAny(TType ttype, FieldInfo *field_info,
                                VALUE actual, VALUE outer_struct,
                                VALUE field_sym) {
  switch (ttype) {
    HANDLE_TYPE(I16, I16, CONVERT_FIXNUM(SHORT_CONVERT, protocol::T_I16))
    HANDLE_TYPE(I32, I32, CONVERT_FIXNUM(FIX2INT, protocol::T_I32))
    HANDLE_TYPE(I64, I64, CONVERT_I64)
    HANDLE_TYPE(BOOL, Bool, convertBool(actual, outer_struct, field_sym))
    HANDLE_TYPE(DOUBLE, Double, CONVERT_FLOAT(NUM2DBL))
    HANDLE_TYPE(BYTE, Byte, CONVERT_FIXNUM(byte_convert, protocol::T_BYTE))

    case protocol::T_STRING: {
      Sparsam_Check_Type(actual, T_STRING, outer_struct, field_sym);

      string data = string(StringValuePtr(actual), RSTRING_LEN(actual));
      if (field_info->isBinaryString) {
        this->tprot->writeBinary(data);
      } else {
        this->tprot->writeString(data);
      }
      break;
    }

    case protocol::T_LIST: {
      Sparsam_Check_Type(actual, T_ARRAY, outer_struct, field_sym);

      long length = RARRAY_LEN(actual);
      TType elem = field_info->elementType->ftype;
      this->tprot->writeListBegin(elem, static_cast<size_t>(length));
      for (long i = 0; i < length; i++) {
        this->writeAny(elem, field_info->elementType, rb_ary_entry(actual, i),
                       outer_struct, field_sym);
      }
      this->tprot->writeListEnd();
      break;
    }

    case protocol::T_SET: {
      if (rb_class_real(CLASS_OF(actual)) != SetClass) {
        raise_type_mismatch(outer_struct, field_sym, protocol::T_SET, actual);
      }

      VALUE ary = rb_funcall(actual, intern_for_to_a, 0);
      long length = RARRAY_LEN(ary);
      TType elem = field_info->elementType->ftype;
      this->tprot->writeListBegin(elem, static_cast<size_t>(length));
      for (long i = 0; i < length; i++) {
        this->writeAny(elem, field_info->elementType, rb_ary_entry(ary, i),
                       outer_struct, field_sym);
      }
      this->tprot->writeListEnd();
      break;
    }

    case protocol::T_MAP: {
      Sparsam_Check_Type(actual, T_HASH, outer_struct, field_sym);

      TType keyTType = field_info->keyType->ftype,
            valueTType = field_info->elementType->ftype;
      this->tprot->writeMapBegin(keyTType, valueTType,
                                 static_cast<size_t>(RHASH_SIZE(actual)));
      HASH_FOREACH_BEGIN(actual, this, field_info, &outer_struct, &field_sym)
      ThriftSerializer *that = (ThriftSerializer *)argv[0];
      FieldInfo *field_info = (FieldInfo *)argv[1];
      VALUE *outer_struct = (VALUE *)argv[2];
      VALUE *field_sym = (VALUE *)argv[3];
      that->writeAny(field_info->keyType->ftype, field_info->keyType, k,
                     *outer_struct, *field_sym);
      that->writeAny(field_info->elementType->ftype, field_info->elementType, v,
                     *outer_struct, *field_sym);
      HASH_FOREACH_END()
      this->tprot->writeMapEnd();
      break;
    }

    case protocol::T_STRUCT: {
      if (rb_class_real(CLASS_OF(actual)) != field_info->klass) {
        raise_type_mismatch(outer_struct, field_sym, field_info->klass, actual);
      }

      static const string cname = "";
      this->tprot->writeStructBegin(cname.c_str());
      this->writeStruct(field_info->klass, actual);
      this->tprot->writeFieldStop();
      this->tprot->writeStructEnd();
      break;
    }

    default: { break; }
  }
}

#undef HANDLE_TYPE

static bool checkRequiredFields(VALUE klass, VALUE data) {
  auto fields = FindOrCreateFieldInfoMap(klass);
  for (auto const &entry : *fields) {
    if (!entry.second->isOptional) {
      VALUE val = rb_ivar_get(data, entry.second->ivarName);
      if (NIL_P(val)) {
        raise_exc_with_struct_and_field_names(
            SparsamMissingMandatory, rb_str_new2("Required field missing"),
            klass, entry.second->symName);
        return false;
      }
    }
  }

  return true;
}

void ThriftSerializer::writeStruct(VALUE klass, VALUE data) {
  static const string cname = "";
  FieldBegin fieldBegin;
  FieldInfo *fieldInfo;
  auto fields = FindOrCreateFieldInfoMap(klass);

  if (!checkRequiredFields(klass, data)) {
    return;
  }

  for (auto const &entry : *fields) {
    fieldBegin.fid = entry.first;
    fieldInfo = entry.second;
    fieldBegin.ftype = fieldInfo->ftype;
    VALUE actual = rb_ivar_get(data, fieldInfo->ivarName);
    if (!NIL_P(actual)) {
      this->tprot->writeFieldBegin(cname.c_str(), fieldBegin.ftype,
                                   fieldBegin.fid);
      this->writeAny(fieldBegin.ftype, entry.second, actual, data,
                     fieldInfo->symName);
      this->tprot->writeFieldEnd();
    }
  }
}

VALUE serializer_writeStruct(VALUE self, VALUE klass, VALUE data) {
  watch_for_texcept() get_ts();

  if (rb_class_real(CLASS_OF(data)) != klass) {
    VALUE expected_name = rb_class_name(klass);
    VALUE actual_name = rb_class_name(CLASS_OF(data));

    raise_exc_with_struct_and_field_names(
        SparsamTypeMismatchError,
        rb_sprintf("Mismatched type passed to serialize (expected: %s got: %s)",
                   RSTRING_PTR(expected_name), RSTRING_PTR(actual_name)),
        CLASS_OF(data), ID2SYM(rb_intern("(root)")));

    RB_GC_GUARD(expected_name);
    RB_GC_GUARD(actual_name);
  }

  static const string cname = "";
  ts->tprot->writeStructBegin(cname.c_str());
  ts->writeStruct(klass, data);
  ts->tprot->writeFieldStop();
  ts->tprot->writeStructEnd();
  string retval = ts->tmb->getBufferAsString();
  return rb_str_new(retval.c_str(), retval.size());
  catch_thrift_and_reraise();
}

VALUE serializer_readStruct(VALUE self, VALUE klass) {
  watch_for_texcept() get_ts();
  string cname;
  VALUE ret;
  ts->tprot->readStructBegin(cname);
  ret = ts->readStruct(klass);
  ts->tprot->readStructEnd();
  return ret;
  catch_thrift_and_reraise();
}

VALUE serializer_readUnion(VALUE self, VALUE klass) {
  watch_for_texcept() get_ts();
  string cname;
  VALUE ret;
  ts->tprot->readStructBegin(cname);
  ret = ts->readUnion(klass);
  ts->tprot->readStructEnd();
  return ret;
  catch_thrift_and_reraise();
}

#define R_FIX_TO_TTYPE(x) (static_cast<TType>(FIX2INT(x)))

FieldInfoMap *FindOrCreateFieldInfoMap(VALUE klass) {
  auto iter = klassCache.find(klass);
  if (iter == klassCache.end()) {
    if (RTEST(rb_class_inherited_p(klass, klass_for_union))) {
      unions.insert(klass);
    }
    auto ret = CreateFieldInfoMap(klass);
    klassCache[klass] = ret;
    return ret;
  } else {
    return iter->second;
  }
}

ID field_name_to_ivar_id(VALUE str_name) {
  if (str_name != Qnil) {
    return rb_intern_str(rb_str_concat(rb_str_new2("@"), str_name));
  } else {
    return 0;
  }
}

VALUE field_name_to_sym(VALUE str_name) {
  if (str_name != Qnil) {
    return ID2SYM(rb_intern_str(str_name));
  } else {
    return 0;
  }
}

// each FieldInfoMap has multiple FieldInfos
FieldInfo *CreateFieldInfo(VALUE field_map_entry) {
  FieldInfo *fieldInfo = new FieldInfo();
  fieldInfo->ftype =
      R_FIX_TO_TTYPE(rb_hash_aref(field_map_entry, sym_for_type));
  fieldInfo->isOptional =
      RTEST(rb_hash_aref(field_map_entry, sym_for_optional));
  fieldInfo->ivarName =
      field_name_to_ivar_id(rb_hash_aref(field_map_entry, sym_for_name));
  fieldInfo->symName =
      field_name_to_sym(rb_hash_aref(field_map_entry, sym_for_name));
  switch (fieldInfo->ftype) {
    case protocol::T_STRING: {
      if (RTEST(rb_hash_aref(field_map_entry, sym_for_binary))) {
        fieldInfo->isBinaryString = true;
      }
      break;
    }
    case protocol::T_STRUCT: {
      fieldInfo->klass = rb_hash_aref(field_map_entry, sym_for_class);
      break;
    }
    case protocol::T_LIST:
    case protocol::T_SET: {
      fieldInfo->elementType =
          CreateFieldInfo(rb_hash_aref(field_map_entry, sym_for_element));
      break;
    }
    case protocol::T_MAP: {
      fieldInfo->keyType =
          CreateFieldInfo(rb_hash_aref(field_map_entry, sym_for_key));
      fieldInfo->elementType =
          CreateFieldInfo(rb_hash_aref(field_map_entry, sym_for_value));
      break;
    }
    default:
      break;
  }
  return fieldInfo;
}

// each klass has a FieldInfoMap
FieldInfoMap *CreateFieldInfoMap(VALUE klass) {
  FieldInfoMap *fieldMap = new FieldInfoMap();
  VALUE field_map = rb_const_get_at(klass, intern_for_FIELDS);

  HASH_FOREACH_BEGIN(field_map, fieldMap)
  FieldInfoMap *fieldMap = (FieldInfoMap *)argv[0];
  (*fieldMap)[FIX2INT(k)] = CreateFieldInfo(v);
  HASH_FOREACH_END()
  return fieldMap;
}

VALUE cache_fields(VALUE self, VALUE klass) {
  FindOrCreateFieldInfoMap(klass);
  return Qnil;
}

#undef HASH_FOREACH_BEGIN
#undef HASH_FOEACH_RET
#undef HASH_FOREACH_END

#undef R_FIX_TO_TTYPE

ext: cpp

max_stars_repo_name: gridgentoo/sparsam

max_issues_repo_name: gridgentoo/sparsam

max_forks_repo_name: gridgentoo/sparsam


text:
/* Copyright 2013 Stanford University.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 vd* 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 copyright holders 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.
 */

/*
 * Author: Chinmayee Shah <chinmayee.shah@stanford.edu>
 */

#include "applications/physbam/water//app_utils.h"
#include "applications/physbam/water//app_data_include.h"
#include "applications/physbam/water//app_data_prototypes.h"
#include "applications/physbam/water//data_include.h"
#include "applications/physbam/water//parameters.h"
#include "applications/physbam/water//job_include.h"
#include "applications/physbam/water//reg_def.h"
#include "applications/physbam/water//water_app.h"
#include "src/data/physbam/translator_physbam.h"
#include "src/data/physbam/translator_physbam_old.h"
#include "src/data/scalar_data.h"
#include "src/data/scratch_data_helper.h"
#include <PhysBAM_Tools/Log/DEBUG_SUBSTEPS.h>
#include <PhysBAM_Tools/Log/LOG.h>
#include <PhysBAM_Tools/Read_Write/Utilities/FILE_UTILITIES.h>
#include "src/shared/dbg.h"
#include "src/shared/geometric_region.h"
#include "src/shared/log.h"
#include "src/shared/nimbus.h"
#include "src/shared/timer.h"
#include "stdio.h"

#include <boost/program_options.hpp>
extern "C" nimbus::Application * ApplicationBuilder(int argc, char *argv[]) {
  namespace po = boost::program_options;

  uint64_t scale;
  uint64_t part_num_x;
  uint64_t part_num_y;
  uint64_t part_num_z;
  uint64_t projection_part_num_x;
  uint64_t projection_part_num_y;
  uint64_t projection_part_num_z;
  uint64_t last_frame;
  uint64_t max_iterations;
  uint64_t iteration_batch;
  uint64_t smart_projection_level;
  float water_level;

  po::options_description desc("Water Application Options");
  desc.add_options()
    ("help,h", "produce help message")

    // Optinal arguments
    ("scale,s", po::value<uint64_t>(&scale), "scale of the simulation") //NOLINT
    ("pnx", po::value<uint64_t>(&part_num_x), "partition number along x") //NOLINT
    ("pny", po::value<uint64_t>(&part_num_y), "partition number along y") //NOLINT
    ("pnz", po::value<uint64_t>(&part_num_z), "partition number along z") //NOLINT
    ("ppnx", po::value<uint64_t>(&projection_part_num_x), "projection partition number along x") //NOLINT
    ("ppny", po::value<uint64_t>(&projection_part_num_y), "projection partition number along y") //NOLINT
    ("ppnz", po::value<uint64_t>(&projection_part_num_z), "projection partition number along z") //NOLINT
    ("last_frame,e", po::value<uint64_t>(&last_frame), "last frame to compute") //NOLINT
    ("maxi", po::value<uint64_t>(&max_iterations), "maximum projection iterations") //NOLINT
    ("ibatch", po::value<uint64_t>(&iteration_batch), "projection iteration batch") //NOLINT
    ("psl", po::value<uint64_t>(&smart_projection_level), "smart projection level") //NOLINT
    ("wl", po::value<float>(&water_level), "initial water level, float between 0 and 1") //NOLINT
    ("dpb", "deactivate projection bottleneck job")
    ("dgw", "deactivate one global write per frame");

  po::variables_map vm;
  try {
    po::store(po::parse_command_line(argc, argv, desc), vm);
  }
  catch(std::exception& e) { // NOLINT
    std::cerr << "ERROR: " << e.what() << "\n";
    return NULL;
  }

  if (vm.count("help")) {
    std::cout << desc << "\n";
    return NULL;
  }

  try {
    po::notify(vm);
  }
  catch(std::exception& e) { // NOLINT
    std::cerr << "ERROR: " << e.what() << "\n";
    return NULL;
  }

  application::WaterApp *app = new application::WaterApp();

  if (vm.count("scale")) {
    app->set_scale(scale);
  }

  if (vm.count("pnx")) {
    app->set_part_num_x(part_num_x);
  }

  if (vm.count("pny")) {
    app->set_part_num_y(part_num_y);
  }

  if (vm.count("pnz")) {
    app->set_part_num_z(part_num_z);
  }

  if (vm.count("ppnx")) {
    app->set_projection_part_num_x(projection_part_num_x);
  }

  if (vm.count("ppny")) {
    app->set_projection_part_num_y(projection_part_num_y);
  }

  if (vm.count("ppnz")) {
    app->set_projection_part_num_z(projection_part_num_z);
  }

  if (vm.count("last_frame")) {
    app->set_last_frame(last_frame);
  }

  if (vm.count("psl")) {
    app->set_smart_projection(smart_projection_level);
  }

  if (vm.count("wl")) {
    app->set_water_level(water_level);
  }

  if (vm.count("maxi")) {
    app->set_max_iterations(max_iterations);
  }

  if (vm.count("ibatch")) {
    app->set_iteration_batch(iteration_batch);
  }

  if (vm.count("dpb")) {
    app->set_spawn_projection_loop_bottleneck(false);
  }

  if (vm.count("dgw")) {
    app->set_global_write(false);
  }

  return app;
}

namespace application {

nimbus::PartitionHandler ph;

uint64_t kScale;
uint64_t kAppPartNum;
uint64_t kAppPartNumX;
uint64_t kAppPartNumY;
uint64_t kAppPartNumZ;
uint64_t kProjAppPartNum;
uint64_t kProjAppPartNumX;
uint64_t kProjAppPartNumY;
uint64_t kProjAppPartNumZ;
nimbus::GeometricRegion kDefaultRegion;
uint64_t kLastFrame;
uint64_t kMaxIterations;
uint64_t kIterationBatch;
bool kSpawnProjectionLoopBottleneck;
bool kUseGlobalWrite;

    WaterApp::WaterApp() {
      scale_ = DEFAULT_SCALE;
      part_num_x_ = DEFAULT_APP_PART_NUM_X;
      part_num_y_ = DEFAULT_APP_PART_NUM_Y;
      part_num_z_ = DEFAULT_APP_PART_NUM_Z;
      projection_part_num_x_ = DEFAULT_APP_PROJ_PART_NUM_X;
      projection_part_num_y_ = DEFAULT_APP_PROJ_PART_NUM_Y;
      projection_part_num_z_ = DEFAULT_APP_PROJ_PART_NUM_Z;
      last_frame_ = DEFAULT_LAST_FRAME;
      max_iterations_ = DEFAULT_MAX_ITERATIONS;
      iteration_batch_ = DEFAULT_ITERATION_BATCH;
      spawn_projection_loop_bottleneck_ =
        DEFAULT_SPAWN_PROJECTION_LOOP_BOTTLENECK;
      global_write_ = DEFAULT_USE_GLOBAL_WRITE;
      translator_log = NULL;
      smart_projection_ = DEFAULT_PROJECTION_SMART;
      water_level_ = DEFAULT_WATER_LEVEL;
    }

    float WaterApp::water_level() {
      return water_level_;
    }

    void WaterApp::set_scale(uint64_t scale) {
      scale_ = scale;
    }

    void WaterApp::set_part_num_x(uint64_t part_num_x) {
      part_num_x_ = part_num_x;
    }

    void WaterApp::set_part_num_y(uint64_t part_num_y) {
      part_num_y_ = part_num_y;
    }

    void WaterApp::set_part_num_z(uint64_t part_num_z) {
      part_num_z_ = part_num_z;
    }

    void WaterApp::set_projection_part_num_x(uint64_t projection_part_num_x) {
      projection_part_num_x_ = projection_part_num_x;
    }

    void WaterApp::set_projection_part_num_y(uint64_t projection_part_num_y) {
      projection_part_num_y_ = projection_part_num_y;
    }

    void WaterApp::set_projection_part_num_z(uint64_t projection_part_num_z) {
      projection_part_num_z_ = projection_part_num_z;
    }

    void WaterApp::set_last_frame(uint64_t last_frame) {
      last_frame_ = last_frame;
    }

    void WaterApp::set_max_iterations(uint64_t max_iterations) {
      max_iterations_ = max_iterations;
    }

    void WaterApp::set_iteration_batch(uint64_t iteration_batch) {
      iteration_batch_ = iteration_batch;
    }

    void WaterApp::set_spawn_projection_loop_bottleneck(bool flag) {
      spawn_projection_loop_bottleneck_ = flag;
    }

    void WaterApp::set_global_write(bool flag) {
      global_write_ = flag;
    }

    void WaterApp::set_smart_projection(uint64_t level) {
      smart_projection_ = level;
    }

    void WaterApp::set_water_level(float water_level) {
      water_level_ = water_level;
    }

    /* Register data and job types and initialize constant quantities used by
     * application jobs. */
    void WaterApp::Load() {
        ScopeTimer::Initialize(true);
        //nimbus::Timer::Initialize();
        PhysBAM::Initialize_Particles();
        PhysBAM::Initialize_Read_Write_General_Structures();

        // initialize application specific parameters and constants
        assert((part_num_x_ % projection_part_num_x_) == 0);
        assert((part_num_y_ % projection_part_num_y_) == 0);
        assert((part_num_z_ % projection_part_num_z_) == 0);
        assert((water_level_ >= 0) && (water_level_ <= 1));
        kScale = scale_;
        kAppPartNumX = part_num_x_;
        kAppPartNumY = part_num_y_;
        kAppPartNumZ = part_num_z_;
        kAppPartNum = part_num_x_ * part_num_y_ * part_num_z_;
        kProjAppPartNumX = projection_part_num_x_;
        kProjAppPartNumY = projection_part_num_y_;
        kProjAppPartNumZ = projection_part_num_z_;
        kProjAppPartNum = projection_part_num_x_ * projection_part_num_y_ * projection_part_num_z_;
        kDefaultRegion.Rebuild(1, 1, 1, scale_, scale_, scale_);
        kLastFrame = last_frame_;
        kMaxIterations = max_iterations_;
        kIterationBatch = iteration_batch_;
        kSpawnProjectionLoopBottleneck = spawn_projection_loop_bottleneck_;
        kUseGlobalWrite = global_write_;

        // Initialize the app data prototypes.
        InitializeAppDataPrototypes(kDefaultRegion);

        // dbg_add_mode(APP_LOG_STR);
        // dbg_add_mode(TRANSLATE_STR);

        dbg(APP_LOG, "Loading water application\n");

        // Region constants
        // Old code, using the code generated by python scripts -omidm
        // InitializeRegions();

        ph.AddPartitions("kRegW3",
                         kScale, kScale, kScale, 3, 3, 3, 1, 1, 1,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegW3",
                         kScale, kScale, kScale, 3, 3, 3, 1, 1, 1,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kRegW3",
                         kScale, kScale, kScale, 3, 3, 3, 1, 1, 1,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kRegW1",
                         kScale, kScale, kScale, 1, 1, 1, 1, 1, 1,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegW1",
                         kScale, kScale, kScale, 1, 1, 1, 1, 1, 1,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kRegW1",
                         kScale, kScale, kScale, 1, 1, 1, 1, 1, 1,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kRegW0",
                         kScale, kScale, kScale, 0, 0, 0, 1, 1, 1,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL_WGB);

        ph.AddPartitions("kRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL_WGB);

        ph.AddPartitions("kRegY2W0",
                         kScale, kScale, kScale, 0, 0, 0,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kRegY2W8",
                         kScale, kScale, kScale, 8, 8, 8,
                         kAppPartNumX, kAppPartNumY, kAppPartNumZ,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kProjRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kProjRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kProjRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kProjRegY2W3",
                         kScale, kScale, kScale, 3, 3, 3,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL_WGB);

        ph.AddPartitions("kProjRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        ph.AddPartitions("kProjRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::OUTER);

        ph.AddPartitions("kProjRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::INNER);

        ph.AddPartitions("kProjRegY2W1",
                         kScale, kScale, kScale, 1, 1, 1,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL_WGB);

        ph.AddPartitions("kProjRegY2W0",
                         kScale, kScale, kScale, 0, 0, 0,
                         kProjAppPartNumX, kProjAppPartNumY, kProjAppPartNumZ,
                         nimbus::PartitionHandler::CENTRAL);

        // PhysBAM logging and R/W
        PhysBAM::LOG::Initialize_Logging(false, false, 1<<30, true, kThreadsNum);
        PhysBAM::FILE_UTILITIES::Create_Directory(kOutputDir+"/common");
        PhysBAM::LOG::Instance()->Copy_Log_To_File(kOutputDir+"/common/log.txt", false);

        // REGISTER DATA
        {
            dbg(APP_LOG, "Registering %s\n", APP_FACE_VEL);
            RegisterData(APP_FACE_VEL, new DataFaceArray<float>(APP_FACE_VEL));
            dbg(APP_LOG, "Registering %s\n", APP_FACE_VEL_GHOST);
            RegisterData(APP_FACE_VEL_GHOST, new DataFaceArray<float>(APP_FACE_VEL_GHOST));
            dbg(APP_LOG, "Registering %s\n", APP_PHI);
            RegisterData(APP_PHI, new DataScalarArray<float>(APP_PHI));
            dbg(APP_LOG, "Registering %s\n", APP_POS_PARTICLES);
            RegisterData(APP_POS_PARTICLES, new DataParticleArray(APP_POS_PARTICLES));
            dbg(APP_LOG, "Registering %s\n", APP_NEG_PARTICLES);
            RegisterData(APP_NEG_PARTICLES, new DataParticleArray(APP_NEG_PARTICLES));
            dbg(APP_LOG, "Registering %s\n", APP_POS_REM_PARTICLES);
            RegisterData(APP_POS_REM_PARTICLES, new DataParticleArray(APP_POS_REM_PARTICLES));
            dbg(APP_LOG, "Registering %s\n", APP_NEG_REM_PARTICLES);
            RegisterData(APP_NEG_REM_PARTICLES, new DataParticleArray(APP_NEG_REM_PARTICLES));
            dbg(APP_LOG, "Registering %s\n", APP_LAST_UNIQUE_PARTICLE_ID);
            RegisterData(APP_LAST_UNIQUE_PARTICLE_ID, new nimbus::ScalarData<int>(APP_LAST_UNIQUE_PARTICLE_ID));
            dbg(APP_LOG, "Registering %s\n", APP_DT);
            RegisterData(APP_DT, new nimbus::ScalarData<float>(APP_DT));

            // These Nimbus data types are used in projection but not used in
            // internal projection loop. They are generally used to generate the
            // matrixes or vectors used in internal projection loop.
            // PSI_D.
            dbg(APP_LOG, "Registering %s\n", APP_PSI_D);
            RegisterData(APP_PSI_D, new DataScalarArray<bool>(APP_PSI_D));
            // PSI_N.
            dbg(APP_LOG, "Registering %s\n", APP_PSI_N);
            RegisterData(APP_PSI_N, new DataFaceArray<bool>(APP_PSI_N));
            // PRESSURE.
            dbg(APP_LOG, "Registering %s\n", APP_PRESSURE);
            RegisterData(APP_PRESSURE, new DataScalarArray<float>(APP_PRESSURE));
            // FILLED_REGION_COLORS.
            dbg(APP_LOG, "Registering %s\n", APP_FILLED_REGION_COLORS);
            RegisterData(APP_FILLED_REGION_COLORS,
                new DataScalarArray<int>(APP_FILLED_REGION_COLORS));
            // DIVERGENCE.
            dbg(APP_LOG, "Registering %s\n", APP_DIVERGENCE);
            RegisterData(APP_DIVERGENCE, new DataScalarArray<float>(APP_DIVERGENCE));

            // These Nimbus data types are used in internal projection loop. They
            // are derived from boundary conditions and act as the linkage between
            // inside projection and outside projection.
            // MATRIX_A.
            dbg(APP_LOG, "Registering %s\n", APP_MATRIX_A);
            RegisterData(APP_MATRIX_A, new DataSparseMatrix(APP_MATRIX_A));
            // VECTOR_B.
            dbg(APP_LOG, "Registering %s\n", APP_VECTOR_B);
            RegisterData(APP_VECTOR_B, new DataRawVectorNd(APP_VECTOR_B));
            // VECTOR_X.
            // dbg(APP_LOG, "Registering %s\n", APP_VECTOR_X);
            // RegisterData(APP_VECTOR_X, new DataRawVectorNd(APP_VECTOR_X));
            // INDEX_C2M.
            dbg(APP_LOG, "Registering %s\n", APP_INDEX_C2M);
            RegisterData(APP_INDEX_C2M, new DataRawGridArray(APP_INDEX_C2M));
            // INDEX_M2C.
            dbg(APP_LOG, "Registering %s\n", APP_INDEX_M2C);
            RegisterData(APP_INDEX_M2C, new DataRawArrayM2C(APP_INDEX_M2C));
            // PROJECTION_LOCAL_N.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_LOCAL_N);
            RegisterData(APP_PROJECTION_LOCAL_N,
                new nimbus::ScalarData<int>(APP_PROJECTION_LOCAL_N));
            // PROJECTION_INTERIOR_N.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_INTERIOR_N);
            RegisterData(APP_PROJECTION_INTERIOR_N,
                new nimbus::ScalarData<int>(APP_PROJECTION_INTERIOR_N));

            // These Nimbus data types are used only in internal projection loop.
            // They are mostly static configurations set for projeciton.
            // PROJECTION_LOCAL_TOLERANCE.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_LOCAL_TOLERANCE);
            RegisterData(APP_PROJECTION_LOCAL_TOLERANCE,
                new nimbus::ScalarData<float>(APP_PROJECTION_LOCAL_TOLERANCE));
            // PROJECTION_GLOBAL_TOLERANCE.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_GLOBAL_TOLERANCE);
            RegisterData(APP_PROJECTION_GLOBAL_TOLERANCE,
                new nimbus::ScalarData<float>(APP_PROJECTION_GLOBAL_TOLERANCE));
            // PROJECTION_GLOBAL_N.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_GLOBAL_N);
            RegisterData(APP_PROJECTION_GLOBAL_N,
                new nimbus::ScalarData<int>(APP_PROJECTION_GLOBAL_N));
            // PROJECTION_DESIRED_N.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_DESIRED_ITERATIONS);
            RegisterData(APP_PROJECTION_DESIRED_ITERATIONS,
                new nimbus::ScalarData<int>(APP_PROJECTION_DESIRED_ITERATIONS));

            // These Nimbus data types are used in internal projection loop.
            // PROJECTION_LOCAL_RESIDUAL.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_LOCAL_RESIDUAL);
            RegisterData(APP_PROJECTION_LOCAL_RESIDUAL,
                new nimbus::ScalarData<double>(APP_PROJECTION_LOCAL_RESIDUAL));
            // PROJECTION_LOCAL_RHO.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_LOCAL_RHO);
            RegisterData(APP_PROJECTION_LOCAL_RHO,
                new nimbus::ScalarData<double>(APP_PROJECTION_LOCAL_RHO));
            // PROJECTION_GLOBAL_RHO.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_GLOBAL_RHO);
            RegisterData(APP_PROJECTION_GLOBAL_RHO,
                new nimbus::ScalarData<double>(APP_PROJECTION_GLOBAL_RHO));
            // PROJECTION_GLOBAL_RHO_OLD.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_GLOBAL_RHO_OLD);
            RegisterData(APP_PROJECTION_GLOBAL_RHO_OLD,
                new nimbus::ScalarData<double>(APP_PROJECTION_GLOBAL_RHO_OLD));
            // PROJECTION_LOCAL_DOT_PRODUCT_FOR_ALPHA.
            dbg(APP_LOG, "Registering %s\n",
                APP_PROJECTION_LOCAL_DOT_PRODUCT_FOR_ALPHA);
            RegisterData(APP_PROJECTION_LOCAL_DOT_PRODUCT_FOR_ALPHA,
                new nimbus::ScalarData<double>(
                    APP_PROJECTION_LOCAL_DOT_PRODUCT_FOR_ALPHA));
            // PROJECTION_ALPHA.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_ALPHA);
            RegisterData(APP_PROJECTION_ALPHA,
                new nimbus::ScalarData<float>(APP_PROJECTION_ALPHA));
            // PROJECTION_BETA.
            dbg(APP_LOG, "Registering %s\n", APP_PROJECTION_BETA);
            RegisterData(APP_PROJECTION_BETA,
                new nimbus::ScalarData<float>(APP_PROJECTION_BETA));
            // MATRIX_C.
            dbg(APP_LOG, "Registering %s\n", APP_MATRIX_C);
            RegisterData(APP_MATRIX_C, new DataSparseMatrix(APP_MATRIX_C));
            // VECTOR_Z.
            dbg(APP_LOG, "Registering %s\n", APP_VECTOR_PRESSURE);
            RegisterData(APP_VECTOR_PRESSURE,
                         new DataRawVectorNd(APP_VECTOR_PRESSURE));
            // VECTOR_Z.
            dbg(APP_LOG, "Registering %s\n", APP_VECTOR_Z);
            RegisterData(APP_VECTOR_Z, new DataRawVectorNd(APP_VECTOR_Z));
            // VECTOR_P_META_FORMAT.
            dbg(APP_LOG, "Registering %s\n", APP_VECTOR_P_META_FORMAT);
            RegisterData(APP_VECTOR_P_META_FORMAT,
                         new DataCompressedScalarArray<float>(APP_VECTOR_P_META_FORMAT));
            // VECTOR_TEMP.
            dbg(APP_LOG, "Registering %s\n", APP_VECTOR_TEMP);
            RegisterData(APP_VECTOR_TEMP, new DataRawVectorNd(APP_VECTOR_TEMP));

            dbg(APP_LOG, "Registering scratch %s\n", APP_POS_PARTICLES);
            kScratchPosParticles.RegisterScratchNames(this, new DataParticleArray(APP_POS_PARTICLES));
            dbg(APP_LOG, "Registering scratch %s\n", APP_NEG_PARTICLES);
            kScratchNegParticles.RegisterScratchNames(this, new DataParticleArray(APP_NEG_PARTICLES));
            dbg(APP_LOG, "Registering scratch %s\n", APP_POS_REM_PARTICLES);
            kScratchPosRemParticles.RegisterScratchNames(this, new DataParticleArray(APP_POS_REM_PARTICLES));
            dbg(APP_LOG, "Registering scratch %s\n", APP_NEG_REM_PARTICLES);
            kScratchNegRemParticles.RegisterScratchNames(this, new DataParticleArray(APP_NEG_REM_PARTICLES));
        }

        // APP DATA INNER REGIONS
        {
            kAppDataFaceVel.set_inner_delta(ph.map()["kRegY2W3Inner"][0].Delta());
            kAppDataFaceVelGhost.set_inner_delta(ph.map()["kRegY2W3Inner"][0].Delta());
            kAppDataPsiN.set_inner_delta(ph.map()["kRegY2W1Inner"][0].Delta());
            kAppDataPsiD.set_inner_delta(ph.map()["kRegY2W1Inner"][0].Delta());
            kAppDataPhi3.set_inner_delta(ph.map()["kRegY2W8Inner"][0].Delta());
            kAppDataPhi7.set_inner_delta(ph.map()["kRegY2W8Inner"][0].Delta());
            kAppDataPhi8.set_inner_delta(ph.map()["kRegY2W8Inner"][0].Delta());
            kAppDataPLE.set_inner_delta(ph.map()["kRegY2W3Inner"][0].Delta());
            kAppDataPressure.set_inner_delta(ph.map()["kRegY2W1Inner"][0].Delta());
            kAppDataColors.set_inner_delta(ph.map()["kProjRegY2W1Inner"][0].Delta());
            kAppDataDivergence.set_inner_delta(ph.map()["kProjRegY2W1Inner"][0].Delta());
            kAppDataIndexC2M.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataArrayM2C.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataMetaP.set_inner_delta(ph.map()["kProjRegY2W1Inner"][0].Delta());
            kAppDataSparseMatrixA.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataSparseMatrixC.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataVectorB.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataVectorPressure.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataVectorZ.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
            kAppDataVectorTemp.set_inner_delta(ph.map()["kProjRegY2W0Central"][0].Delta());
        }

        // REGISTER JOBS
        {
            RegisterJob(MAIN, new JobMain(this));
            RegisterJob(INITIALIZE, new JobInitialize(this));
            RegisterJob(CALCULATE_DT, new JobCalculateDt(this));
            RegisterJob(UPDATE_GHOST_VELOCITIES, new JobAdjustPhiWithObjects(this));
            RegisterJob(EXTRAPOLATE_PHI, new JobExtrapolatePhi(this));
            RegisterJob(ADVECT_PHI, new JobAdvectPhi(this));
            RegisterJob(STEP_PARTICLES, new JobStepParticles(this));
            RegisterJob(ADVECT_REMOVED_PARTICLES, new JobAdvectRemovedParticles(this));
            RegisterJob(ADVECT_V, new JobAdvectV(this));
            RegisterJob(APPLY_FORCES, new JobApplyForces(this));
            RegisterJob(LOOP_ITERATION, new JobLoopIteration(this));
            RegisterJob(LOOP_ITERATION_PART_TWO, new JobLoopIterationPartTwo(this));
            RegisterJob(LOOP_FRAME, new JobLoopFrame(this));
            RegisterJob(WRITE_OUTPUT, new JobWriteOutput(this));
            RegisterJob(RESEED_PARTICLES, new JobReseedParticles(this));
            RegisterJob(MODIFY_LEVELSET_PART_ONE, new JobModifyLevelsetPartOne(this));
            RegisterJob(MODIFY_LEVELSET_PART_TWO, new JobModifyLevelsetPartTwo(this));
            RegisterJob(MAKE_SIGNED_DISTANCE, new JobMakeSignedDistance(this));
            RegisterJob(ADJUST_PHI, new JobAdjustPhi(this));
            RegisterJob(DELETE_PARTICLES, new JobDeleteParticles(this));
            RegisterJob(REINCORPORATE_PARTICLES, new JobReincorporateRemovedParticles(this));
            RegisterJob(PROJECTION_MAIN, new JobProjectionMain(this));
            RegisterJob(PROJECTION_TRANSFORM_PRESSURE,
                        new JobProjectionTransformPressure(this));
            RegisterJob(PROJECTION_CALCULATE_BOUNDARY_CONDITION_PART_ONE,
                        new JobProjectionCalculateBoundaryConditionPartOne(this));
            RegisterJob(PROJECTION_CALCULATE_BOUNDARY_CONDITION_PART_TWO,
                        new JobProjectionCalculateBoundaryConditionPartTwo(this));
            RegisterJob(PROJECTION_CONSTRUCT_MATRIX,
                        new JobProjectionConstructMatrix(this));
            RegisterJob(PROJECTION_WRAPUP, new JobProjectionWrapup(this));
            RegisterJob(EXTRAPOLATION, new JobExtrapolation(this));
            RegisterJob(SYNCHRONIZE_PARTICLES, new JobSynchronizeParticles(this));

            RegisterJob(PROJECTION_GLOBAL_INITIALIZE,
                        new JobProjectionGlobalInitialize(this));
            RegisterJob(PROJECTION_LOCAL_INITIALIZE,
                        new JobProjectionLocalInitialize(this));
            RegisterJob(PROJECTION_LOOP_ITERATION, new JobProjectionLoopIteration(this));
            RegisterJob(PROJECTION_LOOP_BOTTLENECK, new JobProjectionLoopBottleneck(this));
            RegisterJob(PROJECTION_STEP_ONE, new JobProjectionStepOne(this));
            RegisterJob(PROJECTION_REDUCE_RHO, new JobProjectionReduceRho(this));
            RegisterJob(PROJECTION_STEP_TWO, new JobProjectionStepTwo(this));
            RegisterJob(PROJECTION_STEP_THREE, new JobProjectionStepThree(this));
            RegisterJob(PROJECTION_REDUCE_ALPHA,
                        new JobProjectionReduceAlpha(this));
            RegisterJob(PROJECTION_STEP_FOUR, new JobProjectionStepFour(this));
            RegisterJob(BARRIER_JOB, new JobBarrier(this));

            RegisterStaticConfigPrototype(
                STATIC_CONFIG_VALID_MASK, new StaticConfigValidMask(kDefaultRegion));
            RegisterStaticConfigPrototype(
                STATIC_CONFIG_U_INTERFACE, new StaticConfigUInterface(kDefaultRegion));
            RegisterStaticConfigPrototype(
                STATIC_CONFIG_FORCE, new StaticConfigForce(kDefaultRegion));
            RegisterStaticConfigPrototype(
                STATIC_CONFIG_COLLISION_BODY,
                new StaticConfigCollisionBody(kDefaultRegion));
        }

        nimbus::TranslatorPhysBAM<float>::log = translator_log;
        nimbus::TranslatorPhysBAMOld<PhysBAM::VECTOR<float, 3> >::log = translator_log;
        dbg(APP_LOG, "Completed loading water application\n");
    }

} // namespace application

ext: cc

max_stars_repo_name: schinmayee/nimbus

max_issues_repo_name: schinmayee/nimbus

max_forks_repo_name: schinmayee/nimbus