text:
#include "multi_alignment_model.h"
#include <blitz/tinyvec.h>
#include "io.h"

multi_alignment_model::multi_alignment_model( int i ) : nmodels(i), PHMM(kstates+i*SUBSIZE), dirty(true), stop_at_repeat(false) {
	set_default_parameters( blitz::TinyVector<double,4>(0.25,0.25,0.25,0.25) );
}

void multi_alignment_model::reset() { store.clear(); dstore.clear(); }

void multi_alignment_model::start_training() {
	PHMM.start_training();
}

ereal multi_alignment_model::end_training() {
	ereal pr = PHMM.end_training();

	// We need to normalize the model.
	for( int k = 0; k < nmodels; ++k ) {
		int len = get_region_length(k);
		nw_model_parameters_ptr p = get_parameters(k);
		set_parameters( k, *p, len );
	}
	return pr;
}

void multi_alignment_model::set_default_parameters( blitz::TinyVector<double,4> nq ) {
	using namespace std;
	using namespace boost;
	assert(nmodels < 30);

	typedef boost::tuple<double,double,double> double3;
	vector<double3> custom_pts;
	custom_pts.push_back( double3(0.4,0.0625,2.0) );
	custom_pts.push_back( double3(0.2,0.03125,1.333) );
	custom_pts.push_back( double3(0.6,0.0625,3.41424) );
	custom_pts.push_back( double3(0.8,0.03125,2.0) );

	assert( (unsigned int)nmodels <= custom_pts.size() );
	for( int k = 0; k < nmodels; ++k ) {
		nw_model_parameters newp;
		newp.q = nq;
		newp.set_p_hky( custom_pts[k].get<0>(), 1.0 );
		newp.pr_open = custom_pts[k].get<1>();
		newp.mean_gap_length = custom_pts[k].get<2>();
		set_parameters( k, newp, 200 );

		int ms = kstates+k*SUBSIZE;
		// set switch
		PHMM.set_transition(state_switch,ms+state_S,1.0/nmodels);

		// set end to switch
		PHMM.set_transition(ms+state_E,state_switch,1.0);

		// set silent states for model 0
		PHMM.set_emmission(ms+state_S,APHMM::NOTHING,APHMM::NOTHING,1.0);
		PHMM.set_emmission(ms+state_E,APHMM::NOTHING,APHMM::NOTHING,1.0);
	}

	// set silent states for switch and random
	PHMM.set_emmission(state_switch,APHMM::NOTHING,APHMM::NOTHING,1.0);

	// set start and end states
	PHMM.set_start_state(state_switch);
	PHMM.set_end_state(state_switch);

	set_random_distribution(nq);
	PHMM.optimize();
	dirty = false;
}

int multi_alignment_model::num_models() const {
	return nmodels;
}

double multi_alignment_model::get_region_pr( int k ) const {
	return (double)PHMM.get_transition(state_switch,kstates + k*SUBSIZE + state_S);
}

double multi_alignment_model::get_region_length( int k ) const {
	// since we didn't tie the "t" parts, we'll just average all three values
	int ms = kstates + k*SUBSIZE;
	double sum = (double)PHMM.get_transition(ms+state_M,ms+state_E);
	sum += (double)PHMM.get_transition(ms+state_Ia,ms+state_E);
	sum += (double)PHMM.get_transition(ms+state_Ib,ms+state_E);
	double p_1mt = sum/3.0;
	return 1.0/p_1mt;
}

void multi_alignment_model::set_random_distribution( blitz::TinyVector<double,4> nq ) {
	dirty = true;
	for( int i = 0; i < 4; ++i) q(i) = nq(i);
}

void multi_alignment_model::set_stop_at_repeat( bool t ) {
	stop_at_repeat = t;
}

void multi_alignment_model::set_parameters( int k, nw_model_parameters const &nwp, int len ) {
	dirty = true;
	double pt = 1.0-1.0/(double)len;
	int ms = kstates + k*SUBSIZE;

	// set open
    double p_at = nwp.pr_open*pt;
	double p_1m2at = (1.0-2.0*nwp.pr_open)*pt;
	PHMM.set_transition(ms+state_M,ms+state_Ia,p_at);
	PHMM.set_transition(ms+state_M,ms+state_Ib,p_at);
	PHMM.set_transition(ms+state_M,ms+state_M,p_1m2at);

	// set start
	double p_a = nwp.pr_open;
	double p_1m2a = 1.0 - 2.0*nwp.pr_open;
	PHMM.set_transition(ms+state_S,ms+state_Ia,p_a);
	PHMM.set_transition(ms+state_S,ms+state_Ib,p_a);
	PHMM.set_transition(ms+state_S,ms+state_M,p_1m2a);

	// Since a gap is length 1 by the virtue of opening it, we adjust n down by 1.
	// The remaining gap extension probability corresponds to the mean of a
	// geometric distribution allowing zero extensions.
	double p_b = 1.0-1.0/nwp.mean_gap_length;
	double p_bt = p_b*pt;
	double p_1mbt = (1.0 - p_b)*pt;
	PHMM.set_transition(ms+state_Ia,ms+state_Ia,p_bt);
	PHMM.set_transition(ms+state_Ia,ms+state_M,p_1mbt);
	PHMM.set_transition(ms+state_Ib,ms+state_Ib,p_bt);
	PHMM.set_transition(ms+state_Ib,ms+state_M,p_1mbt);

	// set end transitions
	double p_1mt = 1.0 - pt;
	PHMM.set_transition( ms+state_M,ms+state_E,p_1mt);
	PHMM.set_transition( ms+state_Ia,ms+state_E,p_1mt);
	PHMM.set_transition( ms+state_Ib,ms+state_E,p_1mt);

	// set background
	for( int i = 0; i < 4; i++ ) {
      PHMM.set_emmission(ms+state_Ia,i,APHMM::NOTHING,nwp.q(i));
      PHMM.set_emmission(ms+state_Ib,APHMM::NOTHING,i,nwp.q(i));
   }

	// set p
	for( int i = 0; i < 4; i++ ) {
		for( int j = 0; j < 4; j++ ) {
			PHMM.set_emmission(ms+state_M,i,j,nwp.p(i,j));
		}
   }
}

nw_model_parameters_ptr multi_alignment_model::get_parameters( int k ) const {
	nw_model_parameters_ptr nwp( new nw_model_parameters() );
	int ms = kstates + SUBSIZE*k;
	double t = 1.0 - 1.0/get_region_length(k);
	nwp->pr_open = (double) PHMM.get_transition(ms+state_M,ms+state_Ia)/t;
	nwp->pr_open += (double) PHMM.get_transition(ms+state_M,ms+state_Ib)/t;
	nwp->pr_open /= 2.0;

	double gap_a = 1.0/(1.0- (double)PHMM.get_transition(ms+state_Ia,ms+state_Ia)/t );
	double gap_b = 1.0/(1.0- (double)PHMM.get_transition(ms+state_Ib,ms+state_Ib)/t );
	nwp->mean_gap_length = (gap_a + gap_b)/2.0;

	for( int x = 0; x < APHMM::NOTHING; ++x ) {
		double qx = 0.0;
		for( int y = 0; y < APHMM::NOTHING; ++y ) {
			nwp->p(x,y) = (double)PHMM.get_emmission(ms+state_M,x,y);
			qx += nwp->p(x,y);
		}
		nwp->q(x) = qx;
	}

	return nwp;
}

blitz::TinyVector<double,4> multi_alignment_model::get_random_distribution() const {
	return q;
}

void multi_alignment_model::pretty_print_parameters( std::ostream &out ) const {
	using namespace std;

	out  << "Random background: " << get_random_distribution() << endl;

	for( int k = 0; k < num_models(); ++k ) {
		out << endl << "Model " << k << ": " << get_region_pr(k) << endl;
		out << "Length: " << get_region_length(k) << endl;
		nw_model_parameters_ptr nwp = get_parameters(k);
		nwp->pretty_print(out);
		double pr_match = 0;
		for( int i = 0; i < 4; ++i ) pr_match += nwp->p(i,i);
	}
}

void multi_alignment_model::write_parameters( std::ostream &out ) const {
	using namespace std;

	out  << "# Random background: " << get_random_distribution() << endl << endl;

	for( int j = 0; j < num_models(); ++j ) {
		out << "# Region " << j << " with pr " << get_region_pr(j);
		out << " and length " << get_region_length(j) << "." << endl;
		nw_model_parameters_ptr nwp = get_parameters(j);
		nwp->pretty_print( out, "# " );
		out << endl;
	}
	out << endl << "# Machine readable parameters follow. " << "---------------------------------------------" << endl;

	// Write dummy lengths for random background.
	out << "0.0 200 200 ";
	blitz::TinyVector<double,4> q = get_random_distribution();
	for( int i = 0; i < 4; ++i )  out << q(i) << " ";
	out << endl;

	out << num_models() << endl;
	for( int j = 0; j < num_models(); ++j ) {
		out <<  get_region_pr(j) << " " << get_region_length(j) << " " << *get_parameters(j) << endl;
	}

}

// We don't actually do anything here. This is legacy to take care of mm_aligner.
void multi_alignment_model::disable_random_model() {}

void multi_alignment_model::read_parameters( std::istream &in ) {
	using namespace std;
	vector<double> mdist;
	mdist.push_back( read<double>(in) );
	// Skip unused random model parameters.
	read<double>(in);	read<double>(in);
	for( int i = 0; i < 4; ++i ) q(i) = read<double>(in);

	int nmodels = read<int>(in);
	for( int j = 0; j < nmodels; ++j ) {
		mdist.push_back( read<double>(in) );
		int length = (int)read<double>(in);
		nw_model_parameters param;
		in >> param;
		set_parameters( j, param, length );
	}
	set_region_dist( mdist );
	PHMM.optimize();
	dirty = false;
}

max_stars_repo_name: akhudek/feast

max_issues_repo_name: akhudek/feast

max_forks_repo_name: akhudek/feast


text:
/**
 * GeoDa TM, Copyright (C) 2011-2015 by Luc Anselin - all rights reserved
 *
 * This file is part of GeoDa.
 *
 * GeoDa 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.
 *
 * GeoDa 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/>.
 */

#include <algorithm> // std::sort
#include <cfloat>
#include <iomanip>
#include <iostream>
#include <limits>
#include <math.h>
#include <sstream>
#include <boost/foreach.hpp>
#include <wx/dcclient.h>
#include <wx/dcmemory.h>
#include <wx/msgdlg.h>
#include <wx/xrc/xmlres.h>
#include "../DialogTools/HistIntervalDlg.h"
#include "../GdaConst.h"
#include "../GeneralWxUtils.h"
#include "../GenGeomAlgs.h"
#include "../logger.h"
#include "../GeoDa.h"
#include "../Project.h"
#include "../ShapeOperations/ShapeUtils.h"
#include "SimpleBinsHistCanvas.h"

IMPLEMENT_CLASS(SimpleBinsHistCanvas, TemplateCanvas)
BEGIN_EVENT_TABLE(SimpleBinsHistCanvas, TemplateCanvas)
	EVT_PAINT(TemplateCanvas::OnPaint)
	EVT_ERASE_BACKGROUND(TemplateCanvas::OnEraseBackground)
	EVT_MOUSE_EVENTS(TemplateCanvas::OnMouseEvent)
	EVT_MOUSE_CAPTURE_LOST(TemplateCanvas::OnMouseCaptureLostEvent)
END_EVENT_TABLE()

const int SimpleBinsHistCanvas::MAX_INTERVALS = 200;
const double SimpleBinsHistCanvas::left_pad_const = 0;
const double SimpleBinsHistCanvas::right_pad_const = 0;
const double SimpleBinsHistCanvas::interval_width_const = 10;
const double SimpleBinsHistCanvas::interval_gap_const = 0;

SimpleBinsHistCanvas::SimpleBinsHistCanvas(wxWindow *parent,
									TemplateFrame* t_frame,
									Project* project,
									HLStateInt* hl_state_int,
									SimpleBinsHistCanvasCbInt* sbh_canv_cb_,
									const std::vector<SimpleBin>& hist_bins_,
									const wxString& Xname_,
									double Xmin_, double Xmax_,
									const wxString& right_click_menu_id_,
									bool show_axes_,
									const wxPoint& pos,
									const wxSize& size)
: TemplateCanvas(parent, t_frame, project, hl_state_int,
								 pos, size, false, true),
