text:
/*

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

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

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifndef ELECTRODES_HPP_
#define ELECTRODES_HPP_

#include <boost/shared_ptr.hpp>
#include "ChasteSerialization.hpp"
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/shared_ptr.hpp>

#include "AbstractTetrahedralMesh.hpp"
#include "DistributedVector.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "ConstBoundaryCondition.hpp"

/**
 *  A class for setting up the boundary conditions associated with electrodes.
 *  There are two modes: grounding the second electrode, in which case the first
 *  electrode has an input flux (Neumann boundary condition) of the specified
 *  magnitude, the extracellular potential is fixed on the second electrode; or
 *  not grounded, in which case the opposite electrode has an equal and opposite
 *  output flux.
 *
 *  This class assumes the given mesh is cuboid, and the electrodes are taken to be
 *  the specified opposite surfaces.
 *
 *  Note the class now includes a pointer to the mesh as a member variable, as this
 *  is required to archive (reconstruct) itself. Because boost is clever it will
 *  still only archive one copy of the mesh.
 */
template<unsigned DIM>
class Electrodes
{
friend class TestBidomainWithBathProblem;
friend class TestElectrodes;

private:
    /** Whether the second electrode is grounded */
    bool mGroundSecondElectrode;
    /** The created bcc, which BidomainProblem will use */
    boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,2> > mpBoundaryConditionsContainer;
    /** The time the electrodes are switched on */
    double mStartTime;
    /** The time the electrodes are switched off */
    double mEndTime;
    /** Whether the electrodes are currently switched on */
    bool mAreActive;
    /**
     * This is needed for archiving:
     * the boundary conditions refer to nodes and/or elements, so need
     * the mesh to be archived, but don't have a pointer to the mesh itself.
     */
    AbstractTetrahedralMesh<DIM,DIM>* mpMesh;

    /** Left electrode area*/
    double mLeftElectrodeArea;

    /** Right electrode area*/
    double mRightElectrodeArea;

    /**
     * Helper method to compute electrodes area and check they are equal. Throws if they are not.
     *
     *  @param index The value i when applying the electrodes to x_i=a and x_i=b (a<b)
     *  @param lowerValue The value a when applying the electrodes to x_i=a and x_i=b (a<b) (should
     *    be the minimum value of x_i for the given mesh)
     *  @param upperValue The value b when applying the electrodes to x_i=a and x_i=b (a<b) (should
     *    be the maximum value of x_i for the given mesh)
     */
    void ComputeElectrodesAreasAndCheckEquality(unsigned index, double lowerValue, double upperValue);

    /** Needed for serialization. */
    friend class boost::serialization::access;
    /**
     * Save the Electrodes class
     *
     * @param archive the archive
     * @param version the current version of this class
     */
    template<class Archive>
    void save(Archive & archive, const unsigned int version) const
    {
        archive & mGroundSecondElectrode;
        archive & mStartTime;
        archive & mEndTime;
        archive & mAreActive;
        archive & mpMesh;
        (*ProcessSpecificArchive<Archive>::Get()) & mpBoundaryConditionsContainer;
    }
    /**
     * Load the Electrodes class
     *
     * @param archive the archive
     * @param version the current version of this class
     */
    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        archive & mGroundSecondElectrode;
        archive & mStartTime;
        archive & mEndTime;
        archive & mAreActive;
        archive & mpMesh;

        (*ProcessSpecificArchive<Archive>::Get()) & mpBoundaryConditionsContainer;
        assert(mpBoundaryConditionsContainer);
        // Allow the new object to load itself from the archive.
        mpBoundaryConditionsContainer->LoadFromArchive(*ProcessSpecificArchive<Archive>::Get(), mpMesh);
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()

    /**
     * Private default constructor for archiving only.
     */
    Electrodes(){};

public:

    /** Constructor.
     * Needs only a reference to a mesh.
     * All other parameters from the HeartConfig class
     *
     *  @param rMesh The mesh, assumed to be a cuboid.
     */
    Electrodes(AbstractTetrahedralMesh<DIM,DIM>& rMesh); // implemented in cpp

    /**
     *  @return the boundary conditions container in which is set up the Neumann
     *  fluxes for the first electrode, and the opposite fluxes for the second
     *  electrode if the the second electrode isn't grounded
     */
    boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,2> > GetBoundaryConditionsContainer();

    /**
     *  @return whether it is time to switch off the electrodes yet. THIS ONLY RETURNS
     *  TRUE ONCE - the first appropriate time. After that the electrodes assume
     *  they have been switched off and therefore this returns false.
     *
     * @param time  the current time
     */
    bool SwitchOff(double time);

    /**
     *  @return whether it is time to switch on the electrodes yet. THIS ONLY RETURNS
     *  TRUE ONCE - the first appropriate time. After that the electrodes assume
     *  they have been switched on and therefore this returns false.
     *
     * @param time  the current time
     */
    bool SwitchOn(double time);

    /** @return the time the electrodes are switched on */
    double GetSwitchOnTime()
    {
        return mStartTime;
    }

    /** @return the time the electrodes are switched off */
    double GetSwitchOffTime()
    {
        return mEndTime;
    }

    /** @return whether the second electrode is grounded or not */
    bool HasGroundedElectrode()
    {
        return mGroundSecondElectrode;
    }
};

#include "SerializationExportWrapper.hpp"
EXPORT_TEMPLATE_CLASS_SAME_DIMS(Electrodes)

#endif /*ELECTRODES_HPP_*/

max_stars_repo_name: SoftMatterMechanics/ApicalStressFibers

max_issues_repo_name: SoftMatterMechanics/ApicalStressFibers

max_forks_repo_name: SoftMatterMechanics/ApicalStressFibers


text:
#ifndef LEXER_TEST_HPP
#define LEXER_TEST_HPP

#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>

#include <cif/lexer.hpp>
#include "lexer_test_parser.hpp"
#include "lexer_file_keywords.hpp"
#include "lexer_file_iterator.hpp"

#include <boost/spirit/include/support_istream_iterator.hpp>
//#include <boost/spirit/include/support_multi_pass.hpp>

#include <boost/lexical_cast.hpp>

#include <limits>
#include <iostream>
#include <fstream>
#include <vector>
#include <utility>

using namespace boost::unit_test;

typedef std::pair<std::string,std::string> strPair;

//convert function with overload to accomodate non_type
template <typename type_t>
type_t convertVal(const std::string& str){
	return boost::lexical_cast<type_t>(str);
}
template <>
cif::detail::non_type
convertVal<cif::detail::non_type>(const std::string& str){
	return cif::detail::non_type();
}

template <typename lexer_tp>
class	TestFunc {
	public:
		template <typename A, typename B>
		struct result{ typedef void type; };
		//template <typename A> struct result<A,A>{ typedef void type; };

		typedef lexer_tp lexer_t;
		typedef typename lexer_t::non_t non_t;
		typedef typename lexer_t::int_t int_t;
		typedef typename lexer_t::float_t float_t;
		typedef typename lexer_t::string_t string_t;

		TestFunc(std::istream& is) :
			it(is),
			endIt(),
			typeMap(cif::make_lexer_type_map<lexer_t>())
		{}

		template <typename type_t>
		void operator()(const std::string& name, const type_t& val){
			if(is_end()){
				BOOST_CHECK(false);
				return;
			}

			//name
			BOOST_CHECK_EQUAL(name, it->name);

			//type
			BOOST_CHECK_EQUAL(boost::fusion::at_key<type_t>(typeMap), it->type);

			//value
			bool goodCast = true;
			type_t outVal;
			try{
				outVal = convertVal<type_t>(it->value);
			}
			catch(...){
				goodCast = false;
				BOOST_ERROR("Error casting '" + it->value + "' to " +
				            boost::fusion::at_key<type_t>(typeMap));
			}
			if(goodCast){
				comp(val, outVal);
			}

			++it;
		}

		bool is_end(void) const { return it == endIt; }

	private:
		//value comparisons w/ overload template to catch bad cases.
		static void comp(const int_t& v1, const int_t& v2){
			BOOST_CHECK_EQUAL(v1, v2);
		}
		static void comp(const float_t& v1, const float_t& v2){
			BOOST_CHECK_CLOSE(v1, v2, std::numeric_limits<float_t>::epsilon());
		}
		static void comp(const string_t& v1, const string_t& v2){
			BOOST_CHECK_EQUAL(v1, v2);
		}
		static void comp(const non_t&, const non_t&){
			BOOST_CHECK(true);
		}
		template <typename t1, typename t2>
		static void comp(const t1&, const t2&){
			BOOST_CHECK(false);
		}

		cif::lexer_file_iterator it, endIt;
		typename cif::lexer_type_map<lexer_t>::type typeMap;
};

void lexer_test(const strPair& sp){
	typedef cif::lexer<boost::spirit::istream_iterator> lexer_t;
	typedef cif::lexer_test_parser<lexer_t::iterator_type> parser_t;

	std::ifstream outFH(sp.second.c_str());
	if(!outFH.is_open()){
		BOOST_FAIL("Could not open test file '" + sp.second + "'");
		return;
	}
	TestFunc<lexer_t> testFunc(outFH);

	std::ifstream inFH(sp.first.c_str());
	if(!inFH.is_open()){
		BOOST_FAIL("Could not open test file '" + sp.first + "'");
		return;
	}
	inFH.unsetf(std::ios::skipws);

	lexer_t lexer;
	parser_t parser(lexer, testFunc);

	boost::spirit::istream_iterator it(inFH), endIt;
	bool r = tokenize_and_parse(it, endIt, lexer, parser);

	if(!r || it != endIt){
		BOOST_ERROR("Test case '"+sp.first+"' vs. '"+sp.second+
		            "'gave a truncated output");
	}
}

#endif // LEXER_TEST_HPP

ext: hpp

max_stars_repo_name: academicRobot/mmstructlib

max_issues_repo_name: academicRobot/mmstructlib

max_forks_repo_name: academicRobot/mmstructlib


text:
extern "C" {
	#include "photospline/splinetable.h"
	#include "photospline/bspline.h"
}

#include <I3Test.h>
#include <boost/filesystem.hpp>
#include <sys/time.h>
#include <limits>

namespace fs = boost::filesystem;

struct TableSet {
	fs::path abs, prob;
};

static void
splinetable_destructor(struct splinetable *table) {
	if (!table) return;

	splinetable_free(table);
	delete table;

}

static TableSet
get_splinetables()
{
	ENSURE(getenv("I3_TESTDATA") != NULL,
	    "I3_TESTDATA must be defined in the parent shell.");

	const std::string I3_TESTDATA(getenv("I3_TESTDATA"));

	fs::path abs_table(I3_TESTDATA + "/photospline/ems_z0_a0.pt.abs.fits");
	fs::path prob_table(I3_TESTDATA + "/photospline/ems_z0_a0.pt.prob.fits");

	ENSURE(fs::exists(abs_table), "Amplitude table does not exist.");
	ENSURE(fs::exists(prob_table), "Quantile table does not exists.");

	TableSet tabset;
	tabset.abs = abs_table;
	tabset.prob = prob_table;

	return tabset;
}

static boost::shared_ptr<struct splinetable>
load_splinetable(fs::path &fname)
{
	boost::shared_ptr<struct splinetable> table(new struct splinetable, splinetable_destructor);
	ENSURE(readsplinefitstable(fname.string().c_str(), table.get()) == 0, "Table can be read.");

	return table;
}

static void
compare_tables(struct splinetable *oldtable, struct splinetable *newtable)
{
	ENSURE_EQUAL(oldtable->ndim, newtable->ndim, "Number of dimensions match");
	size_t arraysize = 1;
	for (int i=0; i < oldtable->ndim; i++) {
		ENSURE_EQUAL(oldtable->naxes[i], newtable->naxes[i], "Coefficient array sizes match");
		arraysize *= oldtable->naxes[i];
		ENSURE_EQUAL(oldtable->order[i], newtable->order[i], "Spline orders match");
		ENSURE_EQUAL(oldtable->periods[i], newtable->periods[i], "Spline periods match");
		ENSURE_EQUAL(oldtable->nknots[i], newtable->nknots[i], "Size of knot vectors match");
		for (int j=0; j < oldtable->nknots[i]; j++)
			ENSURE_EQUAL(oldtable->knots[i][j], newtable->knots[i][j], "Knot positions match");
		ENSURE_EQUAL(oldtable->extents[i][0], newtable->extents[i][0], "Lower extents match");
		ENSURE_EQUAL(oldtable->extents[i][1], newtable->extents[i][1], "Upper extents match");
	}

	for (size_t i=0; i < arraysize; i++)
		ENSURE_EQUAL(oldtable->coefficients[i], newtable->coefficients[i], "Coefficients match");

	ENSURE_EQUAL(oldtable->naux, newtable->naux, "Number of auxiliary keywords match");
	for (int i=0; i < oldtable->naux; i++) {
		ENSURE_EQUAL(std::string(oldtable->aux[i][0]), std::string(newtable->aux[i][0]), "Keys match");
		ENSURE_EQUAL(std::string(oldtable->aux[i][1]), std::string(newtable->aux[i][1]), "Values match");
	}
}

TEST_GROUP(BoundaryIssues);

TEST(DiskFile)
{
	TableSet tables = get_splinetables();
	boost::shared_ptr<struct splinetable> oldtable = load_splinetable(tables.abs);

	fs::path tmp("photospline-write-test.fits");
	if (fs::exists(tmp))
		fs::remove(tmp);
	ENSURE_EQUAL(writesplinefitstable(tmp.string().c_str(), oldtable.get()), 0, "Table can be written");

	boost::shared_ptr<struct splinetable> newtable = load_splinetable(tmp);

	compare_tables(oldtable.get(), newtable.get());

	fs::remove(tmp);
}

TEST(MemoryFile)
{
	TableSet tables = get_splinetables();
	boost::shared_ptr<struct splinetable> oldtable = load_splinetable(tables.abs);

	struct splinetable_buffer buf;
	buf.mem_alloc = &malloc;
	buf.mem_realloc = &realloc;
	ENSURE_EQUAL(writesplinefitstable_mem(&buf, oldtable.get()), 0, "Table can be written");

	boost::shared_ptr<struct splinetable> newtable(new struct splinetable, splinetable_destructor);
	ENSURE_EQUAL(readsplinefitstable_mem(&buf, newtable.get()), 0, "Table can be read.");
	free(buf.data);

	compare_tables(oldtable.get(), newtable.get());
}

/*
 * Check that analytic convolution works and preserves the monotonicity
 * of the arrival-time CDF.
 */
TEST(Convolution)
{
	TableSet tables = get_splinetables();
	boost::shared_ptr<struct splinetable> table = load_splinetable(tables.prob);
	const unsigned time_dim = 3;
	double sigma = 10.0;
	int n_knots = 3;
	double knots[3] = {-2*sigma, 0, 2*sigma};
	unsigned i;
	unsigned ndim = table->ndim;
	double tablecoords[ndim], base, nudge, q0_raw, q0_conv, q1_raw, q1_conv;
	const double eps = std::numeric_limits<double>::epsilon();
	const double tmax = table->extents[time_dim][1];
	int centers[ndim];

	for (i = 0; i < ndim; i++) {
		double low, high;
		low = table->extents[i][0];
		high = table->extents[i][1];
		tablecoords[i] = (low+high)/2.0;
	}

	tablecoords[time_dim] = table->knots[time_dim][0]+eps;
	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
	q0_raw = ndsplineeval(table.get(), tablecoords, centers, 0);
	tablecoords[time_dim] = table->extents[time_dim][1];
	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
	q1_raw = ndsplineeval(table.get(), tablecoords, centers, 0);

	ENSURE_DISTANCE(q1_raw-q0_raw, 1, 1e-3,
	    "Arrival-time CDF is normalized to within 0.1%%.");

	base = q0_raw;
	for (i = 1; i+table->order[time_dim]+1 < table->nknots[time_dim]; i++) {
		tablecoords[time_dim] = table->knots[time_dim][i];
		ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
		nudge = ndsplineeval(table.get(), tablecoords, centers, 0);
		ENSURE(nudge - base >= 0,
		    "Arrival-time CDF is monotonic.");
		base = nudge;
	}

	int err = splinetable_convolve(table.get(), time_dim, knots, n_knots);
	ENSURE_EQUAL(err, 0, "Convolution succeeds.");

	tablecoords[time_dim] = table->knots[time_dim][0]*(1-eps);
	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
	q0_conv = ndsplineeval(table.get(), tablecoords, centers, 0);

	base = q0_conv;
	for (i = 1; i < table->nknots[time_dim]; i++) {
		tablecoords[time_dim] = table->knots[time_dim][i];
		if (tablecoords[time_dim] > table->extents[time_dim][1])
			break;
		ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
		nudge = ndsplineeval(table.get(), tablecoords, centers, 0);
		ENSURE(nudge - base >= 0,
		    "Arrival-time CDF remains monotonic.");
		base = nudge;
	}

	tablecoords[time_dim] = table->extents[time_dim][1];
	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0);
	q1_conv = ndsplineeval(table.get(), tablecoords, centers, 0);

	ENSURE(table->extents[time_dim][1] < tmax,
	    "t_max is slightly smaller after convolution.");
	ENSURE(q1_conv - q1_raw > -10*std::numeric_limits<float>::epsilon(),
	    "Maximum quantile (at the end of support) is not diminished.");

	ENSURE_DISTANCE(q1_raw-q0_raw, q1_conv-q0_conv, 10*std::numeric_limits<float>::epsilon(),
	    "Arrival-time CDF remains normalized after convolution.");
}

/*
 * The quantiles in the time dimension are continous. This can only be true
 * if the spline evaluation code is capable of handling points near the
 * edges of the knot fields that are supported by fewer than (order+1)
 * splines.
 */
TEST(QuantileContinuity)
{
	const unsigned time_dim = 3;
	TableSet tables = get_splinetables();
	boost::shared_ptr<struct splinetable> table = load_splinetable(tables.prob);

	unsigned i;
	unsigned ndim = table->ndim;
	double tablecoords[ndim], base, nudge;
	const double eps = std::numeric_limits<double>::epsilon();
	int centers[ndim];

	for (i = 0; i < ndim; i++) {
		double low, high;
		low = table->extents[i][0];
		high = table->extents[i][1];
		tablecoords[i] = (low+high)/2.0;
	}

	/* Check the transition into the knot field from the left. */
	tablecoords[time_dim] = table->knots[time_dim][0];
	ENSURE(tablesearchcenters(table.get(), tablecoords, centers) != 0,
	    "tablesearchcenters() fails right at the left edge of the knot field.");

	tablecoords[time_dim] += eps;

	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0,
	    "tablesearchcenters() succeeds just inside the left edge of the knot field.");
	ENSURE_EQUAL(centers[time_dim], table->order[time_dim],
	    "centers[time_dim] holds the first fully-supported knot index.");

	base = ndsplineeval(table.get(), tablecoords, centers, 0);

	ENSURE(base >= 0, "Base quantile is positive.");
	ENSURE_DISTANCE(0, base, std::numeric_limits<float>::epsilon(),
	    "Time quantile is continuous "
	    "(slope is less than FLOAT_EPSILON/DBL_EPSILON)");

	/*
	 * Now, step over the intertior knots, checking for continuity
	 * at every knot crossing.
	 */
	for (i = 1; i+1 < table->nknots[time_dim]; i++) {
		tablecoords[time_dim] = table->knots[time_dim][i]*(1-eps);
		ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0,
		    "tablesearchcenters() succeeds inside the knot field.");

		base = ndsplineeval(table.get(), tablecoords, centers, 0);
		ENSURE(base >= 0);

		tablecoords[time_dim] = table->knots[time_dim][i];
		ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0,
		    "tablesearchcenters() succeeds inside the knot field.");

		nudge = ndsplineeval(table.get(), tablecoords, centers, 0);
		ENSURE(nudge >= 0);

		ENSURE_DISTANCE(base, nudge, std::numeric_limits<float>::epsilon(),
		    "Time quantile is continuous "
		    "(slope is less than FLOAT_EPSILON/DBL_EPSILON)");

		base = nudge;

		tablecoords[time_dim] = table->knots[time_dim][i]*(1+eps);
		ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0,
		    "tablesearchcenters() succeeds inside the knot field.");

		nudge = ndsplineeval(table.get(), tablecoords, centers, 0);
		ENSURE(nudge >= 0);

		ENSURE_DISTANCE(base, nudge, std::numeric_limits<float>::epsilon(),
		    "Time quantile is continuous "
		    "(slope is less than FLOAT_EPSILON/DBL_EPSILON)");
	}

	/* Check the transition into the knot field from the right */
	tablecoords[time_dim] = table->knots[time_dim][table->nknots[time_dim]-1]*(1+eps);
	ENSURE(tablecoords[time_dim] > table->knots[time_dim][table->nknots[time_dim]-1]);
	ENSURE(tablesearchcenters(table.get(), tablecoords, centers) != 0,
	    "tablesearchcenters() fails right at the right edge of the knot field.");

	tablecoords[time_dim] = table->knots[time_dim][table->nknots[time_dim]-1];

	ENSURE_EQUAL(tablesearchcenters(table.get(), tablecoords, centers), 0,
	    "tablesearchcenters() succeeds just inside the right edge of the knot field.");
	ENSURE_EQUAL(centers[time_dim],
	    table->nknots[time_dim]-table->order[time_dim]-2,
	    "centers[time_dim] holds the first fully-supported knot index.");

	base = ndsplineeval(table.get(), tablecoords, centers, 0);

	ENSURE(base >= 0, "Base quantile is positive.");
	ENSURE_DISTANCE(0, base, std::numeric_limits<float>::epsilon(),
	    "Time quantile is continuous "
	    "(slope is less than FLOAT_EPSILON/DBL_EPSILON)");
}

TEST(ndssplineeval_vs_ndssplineeval_gradient)
{
	srand(42);

	TableSet tables = get_splinetables();
	boost::shared_ptr<struct splinetable> table = load_splinetable(tables.prob);

	const int ndim = table->ndim;
	ENSURE(ndim < 6);
	for (int i=0; i < 1000; i++) {
		/* Generate coordinates uniformly in the support of the spline */
		double x[ndim];
		int centers[ndim];
		for (int j=0; j < ndim; j++) {
			double p = double(rand())/double(RAND_MAX);
			x[j] = table->extents[j][0] + p*(table->extents[j][1]-table->extents[j][0]);
		}
		ENSURE_EQUAL(tablesearchcenters(table.get(), x, centers), 0);
		double evaluate;
		double evaluate_with_gradient[ndim+1];
		evaluate = ndsplineeval(table.get(), x, centers, 0);
		ndsplineeval_gradient(table.get(), x, centers, evaluate_with_gradient);
		ENSURE_EQUAL(evaluate, evaluate_with_gradient[0],
		    "ndsplineeval() and ndssplineeval_gradient() yield identical evaluates");
		for (int j=0; j < ndim; j++) {
			evaluate = ndsplineeval(table.get(), x, centers, 1 << j);
			ENSURE_EQUAL(evaluate, evaluate_with_gradient[1+j],
			    "ndsplineeval() and ndssplineeval_gradient() yield identical derivatives");
		}
	}
}

/*
 * bsplvb_simple() can be made to return sensical values anywhere
 * in the knot field.
 */