sbh_canv_cb(sbh_canv_cb_),
hist_bins(hist_bins_), Xname(Xname_), Xmin(Xmin_), Xmax(Xmax_),
x_axis(0), y_axis(0),
right_click_menu_id(right_click_menu_id_),
show_axes(show_axes_)
{
	using namespace Shapefile;
	LOG_MSG("Entering SimpleBinsHistCanvas::SimpleBinsHistCanvas");

	hist_bins_max = 0;
	for (size_t i=0; i<hist_bins.size(); ++i) {
		if (hist_bins[i].count > hist_bins_max) {
			hist_bins_max = hist_bins[i].count;
		}
	}

	highlight_color = GdaConst::highlight_color;
    last_scale_trans.SetFixedAspectRatio(false);
	use_category_brushes = false;
	selectable_shps_type = rectangles;

	PopulateCanvas();

	highlight_state->registerObserver(this);
	SetBackgroundStyle(wxBG_STYLE_CUSTOM);  // default style
	LOG_MSG("Exiting SimpleBinsHistCanvas::SimpleBinsHistCanvas");
}

SimpleBinsHistCanvas::~SimpleBinsHistCanvas()
{
	LOG_MSG("Entering SimpleBinsHistCanvas::~SimpleBinsHistCanvas");
	highlight_state->removeObserver(this);
	LOG_MSG("Exiting SimpleBinsHistCanvas::~SimpleBinsHistCanvas");
}

void SimpleBinsHistCanvas::DisplayRightClickMenu(const wxPoint& pos)
{
	LOG_MSG("Entering SimpleBinsHistCanvas::DisplayRightClickMenu");
	if (right_click_menu_id.IsEmpty()) return;
	// Workaround for right-click not changing window focus in OSX / wxW 3.0
	wxActivateEvent ae(wxEVT_NULL, true, 0, wxActivateEvent::Reason_Mouse);
	template_frame->OnActivate(ae);

	wxMenu* optMenu;
	optMenu = wxXmlResource::Get()->LoadMenu(right_click_menu_id);

	template_frame->UpdateContextMenuItems(optMenu);
	template_frame->PopupMenu(optMenu, pos + GetPosition());
	template_frame->UpdateOptionMenuItems();
	LOG_MSG("Exiting SimpleBinsHistCanvas::DisplayRightClickMenu");
}

/** Override of TemplateCanvas method. */
void SimpleBinsHistCanvas::update(HLStateInt* o)
{
    ResetBrushing();
	layer1_valid = false;
	Refresh();
}

wxString SimpleBinsHistCanvas::GetCanvasTitle()
{
	wxString s("Histogram");
	s << " - " << Xname;
	return s;
}

void SimpleBinsHistCanvas::TimeSyncVariableToggle(int var_index)
{
	LOG_MSG("In SimpleBinsHistCanvas::TimeSyncVariableToggle");
	invalidateBms();
	PopulateCanvas();
	Refresh();
}

void SimpleBinsHistCanvas::FixedScaleVariableToggle(int var_index)
{
	LOG_MSG("In SimpleBinsHistCanvas::FixedScaleVariableToggle");
	invalidateBms();
	PopulateCanvas();
	Refresh();
}

void SimpleBinsHistCanvas::ShowAxes(bool show_axes_s)
{
	show_axes = show_axes_s;
	invalidateBms();
	PopulateCanvas();
	Refresh();
}

void SimpleBinsHistCanvas::PopulateCanvas()
{
	BOOST_FOREACH( GdaShape* shp, background_shps ) { delete shp; }
	background_shps.clear();
	BOOST_FOREACH( GdaShape* shp, selectable_shps ) { delete shp; }
	selectable_shps.clear();
	BOOST_FOREACH( GdaShape* shp, foreground_shps ) { delete shp; }
	foreground_shps.clear();

	double num_ivals_d = (double) hist_bins.size();
	double x_min = 0;
	double x_max = left_pad_const + right_pad_const
	+ interval_width_const * num_ivals_d +
	+ interval_gap_const * (num_ivals_d-1);

	// orig_x_pos is the center of each histogram bar
	std::vector<double> orig_x_pos(hist_bins.size());
	for (size_t i=0; i<hist_bins.size(); ++i) {
		orig_x_pos[i] = left_pad_const + interval_width_const/2.0
		+ i * (interval_width_const + interval_gap_const);
	}

	double y_max = hist_bins_max;

    last_scale_trans.SetData(x_min, 0, x_max, y_max);
	if (show_axes) {
		axis_scale_y = AxisScale(0, y_max, 5, axis_display_precision);
		y_max = axis_scale_y.scale_max;
		y_axis = new GdaAxis("Frequency", axis_scale_y,
							 wxRealPoint(0,0), wxRealPoint(0, y_max),
							 -9, 0);
		foreground_shps.push_back(y_axis);

		axis_scale_x = AxisScale(0, Xmax, 5, axis_display_precision);
		//shps_orig_xmax = axis_scale_x.scale_max;
		axis_scale_x.data_min = Xmin;
		axis_scale_x.data_max = Xmax;
		axis_scale_x.scale_min = axis_scale_x.data_min;
		axis_scale_x.scale_max = axis_scale_x.data_max;
		double range = axis_scale_x.scale_max - axis_scale_x.scale_min;

		axis_scale_x.scale_range = range;
		axis_scale_x.p = floor(log10(range));
		axis_scale_x.ticks = hist_bins.size()+1;
		axis_scale_x.tics.resize(axis_scale_x.ticks);
		axis_scale_x.tics_str.resize(axis_scale_x.ticks);
		axis_scale_x.tics_str_show.resize(axis_scale_x.tics_str.size());
		for (int i=0; i<axis_scale_x.ticks; i++) {
			axis_scale_x.tics[i] =
			axis_scale_x.data_min +
			range*((double) i)/((double) axis_scale_x.ticks-1);
			std::ostringstream ss;
			ss << std::fixed << std::setprecision(3) << axis_scale_x.tics[i];
			axis_scale_x.tics_str[i] = ss.str();
			axis_scale_x.tics_str_show[i] = false;
		}
		int tick_freq = ceil(num_ivals_d/10.0);
		for (int i=0; i<axis_scale_x.ticks; i++) {
			if (i % tick_freq == 0) {
				axis_scale_x.tics_str_show[i] = true;
			}
		}
		axis_scale_x.tic_inc = axis_scale_x.tics[1]-axis_scale_x.tics[0];
		x_axis = new GdaAxis(Xname, axis_scale_x, wxRealPoint(0,0),
							 wxRealPoint(x_max, 0), 0, 9);
		foreground_shps.push_back(x_axis);
	}

	GdaShape* s = 0;
	int table_w=0, table_h=0;

	int virtual_screen_marg_top = 5; //25;
	int virtual_screen_marg_bottom = 5; //25;
	int virtual_screen_marg_left = 5; //25;
	int virtual_screen_marg_right = 5; //25;

	if (show_axes) {
		virtual_screen_marg_bottom += 32;
		virtual_screen_marg_left += 35;
	}
    last_scale_trans.SetMargin(virtual_screen_marg_top,
                               virtual_screen_marg_bottom,
                               virtual_screen_marg_left,
                               virtual_screen_marg_right);

	selectable_shps.resize(hist_bins.size());
	for (size_t i=0; i<hist_bins.size(); ++i) {
		double x0 = orig_x_pos[i] - interval_width_const/2.0;
		double x1 = orig_x_pos[i] + interval_width_const/2.0;
		double y0 = 0;
		double y1 = hist_bins[i].count;
		selectable_shps[i] = new GdaRectangle(wxRealPoint(x0, 0),
											  wxRealPoint(x1, y1));
		int sz = GdaConst::qualitative_colors.size();
		selectable_shps[i]->setPen(GdaConst::qualitative_colors[i%sz]);
		selectable_shps[i]->setBrush(GdaConst::qualitative_colors[i%sz]);
	}

	ResizeSelectableShps();
}

void SimpleBinsHistCanvas::UpdateStatusBar()
{
	if (sbh_canv_cb) {
		sbh_canv_cb->notifyNewHistHover(hover_obs, total_hover_obs);
	}
}

ext: cpp

max_stars_repo_name: chenyoujie/GeoDa

max_issues_repo_name: chenyoujie/GeoDa

max_forks_repo_name: chenyoujie/GeoDa


text:
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
Copyright (C) 2005, 2007, 2009, 2014 Klaus Spanderen
Copyright (C) 2015 CompatibL

This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/

QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license.  You should have received a
copy of the license along with this program; if not, please email
<quantlib-dev@lists.sf.net>. The license is also available online at
<http://quantlib.org/license.shtml>.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#ifndef cl_adjoint_heston_process_impl_hpp
#define cl_adjoint_heston_process_impl_hpp
#pragma once

#include "adjointhestonprocesstest.hpp"
#include "utilities.hpp"
#include "adjointtestutilities.hpp"
#include "adjointtestbase.hpp"
#include <ql/quotes/simplequote.hpp>
#include <ql/math/modifiedbessel.hpp>
#include <ql/processes/hestonprocess.hpp>
#include <ql/time/daycounters/actualactual.hpp>
#include <ql/time/schedule.hpp>
#include <boost/make_shared.hpp>

using namespace QuantLib;
using namespace boost::unit_test_framework;

#define OUTPUT_FOLDER_NAME "AdjointHestonProcess"

/*!
AdjointHestonProcessTest is tested using adjoint differentiation for the class hestonprocess.

Heston model describes the evolution of the volatility of an underlying asset.
It is a stochastic volatility model: such a model assumes that the
volatility of the asset is neither constant, nor even deterministic,
but is a random process.

This class describes the square root stochastic volatility
process governed by
\f[
\begin{array}{rcl}
dS(t, S)  &=& \mu S dt + \sqrt{v} S dW_1 \\
dv(t, S)  &=& \kappa (\theta - v) dt + \sigma \sqrt{v} dW_2 \\
dW_1 dW_2 &=& \rho dt
\end{array}
\f]

Tested functions are expectation, evolve, phi, drift, stdDeviation, covariance.

Phi is the continuous version of a characteristic function
for the exact sampling of the Heston process, s. page 8, formula 13,
M. Broadie, O. Kaya, Exact Simulation of Stochastic Volatility and
other Affine Jump Diffusion Processes
http://finmath.stanford.edu/seminars/documents/Broadie.pdf
Output of function Phi is value of haracteristic function.
Inputs are process, value of variable, start and finish variance distibution, and delta time.
Phi is obtained from fourier transform, so Phi has complex derivative with respect to input variable.

Output of fuction expectation is Expectation
\f$ E(\mathrm{x}_{t_0 + \Delta t}
| \mathrm{x}_{t_0} = \mathrm{x}_0) \f$
of the process after a time interval \f$ \Delta t \f$
according to the given discretization.
Inputs are time, start value, and delta time.
Expectation is differentiated with respect to delta time.

Output of fuction evolve is Evolve,
which returns the asset value after a time interval \f$ \Delta t
\f$ according to the given discretization. By default, it
returns
\f[
E(x_0,t_0,\Delta t) + S(x_0,t_0,\Delta t) \cdot \Delta w
\f]
where \f$ E \f$ is the expectation and \f$ S \f$ the
standard deviation.
Inputs are time, start value, and delta time.
Evolve is differentiated with respect to delta time.

Output of fuction drift is Drift (the change of the average value)
Inputs are time, array of average values.
Drift is differentiated with respect to average values.

Output of fuction stdDeviation is Standart Deviation.
Inputs are time, array of average values, and delta time.
StdDeviation is differentiated with respect to delta time.

Output of fuction covariance is covariances matrix
Inputs are time, array of average values, and delta time.
Covariance is differentiated with respect to delta time.
*/

namespace
{
    enum
    {
#if defined CL_GRAPH_GEN
        // Number of points for dependency plots.
        pointNo = 100,
        // Number of points for performance plot.
        iterNo = 30,
        // Step for portfolio size for performance testing .
        step = 1,
#else
        // Number of points for dependency plots.
        pointNo = 0,
        // Number of points for performance plot.
        iterNo = 0,
        // Step for portfolio size for performance testing .
        step = 1,
#endif
        // Defines performance accuracy. Its value is a minimum number
        // of calling of O(1) complexity methods per one performance test.
        iterNumFactor = 1000,
    };