TEST(bsplvb_simple_vs_bspline)
{
	unsigned i;
	const unsigned n_knots = 10;
	const int order = 2;
	double x, *knots;
	int center, offset;
	float localbasis_bsplvb[order+1], localbasis_bspline[order+1];
	struct timeval thetime;
	std::vector<double> knotvec;

	/* Seed our crappy little PSRNG. */
	gettimeofday(&thetime, NULL);
	srand(thetime.tv_sec);
	/* Generate a random knot field. */
	for (i = 0; i < n_knots; i++)
		knotvec.push_back(double(rand())/double(RAND_MAX));
	std::sort(knotvec.begin(), knotvec.end());
	knots = &knotvec.front();

	/* Before the first fully-supported knot. */
	for (i = 0; i < order+1; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = order; /* First fully-supported spline. */

		ENSURE(int(i) <= center);

		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order+1 /* jhigh */, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset],
			std::numeric_limits<double>::epsilon());
		}
	}

	/* Within the support. */
	for (i = order+1; i < n_knots-order-1; i++) {
		x = (knots[i]+knots[i+1])/2.0;
		center = i;

		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order+1 /* jhigh */, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset],
			std::numeric_limits<double>::epsilon());
		}
	}

	/* After the last first fully-supported knot. */
	for (i = n_knots-order-1; i < n_knots-2; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = n_knots-order-2; /* Last fully-supported spline. */

		ENSURE(int(i) >= center);

		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order+1 /* jhigh */, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset],
			std::numeric_limits<double>::epsilon());
		}
	}
}

/*
 * bspline_deriv_nonzero() can be made to return sensical values anywhere
 * in the knot field.
 */

TEST(bspline_deriv_nonzero_vs_bspline_deriv)
{
	unsigned i;
	const unsigned n_knots = 10;
	const int order = 2;
	double x, *knots;
	int center, offset;
	float localbasis_bsplvb[order+1], localbasis_bspline[order+1];
	struct timeval thetime;
	std::vector<double> knotvec;
	/* This calculation is less stable. */
	double tol = 10*std::numeric_limits<float>::epsilon();

	/* Seed our crappy little PSRNG. */
	gettimeofday(&thetime, NULL);
	srand(thetime.tv_sec);
	/* Generate a random knot field. */
	for (i = 0; i < n_knots; i++)
		knotvec.push_back(double(rand())/double(RAND_MAX));
	std::sort(knotvec.begin(), knotvec.end());
	knots = &knotvec.front();

	/* Before the first fully-supported knot. */
	for (i = 0; i < order+1; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = order; /* First fully-supported spline. */

		ENSURE(int(i) <= center);

		bspline_deriv_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
		}
	}

	/* Within the support. */
	for (i = order+1; i < n_knots-order-1; i++) {
		x = (knots[i]+knots[i+1])/2.0;
		center = i;

		bspline_deriv_nonzero(knots, n_knots, x, center /* left */,
		    order, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
		}
	}

	/* After the last first fully-supported knot. */
	for (i = n_knots-order-1; i < n_knots-2; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = n_knots-order-2; /* Last fully-supported spline. */

		ENSURE(int(i) >= center);

		bspline_deriv_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bsplvb /* biatx */);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
		}
	}
}

/*
 * bspline_nonzero() can be made to return sensical values anywhere
 * in the knot field.
 */

TEST(bspline_nonzero_vs_bspline)
{
	unsigned i;
	const unsigned n_knots = 10;
	const int order = 2;
	double x, *knots;
	int center, offset;
	float localbasis_bsplvb[order+1], localbasis_bspline[order+1],
	    localbasis_bsplvb_deriv[order+1], localbasis_bspline_deriv[order+1];
	struct timeval thetime;
	std::vector<double> knotvec;
	/* This calculation is less stable. */
	double tol = 10*std::numeric_limits<float>::epsilon();

	/* Seed our crappy little PSRNG. */
	gettimeofday(&thetime, NULL);
	srand(thetime.tv_sec);
	/* Generate a random knot field. */
	for (i = 0; i < n_knots; i++)
		knotvec.push_back(double(rand())/double(RAND_MAX));
	std::sort(knotvec.begin(), knotvec.end());
	knots = &knotvec.front();

	/* Before the first fully-supported knot. */
	for (i = 0; i < order+1; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = order; /* First fully-supported spline. */

		ENSURE(int(i) <= center);

		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bsplvb, localbasis_bsplvb_deriv);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
			localbasis_bspline_deriv[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
			ENSURE_DISTANCE(localbasis_bspline_deriv[offset], localbasis_bsplvb_deriv[offset], tol);

		}
	}

	/* Within the support. */
	for (i = order+1; i < n_knots-order-1; i++) {
		x = (knots[i]+knots[i+1])/2.0;
		center = i;

		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bsplvb, localbasis_bsplvb_deriv);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
			localbasis_bspline_deriv[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
			ENSURE_DISTANCE(localbasis_bspline_deriv[offset], localbasis_bsplvb_deriv[offset], tol);

		}
	}

	/* After the last first fully-supported knot. */
	for (i = n_knots-order-1; i < n_knots-2; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = n_knots-order-2; /* Last fully-supported spline. */

		ENSURE(int(i) >= center);

		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bsplvb, localbasis_bsplvb_deriv);

		for (offset = -order; offset <= 0; offset++) {
			ENSURE(offset+order >= 0);
			localbasis_bspline[offset+order] =
			    bspline(knots, x, center + offset, order);
			localbasis_bspline_deriv[offset+order] =
			    bspline_deriv(knots, x, center + offset, order, 1);
		}

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_DISTANCE(localbasis_bspline[offset], localbasis_bsplvb[offset], tol);
			ENSURE_DISTANCE(localbasis_bspline_deriv[offset], localbasis_bsplvb_deriv[offset], tol);

		}
	}
}

/*
 * bspline_nonzero() gives the same result as bsplvb_simple(), ergo the
 * value bases used in ndsplineeval() and ndsplineeval_gradient() are
 * identical. Roundoff error in the derivative basis calculation
 * scales with the spline order, so we use a large number here.
 */

TEST(single_basis_vs_multi)
{
	unsigned i;
	const unsigned n_knots = 100;
	const int order = 5;
	double x, *knots;
	int center, offset;
	float localbasis_bspline_nonzero[order+1], localbasis_bspline_nonzero_deriv[order+1],
	    localbasis_bsplvb_simple[order+1], localbasis_bspline_deriv_nonzero[order+1];
	// bsplvb() may access up to *order* elements off either end
	// Pad accordingly.
	std::vector<double> knotvec(n_knots + 2*order, 0.);

	srand(0);
	/* Generate a random knot field. */
	for (i = 0; i < n_knots; i++)
		knotvec[i+order] = (double(rand())/double(RAND_MAX));
	std::sort(knotvec.begin()+order, knotvec.end()-order);
	knots = &knotvec[order];

	/* Before the first fully-supported knot. */
	for (i = 0; i < order+1; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = order; /* First fully-supported spline. */

		ENSURE(int(i) <= center);

		/* As used in ndssplineeval_gradient() */
		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bspline_nonzero, localbasis_bspline_nonzero_deriv);
		/* As used in ndsplineeval() */
		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order + 1 /* jhigh */, localbasis_bsplvb_simple);
		bspline_deriv_nonzero(knots, n_knots, x, center,
		    order, localbasis_bspline_deriv_nonzero);

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_EQUAL(localbasis_bspline_nonzero[offset], localbasis_bsplvb_simple[offset]);
			ENSURE_EQUAL(localbasis_bspline_nonzero_deriv[offset], localbasis_bspline_deriv_nonzero[offset]);

		}

	}

	/* Within the support. */
	for (i = order+1; i < n_knots-order-1; i++) {
		x = (knots[i]+knots[i+1])/2.0;
		center = i;

		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bspline_nonzero, localbasis_bspline_nonzero_deriv);
		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order + 1 /* jhigh */, localbasis_bsplvb_simple);
		bspline_deriv_nonzero(knots, n_knots, x, center,
		    order, localbasis_bspline_deriv_nonzero);

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_EQUAL(localbasis_bspline_nonzero[offset], localbasis_bsplvb_simple[offset]);
			ENSURE_EQUAL(localbasis_bspline_nonzero_deriv[offset], localbasis_bspline_deriv_nonzero[offset]);
		}
	}

	/* After the last first fully-supported knot. */
	for (i = n_knots-order-1; i < n_knots-2; i++) {
		x = (knots[i]+knots[i+1])/2.0; /* Less than fully-supported */
		center = n_knots-order-2; /* Last fully-supported spline. */

		ENSURE(int(i) >= center);

		bspline_nonzero(knots, n_knots, x, center /* left */,
		    order /* jhigh */, localbasis_bspline_nonzero, localbasis_bspline_nonzero_deriv);
		bsplvb_simple(knots, n_knots, x, center /* left */,
		    order + 1 /* jhigh */, localbasis_bsplvb_simple);
		bspline_deriv_nonzero(knots, n_knots, x, center,
		    order, localbasis_bspline_deriv_nonzero);

		for (offset = 0; offset < order+1; offset++) {
			ENSURE_EQUAL(localbasis_bspline_nonzero[offset], localbasis_bsplvb_simple[offset]);
			ENSURE_EQUAL(localbasis_bspline_nonzero_deriv[offset], localbasis_bspline_deriv_nonzero[offset]);
		}
	}
}

ext: cxx


text:
// Andrew Naplavkov

#ifndef BRIG_DETAIL_PAGE_HPP
#define BRIG_DETAIL_PAGE_HPP

#include <algorithm>
#include <array>
#include <boost/utility.hpp>
#include <brig/global.hpp>
#include <brig/rowset.hpp>
#include <brig/variant.hpp>
#include <vector>

namespace brig { namespace detail {

class page : ::boost::noncopyable { // ::boost::circular_buffer is slow
  std::array<std::vector<variant>, PageSize + 1> m_rows;
  size_t m_beg, m_end;

  static size_t next(size_t pos)  { return pos < PageSize? pos + 1: 0; }

public:
  page() : m_beg(0), m_end(0)  {}
  void swap(page& r);
  bool full() const  { return m_beg == next(m_end); }

  void clear()  { m_beg = m_end = 0; }
  bool empty() const  { return m_beg == m_end; }
  bool fetch(std::vector<variant>& row);
  void fill(rowset* rs);
}; // page

inline void page::swap(page& r)
{
  m_rows.swap(r.m_rows);
  std::swap(m_beg, r.m_beg);
  std::swap(m_end, r.m_end);
}

inline bool page::fetch(std::vector<variant>& row)
{
  if (m_beg == m_end) return false;
  row.swap(m_rows[m_beg]);
  m_beg = next(m_beg);
  return true;
}

inline void page::fill(rowset* rs)
{
  while (true)
  {
    const size_t end(next(m_end));
    if (m_beg == end || !rs || !rs->fetch(m_rows[m_end])) return;
    m_end = end;
  }
} // page::

} } // brig::detail

#endif // BRIG_DETAIL_PAGE_HPP

max_stars_repo_name: isadchenko/brig

max_issues_repo_name: isadchenko/brig

max_forks_repo_name: isadchenko/brig


text: #include <boost/spirit/home/x3/support/context.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text: #include <boost/mpl/aux_/config/static_constant.hpp>


text:
/*******************************************************************************
 *
 * MIT License
 *
 * Copyright (c) 2017 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 *******************************************************************************/

#include <miopen/solver.hpp>

#include <miopen/algorithm.hpp>
#include <miopen/conv/data_invoke_params.hpp>
#include <miopen/conv_solution.hpp>
#include <miopen/convolution_fft.hpp>
#include <miopen/env.hpp>
#include <miopen/tensor.hpp>

#include <boost/any.hpp>

namespace miopen {
namespace solver {

MIOPEN_DECLARE_ENV_VAR(MIOPEN_DEBUG_CONV_FFT)

static void cgemm_grid(size_t* global_work_size,
                       size_t* local_work_size,
                       int cgemm_choice,
                       const int N,
                       const int out_c,
                       const int out_n)
{
    unsigned int threadTile[2];
    unsigned int groupSize[2];

    // grid for cgemm
    if(cgemm_choice == 1)
    {
        threadTile[0] = 4;
        threadTile[1] = 4;

        groupSize[0] = 16;
        groupSize[1] = 16;

        local_work_size[0] = 16;
        local_work_size[1] = 16;
    }
    else if(cgemm_choice == 2)
    {
        threadTile[0] = 4;
        threadTile[1] = 4;

        groupSize[0] = 4;
        groupSize[1] = 4;

        local_work_size[0] = 64;
        local_work_size[1] = 1;
    }
    else
    {
        threadTile[0] = 2;
        threadTile[1] = 2;

        groupSize[0] = 4;
        groupSize[1] = 4;

        local_work_size[0] = 64;
        local_work_size[1] = 1;
    }

    global_work_size[2] = 1;
    global_work_size[2] *= N;

    unsigned int sizeOfC0         = out_c;
    unsigned int sizeOfC1         = out_n;
    auto macroTile0               = static_cast<unsigned int>(groupSize[0] * threadTile[0]);
    auto macroTile1               = static_cast<unsigned int>(groupSize[1] * threadTile[1]);
    unsigned int totalWorkGroups0 = sizeOfC0 / macroTile0;
    unsigned int totalWorkGroups1 = sizeOfC1 / macroTile1;
    // b/c single kernel, add extra work-group here if edge needed
    if(totalWorkGroups0 * macroTile0 < sizeOfC0)
    {
        totalWorkGroups0++;
    }
    if(totalWorkGroups1 * macroTile1 < sizeOfC1)
    {
        totalWorkGroups1++;
    }
    global_work_size[0] = totalWorkGroups0 * local_work_size[0];
    global_work_size[1] = totalWorkGroups1 * local_work_size[1];
}

bool fft::IsApplicable(const ConvolutionContext& ctx) const
{
    // disable running any FFT based convolutions by checking this env variable
    if(ctx.direction.IsBackwardWrW() || !ctx.conv_problem.IsFp32())
        return false;

    if(!ctx.IsLayoutDefault())
    {
        return false;
    }

    const auto is_fwd    = ctx.direction.IsForward();
    decltype(auto) conv  = ctx.conv_problem.GetConv();
    decltype(auto) xDesc = is_fwd ? ctx.conv_problem.GetIn() : ctx.conv_problem.GetOut();
    decltype(auto) yDesc = is_fwd ? ctx.conv_problem.GetOut() : ctx.conv_problem.GetIn();
    decltype(auto) wDesc = ctx.conv_problem.GetWeights();

    if(conv.GetSpatialDimension() != 2 || conv.group_count != 1 ||
       !miopen::all_of(conv.GetConvDilations(), [](auto v) { return v == 1; }))
        return false;

    int in_n, in_c, in_h, in_w;
    int out_n, out_c, out_h, out_w;
    int wei_k, wei_c, wei_h, wei_w;
    std::tie(in_n, in_c, in_h, in_w)     = miopen::tien<4>(xDesc.GetLengths());
    std::tie(out_n, out_c, out_h, out_w) = miopen::tien<4>(yDesc.GetLengths());
    std::tie(wei_k, wei_c, wei_h, wei_w) = miopen::tien<4>(wDesc.GetLengths());

    // FFT convolutions only works for specific config(s)
    // coverage to expand gradually

    if((in_n < 1) || (in_n > 512) || (wei_k < 1) || (wei_k > 512) || ((in_c * in_n) % 16 != 0) ||
       (wei_c * wei_k) % 16 != 0 || (out_c * out_n) % 16 != 0)
        return false;

    if((std::tie(in_h, in_w) != std::make_tuple(28, 28)) &&
       (std::tie(in_h, in_w) != std::make_tuple(27, 27)) &&
       (std::tie(in_h, in_w) != std::make_tuple(14, 14)) &&
       (std::tie(in_h, in_w) != std::make_tuple(7, 7)))
        return false;

    const auto cparam = std::make_tuple(conv.GetConvPads()[0],
                                        conv.GetConvPads()[1],
                                        conv.GetConvStrides()[0],
                                        conv.GetConvStrides()[1]);

    return std::tie(wei_h, wei_w) == std::make_tuple(5, 5) && cparam == std::make_tuple(2, 2, 1, 1);
}

size_t fft::GetWorkspaceSize(const ConvolutionContext& ctx) const
{
    const auto fwd       = ctx.direction.IsForward();
    decltype(auto) xDesc = fwd ? ctx.conv_problem.GetIn() : ctx.conv_problem.GetOut();
    decltype(auto) yDesc = fwd ? ctx.conv_problem.GetOut() : ctx.conv_problem.GetIn();
    decltype(auto) wDesc = ctx.conv_problem.GetWeights();

    int in_n, in_c, in_h, in_w;
    std::tie(in_n, in_c, in_h, in_w) = miopen::tien<4>(xDesc.GetLengths());

    int out_n, out_c, out_h, out_w;
    std::tie(out_n, out_c, out_h, out_w) = miopen::tien<4>(yDesc.GetLengths());

    int wei_k, wei_c, wei_h, wei_w;
    std::tie(wei_k, wei_c, wei_h, wei_w) = miopen::tien<4>(wDesc.GetLengths());

    // FFT convolutions only works for specific config(s)
    // coverage to expand gradually

    const int N       = FFTConvParams::TileSize(in_h, in_w);
    const int Padding = FFTConvParams::TransposePadding;

    int temp_size = 0;

    if(fwd)
    {
        int temp_size1 = (in_c * in_n + Padding) + (wei_k * wei_c + Padding);
        int temp_size2 = (out_n * out_c + Padding);
        temp_size      = std::max(temp_size1, temp_size2);
    }
    else
    {
        int temp_size1 = (out_n * out_c + Padding) + (wei_k * wei_c + Padding);
        int temp_size2 = (in_c * in_n + Padding);
        temp_size      = std::max(temp_size1, temp_size2);
    }

    return 2 * 2 * N * temp_size * sizeof(float);
}

ConvSolution fft::GetSolution(const ConvolutionContext& ctx) const
{
    int in_n = ctx.batch_sz, in_c = ctx.n_inputs, in_h = ctx.in_height, in_w = ctx.in_width;
    int out_n = ctx.batch_sz, out_c = ctx.n_outputs;

    const int N          = FFTConvParams::TileSize(in_h, in_w);
    const int NumKernels = FFTConvParams::NumKernels;

    size_t global_work_size[NumKernels][3];
    size_t local_work_size[NumKernels][3];

    for(int ik = 0; ik < NumKernels; ik++)
    {
        global_work_size[ik][0] = local_work_size[ik][0] = 1;
        global_work_size[ik][1] = local_work_size[ik][1] = 1;
        global_work_size[ik][2] = local_work_size[ik][2] = 1;
    }

    // grid for FFT kernels
    if((in_h == 7) && (in_w == 7))
    {
        local_work_size[0][0]  = 192;
        global_work_size[0][0] = ((in_c * out_n) / 16) * local_work_size[0][0];

        local_work_size[1][0]  = 192;
        global_work_size[1][0] = ((in_c * out_c) / 16) * local_work_size[1][0];

        local_work_size[6][0]  = 192;
        global_work_size[6][0] = ((out_n * out_c) / 16) * local_work_size[6][0];
    }
    else if((in_h == 14) && (in_w == 14))
    {
        local_work_size[0][0]  = 128;
        global_work_size[0][0] = ((in_c * out_n) / 4) * local_work_size[0][0];

        local_work_size[1][0]  = 128;
        global_work_size[1][0] = ((in_c * out_c) / 4) * local_work_size[1][0];

        local_work_size[6][0]  = 128;
        global_work_size[6][0] = ((out_n * out_c) / 4) * local_work_size[6][0];
    }
    else
    {
        local_work_size[0][0]  = 64;
        global_work_size[0][0] = in_c * out_n * local_work_size[0][0];

        local_work_size[1][0]  = 64;
        global_work_size[1][0] = in_c * out_c * local_work_size[1][0];

        local_work_size[6][0]  = 64;
        global_work_size[6][0] = out_n * out_c * local_work_size[6][0];
    }

    // decide tranpose kernel options based on params
    int in_tranpose_choice = 0;
    int wt_tranpose_choice = 0;
    int ot_tranpose_choice = 0;

    // grid for transpose kernels
    if((in_h == 7) && (in_w == 7))
    {
        local_work_size[5][0]  = 256;
        global_work_size[5][0] = (1 + N / 16) * (out_n * out_c / 16) * local_work_size[5][0];
    }
    else if((in_h == 14) && (in_w == 14))
    {
        local_work_size[2][0]  = 256;
        global_work_size[2][0] = (1 + N / 16) * (in_c * out_n / 16) * local_work_size[2][0];

        local_work_size[3][0]  = 256;
        global_work_size[3][0] = (1 + N / 16) * (in_c * out_c / 16) * local_work_size[3][0];

        local_work_size[5][0]  = 256;
        global_work_size[5][0] = (1 + N / 16) * (out_n * out_c / 16) * local_work_size[5][0];
    }
    else
    {
        if((in_n * in_c >= 64) && ((in_n * in_c) % 32 == 0))
            in_tranpose_choice = 1;
        if((out_c * in_c >= 64) && ((out_c * in_c) % 32 == 0))
            wt_tranpose_choice = 1;
        if((out_n * out_c >= 64) && ((out_n * out_c) % 32 == 0))
            ot_tranpose_choice = 1;

        int in_tranpose_bwidth = in_tranpose_choice != 0 ? 32 : 16;
        int wt_tranpose_bwidth = wt_tranpose_choice != 0 ? 32 : 16;
        int ot_tranpose_bwidth = ot_tranpose_choice != 0 ? 32 : 16;

        local_work_size[2][0] = 256;
        global_work_size[2][0] =
            (N / in_tranpose_bwidth) * (in_c * out_n / in_tranpose_bwidth) * local_work_size[2][0];

        local_work_size[3][0] = 256;
        global_work_size[3][0] =
            (N / wt_tranpose_bwidth) * (in_c * out_c / wt_tranpose_bwidth) * local_work_size[3][0];

        local_work_size[5][0] = 256;
        global_work_size[5][0] =
            (N / ot_tranpose_bwidth) * (out_n * out_c / ot_tranpose_bwidth) * local_work_size[5][0];
    }

    // cgemm kernel options
    int cgemm_choice = 0;

    if(((in_h == 28) && (in_w == 28)) || ((in_h == 14) && (in_w == 14)) ||
       ((in_h == 7) && (in_w == 7)))
        cgemm_choice = 2;
    else if((in_h == 27) && (in_w == 27))
        cgemm_choice = 1;

    if((in_n < 16) || (in_c < 16) || (out_c < 16))
        cgemm_choice = 0;

    cgemm_grid(global_work_size[4], local_work_size[4], cgemm_choice, N, out_c, out_n);

    std::string parms;

    if(in_tranpose_choice == 0)
        parms += " -DCFF_TRANSP_IN_MOD16=1";
    if(wt_tranpose_choice == 0)
        parms += " -DCFF_TRANSP_WT_MOD16=1";
    if(ot_tranpose_choice == 0)
        parms += " -DCFF_TRANSP_OT_MOD16=1";

    switch(cgemm_choice)
    {
    case 1: parms += " -DCFF_CGEMM_CHOICE_1=1"; break;
    case 2: parms += " -DCFF_CGEMM_CHOICE_2=1"; break;
    default: parms += " -DCFF_CGEMM_CHOICE_0=1"; break;
    }

    if((in_h == 28) && (in_w == 28))
        parms += " -DCFF_IMG_SZ_28_28";
    else if((in_h == 27) && (in_w == 27))
        parms += " -DCFF_IMG_SZ_27_27";
    else if((in_h == 14) && (in_w == 14))
        parms += " -DCFF_IMG_SZ_14_14";
    else if((in_h == 7) && (in_w == 7))
        parms += " -DCFF_IMG_SZ_7_7";

    const auto workSpaceSize = GetWorkspaceSize(ctx);

    parms += " -DCFF_IMG_H=";
    parms += std::to_string(in_h);
    parms += " -DCFF_IMG_W=";
    parms += std::to_string(in_w);
    parms += " -DCFF_BATCH=";
    parms += std::to_string(in_n);
    parms += " -DCFF_NFILTER=";
    parms += std::to_string(out_c);
    parms += " -DCFF_CHANNELS=";
    parms += std::to_string(in_c);
    parms += " -DCFF_HALFW=";
    parms += std::to_string(workSpaceSize / (2 * 2 * sizeof(float)));

    if(!ctx.direction.IsForward())
    {
        parms += " -DCFF_BACKWARD";
    }

    const std::string algorithm    = "miopenConvolutionFwdAlgoFFT";
    const std::string program_name = "MIOpenConvFFT.cl";

    auto sol        = ConvSolution{miopenStatusSuccess};
    sol.workspce_sz = workSpaceSize;

    // skip front transposes for 7x7
    const auto skip_front_transposes = (in_h == 7) && (in_w == 7);

    for(int ik = 0; ik < NumKernels; ik++)
    {
        if(skip_front_transposes && ((ik == 2) || (ik == 3)))
            continue;

        const auto kernel_name = [=]() {
            switch(ik)
            {
            case 0: return "MIOpenConvFFT_fwd_in";
            case 1: return "MIOpenConvFFT_fwd_we";
            case 2: return "MIOpenConvFFT_transpose_in";
            case 3: return "MIOpenConvFFT_transpose_we";
            case 4: return "MIOpenConvFFT_cgemm";
            case 5: return "MIOpenConvFFT_transpose_out";
            case 6: return "MIOpenConvFFT_inv_out";
            default: assert(false); return "";
            }
        }();

        std::vector<size_t> vld(3);
        std::vector<size_t> vgd(3);

        vld[0] = local_work_size[ik][0];
        vld[1] = local_work_size[ik][1];
        vld[2] = local_work_size[ik][2];

        vgd[0] = global_work_size[ik][0];
        vgd[1] = global_work_size[ik][1];
        vgd[2] = global_work_size[ik][2];

        auto kernel         = KernelInfo{};
        kernel.kernel_file  = program_name;
        kernel.kernel_name  = kernel_name;
        kernel.comp_options = parms;
        kernel.g_wk         = vgd;
        kernel.l_wk         = vld;
        sol.construction_params.push_back(kernel);
    }

    sol.invoker_factory = [=](const std::vector<Kernel>& kernels) {
        int halfw = static_cast<int>(workSpaceSize) / (2 * 2 * static_cast<int>(sizeof(float)));
        const int padding = FFTConvParams::TransposePadding;

        return [=](const Handle& handle, const AnyInvokeParams& primitive_params) {
            const auto& params  = primitive_params.CastTo<conv::DataInvokeParams>();
            const auto& tensors = params.tensors;

            if(params.workSpaceSize < workSpaceSize)
                MIOPEN_THROW("Not enough workspace for FFT: expected " +
                             std::to_string(workSpaceSize) + ", got " +
                             std::to_string(params.workSpaceSize));

            float time_fft = 0;
            int kernel_id  = 0;
            for(int ik = 0; ik < NumKernels; ik++)
            {
                if(skip_front_transposes && ((ik == 2) || (ik == 3)))
                    continue;

                const auto& k = handle.Run(kernels[kernel_id++]);

                switch(ik)
                {
                case 0: k(tensors.in, params.workSpace); break;
                case 1: k(tensors.w, params.workSpace); break;
                case 4: {
                    k(params.workSpace,
                      0,
                      halfw + N * (in_n * in_c + padding),
                      halfw + 0,
                      out_c,
                      out_n * out_c + padding,
                      in_c,
                      in_c * out_c + padding,
                      in_c,
                      in_n * in_c + padding,
                      out_c,
                      in_n,
                      N,
                      in_c);
                    break;
                }
                case 6: k(params.workSpace, tensors.out); break;
                case 2:
                case 3:
                case 5: k(params.workSpace); break;
                default: assert(false);
                }

                if(handle.IsProfilingEnabled())
                    time_fft += handle.GetKernelTime();
            }

            if(handle.IsProfilingEnabled())
            {
                handle.ResetKernelTime();
                handle.AccumKernelTime(time_fft);
            }
        };
    };

    return sol;
}

} // namespace solver
} // namespace miopen