    // copy of function Phi from hestonprocess.hpp
    // This is the continuous version of a characteristic function
    // for the exact sampling of the Heston process, s. page 8, formula 13,
    // M. Broadie, O. Kaya, Exact Simulation of Stochastic Volatility and
    // other Affine Jump Diffusion Processes
    // http://finmath.stanford.edu/seminars/documents/Broadie.pdf
    std::complex<Real> Phi1(const HestonProcess& process,
        const std::complex<Real>& a,
        Real nu_0, Real nu_t, Time dt) {
        const Real theta = process.theta();
        const Real kappa = process.kappa();
        const Real sigma = process.sigma();

        const Volatility sigma2 = sigma*sigma;
        const std::complex<Real> ga = std::sqrt(
            kappa*kappa - 2 * sigma2*a*std::complex<Real>(0.0, 1.0));
        const Real d = 4 * theta*kappa / sigma2;

        const Real nu = 0.5*d - 1;
        const std::complex<Real> z
            = ga*std::exp(-0.5*ga*dt) / (1.0 - std::exp(-ga*dt));
        const std::complex<Real> log_z
            = -0.5*ga*dt + std::log(ga / (1.0 - std::exp(-ga*dt)));

        const std::complex<Real> alpha
            = 4.0*ga*std::exp(-0.5*ga*dt) / (sigma2*(1.0 - std::exp(-ga*dt)));

        const std::complex<Real> beta = 4.0*kappa*std::exp(-0.5*kappa*dt)
            / (sigma2*(1.0 - std::exp(-kappa*dt)));

        return ga*std::exp(-0.5*(ga - kappa)*dt)*(1 - std::exp(-kappa*dt))
            / (kappa*(1.0 - std::exp(-ga*dt)))
            *std::exp((nu_0 + nu_t) / sigma2 * (
            kappa*(1.0 + std::exp(-kappa*dt)) / (1.0 - std::exp(-kappa*dt))
            - ga*(1.0 + std::exp(-ga*dt)) / (1.0 - std::exp(-ga*dt))))
            *std::exp(nu*log_z) / std::pow(z, nu)
            *((nu_t > 1e-8)
            ? modifiedBesselFunction_i(
            nu, std::sqrt(nu_0*nu_t)*alpha)
            / modifiedBesselFunction_i(
            nu, std::sqrt(nu_0*nu_t)*beta)
            : std::pow(alpha / beta, nu)
            );
    }

    struct CommonVars0
    {
        CommonVars0()
        : dayCounter_()
        {
            dayCounter_ = ActualActual();
        }
        DayCounter dayCounter_;

    };

    struct CommonVars : public CommonVars0
    {
        explicit CommonVars()
        : CommonVars0()
        , s0_(boost::make_shared<SimpleQuote>(SimpleQuote(1.05)))
        , process_(Handle<YieldTermStructure>(flatRate(0.04, dayCounter_)), Handle<YieldTermStructure>(flatRate(0.5, dayCounter_)), s0_, 0.3, 1.16, 0.2, 0.8, 0.8,
        HestonProcess::QuadraticExponentialMartingale)
        {
        }

        Handle<Quote> s0_;
        HestonProcess process_;
    };

    enum VariableName { expectation, evolve, phi, drift, stdDeviation, covariance };

    std::string toString(VariableName name)
    {
        switch (name)
        {
        case expectation:
            return "Expectation";
        case evolve:
            return "Evolve";
        case phi:
            return "Phi";
        case drift:
            return "Drift";
        case stdDeviation:
            return "stdDeviation";
        case covariance:
            return "covariance";
        default:
            return "";
        }
    }

    struct HestonProcessTestData : public CommonVars
    {
        explicit HestonProcessTestData(VariableName type)
        : CommonVars()
        , type_(type)
        , t(0.5)
        , x(2, 0.3, 0.2)
        , dw(2, 0.1, 0.4)
        {}

        // Calculate the an expectation in case of test of function Expectation
        // and an evolve in case of test of function Evolve.
        Real calculateRealFromDeltaTime(Real deltaTimes)
        {
            if (type_ == expectation)
            {
                return process_.expectation(t, x, deltaTimes)[0];
            }

            if (type_ == evolve)
            {
                return process_.evolve(t, x, deltaTimes, dw)[0];
            }

            return Real(0.0);
        }

        Matrix calculateMatrixFromDeltaTime(Real deltaTimes)
        {
            if (type_ == stdDeviation)
            {
                return process_.stdDeviation(t, x, deltaTimes);
            }

            if (type_ == covariance)
            {
                return process_.covariance(t, x, deltaTimes);
            }

            return Matrix();
        }

        VariableName type_;
        Real t;
        Array x;
        Array dw;
    };

    // Test of functions expectation and evolve.
    // Testing derivative with respect to delta time.
    struct HestonProcessRespectToDeltaTimeTest
        : public cl::AdjointTest<HestonProcessRespectToDeltaTimeTest>
    {
        explicit HestonProcessRespectToDeltaTimeTest(HestonProcessTestData* data, Size size, cl::tape_empty_test_output* logger = nullptr, double shift = 1e-10)
        : AdjointTest()
        , data_(data)
        , size_(size)
        , shift_(shift)
        , deltaTimes_(size)
        , Variables_(size)
        , SumVariables_(1)
        {
            setLogger(logger);
            for (int i = 0; i < size_; i++)
            {
                deltaTimes_[i] = (i + 1) * 0.2 / double(size);
            }
        }

        Size indepVarNumber() { return size_; }

        Size depVarNumber() { return 1; }

        Size minPerfIteration() { return iterNumFactor; }

        double absTol() const { return 1e-5; }

        void calcAnalytical()
        {
            analyticalResults_.resize(size_);
            for (int i = 0; i < size_; i++)
            {
                double deltaTimeUp = double(deltaTimes_[i]) + shift_;
                double deltaTimeDown = double(deltaTimes_[i]) - shift_;
                double varUp = double(data_->calculateRealFromDeltaTime(deltaTimeUp));
                double varDown = double(data_->calculateRealFromDeltaTime(deltaTimeDown));
                analyticalResults_[i] = (varUp - varDown) / 2 / shift_;
            }
        }

        void calculateVariables()
        {
            SumVariables_[0] = 0.0;
            for (int i = 0; i < size_; i++)
            {
                Variables_[i] = data_->calculateRealFromDeltaTime(deltaTimes_[i]);
                SumVariables_[0] += Variables_[i];
            }
        }

        void recordTape()
        {
            Independent(deltaTimes_);
            calculateVariables();
            f_ = std::make_unique<cl::tape_function<double>>(deltaTimes_, SumVariables_);
        }

        HestonProcessTestData* data_;
        Size size_;
        double shift_;
        std::vector<Real> deltaTimes_;
        std::vector<Real> Variables_;
        std::vector<Real> SumVariables_;
    };

    // Test of matrix functions stdDeviation and covariance.
    struct HestonProcessMatrixFunctionsTest :
        public cl::AdjointTest<HestonProcessMatrixFunctionsTest>
    {
        static const CalcMethod default_method = other;

        explicit HestonProcessMatrixFunctionsTest(HestonProcessTestData* data, cl::tape_empty_test_output* logger = nullptr, double shift = 1e-10)
        : AdjointTest()
        , data_(data)
        , shift_(shift)
        , deltaTimes_(1, 0.1)
        , size_(4)
        {}

        Size indepVarNumber() { return 1; }

        Size depVarNumber() { return size_; }

        Size minPerfIteration() { return iterNumFactor; }

        double absTol() const { return 1e-6; }

        void calculateMatrixRes()
        {
            Matrix StdDeviation = data_->calculateMatrixFromDeltaTime(deltaTimes_[0]);
            Results_.assign(StdDeviation.begin(), StdDeviation.end());
        }

        void recordTape()
        {
            Independent(deltaTimes_);
            calculateMatrixRes();
            f_ = std::make_unique<cl::tape_function<double>>(deltaTimes_, Results_);
        }

        void calcAdjoint()
        {
            int num = Results_.size();
            std::vector<double> w(1, 1);
            adjointResults_ = f_->Forward(1, w);
        }

        void calcAnalitical()
        {
            Time dtUp = deltaTimes_[0] + shift_;
            Time dtDown = deltaTimes_[0] - shift_;

            Matrix ResUp = data_->calculateMatrixFromDeltaTime(dtUp);
            Matrix ResDown = data_->calculateMatrixFromDeltaTime(dtDown);

            for (int i = 0; i < ResUp.columns(); i++)
            {
                for (int j = 0; j < ResUp.rows(); j++)
                {
                    double analyticalRes = double((ResUp[i][j] - ResDown[i][j]) / 2 / shift_);
                    analyticalResults_.push_back(analyticalRes);
                }
            }
        }

        HestonProcessTestData* data_;
        double shift_;
        std::vector<Real> deltaTimes_;
        std::vector<Real> Results_;
        Size size_;
    };

    struct HestonProcessDriftTest : public cl::AdjointTest<HestonProcessDriftTest>
    {
        static const CalcMethod default_method = other;

        explicit HestonProcessDriftTest(CommonVars* data, cl::tape_empty_test_output* logger = nullptr, double shift = 1e-10)
            : AdjointTest()
            , data_(data)
            , shift_(shift)
            , size_(2)
            , values_(size_)
            , x_(size_)
            , drift_(2)
            , t_(0.5)
            , indepDouble_(2)
        {
            for (int i = 0; i < size_; i++)
            {
                indepDouble_[i] = 0.3 + 0.2 * i;
                values_[i] = indepDouble_[i];
            }
        }

        Size indepVarNumber() { return size_; }

        Size depVarNumber() { return size_; }

        Size minPerfIteration() { return iterNumFactor; }

        double absTol() const { return 1e-6; }

        void calculateDrift()
        {
            for (int i = 0; i < size_; i++)
            {
                x_[i] = values_[i];
            }
            Array drift = data_->process_.drift(t_, x_);
            for (int i = 0; i < size_; i++)
            {
                drift_[i] = drift[i];
            }
        }

        void recordTape()
        {
            cl::Independent(values_);
            calculateDrift();
            f_ = std::make_unique<cl::tape_function<double>>(values_, drift_);
        }

        void calcAdjoint()
        {
            adjointResults_ = f_->Jacobian(indepDouble_);
        }

        void calcAnalitical()
        {
            analyticalResults_.resize(size_ * size_);
            Array xUp = x_;
            Array xDown = x_;
            for (int i = 0; i < size_; i++)
            {
                xUp[i] += shift_;
                Array DriftUp = data_->process_.drift(t_, xUp);
                xUp[i] -= shift_;

                xDown[i] -= shift_;
                Array DriftDown = data_->process_.drift(t_, xDown);
                xDown[i] += shift_;

                for (int j = 0; j < size_; j++)
                {
                    analyticalResults_[i + size_ * j] = (DriftUp[j] - DriftDown[j]) / 2 / shift_;
                }
            }
        }

        CommonVars* data_;
        double shift_;
        Size size_;
        std::vector<Real> values_;
        Array x_;
        std::vector<Real> drift_;
        Time t_;
        std::vector<double> indepDouble_;
    };

    // Struct for plots recording.
    struct DeltaTimeResult
    {
        static std::deque<std::string > get_columns()
        {
            static std::deque<std::string > columns =
            {
                "Delta time", ""
            };

            return columns;
        }

        template <typename stream_type>
        friend inline stream_type&
            operator << (stream_type& stm, DeltaTimeResult& v)
        {
                stm << v.deltaTime_
                    << ";" << v.result_
                    << std::endl;
                return stm;
            }

        Real deltaTime_;
        Real result_;
    };