ext: cpp


text:
#ifndef BART_SRC_DOMAIN_DOMAIN_HPP_
#define BART_SRC_DOMAIN_DOMAIN_HPP_

#include <memory>
#include <vector>

#include <deal.II/base/index_set.h>
#include <deal.II/base/iterator_range.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>

#include "domain/domain_i.hpp"
#include "domain/finite_element/finite_element_i.hpp"
#include "domain/mesh/mesh_i.hpp"
#include "problem/parameter_types.hpp"
#include "system/system_types.h"
#include "utility/has_dependencies.h"

namespace bart::domain {

/*! This struct provides the type of the triangulation to use. For 2D and 3D we
 * will use the built-in distributed triangulation, which handles distribution
 * of the triangulation across processors. For the 1D case, we need to make a
 * normal triangulation and use a separate function to distribute it. This is
 * a product of the way dealii is coded.
 */
template <int dim> struct TriangulationType { using type = dealii::parallel::distributed::Triangulation<dim>; };
template <> struct TriangulationType<1> { using type = dealii::Triangulation<1>; };
template struct TriangulationType<2>;
template struct TriangulationType<3>;

/*! \brief Default implementation of a domain definition.
 *
 * See the base class for general information about domain definitions, an
 * example of use is shown below:
 *
 * \code{.cpp}
 *
 * domain::definition<3> domain_definition(std::move(mesh_ptr),
 *                                         finite_element_ptr,
 *                                         problem::DiscretizationType::kContinuousFEM);
 * domain_definition.SetUpMesh().SetUpDOF();
 *
 * auto cell_matrix = domain_definition.GetCellMatrix();
 *
 * \endcode
 *
 * \author Joshua Rehak
 * \date 2019/02
 */
template <int dim>
class Domain : public DomainI<dim>, public utility::HasDependencies {
 public:
  typedef std::vector<typename dealii::DoFHandler<dim>::active_cell_iterator> CellRange;

  /*! \brief Constructor.
   * Takes ownership of injected dependencies (MeshI and FiniteElementI) and
   * sets the type of discretization (default: continuous FEM).
   */
  Domain(std::unique_ptr<domain::mesh::MeshI<dim>> mesh,
         std::shared_ptr<domain::finite_element::FiniteElementI<dim>> finite_element,
         problem::DiscretizationType discretization = problem::DiscretizationType::kContinuousFEM);
  ~Domain() = default;

  auto SetUpDOF() -> Domain<dim>& override;
  auto SetUpMesh() -> Domain<dim>& override;
  auto SetUpMesh(const int global_refinements) -> Domain<dim>& override;

  auto GetCellMatrix() const -> dealii::FullMatrix<double> override;
  auto GetCellVector() const -> dealii::Vector<double> override;
  auto MakeSystemMatrix() const -> std::shared_ptr<system::MPISparseMatrix> override;
  auto MakeSystemVector() const -> std::shared_ptr<system::MPIVector> override;

  auto Cells() const -> CellRange override { return local_cells_; };
  auto discretization_type() const -> problem::DiscretizationType override { return discretization_type_; }
  auto total_degrees_of_freedom() const -> int override ;
  auto dof_handler() const -> const dealii::DoFHandler<dim>& override { return dof_handler_; }
  auto locally_owned_dofs() const -> dealii::IndexSet override { return locally_owned_dofs_; }

 private:

  //! Internal owned mesh object.
  std::unique_ptr<domain::mesh::MeshI<dim>> mesh_;

  //! Internal owned finite element object
  std::shared_ptr<domain::finite_element::FiniteElementI<dim>> finite_element_;

  //! Internal distributed triangulation object
  typename TriangulationType<dim>::type triangulation_;

  //! Internal DoFHandler object
  dealii::DoFHandler<dim> dof_handler_;

  //! Total degrees of freedom
  mutable int total_degrees_of_freedom_{ 0 };

  //! Index of locally owned dofs owned by the current processor
  dealii::IndexSet locally_owned_dofs_;

  /*! Index of locally relevant dofs to the current processor
   * This includes indices of dofs owned by the current processor but also
   * the indices of neighboring cells (ghost cells) that may be relevant.
   */
  dealii::IndexSet locally_relevant_dofs_;

  /*! Constraint matrix */
  dealii::AffineConstraints<double> constraint_matrix_;

  /*! Dynamic sparsity pattern for MPI matrices */
  dealii::DynamicSparsityPattern dynamic_sparsity_pattern_;

  /*! local cells */
  CellRange local_cells_;

  /*! Discretization type */
  const problem::DiscretizationType discretization_type_;
};

} // namespace bart::domain

#endif // BART_SRC_DOMAIN_DOMAIN_HPP_

max_stars_repo_name: SlaybaughLab/Transport

max_issues_repo_name: SlaybaughLab/Transport

max_forks_repo_name: SlaybaughLab/Transport