    struct TestDataRespectDeltaTime
        : public HestonProcessTestData
    {
        explicit TestDataRespectDeltaTime(VariableName name)
        : HestonProcessTestData(name)
        , outPerform_title_("Swaption " + toString(name) + " differentiation performance with respect to delta time")
        , outAdjoint_title_("Swaption " + toString(name) + " adjoint differentiation with respect to fixed delta time")
        , out_title_(toString(name) + " dependence on delta time")
        , out_ylabel_(toString(name))
        , outPerform_(OUTPUT_FOLDER_NAME "//" + toString(name), {
            { "title", outPerform_title_ }
            , { "not_clear", "Not" }
            , { "ylabel", "Time (s)" }
            , { "xlabel", "Number of delta times" }
            , { "line_box_width", "-5" }
            , { "cleanlog", "true" }
            , { "smooth", "default" }
        })
            , outAdjoint_(OUTPUT_FOLDER_NAME "//" + toString(name), {
                { "title", outAdjoint_title_ }
            , { "filename", "Adjoint" }
            , { "not_clear", "Not" }
            , { "ylabel", "Time (s)" }
            , { "xlabel", "Number of delta times" }
            , { "cleanlog", "false" }
            , { "smooth", "default" }
        })
            , outSize_(OUTPUT_FOLDER_NAME "//" + toString(name), {
                { "title", "Tape size dependence on number of delta time" }
            , { "filename", "TapeSize" }
            , { "not_clear", "Not" }
            , { "ylabel", "Memory (MB)" }
            , { "cleanlog", "false" }
        })

            , out_(OUTPUT_FOLDER_NAME "//" + toString(name) + "//output", {
                { "filename", "DeltaTimeDependence" }
            , { "not_clear", "Not" }
            , { "title", out_title_ }
            , { "ylabel", out_ylabel_ }
            , { "xlabel", "delta time" }
            , { "cleanlog", "false" }
        })
        { }

        bool makeOutput()
        {
            bool ok = true;

            if (pointNo > 0)
            {
                ok &= recordDependencePlot();
            }
            ok &= cl::recordPerformance(*this, iterNo, step);

            return ok;
        }

        std::shared_ptr<HestonProcessRespectToDeltaTimeTest> getTest(size_t size)
        {
            return std::make_shared<HestonProcessRespectToDeltaTimeTest>(this, size, &outPerform_);
        }

        // Makes plots for derivative of blackVariance dependence.
        bool recordDependencePlot()
        {
            std::vector<DeltaTimeResult> outData(pointNo);
            auto test = getTest(pointNo);
            bool ok = test->testReverse();
            for (Size i = 0; i < pointNo; i++)
            {
                outData[i] = { test->deltaTimes_[i], test->Variables_[i] };
            }

            out_ << outData;
            return ok;

            return true;
        }

        std::string outPerform_title_;
        std::string outAdjoint_title_;
        std::string out_title_;
        std::string out_ylabel_;
        cl::tape_empty_test_output outPerform_;
        cl::tape_empty_test_output outAdjoint_;
        cl::tape_empty_test_output outSize_;
        cl::tape_empty_test_output out_;
    };

}

#endif

ext: hpp

max_stars_repo_name: fduffy/QuantLibAdjoint

max_issues_repo_name: fduffy/QuantLibAdjoint

max_forks_repo_name: fduffy/QuantLibAdjoint


text:
#pragma once
#include <boost/utility/string_view.hpp>
#include "sobject_cloud_vm_base.h"
#pragma GCC visibility push(default)
template<typename Device,typename DEBUG_LOGGER>
struct device_capabilities_impl{
     //assumes copy constructable logger
     static DEBUG_LOGGER logger() {return DEBUG_LOGGER();}
     using VM=typename Device::VM;
     using ScreenManager=typename Device::ScreenManager;
     static std::string getModel(){return Device::getModel();}
     static std::string getCountry(){return Device::getCountry();}
     static std::string getLanguage(){return Device::getLanguage();}
     static std::string getAuthToken(std::string name){return Device::getAuthToken(name);}
     static std::tuple<std::string,std::string,std::string> get_auth_data(){return Device::get_auth_data();}
     static std::string getCloudId(){ return Device::getCloudId(); }
     static std::string getCloudVersion(){ return Device::getCloudVersion(); }
     static std::string expType(){ return Device::expType(); }
     static std::string get_proxy(){return Device::get_proxy();}
     static bool getSmartTasksStatus(){return Device::getSmartTasksStatus();}
     static std::string getRootUrl(){return  Device::getRootUrl();}
     static std::chrono::milliseconds timeout(){return Device::timeout();}
     static std::unique_ptr<VM> createVM(){ return Device::createVM();}
     static std::unique_ptr<ScreenManager> createScreenManager(){return Device::createScreenManager();}
     template<typename Handler>
     static auto make_app_handler(Handler&& handler)
     {
      return  Device::make_app_handler(std::forward<Handler>(handler));
     }

};
#pragma GCC visibility pop


text:
#include <armadillo>
#include <iostream>
#include "gtest/gtest.h"

#include "svm_helpers.h"
#include "svm_client.h"
#include "utils/logger.h"

namespace {

int log_level = LogLevel::DEBUG_LEVEL;

// The fixture
class SVMClientTest: public ::testing::Test {

protected:

    SVMClientTest() {
        std::cout << "Creating SVMConfiguration..." << std::endl << std::flush;
        svm_config = SVMConfiguration();
        std::cout << "Setting logger..." << std::endl << std::flush;
        svm_config.log.verbosity = log_level;

        // Exemplary sparse matrix
        sparse_matrix_csc_01_row_indices =
            helper_sparse_matrix_csc_01_row_indices();
        sparse_matrix_csc_01_column_pointers =
            helper_sparse_matrix_csc_01_column_pointers();
        sparse_matrix_csc_01_values =
            helper_sparse_matrix_csc_01_values();
        sparse_matrix_csc_01_nrow =
            helper_sparse_matrix_csc_01_nrow();
        sparse_matrix_csc_01_ncol =
            helper_sparse_matrix_csc_01_ncol();
        sparse_matrix_csc_01_sp_mat =
            helper_sparse_matrix_csc_01_sp_mat();
        sparse_matrix_csc_01_learning_target =
            helper_sparse_matrix_csc_01_learning_target();

        std::cout << "Starting test..." << std::endl << std::flush;
    }

    virtual ~SVMClientTest() {}

    /// Called immediately after the constructor (righ before each test)
    virtual void SetUp() {}

    /// Called immediately after each test (right before the destructor)
    virtual void TearDown() {}

    /* Objects declared here can be used by all tests in the test case */

    SVMConfiguration svm_config;

	arma::uvec sparse_matrix_csc_01_row_indices;
    arma::uvec sparse_matrix_csc_01_column_pointers;
    arma::vec sparse_matrix_csc_01_values;
    size_t sparse_matrix_csc_01_nrow;
    size_t sparse_matrix_csc_01_ncol;
    arma::sp_mat sparse_matrix_csc_01_sp_mat;
    arma::vec sparse_matrix_csc_01_learning_target;
};

} // end namespace

/* Fixture tests */

TEST_F(SVMClientTest, DISABLED_sparse_data_test) {
    std::cout << "Testing learning..." << std::endl << std::flush;
    svm_config.setSparseData(
        sparse_matrix_csc_01_row_indices,
        sparse_matrix_csc_01_column_pointers,
        sparse_matrix_csc_01_values,
        sparse_matrix_csc_01_nrow,
        sparse_matrix_csc_01_ncol,
        true
    );
    //svm_config.sparse_data = sparse_matrix_csc_01_sp_mat;
    svm_config.target = sparse_matrix_csc_01_learning_target;
    svm_config.setLibrary("libsvm");
    SVMClient *svm_client = new SVMClient(&svm_config);
    // XXX: Why access violation here?
    svm_client->train();

    std::cout << "Testing prediction..." << std::endl << std::flush;
    svm_client->sparse_predict(
        sparse_matrix_csc_01_row_indices,
        sparse_matrix_csc_01_column_pointers,
        sparse_matrix_csc_01_values,
        sparse_matrix_csc_01_nrow,
        sparse_matrix_csc_01_ncol
    );
    SVMConfiguration client_config = svm_client->getConfiguration();

    for (int i = 0; i < sparse_matrix_csc_01_learning_target.n_rows; ++i) {
        ASSERT_DOUBLE_EQ(
            client_config.result[i], sparse_matrix_csc_01_learning_target[i]);
    }
}


text:
  #include <pcl/io/pcd_io.h>
  #include <pcl/point_cloud.h>
  #include <pcl/correspondence.h>

  #include <pcl/features/normal_3d_omp.h>
  //#include <pcl/features/spin_image.h>
  #include <pcl/features/pfh.h>
  #include <pcl/features/shot_omp.h>
  #include <pcl/features/fpfh.h>
  #include <pcl/features/board.h>

  //#include <pcl/filters/uniform_sampling.h>
  #include <pcl/keypoints/uniform_sampling.h>
  #include <pcl/keypoints/iss_3d.h>

  #include <pcl/recognition/cg/hough_3d.h>
  #include <pcl/recognition/cg/geometric_consistency.h>

  #include <pcl/visualization/pcl_visualizer.h>

  #include <pcl/kdtree/kdtree_flann.h>
  #include <pcl/kdtree/impl/kdtree_flann.hpp>

  #include <pcl/common/transforms.h>
  #include <pcl/console/parse.h>

  // ensaio de acerto de resolução
  #include <iostream>

  #include <string.h>
  #include <stdlib.h>
  #include <boost/lexical_cast.hpp>

  #include <pcl/point_types.h>
  #include <pcl/filters/voxel_grid.h>

  #include <sstream>

  // tempo
  #include <ctime>

  typedef pcl::PointXYZ PointType;
  typedef pcl::Normal NormalType;
  typedef pcl::ReferenceFrame RFType;

  typedef pcl::PointWithViewpoint PointTypeViewPoint;

  //typedef pcl::Histogram<153> DescriptorType;

  std::string model_filename_;
  std::string scene_filename_;
  std::vector<std::string> model_filenames_;

  //Algorithm params
  bool show_keypoints_ (false);
  bool show_correspondences_ (false);
  bool show_normals_(false);
  bool use_cloud_resolution_ (false);
  bool use_hough_ (true);// default

  //Keypoints parameters
  bool use_uniformsampling(true);// default
  bool use_iss(false);
  bool use_harris3d(false);

  //Features parameters
  bool use_shot(true);//default
  bool use_pfh(false);
  bool use_usd(false);

  float model_ss_ (0.005f); //Amostragem do modelo para cálculo dos keyoints
  float scene_ss_ (0.005f); //Amostragem da cena para cálculo dos keyoints

  float rf_rad_ (0.015f);
  float descr_rad_ (0.04f); //Raio para o cálculo do descritor

  float cg_size_ (0.08f); //Cluster size
  float cg_thresh_ (40.0f);

  std::string model_directory_;

  // Multi-view data storage
  struct chave
  {
  std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f> > rototranslations;
  std::vector<pcl::Correspondences> clustered_corrs;
  };

  std::vector<chave> model_view;

  // Cronometro
  std::clock_t start;
  std::clock_t end;

  double diffclock(std::clock_t clock1,std::clock_t clock2)
  {
    double diffticks = clock2 - clock1;
    double diffms = diffticks / (CLOCKS_PER_SEC / 1000 );
    return diffms;
  }

  //***********************************************************************************************
  //                                  Função ajuda a configuração
  //***********************************************************************************************
  void
  showHelp (char *filename)
  {
    std::cout << std::endl;
    std::cout << "***************************************************************************" << std::endl;
    std::cout << "*                            Usage Guide                                  *" << std::endl;
    std::cout << "***************************************************************************" << std::endl << std::endl;
    std::cout << "Usage: " << filename << " model_filename.pcd scene_filename.pcd [Options]" << std::endl << std::endl;
    std::cout << "Options:" << std::endl;
    std::cout << "     -h:                     Show this help." << std::endl;
    std::cout << "     -k:                     Show used keypoints." << std::endl;
    std::cout << "     -c:                     Show used correspondences." << std::endl;
    std::cout << "     -n:                     Show used Normals." << std::endl;
    std::cout << "     -r:                     Compute the model cloud resolution and multiply" << std::endl;
    std::cout << "                             each radius given by that value." << std::endl;
    std::cout << "     --Keypoints: keypoints downsample algorithm used (default UniformSampling)." << std::endl;
    std::cout << "                  ISS (Intrinsic Shape Signatures)" << std::endl;
    //std::cout << "                  Harris3D                        " << std::endl;
    std::cout << "     --Features:  descritor features algorithm used (default SHOT)." << std::endl;
    //std::cout << "                  PFH (Point Features Histogram)" << std::endl;
    //std::cout << "                  ISS (Intrinsic Shape Signatures)" << std::endl;
    //std::cout << "                  USD (Unique Shape Context)" << std::endl;
    std::cout << "     --algorithm (Hough|GC): Clustering algorithm used (default Hough)." << std::endl;
    std::cout << "     --model_ss val:         Model uniform sampling radius (default 0.01)" << std::endl;
    std::cout << "     --scene_ss val:         Scene uniform sampling radius (default 0.03)" << std::endl;
    std::cout << "     --rf_rad val:           Reference frame radius (default 0.015)" << std::endl;
    std::cout << "     --descr_rad val:        Descriptor radius (default 0.02)" << std::endl;
    std::cout << "     --cg_size val:          Cluster size (default 0.01)" << std::endl;
    std::cout << "     --cg_thresh val:        Clustering threshold (default 5)" << std::endl << std::endl;
  }

//

  //***********************************************************************************************
  //                                  Função parametros de entrada
  //***********************************************************************************************
  void
  parseCommandLine (int argc, char *argv[])
  {
    //Show help
    if (pcl::console::find_switch (argc, argv, "-h"))
    {
      showHelp (argv[0]);
      exit (0);
    }

    //Model & scene filenames
    std::vector<int> filenames;
    filenames = pcl::console::parse_file_extension_argument (argc, argv, ".pcd");
    if (filenames.size () != 1)
    {
      std::cout << "Filenames missing.\n";
      showHelp (argv[0]);
      exit (-1);
      }

    //model_filename_ = argv[filenames[0]];
    scene_filename_ = argv[filenames[0]];

    //Program visualization behavior

    if (pcl::console::find_switch (argc, argv, "-k"))
    {
      show_keypoints_ = true;
    }
    if (pcl::console::find_switch (argc, argv, "-c"))
    {
      show_correspondences_ = true;
    }
    if (pcl::console::find_switch (argc, argv, "-r"))
    {
      use_cloud_resolution_ = true;
    }

    if (pcl::console::find_switch (argc, argv, "-n"))
    {
      show_normals_ = true;
    }

    // selection of algorithm for clustering
    std::string used_algorithm;
    if (pcl::console::parse_argument (argc, argv, "--algorithm", used_algorithm) != -1)
    {
      if (used_algorithm.compare ("Hough") == 0)
      {
        use_hough_ = true;
      }
      else if (used_algorithm.compare ("GC") == 0)
      {
        use_hough_ = false;
      }
      else
      {
        std::cout << "Wrong algorithm name.\n";
        showHelp (argv[0]);
        exit (-1);
      }
    }

    // Selection of algorithm for Keypoint extration
    std::string used_keypoints;
    if (pcl::console::parse_argument (argc, argv, "--keypoint", used_keypoints) != -1)
    {
      if (used_keypoints.compare ("UniformSample") == 0)
      {

      }
      else if (used_keypoints.compare ("ISS") == 0)
      {
        use_uniformsampling = false;
        use_iss = true;
      }
      else if (used_keypoints.compare ("Harris3D") == 0)
      {
        use_uniformsampling = false;
        use_harris3d = true;
      }
      else
      {
        std::cout << "Wrong Keypoint name.\n";
        showHelp (argv[0]);
        exit (-1);
      }
    }

start = std::clock();                                    // Debug time

  ///std::string model_directory_;
    if (pcl::console::parse_argument (argc, argv, "--md", model_directory_) != -1)
    {
      std::cout << "model_directory_=" << model_directory_ << std::endl;
  //TODO criar a list de strings usando uma analise da diretoria (dir)

  //model_filenames_.push_back("../models/2.pcd");
 // model_filenames_.push_back("../models/8.pcd");
 // model_filenames_.push_back("../models/9.pcd");

 // model_filenames_.push_back("../models/model062.pcd");
  //model_filenames_.push_back("../models/16.pcd");

      model_filenames_.push_back("../models/0.pcd");
      model_filenames_.push_back("../models/1.pcd");
      model_filenames_.push_back("../models/2.pcd");
      model_filenames_.push_back("../models/3.pcd");
      model_filenames_.push_back("../models/4.pcd");
      model_filenames_.push_back("../models/5.pcd");
      model_filenames_.push_back("../models/6.pcd");
      model_filenames_.push_back("../models/7.pcd");
      model_filenames_.push_back("../models/8.pcd");
      model_filenames_.push_back("../models/9.pcd");
      model_filenames_.push_back("../models/10.pcd");
      model_filenames_.push_back("../models/11.pcd");
      model_filenames_.push_back("../models/12.pcd");
      model_filenames_.push_back("../models/13.pcd");
      model_filenames_.push_back("../models/14.pcd");
      model_filenames_.push_back("../models/15.pcd");
      model_filenames_.push_back("../models/16.pcd");
      model_filenames_.push_back("../models/17.pcd");
      model_filenames_.push_back("../models/18.pcd");
      model_filenames_.push_back("../models/19.pcd");
      model_filenames_.push_back("../models/20.pcd");
      model_filenames_.push_back("../models/21.pcd");
      model_filenames_.push_back("../models/22.pcd");
      model_filenames_.push_back("../models/23.pcd");
      model_filenames_.push_back("../models/24.pcd");
      model_filenames_.push_back("../models/25.pcd");
      model_filenames_.push_back("../models/26.pcd");
      model_filenames_.push_back("../models/27.pcd");
      model_filenames_.push_back("../models/28.pcd");
      model_filenames_.push_back("../models/29.pcd");
      model_filenames_.push_back("../models/30.pcd");
      model_filenames_.push_back("../models/31.pcd");
      model_filenames_.push_back("../models/32.pcd");
      model_filenames_.push_back("../models/33.pcd");
      model_filenames_.push_back("../models/34.pcd");
      model_filenames_.push_back("../models/35.pcd");
      model_filenames_.push_back("../models/36.pcd");
      model_filenames_.push_back("../models/37.pcd");
      model_filenames_.push_back("../models/38.pcd");
      model_filenames_.push_back("../models/39.pcd");
      model_filenames_.push_back("../models/40.pcd");

    }
 end = std::clock();                                       // Debug time
 std::cout << "Time loading Files:  " << diffclock(start, end) << " ms" << endl << endl;

    // Selection of features
  /*  std::string used_features;
    if (pcl::console::parse_argument (argc, argv, "--features", used_features) != -1)
    {
      if (used_algorithm.compare ("SHOT") == 0)
      {

      }
      else if (used_algorithm.compare ("PFH") == 0)
      {
          use_shot = false;
          use_pfh = true;
      }
      else if (used_algorithm.compare ("USD") == 0)
      {
          use_shot = false;
          use_usd = true;
      }
      else
      {
        std::cout << "Wrong Feature name.\n";
        showHelp (argv[0]);
        exit (-1);
      }
    }*/

    //General parameters
    pcl::console::parse_argument (argc, argv, "--model_ss", model_ss_);
    pcl::console::parse_argument (argc, argv, "--scene_ss", scene_ss_);
    pcl::console::parse_argument (argc, argv, "--rf_rad", rf_rad_);
    pcl::console::parse_argument (argc, argv, "--descr_rad", descr_rad_);
    pcl::console::parse_argument (argc, argv, "--cg_size", cg_size_);
    pcl::console::parse_argument (argc, argv, "--cg_thresh", cg_thresh_);
  }

  //***********************************************************************************************
  //                                  Função Cálculo resolução nuvem
  //***********************************************************************************************
  double
  computeCloudResolution (const pcl::PointCloud<PointType>::ConstPtr &cloud)
  {
    double res = 0.0;
    int n_points = 0;
    int nres;

    std::vector<int> indices (2);
    std::vector<float> sqr_distances (2);
    pcl::search::KdTree<PointType> tree;
    tree.setInputCloud (cloud);

    for (size_t i = 0; i < cloud->size (); ++i)
    {
      if (! pcl_isfinite ((*cloud)[i].x))
      {
        continue;
      }
      //Considering the second neighbor since the first is the point itself.
      nres = tree.nearestKSearch (i, 2, indices, sqr_distances);
      if (nres == 2)
      {
        res += sqrt (sqr_distances[1]);
        ++n_points;
      }
    }
    if (n_points != 0)
    {
      res /= n_points;
    }
    return res;
  }

  //***********************************************************************************************
  //                                          Função PRINCIPAL
  //***********************************************************************************************
  int
  main (int argc, char *argv[])
  {
    parseCommandLine (argc, argv);

    //if(use_shot)
    //{
          //pcl::SHOT352 DescriptorType;
    //}
    //else if(use_pfh)
    //{
          typedef pcl::FPFHSignature33 DescriptorType;
    //}
    //else if(use_usd)
    //{

    //}

    pcl::PointCloud<PointTypeViewPoint>::Ptr modelviewpoint (new pcl::PointCloud<PointTypeViewPoint> ());

    pcl::PointCloud<PointType>::Ptr model (new pcl::PointCloud<PointType> ());
    pcl::PointCloud<PointType>::Ptr scene (new pcl::PointCloud<PointType> ());

    pcl::PointCloud<PointType>::Ptr model_keypoints (new pcl::PointCloud<PointType> ());
    pcl::PointCloud<PointType>::Ptr scene_keypoints (new pcl::PointCloud<PointType> ());

    pcl::PointCloud<NormalType>::Ptr model_normals (new pcl::PointCloud<NormalType> ());
    pcl::PointCloud<NormalType>::Ptr scene_normals (new pcl::PointCloud<NormalType> ());

    pcl::PointCloud<NormalType>::Ptr model_keypoints_normals (new pcl::PointCloud<NormalType> ());
    pcl::PointCloud<NormalType>::Ptr scene_keypoints_normals (new pcl::PointCloud<NormalType> ());

    pcl::PointCloud<DescriptorType>::Ptr model_descriptors (new pcl::PointCloud<DescriptorType> ());
    pcl::PointCloud<DescriptorType>::Ptr scene_descriptors (new pcl::PointCloud<DescriptorType> ());

    //
    //********************************  Load clouds model and scene  ******************************
    //

    if (pcl::io::loadPCDFile (scene_filename_, *scene) < 0)
    {
      std::cout << "Error loading scene cloud." << std::endl;
      showHelp (argv[0]);
      return (-1);
    }

  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  printf("Loading %ld models\n", model_filenames_.size());
  for(size_t v=0; v<model_filenames_.size(); ++v)///////////////////////////////////////////////////////////////////////////////
  {

    printf("Analysing for model view %ld\n", v);
          if (pcl::io::loadPCDFile (model_filenames_[v], *modelviewpoint) < 0)
          {
          std::cout << "Error loading model cloud." << std::endl;
          showHelp (argv[0]);
          return (-1);
          }
          else
          {
          std::cout << "loading model cloud " << model_filenames_[v] << std::endl;
          }
  // Recolha dos view points para calculo das normais
  float vpx = modelviewpoint->points[0].vp_x;
  float vpy = modelviewpoint->points[0].vp_y;
  float vpz = modelviewpoint->points[0].vp_z;

  // Copia da nuvem para tipologia sem view points
  pcl::copyPointCloud<PointTypeViewPoint,PointType>(*modelviewpoint, *model);

  // see the resolution of the clouds
start = std::clock();                                    // Debug time

  double res_scene = computeCloudResolution(scene);
  double res_model = computeCloudResolution(model);

end = std::clock();                                       // Debug time
std::cout << "Time computing clouds resolution:  " << diffclock(start, end) << " ms" << endl << endl;

  std::cout << "Scene resolution:  " << res_scene << std::endl;
  std::cout << "Model resolution:  " << res_model << std::endl;

    //
    //  Set up resolution invariance
    //
    if (use_cloud_resolution_)
    {
      float resolution = static_cast<float> (computeCloudResolution (model));
      if (resolution != 0.0f)
      {
        model_ss_   *= resolution;
        scene_ss_   *= resolution;
        rf_rad_     *= resolution;
        descr_rad_  *= resolution;
        cg_size_    *= resolution;
      }

      std::cout << "Model resolution:       " << resolution << std::endl;
      std::cout << "Model sampling size:    " << model_ss_ << std::endl;
      std::cout << "Scene sampling size:    " << scene_ss_ << std::endl;
      std::cout << "LRF support radius:     " << rf_rad_ << std::endl;
      std::cout << "SHOT descriptor radius: " << descr_rad_ << std::endl;
      std::cout << "Clustering bin size:    " << cg_size_ << std::endl << std::endl;
    }

    //
    //***********************************  Compute Normals  ************************************
    //
    start = std::clock();                                    // Debug time

    pcl::NormalEstimationOMP<PointType, NormalType> norm_est;
    norm_est.setKSearch (5);

    norm_est.setInputCloud (scene);
    norm_est.compute (*scene_normals);

    norm_est.setInputCloud (model);
    norm_est.setViewPoint(vpx,vpy,vpz);
    norm_est.compute (*model_normals);

    end = std::clock();                               // Debug time
    std::cout << "Time computing normals:  " << diffclock(start, end) << " ms" << endl << endl;

    //
    //************************  Downsample Clouds to Extract keypoints  *************************
    //

          //****************************  UniformSampling
          if(use_uniformsampling)
          {
            start = std::clock();                                    // Debug time
          pcl::UniformSampling<PointType> uniform_sampling;
          uniform_sampling.setInputCloud (model);
          uniform_sampling.setRadiusSearch (model_ss_);
          //uniform_sampling.filter (*model_keypoints);
          pcl::PointCloud<int> keypointIndices1;
          uniform_sampling.compute(keypointIndices1);

          pcl::copyPointCloud(*model, keypointIndices1.points, *model_keypoints);
          pcl::copyPointCloud(*model_normals, keypointIndices1.points, *model_keypoints_normals);

            end = std::clock();                               // Debug time
            std::cout << "Time computing Model Keypoints:  " << diffclock(start, end) << " ms" << endl << endl;

          std::cout << "Model total points: " << model->size () << "; Selected Keypoints: " << model_keypoints->size () << std::endl;

            start = std::clock();                                    // Debug time

          uniform_sampling.setInputCloud (scene);
          uniform_sampling.setRadiusSearch (scene_ss_);
          //uniform_sampling.filter (*scene_keypoints);
          pcl::PointCloud<int> keypointIndices2;
          uniform_sampling.compute(keypointIndices2);

          pcl::copyPointCloud(*scene, keypointIndices2.points, *scene_keypoints);
          pcl::copyPointCloud(*scene_normals, keypointIndices2.points, *scene_keypoints_normals);

            end = std::clock();                               // Debug time
            std::cout << "Time computing Model Keypoints:  " << diffclock(start, end) << " ms" << endl << endl;

          std::cout << "Scene total points: " << scene->size () << "; Selected Keypoints: " << scene_keypoints->size () << std::endl;
          }

          //****************************  ISSKeypoint3D
          if(use_iss)
          {
          // ISS keypoint detector model.
          pcl::ISSKeypoint3D<pcl::PointXYZ, pcl::PointXYZ> detector;
          detector.setInputCloud(model);
          pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_model(new pcl::search::KdTree<pcl::PointXYZ>);
          detector.setSearchMethod(kdtree_model);
          double resolution_model = computeCloudResolution(model);
          // Set the radius of the spherical neighborhood used to compute the scatter matrix.
          detector.setSalientRadius(6 * resolution_model);
          // Set the radius for the application of the non maxima supression algorithm.
          detector.setNonMaxRadius(4 * resolution_model);
          // Set the minimum number of neighbors that has to be found while applying the non maxima suppression algorithm.
          detector.setMinNeighbors(5);
          // Set the upper bound on the ratio between the second and the first eigenvalue.
          detector.setThreshold21(0.975);
          // Set the upper bound on the ratio between the third and the second eigenvalue.
          detector.setThreshold32(0.975);
          // Set the number of prpcessing threads to use. 0 sets it to automatic.
          detector.setNumberOfThreads(4);

          detector.setNormalRadius (4 * resolution_model);
          detector.setBorderRadius (1 * resolution_model);

          detector.compute(*model_keypoints);
          std::cout << "model_keypoints: " << model_keypoints->size () << std::endl;//debug

          // ISS keypoint detector scene.
          //pcl::ISSKeypoint3D<pcl::PointXYZ, pcl::PointXYZ> detector;
          detector.setInputCloud(scene);
          pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_scene(new pcl::search::KdTree<pcl::PointXYZ>);
          detector.setSearchMethod(kdtree_scene);
          double resolution_scene = computeCloudResolution(scene);
          // Set the radius of the spherical neighborhood used to compute the scatter matrix.
          detector.setSalientRadius(6 * resolution_scene);
          // Set the radius for the application of the non maxima supression algorithm.
          detector.setNonMaxRadius(4 * resolution_scene);
          // Set the minimum number of neighbors that has to be found while applying the non maxima suppression algorithm.
          detector.setMinNeighbors(5);
          // Set the upper bound on the ratio between the second and the first eigenvalue.
          detector.setThreshold21(0.975);
          // Set the upper bound on the ratio between the third and the second eigenvalue.
          detector.setThreshold32(0.975);
          // Set the number of prpcessing threads to use. 0 sets it to automatic.
          detector.setNumberOfThreads(4);

          detector.setNormalRadius (4 * resolution_scene);
          detector.setBorderRadius (1 * resolution_scene);

          detector.compute(*scene_keypoints);
          std::cout << "scene_keypoints: " << scene_keypoints->size () << std::endl;//debug
          }
    //
    //  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //

    //
    //  Compute Descriptor for keypoints SHOT
    //
  /*	if(use_shot)
          {
          pcl::SHOTEstimationOMP<PointType, NormalType, DescriptorType> descr_est;
          descr_est.setRadiusSearch (descr_rad_);

          descr_est.setInputCloud (model_keypoints);
          descr_est.setInputNormals (model_normals);
          descr_est.setSearchSurface (model);
          descr_est.compute (*model_descriptors);

          descr_est.setInputCloud (scene_keypoints);
          descr_est.setInputNormals (scene_normals);
          descr_est.setSearchSurface (scene);
          descr_est.compute (*scene_descriptors);
          }
  */

    //
    //  Compute Descriptor for keypoints PFH
    //
  /*	if(use_pfh)
          {
          pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>);

          pcl::PFHEstimation<PointType, NormalType, DescriptorType> descr_est;
          descr_est.setRadiusSearch (descr_rad_);

          descr_est.setInputCloud (model_keypoints);
          descr_est.setInputNormals (model_keypoints_normals);
          descr_est.setSearchMethod(kdtree);
          // Search radius, to look for neighbors. Note: the value given here has to be
          // larger than the radius used to estimate the normals.
          descr_est.compute (*model_descriptors);

          descr_est.setInputCloud (scene_keypoints);
          descr_est.setInputNormals (scene_keypoints_normals);
          descr_est.setSearchMethod(kdtree);
          // Search radius, to look for neighbors. Note: the value given here has to be
          // larger than the radius used to estimate the normals.
          descr_est.compute (*scene_descriptors);
          }
  */

  ///////////////////////////////////////////////////////////////////////////////////////////////////////
  //  Compute Descriptor for keypoints FPFH
  //////////////////////////////////////////////////////////////////////////////////////////////////////

          //if(use_pfh)
          //{

          pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>);

          pcl::FPFHEstimation<PointType, NormalType, DescriptorType> descr_est;
          descr_est.setRadiusSearch (descr_rad_);
            start = std::clock();                                    // Debug time
          descr_est.setInputCloud (model_keypoints);
          descr_est.setInputNormals (model_keypoints_normals);
          descr_est.setSearchMethod(kdtree);
          // Search radius, to look for neighbors. Note: the value given here has to be
          // larger than the radius used to estimate the normals.
          descr_est.compute (*model_descriptors);

            end = std::clock();                               // Debug time
            std::cout << "Time computing Model Descriptors:  " << diffclock(start, end) << " ms" << endl << endl;

            start = std::clock();                                    // Debug time
          descr_est.setInputCloud (scene_keypoints);
          descr_est.setInputNormals (scene_keypoints_normals);
          descr_est.setSearchMethod(kdtree);
          // Search radius, to look for neighbors. Note: the value given here has to be
          // larger than the radius used to estimate the normals.
          descr_est.compute (*scene_descriptors);

            end = std::clock();                               // Debug time
            std::cout << "Time computing Scene Descriptors:  " << diffclock(start, end) << " ms" << endl << endl;
          //}

  /*
    //
    //  Compute Descriptor for keypoints Spin
    //
    pcl::SpinImageEstimation<PointType, NormalType, DescriptorType> descr_est;
    // Radius of the support cylinder.
    descr_est.setRadiusSearch(descr_rad_);
          // Set the resolution of the spin image (the number of bins along one dimension).
          // Note: you must change the output histogram size to reflect this.
    descr_est.setImageWidth(8);

    descr_est.setInputCloud (model_keypoints);
    descr_est.setInputNormals (model_normals);
    descr_est.compute (*model_descriptors);

    descr_est.setInputCloud (scene_keypoints);
    descr_est.setInputNormals (scene_normals);
    descr_est.compute (*scene_descriptors);
  */

    //
    //  Find Model-Scene Correspondences with KdTree for SHOT
    //
  /*if(use_shot)
  {
    pcl::CorrespondencesPtr model_scene_corrs (new pcl::Correspondences ());

    pcl::KdTreeFLANN<DescriptorType> match_search;
    match_search.setInputCloud (model_descriptors);

    //  For each scene keypoint descriptor, find nearest neighbor into the model keypoints descriptor cloud and add it to the correspondences vector.
    //  Para cada descritor de ponto chave da cena, encontra o mais proximo vizinho dentro dos descritores de pontos chave do modelo e junta os ao vector de correspondecias.
    for (size_t i = 0; i < scene_descriptors->size (); ++i)
    {
      std::vector<int> neigh_indices (1);
      std::vector<float> neigh_sqr_dists (1);

      if (!pcl_isfinite (scene_descriptors->at (i).descriptor[0])) //skipping NaNs
      {
        continue;
      }

      int found_neighs = match_search.nearestKSearch (scene_descriptors->at (i), 1, neigh_indices, neigh_sqr_dists);
      if(found_neighs == 1 && neigh_sqr_dists[0] < 0.30f) //  add match only if the squared descriptor distance is less than 0.25 (SHOT descriptor distances are between 0 and 1 by design)
      {
        pcl::Correspondence corr (neigh_indices[0], static_cast<int> (i), neigh_sqr_dists[0]);
        model_scene_corrs->push_back (corr);
      }
    }
    std::cout << "Correspondences found: " << model_scene_corrs->size () << std::endl;
  }*/

    //
    //  Find Model-Scene Correspondences with KdTree for PFH
    //

  //if(use_pfh)
  //{
    pcl::CorrespondencesPtr model_scene_corrs (new pcl::Correspondences ());

    pcl::KdTreeFLANN<DescriptorType> match_search;
    match_search.setInputCloud (model_descriptors);

    //  For each scene keypoint descriptor, find nearest neighbor into the model keypoints descriptor cloud and add it to the correspondences vector.
    //  Para cada descritor de ponto chave da cena, encontra o mais proximo vizinho dentro dos descritores de pontos chave do modelo e junta os ao vector de correspondecias.

      start = std::clock();                                    // Debug time

    for (size_t i = 0; i < scene_descriptors->size (); ++i)
    {
      std::vector<int> neigh_indices (1);
      std::vector<float> neigh_sqr_dists (1);

      if (!pcl_isfinite (scene_descriptors->at (i).histogram[0])) //skipping NaNs
      {
        continue;
      }

      int found_neighs = match_search.nearestKSearch (scene_descriptors->at (i), 1, neigh_indices, neigh_sqr_dists);
          //std::cout << "neigh_sqr_dists: " << neigh_sqr_dists[0] << std::endl; //debugg
      if(found_neighs == 1 && neigh_sqr_dists[0] < 700.0f) //  add match only if the squared descriptor distance is less than 0.25 (SHOT descriptor distances are between 0 and 1 by design)
      {
        pcl::Correspondence corr (neigh_indices[0], static_cast<int> (i), neigh_sqr_dists[0]);
        model_scene_corrs->push_back (corr);
      }
    }
    std::cout << "Correspondences found: " << model_scene_corrs->size () << std::endl;

      end = std::clock();                               // Debug time
      std::cout << "Time to find scene-model correspondences:  " << diffclock(start, end) << " ms" << endl << endl;
  //}

model_view.push_back(chave());

    //
    //  Actual Clustering
    //
    //std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f> > rototranslations;
    //std::vector<pcl::Correspondences> clustered_corrs;

    //**************  Using Hough3D
          if (use_hough_)
          {
          //  Compute (Keypoints) Reference Frames only for Hough
          pcl::PointCloud<RFType>::Ptr model_rf (new pcl::PointCloud<RFType> ());
          pcl::PointCloud<RFType>::Ptr scene_rf (new pcl::PointCloud<RFType> ());

          pcl::BOARDLocalReferenceFrameEstimation<PointType, NormalType, RFType> rf_est;
          rf_est.setFindHoles (true);
          rf_est.setRadiusSearch (rf_rad_);

          rf_est.setInputCloud (model_keypoints);
          rf_est.setInputNormals (model_normals);
          rf_est.setSearchSurface (model);
          rf_est.compute (*model_rf);

          rf_est.setInputCloud (scene_keypoints);
          rf_est.setInputNormals (scene_normals);
          rf_est.setSearchSurface (scene);
          rf_est.compute (*scene_rf);

          //  Clustering
          pcl::Hough3DGrouping<PointType, PointType, RFType, RFType> clusterer;
          clusterer.setHoughBinSize (cg_size_);
          clusterer.setHoughThreshold (cg_thresh_);
          clusterer.setUseInterpolation (true);
          clusterer.setUseDistanceWeight (false);

          clusterer.setInputCloud (model_keypoints);
          clusterer.setInputRf (model_rf);
          clusterer.setSceneCloud (scene_keypoints);
          clusterer.setSceneRf (scene_rf);
          clusterer.setModelSceneCorrespondences (model_scene_corrs);

          //clusterer.cluster (clustered_corrs);
          clusterer.recognize (model_view.at(v).rototranslations, model_view.at(v).clustered_corrs);
          }
          else // Using GeometricConsistency
          {
            start = std::clock();                                    // Debug time
          pcl::GeometricConsistencyGrouping<PointType, PointType> gc_clusterer;
          gc_clusterer.setGCSize (cg_size_);
          gc_clusterer.setGCThreshold (cg_thresh_);

          gc_clusterer.setInputCloud (model_keypoints);
          gc_clusterer.setSceneCloud (scene_keypoints);
          gc_clusterer.setModelSceneCorrespondences (model_scene_corrs);

          //gc_clusterer.cluster (clustered_corrs);

          //printf("Adding %ld rototranslations")
          gc_clusterer.recognize (model_view[v].rototranslations, model_view[v].clustered_corrs);
            end = std::clock();                               // Debug time
            std::cout << "Time clustering:  " << diffclock(start, end) << " ms" << endl << endl;
          }

    //
    //  Output results
    //
    std::cout << "Model instances found: " << model_view.at(v).rototranslations.size () << std::endl;
    for (size_t i = 0; i < model_view.at(v).rototranslations.size (); ++i)
    {
      std::cout << "\n    Instance " << i + 1 << ":" << std::endl;
      std::cout << "        Correspondences belonging to this instance: " << model_view.at(v).clustered_corrs[i].size () << std::endl;

      // Print the rotation matrix and translation vector
      Eigen::Matrix3f rotation = model_view.at(v).rototranslations[i].block<3,3>(0, 0);
      Eigen::Vector3f translation = model_view.at(v).rototranslations[i].block<3,1>(0, 3);

      printf ("\n");
      printf ("            | %6.3f %6.3f %6.3f | \n", rotation (0,0), rotation (0,1), rotation (0,2));
      printf ("        R = | %6.3f %6.3f %6.3f | \n", rotation (1,0), rotation (1,1), rotation (1,2));
      printf ("            | %6.3f %6.3f %6.3f | \n", rotation (2,0), rotation (2,1), rotation (2,2));
      printf ("\n");
      printf ("        t = < %0.3f, %0.3f, %0.3f >\n", translation (0), translation (1), translation (2));
    }

  }                                                         // end of Multi-view For cycle

    //
    //  Visualization
    //
    pcl::visualization::PCLVisualizer viewer ("Correspondence Grouping");

    pcl::visualization::PointCloudColorHandlerCustom<PointType> scene_color_handler (scene, 255, 255, 255);
    viewer.addPointCloud (scene,scene_color_handler, "scene_cloud");

    std::string  	id = "reference";
    viewer.addCoordinateSystem(0.1, id, 0);

    pcl::PointCloud<PointType>::Ptr off_scene_model (new pcl::PointCloud<PointType> ());
    pcl::PointCloud<PointType>::Ptr off_scene_model_keypoints (new pcl::PointCloud<PointType> ());

  for (size_t v = 0; v < model_view.size (); ++v)
  {
  /*
    if (show_correspondences_ || show_keypoints_)
    {
      //  We are translating the model so that it doesn't end in the middle of the scene representation
      pcl::transformPointCloud (*model, *off_scene_model, Eigen::Vector3f (-1,0,0), Eigen::Quaternionf (1, 0, 0, 0));
      pcl::transformPointCloud (*model_keypoints, *off_scene_model_keypoints, Eigen::Vector3f (-1,0,0), Eigen::Quaternionf (1, 0, 0, 0));

      pcl::visualization::PointCloudColorHandlerCustom<PointType> off_scene_model_color_handler (off_scene_model, 255, 255, 0);
      viewer.addPointCloud (off_scene_model, off_scene_model_color_handler, "off_scene_model");
          //viewer.addCoordinateSystem(0.1,0);
    }
  */

  /*
    if (show_keypoints_)
    {
      pcl::visualization::PointCloudColorHandlerCustom<PointType> scene_keypoints_color_handler (scene_keypoints, 0, 128, 255);
      viewer.addPointCloud (scene_keypoints, scene_keypoints_color_handler, "scene_keypoints");
      viewer.setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 5, "scene_keypoints");

      pcl::visualization::PointCloudColorHandlerCustom<PointType> off_scene_model_keypoints_color_handler (off_scene_model_keypoints, 0, 128, 255);
      viewer.addPointCloud (off_scene_model_keypoints, off_scene_model_keypoints_color_handler, "off_scene_model_keypoints");
      viewer.setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 5, "off_scene_model_keypoints");
    }
  */

  /*
    if (show_normals_)
    {

    viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal> (scene, scene_normals, 1, 0.005, "normals");
    viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal> (off_scene_model, model_normals, 1, 0.005, "normal");

    }
  */

    for (size_t i = 0; i < model_view.at(v).rototranslations.size (); ++i)
    {
      pcl::PointCloud<PointType>::Ptr rotated_model (new pcl::PointCloud<PointType> ());
      pcl::transformPointCloud (*model, *rotated_model, model_view.at(v).rototranslations[i]);

      int inst = i + ((v + 1) * 10);

      std::stringstream ss_cloud;
      ss_cloud << "instance" << inst;

      pcl::visualization::PointCloudColorHandlerCustom<PointType> rotated_model_color_handler (rotated_model, 255, 0, 0);
      viewer.addPointCloud (rotated_model, rotated_model_color_handler, ss_cloud.str ());

          // geração da matriz de rotação para os referenciais dos objectos detectados
          Eigen::Affine3f t;

          Eigen::Matrix3f rotation = model_view.at(v).rototranslations[i].block<3,3>(0, 0);
          Eigen::Vector3f translation = model_view.at(v).rototranslations[i].block<3,1>(0, 3);

          t.linear() = rotation;
          t.translation() = translation;

          //std::string car = i;

          std::string number;
          std::stringstream strstream;

          strstream << inst;
          strstream >> number;
          //std::string &id = i;

          viewer.addCoordinateSystem(0.05, t, number, 0);
      /*
      if (show_correspondences_)
      {
        for (size_t j = 0; j < clustered_corrs[i].size (); ++j)
        {
          std::stringstream ss_line;
          ss_line << "correspondence_line" << i << "_" << j;
          PointType& model_point = off_scene_model_keypoints->at (clustered_corrs[i][j].index_query);
          PointType& scene_point = scene_keypoints->at (clustered_corrs[i][j].index_match);

          //  We are drawing a line for each pair of clustered correspondences found between the model and the scene
          viewer.addLine<PointType, PointType> (model_point, scene_point, 0, 255, 0, ss_line.str ());
        }
      }*/
    }
  }

    while (!viewer.wasStopped ())
    {
      viewer.spinOnce ();
    }

    return (0);

  }