text:
/*
 * Copyright (c) 2011, Peter Thorson. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of the WebSocket++ Project nor the
 *       names of its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL PETER THORSON BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef WEBSOCKET_CLIENT_HPP
#define WEBSOCKET_CLIENT_HPP

#include <boost/asio.hpp>
#include <boost/shared_ptr.hpp>

#include <set>

namespace websocketpp {
	class client;
	typedef boost::shared_ptr<client> client_ptr;
}

#include "websocketpp.hpp"
#include "websocket_client_session.hpp"
#include "websocket_connection_handler.hpp"

using boost::asio::ip::tcp;

namespace websocketpp {

class client_error : public std::exception {
public:
	client_error(const std::string& msg)
		: m_msg(msg) {}
	~client_error() throw() {}

	virtual const char* what() const throw() {
		return m_msg.c_str();
	}
private:
	std::string m_msg;
};

class client : public boost::enable_shared_from_this<client> {
	public:
		static const uint16_t CLIENT_STATE_NULL = 0;
		static const uint16_t CLIENT_STATE_INITIALIZED = 1;
		static const uint16_t CLIENT_STATE_CONNECTING = 2;
		static const uint16_t CLIENT_STATE_CONNECTED = 3;

		client(boost::asio::io_service& io_service,
			   connection_handler_ptr defc);

		// INTERFACE FOR LOCAL APPLICATIONS

		// initializes the session. Methods that affect the opening handshake
		// such as add_protocol and set_header must be called after init and
		// before connect.
		void init();

		// starts the connection process. Should be called before
		// io_service.run(), connection process will not start until run() has
		// been called.
		void connect(const std::string& url);

		// Adds a protocol to the opening handshake.
		// Must be called before connect
		void add_subprotocol(const std::string& p);

		// Sets the value of the given HTTP header to be sent during the
		// opening handshake. Must be called before connect
		void set_header(const std::string& key,const std::string& val);

		void set_origin(const std::string& val);

		void set_max_message_size(uint64_t val);

		// Test methods determine if a message of the given level should be
		// written. elog shows all values above the level set. alog shows only
		// the values explicitly set.
		bool test_elog_level(uint16_t level);
		void set_elog_level(uint16_t level);

		bool test_alog_level(uint16_t level);
		void set_alog_level(uint16_t level);
		void unset_alog_level(uint16_t level);

		// INTERFACE FOR SESSIONS

		// Check if message size is within server's acceptable parameters
		bool validate_message_size(uint64_t val);

		// write to the server's logs
		void log(std::string msg,uint16_t level = LOG_ERROR);
		void access_log(std::string msg,uint16_t level);
	private:
		// if no errors starts the session's read loop and returns to the
		// start_accept phase.
		void handle_connect(const boost::system::error_code& error);

	private:
		uint16_t					m_elog_level;
		uint16_t					m_alog_level;

		uint16_t					m_state;

		std::set<std::string>		m_hosts;
		uint64_t					m_max_message_size;
		boost::asio::io_service&	m_io_service;
		tcp::resolver				m_resolver;
		client_session_ptr			m_client_session;
		connection_handler_ptr		m_def_con_handler;
};

}

#endif // WEBSOCKET_CLIENT_HPP

max_stars_repo_name: urbenlegend/WebStreamer

max_issues_repo_name: cqzhanghy/WebStreamer

max_forks_repo_name: cqzhanghy/WebStreamer


text:
// Copyright (c) 2019-2021 The Bigbang developers
// Distributed under the MIT/X11 software license, see the accompanying
// file COPYING or http://www.opensource.org/licenses/mit-license.php.

#include "ctsdb.h"

#include <boost/bind.hpp>

#include "leveldbeng.h"

using namespace std;
using namespace xengine;

namespace bigbang
{
namespace storage
{

//////////////////////////////
// CCTSIndex

CCTSIndex::CCTSIndex()
{
}

CCTSIndex::~CCTSIndex()
{
}

bool CCTSIndex::Initialize(const boost::filesystem::path& pathCTSDB)
{
    CLevelDBArguments args;
    args.path = (pathCTSDB / "index").string();
    args.syncwrite = false;
    CLevelDBEngine* engine = new CLevelDBEngine(args);

    if (!Open(engine))
    {
        delete engine;
        return false;
    }

    return true;
}

void CCTSIndex::Deinitialize()
{
    Close();
}

bool CCTSIndex::Update(const vector<int64>& vTime, const vector<CDiskPos>& vPos, const vector<int64>& vDel)
{
    if (vTime.size() != vPos.size())
    {
        return false;
    }

    if (!TxnBegin())
    {
        return false;
    }

    for (int i = 0; i < vTime.size(); i++)
    {
        Write(vTime[i], vPos[i]);
    }

    for (int i = 0; i < vDel.size(); i++)
    {
        Erase(vDel[i]);
    }

    if (!TxnCommit())
    {
        return false;
    }

    return true;
}

bool CCTSIndex::Retrieve(const int64 nTime, CDiskPos& pos)
{
    return Read(nTime, pos);
}

} // namespace storage
} // namespace bigbang

max_stars_repo_name: BigBang-Foundation/BigBang

max_issues_repo_name: BigBang-Foundation/BigBang

max_forks_repo_name: BigBang-Foundation/BigBang


text:
/*

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

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

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifndef TESTVERTEXBASEDCELLPOPULATIONFORCRYPT_HPP_
#define TESTVERTEXBASEDCELLPOPULATIONFORCRYPT_HPP_

#include <cxxtest/TestSuite.h>

#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "VertexBasedCellPopulation.hpp"
#include "HoneycombVertexMeshGenerator.hpp"
#include "WntCellCycleModel.hpp"
#include "AbstractCellBasedTestSuite.hpp"
#include "WildTypeCellMutationState.hpp"
#include "ApcOneHitCellMutationState.hpp"
#include "ApcTwoHitCellMutationState.hpp"
#include "BetaCateninOneHitCellMutationState.hpp"
#include "DifferentiatedCellProliferativeType.hpp"
#include "CellLabel.hpp"
#include "SmartPointers.hpp"

//This test is always run sequentially (never in parallel)
#include "FakePetscSetup.hpp"

class TestVertexBasedCellPopulationForCrypt : public AbstractCellBasedTestSuite
{
public:

    /**
     * Test that post-#878, WntConcentration copes with a VertexBasedCellPopulation.
     * \todo When vertex-based cell population code is added to cell_based folder, move this
     *       test to TestWntConcentration.hpp
     */
    void TestWntConcentrationWithVertexBasedCellPopulation()
    {
        // Make some nodes
        std::vector<Node<2>*> nodes;
        nodes.push_back(new Node<2>(0, true, 2.0, -1.0));
        nodes.push_back(new Node<2>(1, true, 2.0, 1.0));
        nodes.push_back(new Node<2>(2, true, -2.0, 1.0));
        nodes.push_back(new Node<2>(3, true, -2.0, -1.0));
        nodes.push_back(new Node<2>(4, true, 0.0, 2.0));

        // Make a rectangular element out of nodes 0,1,2,3
        std::vector<Node<2>*> nodes_elem_1;
        nodes_elem_1.push_back(nodes[0]);
        nodes_elem_1.push_back(nodes[1]);
        nodes_elem_1.push_back(nodes[2]);
        nodes_elem_1.push_back(nodes[3]);

        // Make a triangular element out of nodes 1,4,2
        std::vector<Node<2>*> nodes_elem_2;
        nodes_elem_2.push_back(nodes[1]);
        nodes_elem_2.push_back(nodes[4]);
        nodes_elem_2.push_back(nodes[2]);

        std::vector<VertexElement<2,2>*> vertex_elements;
        vertex_elements.push_back(new VertexElement<2,2>(0, nodes_elem_1));
        vertex_elements.push_back(new VertexElement<2,2>(1, nodes_elem_2));

        // Make a vertex mesh
        MutableVertexMesh<2,2> vertex_mesh(nodes, vertex_elements);

        // Create cells
        std::vector<CellPtr> cells;
        MAKE_PTR(WildTypeCellMutationState, p_state);
        MAKE_PTR(DifferentiatedCellProliferativeType, p_diff_type);
        for (unsigned i=0; i<vertex_mesh.GetNumElements(); i++)
        {
            WntCellCycleModel* p_cell_cycle_model = new WntCellCycleModel;
            p_cell_cycle_model->SetDimension(2);

            CellPtr p_cell(new Cell(p_state, p_cell_cycle_model));
            p_cell->SetCellProliferativeType(p_diff_type);
            double birth_time = 0.0 - i;
            p_cell->SetBirthTime(birth_time);
            cells.push_back(p_cell);
        }

        // Create cell population
        VertexBasedCellPopulation<2> cell_population(vertex_mesh, cells);

        // Set the top of this cell_population, for the purposes of computing the WntConcentration
        double crypt_length = 4.0;

        // Set up an instance of the WntConcentration singleton object
        WntConcentration<2>* p_wnt = WntConcentration<2>::Instance();
        TS_ASSERT_EQUALS(p_wnt->IsWntSetUp(), false);

        // Check that the singleton can be set up
        p_wnt->SetType(LINEAR);
        p_wnt->SetCellPopulation(cell_population);
        p_wnt->SetCryptLength(crypt_length);

        TS_ASSERT_EQUALS(p_wnt->IsWntSetUp(), true);

        // Check that the singleton can be destroyed then recreated
        WntConcentration<2>::Destroy();
        WntConcentration<2>::Instance()->SetType(NONE);
        WntConcentration<2>::Instance()->SetCellPopulation(cell_population);
        WntConcentration<2>::Instance()->SetCryptLength(crypt_length);

        TS_ASSERT_EQUALS(WntConcentration<2>::Instance()->IsWntSetUp(), false); // not fully set up now it is a NONE type

        WntConcentration<2>::Destroy();
        WntConcentration<2>::Instance()->SetType(LINEAR);
        WntConcentration<2>::Instance()->SetCellPopulation(cell_population);
        WntConcentration<2>::Instance()->SetCryptLength(crypt_length);

        p_wnt = WntConcentration<2>::Instance();
        TS_ASSERT_EQUALS(p_wnt->IsWntSetUp(), true); // set up again

        double wnt_at_cell0 = p_wnt->GetWntLevel(cell_population.GetCellUsingLocationIndex(0));
        double wnt_at_cell1 = p_wnt->GetWntLevel(cell_population.GetCellUsingLocationIndex(1));

        // We have set the top of the cell population to be 4, so the WntConcentration should decrease linearly
        // up the cell_population, from one at height 0 to zero at height 4.

        // Cell 0 has centre of mass (0,0)
        TS_ASSERT_DELTA(wnt_at_cell0, 1.0, 1e-4);

        // Cell 1 has centre of mass (0, 4/3)
        TS_ASSERT_DELTA(wnt_at_cell1, 2.0/3.0, 1e-4);
    }
};

#endif /*TESTVERTEXBASEDCELLPOPULATIONFORCRYPT_HPP_*/


text:
#include <boost/container/flat_map.hpp>
#include <boost/preprocessor/seq/for_each.hpp>

#include <deip/manifest/plugins.hpp>

namespace deip { namespace plugin {

std::shared_ptr< deip::app::abstract_plugin > create_account_by_key_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_account_history_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_account_statistics_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_auth_util_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_block_info_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_blockchain_statistics_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_debug_node_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_delayed_node_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_raw_block_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_witness_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_tsc_history_plugin( deip::app::application* app );

std::shared_ptr< deip::app::abstract_plugin > create_cr_history_plugin( deip::app::application* app );

std::shared_ptr<deip::app::abstract_plugin> create_eci_history_plugin(deip::app::application* app);

std::shared_ptr< deip::app::abstract_plugin > create_funding_opportunity_history_plugin( deip::app::application* app );

boost::container::flat_map< std::string, std::function< std::shared_ptr< deip::app::abstract_plugin >( deip::app::application* app ) > > plugin_factories_by_name;

void initialize_plugin_factories()
{

   plugin_factories_by_name[ "account_by_key" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_account_by_key_plugin( app );
   };

   plugin_factories_by_name[ "account_history" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_account_history_plugin( app );
   };

   plugin_factories_by_name[ "account_statistics" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_account_statistics_plugin( app );
   };

   plugin_factories_by_name[ "auth_util" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_auth_util_plugin( app );
   };

   plugin_factories_by_name[ "block_info" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_block_info_plugin( app );
   };

   plugin_factories_by_name[ "blockchain_statistics" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_blockchain_statistics_plugin( app );
   };

   plugin_factories_by_name[ "debug_node" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_debug_node_plugin( app );
   };

   plugin_factories_by_name[ "delayed_node" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_delayed_node_plugin( app );
   };

   plugin_factories_by_name[ "raw_block" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_raw_block_plugin( app );
   };

   plugin_factories_by_name[ "witness" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_witness_plugin( app );
   };

   plugin_factories_by_name[ "tsc_history" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_tsc_history_plugin( app );
   };

   plugin_factories_by_name[ "research_content_references_history" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_cr_history_plugin( app );
   };

   plugin_factories_by_name[ "eci_history" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_eci_history_plugin( app );
   };

   plugin_factories_by_name[ "fo_history" ] = []( deip::app::application* app ) -> std::shared_ptr< deip::app::abstract_plugin >
   {
      return create_funding_opportunity_history_plugin( app );
   };

}

std::shared_ptr< deip::app::abstract_plugin > create_plugin( const std::string& name, deip::app::application* app )
{
   auto it = plugin_factories_by_name.find( name );
   if( it == plugin_factories_by_name.end() )
      return std::shared_ptr< deip::app::abstract_plugin >();
   return it->second( app );
}

std::vector< std::string > get_available_plugins()
{
   std::vector< std::string > result;
   for( const auto& e : plugin_factories_by_name )
      result.push_back( e.first );
   return result;
}

} }

max_stars_repo_name: DEIPworld/deip-chain

max_issues_repo_name: DEIPworld/deip-chain

max_forks_repo_name: DEIPworld/deip-chain


text:
#include <iostream>
#include <boost/log/trivial.hpp>
#include "server.h"
#include "database.h"
#include "mailer.h"

int main() {
    BOOST_LOG_TRIVIAL(info) << "Server app startup";
    try {
        auto srv = server();
        auto db = database();
        while (1)
        {
            clientInfo ci;
            auto valid = srv.Listen(ci); //TODO spawn each connection in its own thread
            if (valid)
                db.add(ci);
            //TODO replace hardcoded with XML file
            if ((ci.mem < 50) ||
                (ci.cpu > 80) ||
                (ci.processes > 50))
            {
                printf("\nSending mail\n");
                mailer::sendAlert(ci);
            }
        }
    }
    catch (std::exception& e) {
        BOOST_LOG_TRIVIAL(fatal) << "Caught exception:";
        BOOST_LOG_TRIVIAL(fatal) << e.what() << std::endl;
    }
    return 0;
}