text:
/**
 * @file windvector.cpp
 */

#include "windvector.h"
#include "forecast_time.h"
#include "lambert_conformal_grid.h"
#include "latitude_longitude_grid.h"
#include "level.h"
#include "logger.h"
#include "plugin_factory.h"
#include "stereographic_grid.h"
#include "util.h"
#include <boost/thread.hpp>
#include <iostream>
#include <math.h>

#include "cache.h"
#include "fetcher.h"

using namespace std;
using namespace himan::plugin;

#include "cuda_helper.h"

typedef tuple<double, double, double, double> coefficients;

boost::thread_specific_ptr<map<size_t, coefficients>> myCoefficientCache;

windvector::windvector() : itsCalculationTarget(kUnknownElement), itsVectorCalculation(false)
{
	itsCudaEnabledCalculation = true;

	itsLogger = logger("windvector");
}

void windvector::Process(const std::shared_ptr<const plugin_configuration> conf)
{
	Init(conf);

	/*
	 * Set target parameter to windvector
	 */

	vector<param> theParams;

	param requestedDirParam;
	param requestedSpeedParam;
	param requestedVectorParam;

	if (itsConfiguration->Exists("do_vector") && itsConfiguration->GetValue("do_vector") == "true")
	{
		itsVectorCalculation = true;
	}

	if (itsConfiguration->Exists("for_ice") && itsConfiguration->GetValue("for_ice") == "true")
	{
		requestedSpeedParam = param("IFF-MS", 389, 10, 2, 3);
		requestedDirParam = param("IDD-D", 390, 10, 2, 2);

		itsCalculationTarget = kIce;

		if (itsVectorCalculation)
		{
			itsLogger.Warning("Unable to calculate vector for ice");
		}
	}
	else if (itsConfiguration->Exists("for_sea") && itsConfiguration->GetValue("for_sea") == "true")
	{
		requestedSpeedParam = param("SFF-MS", 163, 10, 1, 1);
		requestedDirParam = param("SDD-D", 164, 10, 1, 0);

		itsCalculationTarget = kSea;

		if (itsVectorCalculation)
		{
			itsLogger.Warning("Unable to calculate vector for sea");
		}
	}
	else if (itsConfiguration->Exists("for_gust") && itsConfiguration->GetValue("for_gust") == "true")
	{
		requestedSpeedParam = param("FFG-MS", 417, 0, 2, 22);

		itsCalculationTarget = kGust;

		if (itsVectorCalculation)
		{
			itsLogger.Warning("Unable to calculate vector for wind gust");
		}
	}
	else
	{
		// By default assume we'll calculate for wind

		requestedDirParam = param("DD-D", 20, 0, 2, 0);
		requestedSpeedParam = param("FF-MS", 21, 0, 2, 1);

		if (itsVectorCalculation)
		{
			requestedVectorParam = param("DF-MS", 22);
		}

		itsCalculationTarget = kWind;
	}

	theParams.push_back(requestedSpeedParam);

	if (itsCalculationTarget != kGust)
	{
		theParams.push_back(requestedDirParam);
	}

	if (itsVectorCalculation)
	{
		theParams.push_back(requestedVectorParam);
	}

	SetParams(theParams);

	Start();
}

/*
 * Calculate()
 *
 * This function does the actual calculation.
 */

void windvector::Calculate(shared_ptr<info> myTargetInfo, unsigned short threadIndex)
{
	if (!myCoefficientCache.get())
	{
		myCoefficientCache.reset(new map<size_t, coefficients>());
	}

	// Required source parameters

	param UParam;
	param VParam;

	double directionOffset = 180;  // For wind direction add this

	switch (itsCalculationTarget)
	{
		case kSea:
			UParam = param("WVELU-MS");
			VParam = param("WVELV-MS");
			directionOffset = 0;
			break;

		case kIce:
			UParam = param("IVELU-MS");
			VParam = param("IVELV-MS");
			directionOffset = 0;
			break;

		case kGust:
			UParam = param("WGU-MS");
			VParam = param("WGV-MS");
			break;

		case kWind:
			UParam = param("U-MS");
			VParam = param("V-MS");
			break;

		default:
			throw runtime_error("Invalid calculation target element: " +
			                    to_string(static_cast<int>(itsCalculationTarget)));
			break;
	}

	auto myThreadedLogger = logger("windvectorThread #" + to_string(threadIndex));

	forecast_time forecastTime = myTargetInfo->Time();
	level forecastLevel = myTargetInfo->Level();
	forecast_type forecastType = myTargetInfo->ForecastType();

	myThreadedLogger.Info("Calculating time " + static_cast<string>(forecastTime.ValidDateTime()) + " level " +
	                      static_cast<string>(forecastLevel));

	info_t UInfo = Fetch(forecastTime, forecastLevel, UParam, forecastType, itsConfiguration->UseCudaForPacking());
	info_t VInfo = Fetch(forecastTime, forecastLevel, VParam, forecastType, itsConfiguration->UseCudaForPacking());

	if (!UInfo || !VInfo)
	{
		myThreadedLogger.Warning("Skipping step " + to_string(forecastTime.Step()) + ", level " +
		                         static_cast<string>(forecastLevel));
		return;
	}

	assert(UInfo->Grid()->AB() == VInfo->Grid()->AB());

	for (myTargetInfo->ResetParam(); myTargetInfo->NextParam();)
	{
		SetAB(myTargetInfo, UInfo);
	}

	assert(UInfo->Grid()->Type() == VInfo->Grid()->Type());

	string deviceType;

#ifdef HAVE_CUDA
	if (itsConfiguration->UseCuda() &&
	    (UInfo->Grid()->Type() == kLatitudeLongitude || UInfo->Grid()->Type() == kRotatedLatitudeLongitude ||
	     UInfo->Grid()->Type() == kLambertConformalConic))
	{
		deviceType = "GPU";

		auto opts = CudaPrepare(myTargetInfo, UInfo, VInfo);

		windvector_cuda::Process(*opts);
	}
	else
#endif
	{
		deviceType = "CPU";

#ifdef HAVE_CUDA
		Unpack({UInfo, VInfo});
#endif
		myTargetInfo->ParamIndex(0);

		auto& FFVec = VEC(myTargetInfo);
		vector<double> DDVec(FFVec.size(), kFloatMissing);

		for (auto&& tup : zip_range(FFVec, DDVec, VEC(UInfo), VEC(VInfo)))
		{
			double& speed = tup.get<0>();
			double& dir = tup.get<1>();
			double U = tup.get<2>();
			double V = tup.get<3>();

			if (U == kFloatMissing || V == kFloatMissing)
			{
				continue;
			}

			speed = sqrt(U * U + V * V);

			if (itsCalculationTarget == kGust)
			{
				continue;
			}

			if (speed > 0)
			{
				dir = himan::constants::kRad * atan2(U, V) + directionOffset;

				// reduce the angle
				dir = fmod(dir, 360);

				// force it to be the positive remainder, so that 0 <= dir < 360
				dir = round(fmod((dir + 360), 360));
			}
		}

		if (myTargetInfo->SizeParams() > 1)
		{
			myTargetInfo->ParamIndex(1);
			myTargetInfo->Data().Set(DDVec);
		}
	}

	myThreadedLogger.Info("[" + deviceType + "] Missing values: " + to_string(myTargetInfo->Data().MissingCount()) +
	                      "/" + to_string(myTargetInfo->Data().Size()));
}