text:
#include <iostream>
#include <vector>
#include <map>
#include <string>
#include <sstream>

#pragma warning( push )
# pragma warning (disable:4800)
#include <Eigen/Dense>
#include <Eigen/StdVector> //NodeArray, EdgeArray
#pragma warning( pop )

namespace cmg {
	typedef double Precision;
	typedef Eigen::Matrix<Precision, Eigen::Dynamic, Eigen::Dynamic> MatT;
	typedef Eigen::Matrix<Precision, Eigen::Dynamic, 1> VecT;
	typedef Eigen::Matrix<Precision, 2, 2> Mat2T;
	typedef Eigen::Matrix<Precision, 2, 1> Vec2T;
	typedef Eigen::Matrix<Precision, 3, 3> Mat3T;
	typedef Eigen::Matrix<Precision, 3, 1> Vec3T;
	typedef Eigen::Matrix<Precision, 4, 4> Mat4T;
	typedef Eigen::Matrix<Precision, 4, 1> Vec4T;
	typedef Eigen::Matrix<Precision, 6, 6> Mat6T;
	typedef Eigen::Matrix<Precision, 6, 1> Vec6T;
	typedef Eigen::Matrix<Precision, 2, 4> Mat2x4T;
	typedef Eigen::Matrix<Precision, 4, 2> Mat4x2T;
	typedef Eigen::Matrix<Precision, 3, 4> Mat3x4T;
	typedef Eigen::Matrix<Precision, 4, 3> Mat4x3T;

	inline static Precision eps() {
		return 10*std::numeric_limits<Precision>::epsilon();
	}

	typedef std::vector<int> Veci;
	inline static void range(const int start, const int end, const int step, Veci& ret)
	{
		ret.clear();
		ret.reserve((end-start+1)/step);
		for(int v=start; v<=end; v+=step) ret.push_back(v);
		ret.resize(ret.size());
	}

	inline static Precision perimeter(const Mat2x4T& u)
	{
		Precision ret=0;
		for(int i=0; i<4; ++i) {
			ret += (u.col(i) - u.col((i+1)%4)).norm();
		}
		return ret;
	}

	struct Calibration {
		Vec4T k;		//[fx, fy, cx, cy]
		Vec2T d;		//[k1, k2]
		Mat4T Ck;		//Cov[k]
		Mat2T Cd;		//Cov[d]
	public:
		EIGEN_MAKE_ALIGNED_OPERATOR_NEW

		inline Mat3T K() const
		{
			Mat3T ret;
			ret.setIdentity();
			ret(0,0)=k(0);
			ret(1,1)=k(1);
			ret(0,2)=k(2);
			ret(1,2)=k(3);
			return ret;
		}

		inline void print() const
		{
			std::cout<<"k="<<k.transpose()<<std::endl;
			std::cout<<"d="<<d.transpose()<<std::endl;
			//std::cout<<"Ck=\n"<<Ck<<std::endl; //TODO: Ck, Cd not used yet
			//std::cout<<"Cd=\n"<<Cd<<std::endl;
		}

		template<typename T>
		void project(const T X[3], T U[2]) const
		{
			T fx(k(0)), fy(k(1)), cx(k(2)), cy(k(3));
			T k1(d(0)), k2(d(1));

			T xn( X[0]/X[2] ), yn( X[1]/X[2] );
			T r2 = xn*xn + yn*yn;
			T factor=T(1)+(k2*r2+k1)*r2;
			T xnp=xn*factor;
			T ynp=yn*factor;
			U[0] = fx*xnp+cx;
			U[1] = fy*ynp+cy;
		}

		Vec2T project(const Vec3T& X) const
		{
			Vec2T ret;
			project(X.data(), ret.data());
			return ret;
		}
	};

	//6D pose
	struct Pose {
		Vec6T p;		//[ra*rx,ra*ry,ra*rz,tx,ty,tz]
		Mat6T Cp;		//Cov[p]
	public:
		EIGEN_MAKE_ALIGNED_OPERATOR_NEW

		static int nParams() { return 6; }

		inline Vec3T transform(const Vec3T& X) const
		{
			return (toT() * X.homogeneous()).hnormalized();
		}

		//return rotation matrix R
		inline Mat3T R() const
		{
			Precision ang = p.head<3>().norm();
			Vec3T axis(1,0,0);
			if(ang>eps())
				axis = p.head<3>().normalized();
			else
				ang = 0;
			Eigen::AngleAxis<Precision> aa(ang, axis);
			return aa.toRotationMatrix();
		}

		//return translation vector t
		inline Vec3T t() const
		{
			return p.tail<3>();
		}

		//return T=[R,t;0,1]
		inline Mat4T toT() const
		{
			Mat4T ret;
			ret.setIdentity();
			ret.topLeftCorner<3,3>() = R();
			ret.topRightCorner<3,1>() = t();
			return ret;
		}

		inline Mat4T invT() const
		{
			Mat3T rot = R().transpose();
			Mat4T T;
			T.setIdentity();
			T.topLeftCorner<3,3>() = rot;
			T.topRightCorner<3,1>() = -rot * t();
			return T;
		}

		inline Vec6T invp() const
		{
			return T2p(invT());
		}

		inline void fromR(const Mat3T& R)
		{
			Eigen::AngleAxis<Precision> aa;
			aa.fromRotationMatrix(R);
			p.head<3>() = aa.axis() * aa.angle();
		}

		inline void fromt(const Vec3T& t)
		{
			p.tail<3>() = t;
		}

		inline void fromT(const Mat4T& T)
		{
			fromR(T.topLeftCorner<3,3>());
			fromt(T.topRightCorner<3,1>());
		}

		inline static Vec6T T2p(const Mat4T& T)
		{
			return Pose(T).p;
		}

		inline static Mat4T p2T(const Vec6T& p)
		{
			return Pose(p).toT();
		}

		Pose() { p.setZero(); Cp.setZero(); }
		Pose(const Vec6T& p_) : p(p_) {
			Cp.setZero();
		}
		Pose(const Mat4T& T) {
			fromT(T);
			Cp.setZero();
		}
		Pose(const Mat4T& T, const Mat6T& Covp) : Cp(Covp)
		{
			fromT(T);
		}
	};
	typedef std::vector< Vec6T, Eigen::aligned_allocator<Vec6T> > Vec6TArray;

	struct Node : public Pose {
		EIGEN_MAKE_ALIGNED_OPERATOR_NEW

		std::string name;	//name of the marker
		Veci vmeids;		//vmEdges' ids linked to this marker (not used now)
	};
	typedef std::vector< Node, Eigen::aligned_allocator<Node> > NodeArray;
	typedef int NID;
	const NID INVALID_NID = -1;

	struct Edge {
		EIGEN_MAKE_ALIGNED_OPERATOR_NEW

		int vid;	//id of view
		int mid;	//id of marker
		Mat2x4T u;	//observations, 2x4, TODO: allow different observations
	};
	typedef std::vector< Edge, Eigen::aligned_allocator<Edge> > EdgeArray;
	typedef int EID;
	const EID INVALID_EID = -1;

	struct Observation {
		EIGEN_MAKE_ALIGNED_OPERATOR_NEW

		std::string name;		//marker's name
		Pose init_marker_pose;	//marker pose in view, i.e., Tmc
		Mat2x4T u;				//observed marker corners in view

		bool operator<(const Observation& other) const {
			return perimeter(u) > perimeter(other.u);
		}
	};
	typedef std::vector< Observation,  Eigen::aligned_allocator<Observation> > ObsArray;
	typedef std::vector<ObsArray> VecObsArray;

	typedef std::map<std::string, int> Str2Int;

	class CMGraph {
	public: //member variables
		NodeArray markers;
		NodeArray views;
		EdgeArray edges;

		Pose fixed_marker_pose;
		NID fixed_marker_id;

		Calibration calib;
		Str2Int name2mid; //name -> makrer id

		Precision marker_half_size;
		bool verbose;

		struct Callback {
			virtual void operator()(const CMGraph& G) = 0;
		};
		typedef Callback* CallbackPtr;
		CallbackPtr cb_addObsFromNewView;
		CallbackPtr cb_optimizeNewViewPose;
		CallbackPtr cb_optimizeOneViewPose;
		CallbackPtr cb_optimizePose;

	public: //member functions
		CMGraph() : fixed_marker_id(INVALID_NID), marker_half_size(1), verbose(true),
			cb_addObsFromNewView(0), cb_optimizeNewViewPose(0), cb_optimizeOneViewPose(0), cb_optimizePose(0)
		{}

		//print input information
		inline void print() const
		{
			std::cout<<"-------------------"<<std::endl;
			std::cout<<"calib:"<<std::endl;
			calib.print();
			std::cout<<std::endl;

			std::cout<<"marker_half_size="<<marker_half_size<<std::endl;
			std::cout<<"verbose="<<verbose<<std::endl;
			std::cout<<"-------------------"<<std::endl;
		}

		//print CMGraph state
		void report(std::ofstream& out) const;

		inline int nParams() const
		{
			return Pose::nParams()*static_cast<int>(markers.size()+views.size());
		}

		inline int nResiduals() const
		{
			return nObsResiduals() + nCstResiduals();
		}

		inline int nObsResiduals() const
		{
			return 8*static_cast<int>(edges.size()); //TODO: allow each marker to have more than 4 point observations
		}

		inline int nCstResiduals() const
		{
			return Pose::nParams();
		}

		inline bool empty() const
		{
			return nParams()<=0;
		}

		NID setFixedMarker(const std::string &fixed_marker_name,
			const Precision p_ang=1e-4,
			const Precision p_pos=1e-2);

		//add all observed markers in this view to the graph, and also add view and edges
		//oa will be sorted by descending order of perimeters of the observed markers in the image
		NID addObsFromNewView(
			ObsArray& oa,
			const std::string &view_name,
			const bool addNewMarker=true);

		bool optimizeNewViewPose(
			ObsArray& oa,
			Node& newView,
			const Precision sigma_u,
			const int max_iter,
			const Precision huber_loss_bandwidth=10, // +/- 10 pixels
			const Precision error_rel_tol=1e-2,
			const bool computeCovariance=false);

		bool optimizeOneViewPose(
			const NID vid,
			const Precision sigma_u,
			const int max_iter,
			const Precision huber_loss_bandwidth=10, // +/- 10 pixels
			const Precision error_rel_tol=1e-2,
			const bool computeCovariance=false);

		//bundle adjustment to optimize all markers' and views' poses
		bool optimizePose(
			const Precision sigma_u,
			const int max_iter,
			const Precision huber_loss_bandwidth=10, // +/- 10 pixels
			const Precision error_rel_tol=1e-2,
			const bool computeCovariance=false);

	public: //static functions
		//return 2D coordinates of a marker's 4 corners
		inline static Mat2x4T marker_x(const Precision half_size=1)
		{
			Precision ret_[]={
				-1,-1,
				 1,-1,
				 1, 1,
				-1, 1
			};

			Mat2x4T ret(ret_);
			return ret * half_size;
		}

		//return 3D coordinates of a marker's 4 corners
		inline static Mat3x4T marker_X(const Precision half_size=1)
		{
			Precision ret_[]={
				-1,-1, 0,
				 1,-1, 0,
				 1, 1, 0,
				-1, 1, 0
			};

			Mat3x4T ret(ret_);
			return ret * half_size;
		}

		static void BatchProcess(
			const VecObsArray &frames, //each frame's Observations' order would be sorted
			const Calibration &calib,
			CMGraph& G,
			const std::string &fixed_marker_name="",
			const Precision p_ang=1e-4,
			const Precision p_pos=1e-2,
			const Precision sigma_u=0.2,
			const int max_iter_per_opt=20,
			const bool do_covariance_estimation=true);