shared_ptr<himan::info> windvector::Fetch(const forecast_time& theTime, const level& theLevel, const param& theParam,
                                          const forecast_type& theType, bool returnPacked) const
{
	auto f = GET_PLUGIN(fetcher);
	f->DoVectorComponentRotation(true);

	info_t ret;

	try
	{
		ret = f->Fetch(itsConfiguration, theTime, theLevel, theParam, theType, itsConfiguration->UseCudaForPacking());

#ifdef HAVE_CUDA
		if (!returnPacked && ret->Grid()->IsPackedData())
		{
			assert(ret->Grid()->PackedData().ClassName() == "simple_packed");

			util::Unpack({ret->Grid()});

			auto c = GET_PLUGIN(cache);

			c->Insert(*ret);
		}
#endif
	}
	catch (HPExceptionType& e)
	{
		if (e != kFileDataNotFound)
		{
			throw runtime_error(ClassName() + ": Unable to proceed");
		}
	}

	return ret;
}

#ifdef HAVE_CUDA

unique_ptr<windvector_cuda::options> windvector::CudaPrepare(shared_ptr<info> myTargetInfo, shared_ptr<info> UInfo,
                                                             shared_ptr<info> VInfo)
{
	unique_ptr<windvector_cuda::options> opts(new windvector_cuda::options);

	opts->N = UInfo->Grid()->Size();

	opts->target_type = itsCalculationTarget;

	opts->u = UInfo->ToSimple();
	opts->v = VInfo->ToSimple();

	myTargetInfo->ParamIndex(0);
	opts->speed = myTargetInfo->ToSimple();

	if (opts->target_type != kGust)
	{
		myTargetInfo->ParamIndex(1);
		opts->dir = myTargetInfo->ToSimple();
	}

	return opts;
}

#endif

max_stars_repo_name: jrintala/fmi-data

max_issues_repo_name: jrintala/fmi-data

max_forks_repo_name: jrintala/fmi-data


text:
// Copyright (c) 2015 - 2016 Anyar, Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "ascent/io/JsonAPI.h"

#include "ascent/algorithms/Interpolation.h"
#include "ascent/core/ModuleCore.h"

#include <Eigen/Dense>

#include <cmath>
#include <limits>

using namespace asc;
using namespace std;

jsoncons::json JsonAPI::jsonModules()
{
   jsoncons::json modules = jsoncons::json::array();

   for (auto& p : ModuleCore::external)
      modules.add(p.first);

   jsoncons::json output;
   output["modules"] = modules;

   return output;
}

jsoncons::json JsonAPI::jsonVariables()
{
   jsoncons::json output = jsoncons::json::array();

   for (auto& p : ModuleCore::external)
   {
      string module = p.first;
      const auto& names = p.second->vars.getNames();
      for (auto& var : names)
      {
         jsoncons::json obj;
         obj["module"] = module;
         obj["type"] = var.first;
         obj["var"] = var.second;
         output.add(std::move(obj));
      }
   }

   return output;
}

jsoncons::json JsonAPI::io(const std::string& input)
{
   jsoncons::json in;

   try
   {
      in = jsoncons::json::parse_file(input);
   }
   catch (const jsoncons::json_exception& e)
   {
      cerr << e.what() << '\n';
   }

   return io(in);
}

jsoncons::json JsonAPI::io(jsoncons::json& input)
{
   jsoncons::json output = jsoncons::json::array();

   for (size_t i = 0; i < input.size(); ++i) // for all modules
   {
      try
      {
         jsoncons::json& module = input[i];
         if (module.count("module")) // this object is a module
         {
            const string name = module["module"].as<string>();
            Module& base = ModuleCore::getExternal(name);

            jsoncons::json module_out;
            module_out["module"] = name;

            if (module.count("chaiscript"))
            {
               base.chai.eval(module["chaiscript"].as<string>());
            }
            else if (module.count("modules"))
            {
               // Want recursion for multiple modules deep.
               module_out["modules"] = io(module["modules"]);
            }
            else if (module.count("vars"))
            {
               jsoncons::json& variables = module["vars"]; // array of variables
               jsoncons::json variables_out = jsoncons::json::array();

               for (size_t v = 0; v < variables.size(); ++v) // for all variables
               {
                  jsoncons::json& variable = variables[v];
                  jsoncons::json variable_out;

                  if (variable.count("var"))
                  {
                     bool success = access<std::string, double, bool, int, unsigned, size_t, vector<double>, vector<bool>, vector<int>, vector<unsigned>, vector<size_t>>(variable, variable_out, base);

                     // IMPROVEMENT: Eigen handling should be inserted into access function, but this will take some thought. I'm not sure how it could be handled generically because Derived is unknown for Eigen::MatrixBase<Derived>.
                     // Important Eigen allignment note: In order to avoid allignment issues, static Eigen parameters are passed via json objects as dynamic Eigen vectors or matrices.
                     if (!success)
                     {
                        std::string type = variable["type"].as<string>();
                        std::string var = variable["var"].as<string>();

                        if (variable.count("value")) // If a "value" member is specified, then we assume the user wants to set the variable.
                        {
                           jsoncons::json& value = variable["value"];

                           if (typeid(Eigen::Vector2d).name() == type)
                              success = base.vars.set(var, static_cast<Eigen::Vector2d>(value.any_cast<Eigen::VectorXd>()));
                           else if (typeid(Eigen::Vector3d).name() == type)
                              success = base.vars.set(var, value.any_cast<Eigen::Vector3d>());
                           else if (typeid(Eigen::Vector4d).name() == type)
                              success = base.vars.set(var, static_cast<Eigen::Vector4d>(value.any_cast<Eigen::VectorXd>()));
                           else if (typeid(Eigen::VectorXd).name() == type)
                              success = base.vars.set(var, value.any_cast<Eigen::VectorXd>());

                           else if (typeid(Eigen::Matrix2d).name() == type)
                              success = base.vars.set(var, static_cast<Eigen::Matrix2d>(value.any_cast<Eigen::MatrixXd>()));
                           else if (typeid(Eigen::Matrix3d).name() == type)
                              success = base.vars.set(var, value.any_cast<Eigen::Matrix3d>());
                           else if (typeid(Eigen::Matrix4d).name() == type)
                              success = base.vars.set(var, static_cast<Eigen::Matrix4d>(value.any_cast<Eigen::MatrixXd>()));
                           else if (typeid(Eigen::MatrixXd).name() == type)
                              success = base.vars.set(var, value.any_cast<Eigen::MatrixXd>());
                        }

                        if (typeid(Eigen::Vector2d).name() == type)
                           variable_out.set("value", jsoncons::json::any(static_cast<Eigen::VectorXd>(base.vars.get<Eigen::Vector2d>(var))));
                        else if (typeid(Eigen::Vector3d).name() == type)
                           variable_out.set("value", jsoncons::json::any(base.vars.get<Eigen::Vector3d>(var)));
                        else if (typeid(Eigen::Vector4d).name() == type)
                           variable_out.set("value", jsoncons::json::any(static_cast<Eigen::VectorXd>(base.vars.get<Eigen::Vector4d>(var))));
                        else if (typeid(Eigen::VectorXd).name() == type)
                           variable_out.set("value", jsoncons::json::any(base.vars.get<Eigen::VectorXd>(var)));

                        else if (typeid(Eigen::Matrix2d).name() == type)
                           variable_out.set("value", jsoncons::json::any(static_cast<Eigen::MatrixXd>(base.vars.get<Eigen::Matrix2d>(var))));
                        else if (typeid(Eigen::Matrix3d).name() == type)
                           variable_out.set("value", jsoncons::json::any(base.vars.get<Eigen::Matrix3d>(var)));
                        else if (typeid(Eigen::Matrix4d).name() == type)
                           variable_out.set("value", jsoncons::json::any(static_cast<Eigen::MatrixXd>(base.vars.get<Eigen::Matrix4d>(var))));
                        else if (typeid(Eigen::MatrixXd).name() == type)
                           variable_out.set("value", jsoncons::json::any(base.vars.get<Eigen::MatrixXd>(var)));
                     }

                     variables_out.add(std::move(variable_out));
                  }
               }

               module_out["vars"] = variables_out;
            }

            output.add(std::move(module_out));
         }
      }
      catch (const jsoncons::json_exception& e)
      {
         cerr << e.what() << '\n';
      }
   }

   return output;
}

jsoncons::json JsonAPI::ioFile(const std::string& filename)
{
   jsoncons::json input;

   try
   {
      input = jsoncons::json::parse_file(filename);
   }
   catch (const jsoncons::json_exception& e)
   {
      cerr << e.what() << '\n';
   }

   return io(input);
}

max_stars_repo_name: AnyarInc/ascent-beta

max_issues_repo_name: AnyarInc/ascent-beta

max_forks_repo_name: AnyarInc/ascent-beta


text:
#ifndef BOOST_MPL_AUX_MSVC_IS_CLASS_HPP_INCLUDED
#define BOOST_MPL_AUX_MSVC_IS_CLASS_HPP_INCLUDED

// Copyright Aleksey Gurtovoy 2002-2004
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
// See http://www.boost.org/libs/mpl for documentation.

// $Id$
// $Date$
// $Revision$

#include <boost/mpl/if.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/mpl/aux_/type_wrapper.hpp>
#include <boost/mpl/aux_/yes_no.hpp>

#include <boost/type_traits/is_reference.hpp>

namespace mars_boost {} namespace boost = mars_boost; namespace mars_boost { namespace mpl { namespace aux {

template< typename T > struct is_class_helper
{
    typedef int (T::* type)();
};

// MSVC 6.x-specific lightweight 'is_class' implementation;
// Distinguishing feature: does not instantiate the type being tested.
template< typename T >
struct msvc_is_class_impl
{
    template< typename U>
    static yes_tag  test(type_wrapper<U>*, /*typename*/ is_class_helper<U>::type = 0);
    static no_tag   test(void const volatile*, ...);

    enum { value = sizeof(test((type_wrapper<T>*)0)) == sizeof(yes_tag) };
    typedef bool_<value> type;
};

// agurt, 17/sep/04: have to check for 'is_reference' upfront to avoid ICEs in
// complex metaprograms
template< typename T >
struct msvc_is_class
    : if_<
          is_reference<T>
        , false_
        , msvc_is_class_impl<T>
        >::type
{
};

}}}

#endif // BOOST_MPL_AUX_MSVC_IS_CLASS_HPP_INCLUDED

max_stars_repo_name: jonetomtom/mars

max_issues_repo_name: jonetomtom/mars

max_forks_repo_name: pengjinning/mars