	protected:
		NID newMarker(const Vec6T& p, const Mat6T& Cp, const std::string& name);

		NID newView(const Vec6T& p, const Mat6T& Cp, const std::string& name);

		EID newEdge(const NID vid, const NID mid, const MatT& u);
	};
}

max_stars_repo_name: simbaforrest/masfm

max_issues_repo_name: simbaforrest/masfm

max_forks_repo_name: simbaforrest/masfm


text:
// SPDX-FileCopyrightText: 2015 - 2021 Marcin Łoś <marcin.los.91@gmail.com>
// SPDX-License-Identifier: MIT

#ifndef TUMOR_TUMOR_HPP
#define TUMOR_TUMOR_HPP

#include <cmath>

#include <boost/format.hpp>

#include "ads/executor/galois.hpp"
#include "ads/output_manager.hpp"
#include "ads/simulation.hpp"
#include "params.hpp"
#include "skin.hpp"
#include "state.hpp"
#include "vasculature.hpp"

namespace bsp = ads::bspline;

namespace tumor {

class tumor_2d : public ads::simulation_2d {
private:
    static constexpr std::size_t Dim = 2;
    using Base = ads::simulation_2d;

    state<Dim> now, prev;
    params p;

    vasc::vasculature vasculature;

    ads::output_manager<2> output;

    int save_every;

    int vasc_update_every = 10;

    ads::galois_executor executor{4};

    ads::bspline::eval_ctx xctx;
    ads::bspline::eval_ctx yctx;
    ads::bspline::eval_ders_ctx xdctx;
    ads::bspline::eval_ders_ctx ydctx;

public:
    tumor_2d(const ads::config_2d& config, const params& params, int save_every,
             vasc::vasculature&& vasculature);

private:
    struct constant {
        double value;

        explicit constant(double v)
        : value{v} { }

        double operator()(double, double) const { return value; }
    };

    double init_tumor(double x, double y) {
        double dx = (x - 1500) / 600;
        double dy = (y - 1500) / 600;
        double r2 = std::min(dx * dx + dy * dy, 1.0);
        return 0.8 * (r2 - 1) * (r2 - 1) * (r2 + 1) * (r2 + 1);
    }

    double init_M(double x, double y) const {
        auto lay = p.skin.layer_at(x, y, x);
        return lay == skin_model::layer::dermis ? 1 : p.init_M;
    }

    void dirichlet() {
        x.fix_left();
        x.fix_right();
        y.fix_left();
        y.fix_right();
    }

    void apply_boundary_conditions(vector_type& v) {
        for (int i = 0; i < y.dofs(); ++i) {
            v(0, i) = 0;
            v(x.dofs() - 1, i) = 0;
        }
        for (int i = 0; i < x.dofs(); ++i) {
            v(i, 0) = 0;
            v(i, y.dofs() - 1) = 0;
        }
    }

    point_type normalize(point_type p) const {
        return {(p[0] - x.a) / (x.b - x.a), (p[1] - y.a) / (y.b - y.a)};
    }

    void prepare_matrices() {
        dirichlet();
        Base::prepare_matrices();
    }

    void plot_vasculature(int iter) {
        using boost::format;
        vasculature.plot_veins(str(format("vasculature_%d.data") % iter));
        vasculature.plot_oxygen(str(format("oxygen_%d.data") % iter));
    }

    void before() override;

    void before_step(int /*iter*/, double /*t*/) override {
        using std::swap;
        swap(now, prev);

        now.clear();
    }

    void step(int iter, double /*t*/) override {
        compute_rhs();
        solve_all();
        update_vasculature(iter);
    }

    void save_to_file(int iter);

    void after_step(int iter, double /*t*/) override {
        int next_iter = iter + 1;
        if (next_iter % save_every == 0) {
            save_to_file(next_iter);
        }
    }

    void solve_all();

    value_type ensure_positive(value_type v) const {
        if (v.val < 0) {
            v.val = 0;
        }
        return v;
    }

    void compute_rhs();

    void update_vasculature(int iter);

    double& val(vector_type& v, index_type idx) { return v(idx[0], idx[1]); }
};

}  // namespace tumor

#endif  // TUMOR_TUMOR_HPP

max_stars_repo_name: Pan-Maciek/iga-ads

max_issues_repo_name: Pan-Maciek/iga-ads

max_forks_repo_name: Pan-Maciek/iga-ads


text:
/******************************************
Copyright (c) 2020, Jeffrey Dudek
******************************************/

#include "util/graded_clauses.h"

#include <algorithm>
#include <iostream>
#include <unordered_map>
#include <utility>
#include <vector>

#include <boost/pending/disjoint_sets.hpp>

namespace util {
GradedClauses::GradedClauses(
  const std::vector<std::vector<size_t>> &clauses
) {
  for (size_t i = 0; i < clauses.size(); i++) {
    components_.push_back(GradedClauses(clauses[i], i));
  }
}

void GradedClauses::write_line_graph(
  std::ostream *output,
  size_t num_variables
) const {
  size_t num_edges = count_line_graph_edges();
  *output << "p tw " << num_variables << " " << num_edges << "\n";
  write_line_graph_edges(output);
}

size_t GradedClauses::count_line_graph_edges() const {
  size_t total = (variables_.size() * (variables_.size()-1)) / 2;
  for (const GradedClauses &clause : components_) {
    total += clause.count_line_graph_edges();
  }
  return total;
}

void GradedClauses::write_line_graph_edges(std::ostream *output) const {
  for (size_t i = 0; i < variables_.size(); i++) {
    for (size_t j = i+1; j < variables_.size(); j++) {
      *output << variables_[i] << " ";
      *output << variables_[j] << "\n";
    }
  }
  for (const GradedClauses &clause : components_) {
    clause.write_line_graph_edges(output);
  }
}

void GradedClauses::group_by(
  const std::vector<size_t> &kept_variables,
  size_t max_var_id
) {
  /*
  Store variables to project in easily accessible way
  */
  std::vector<bool> is_projected_var(max_var_id+1, true);
  for (size_t var : kept_variables) {
    is_projected_var[var] = false;
  }

  /*
  Prepare a disjoint set data structure, whose entries are projected variables
  */
  std::vector<int> rank(is_projected_var.size());
  std::vector<int> parent(is_projected_var.size());
  boost::disjoint_sets<int *, int*> classes(&rank[0], &parent[0]);
  for (size_t var_id = 0; var_id < is_projected_var.size(); var_id++) {
    if (is_projected_var[var_id]) {
      classes.make_set(var_id);
    }
  }

  /*
  Build equivalence classes of variables for the (transitive closure of):
    x ~ y if x=y or (x and y are projected and some clause has both x and y)
  */
  for (const GradedClauses &clause : components_) {
    size_t last = SIZE_MAX;
    for (size_t var_id : clause.variables_) {
      if (is_projected_var[var_id]) {
        if (last != SIZE_MAX) {
          // var_id and last appear together in a clause, so equate them
          classes.union_set(var_id, last);
        }
        last = var_id;
      }
    }
  }

  /*
  Group clauses together if they have equivalent, projected variables
  */
  std::vector<GradedClauses> result;
  std::unordered_map<size_t, size_t> grouped_components;
  for (int i = 0; i < components_.size(); i++) {
    // Choose a GradedClauses for component i
    size_t which_group = SIZE_MAX;
    for (size_t var_id : components_[i].variables_) {
      // Clauses with a projected variable joins a group cooresponding to
      // their equivalence class
      if (is_projected_var[var_id]) {
        size_t rep = classes.find_set(var_id);
        auto group = grouped_components.find(rep);
        if (group != grouped_components.end()) {
          which_group = group->second;
        } else {
          which_group = result.size();
          result.push_back(GradedClauses());
          grouped_components.insert({rep, which_group});
        }
        break;
      }
    }

    if (which_group == SIZE_MAX) {
      // Clauses with no projected variables are placed alone
      result.push_back(std::move(components_[i]));
    } else {
      result[which_group].components_.push_back(std::move(components_[i]));
    }
  }

  /*
  Compute the sets of free variables for the new groups
  */
  for (GradedClauses &new_group : result) {
    // Clauses with no projected variables already have free variables set
    if (new_group.components_.size() == 0) {
      continue;
    }

    for (const GradedClauses &clause : new_group.components_) {
      for (size_t var_id : clause.variables_) {
        if (is_projected_var[var_id]) {
          continue;
        } else {
          new_group.variables_.push_back(var_id);
        }
      }
    }

    std::sort(new_group.variables_.begin(), new_group.variables_.end());
    new_group.variables_.erase(unique(new_group.variables_.begin(),
                                      new_group.variables_.end()),
                                      new_group.variables_.end());
  }

  components_ = std::move(result);
}
}  // namespace util

ext: cc


text:
#ifndef CPP_MICROBENCH_ALIGNED_VECTOR_HPP
#define CPP_MICROBENCH_ALIGNED_VECTOR_HPP

#include <vector>

#include <boost/align/aligned_allocator.hpp>

constexpr std::size_t alignment = 4096;

template <typename T>
using aligned_vector =
std::vector<T,
    boost::alignment::aligned_allocator<T, alignment>>;

#endif//CPP_MICROBENCH_ALIGNED_VECTOR_HPP


text:
#include <boost/algorithm/string.hpp>
#include "URIConfig.h"

namespace paradigm4 {
namespace pico {
namespace core {

bool URIConfig::set_uri(const std::string& uri) {
    clear();
    std::string tmpuri = uri;
    boost::algorithm::trim_if(tmpuri, boost::is_any_of(" \t\r\n"));
    if (tmpuri.length() == 0)
        return true;
    //params
    auto pos = uri.find("?");
    if (pos != uri.npos) {
        tmpuri = URICode::decode(uri.substr(0, pos));
        std::string params = uri.substr(pos + 1);
        const char* p = params.c_str();
        const char* key = p;
        int key_len = 0;
        const char* val = nullptr;
        int val_len = 0;
        while(true) {
            if (*p == '=') {
                key_len = p - key;
                RCHECK(key_len > 0) << uri << " parse key val fail, key_len == 0.";
                val = ++p;
            } else if (*p == '&' || *p =='\0') {
                RCHECK(key_len > 0) << uri << " parse key val fail, no key.";
                val_len = p - val;

                std::string skey = URICode::decode(std::string(key, key_len));
                std::string sval = URICode::decode(std::string(val, val_len));
                SCHECK(_config.set_val(skey, sval, URILVL::URICFG)) << uri << " param "
                                                                    << skey << " set failed.";

                if (*p == '\0')
                    break;
                key_len = 0;
                val_len = 0;
                key = ++p;
            } else {
                ++p;
            }
        }
    } else {
        tmpuri = URICode::decode(uri);
    }

    //storage
    _store = ShellUtility::fs_type(tmpuri, _prefix, _name);
    if (_store == FileSystemType::UNKNOWN) {
        RLOG(FATAL) << uri << " unknown storage";
    }

    return true;
}

void URIConfig::replace_param(const std::map<std::string,
                                             std::pair<std::string, int>>& params) {
    std::string src = uri();
    for (auto it = params.begin(); it != params.end(); ++it) {
        _config.set_val(it->first, it->second.first, it->second.second);
    }
    SLOG(INFO) << "replace uri [" << src << "] -> [" << uri() << "]";
}

std::string URIConfig::value() const {
    std::string str_param = _config.to_string(URILVL::DEFAULT + 1);
    return _prefix + _name + str_param;
}

std::string uri_config_t::to_string(int low_lvl) const {
    std::ostringstream oss;
    int cnt = 0;
    for (auto it = params.begin(); it != params.end(); ++it) {
        if (it->second.second < low_lvl) {
            continue;
        }
        oss << (cnt++ == 0 ? "?" : "&")
            << it->first << "=" << URICode::encode(it->second.first);
    }
    return oss.str();
}

}
} // namespace pico
} // namespace paradigm4