text:
/*

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

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

This file is part of Chaste.

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

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

*/

#ifndef TESTBIDOMAINPROBLEM_HPP_
#define TESTBIDOMAINPROBLEM_HPP_

#include <cxxtest/TestSuite.h>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <vector>

#include "LuoRudy1991.hpp"
#include "FaberRudy2000.hpp"
#include "BidomainProblem.hpp"
#include "MonodomainProblem.hpp"
#include "Hdf5DataReader.hpp"
#include "PlaneStimulusCellFactory.hpp"
#include "HeartEventHandler.hpp"
#include "PetscTools.hpp"
#include "PetscSetupAndFinalize.hpp"
#include "ArchiveOpener.hpp"
#include "CmguiMeshWriter.hpp"
#include "FileFinder.hpp"
#include "OutputFileHandler.hpp"
#include "MemfemMeshReader.hpp"
#include "BidomainSolver.hpp"
#include "CompareHdf5ResultsFiles.hpp"
#include "NumericFileComparison.hpp"
#include "Electrodes.hpp"
#include "SimpleBathProblemSetup.hpp"
#include "FileComparison.hpp"
#include "SingleTraceOutputModifier.hpp"

#ifdef CHASTE_VTK
#define _BACKWARD_BACKWARD_WARNING_H 1 //Cut out the strstream deprecated warning for now (gcc4.3)
#include <vtkVersion.h>
#endif

/**
 * This stimulus causes an exception to be thrown when normally a simple stimulus would activate.
 *
 * Used for testing exception handling in TestBidomain3dErrorHandling
 */
class ExceptionStimulus : public AbstractStimulusFunction
{
private:
    /** Duration of initial stimulus, typically in milliseconds */
    double mDuration;
    /** The time at which the stimulus starts, typically in milliseconds */
    double mTimeOfStimulus;

public:

    /**
     * Constructor.
     *
     * @param duration  Duration of initial stimulus milliseconds
     * @param timeOfStimulus  The time at which the stimulus starts (defaults to 0.0) milliseconds
     */
    ExceptionStimulus(double duration, double timeOfStimulus)
    {
        mDuration = duration;
        mTimeOfStimulus = timeOfStimulus;

        mDuration += (mDuration+mTimeOfStimulus)*DBL_EPSILON;
    }

    /**
     * @param time  time at which to return the stimulus
     * @return zero or throws an Exception at the time specified by constructor.
     */
    double GetStimulus(double time)
    {
        if (mTimeOfStimulus < time && time <= mDuration+mTimeOfStimulus)
        {
            EXCEPTION("Stimulus has caused a bad thing to happen in a cell model.");
        }
        return 0.0;
    }
};

/**
 * ExceptionStimulusCellFactory provides all cells with an ExceptionStimulus.
 *
 * Used for testing exception handling in TestBidomain3dErrorHandling
 */
template<class CELL, unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
class ExceptionStimulusCellFactory : public AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>
{
protected:
    /** The stimulus to apply at stimulated nodes */
    boost::shared_ptr<ExceptionStimulus> mpStimulus;

public:
    /**
     * Constructor
     * @param stimulusDuration  The duration of time the exception should be applied (defaults to 1ms).
     * @param startTime When the exception should be introduced (defaults to 0.0).
     */
    ExceptionStimulusCellFactory(double stimulusDuration=1, double startTime = 0.0)
        : AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>()
    {
        mpStimulus.reset(new ExceptionStimulus(stimulusDuration, startTime));
    }

    /**
     * @param node  The global index of a node
     * @return  A cardiac cell which corresponds to this node.
     */
    AbstractCardiacCellInterface* CreateCardiacCellForTissueNode(Node<SPACE_DIM>* pNode)
    {
        return new CELL(this->mpSolver, mpStimulus);
    }
};

class DelayedTotalStimCellFactory : public AbstractCardiacCellFactory<1>
{
private:
    // define a new stimulus
    boost::shared_ptr<SimpleStimulus> mpIntraStimulus;

public:
    DelayedTotalStimCellFactory(double mag)
        : AbstractCardiacCellFactory<1>(),
          mpIntraStimulus(new SimpleStimulus(  mag, 0.1, 0.1))
    {
    }

    AbstractCardiacCell* CreateCardiacCellForTissueNode(Node<1>* pNode)
    {
        return new CellLuoRudy1991FromCellML(mpSolver, mpIntraStimulus);
    }
};

class TestBidomainProblem : public CxxTest::TestSuite
{
private:
    std::vector<double> mSolutionReplicated1d2ms;///<Used to test differences between tests

public:
    void tearDown()
    {
        HeartConfig::Reset();
    }

    void TestBidomainErrorHandling() throw (Exception)
    {
        if (PetscTools::GetNumProcs() > 2u)
        {
            // There are only 2 nodes in this simulation
            TS_TRACE("This test is not suitable for more than 2 processes.");
            return;
        }
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.75, 1.75, 1.75));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(7.0, 7.0, 7.0));
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1_1_element");
        HeartConfig::Instance()->SetOutputDirectory("BidomainWithErrors");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain");

        // Check the linear system can be solved to a low tolerance (in particular, checks the null space
        // stuff was implemented correctly
        HeartConfig::Instance()->SetUseAbsoluteTolerance(1e-14);

        // We put in a massive stimulus, and alter the printing timestep from default, to cover
        // some memory handling in the case of exceptions being thrown (#1222).
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.1);

        // 100x bigger stimulus than normal.
        ExceptionStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory(1, 0.02);
        BidomainProblem<1> bidomain_problem( &bidomain_cell_factory );
        bidomain_problem.Initialise();

        TS_ASSERT_THROWS_THIS(bidomain_problem.Solve(),
           "Stimulus has caused a bad thing to happen in a cell model.");
    }

    //This test reproduces the problem size of heart/test/data/xml/bidomain1d_small.xml
    void TestBidomainSmallestMesh()
    {
        if (PetscTools::GetNumProcs() > 3u)
        {
            TS_TRACE("This test is not suitable for more than 3 processes.");
            return;
        }
       HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetFibreLength(0.1, 0.05);
        HeartConfig::Instance()->SetSpaceDimension(1u);
        HeartConfig::Instance()->SetOutputDirectory("Bidomain1dSmall");
        HeartConfig::Instance()->SetOutputFilenamePrefix("small");
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.1, 0.1);

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory;
        BidomainProblem<1> bidomain_problem( &bidomain_cell_factory );

        bidomain_problem.Initialise();
        bidomain_problem.Solve();
    }

    // NOTE: This test uses NON-PHYSIOLOGICAL parameters values (conductivities,
    // surface-area-to-volume ratio, capacitance, stimulus amplitude). Essentially,
    // the equations have been divided through by the surface-area-to-volume ratio.
    // (Historical reasons...)
    void TestBidomainDg01DPinned()
    {
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1_100_elements");
        HeartConfig::Instance()->SetOutputDirectory("bidomainDg01d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory;
        BidomainProblem<1> bidomain_problem( &bidomain_cell_factory );

        bidomain_problem.Initialise();

        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);

        std::vector<unsigned> pinned_nodes;

        // check throws if the fixed node num isn't valid
        pinned_nodes.push_back(1000);
        bidomain_problem.SetFixedExtracellularPotentialNodes(pinned_nodes);
        TS_ASSERT_THROWS_THIS( bidomain_problem.Solve(), "Fixed node number must be less than total number nodes" );

        /* HOW_TO_TAG Cardiac/Problem definition
         * Fix phi_e at particular nodes (note: this is not required)
         */
        // Pin extracellular potential of node 100 to 0
        pinned_nodes.clear();
        pinned_nodes.push_back(100);
        bidomain_problem.SetFixedExtracellularPotentialNodes(pinned_nodes);

        try
        {
            bidomain_problem.Solve();
        }
        catch (Exception e)
        {
            TS_FAIL(e.GetMessage());
        }

        DistributedVector striped_voltage = bidomain_problem.GetSolutionDistributedVector();
        DistributedVector::Stripe voltage(striped_voltage,0);

        for (DistributedVector::Iterator index = striped_voltage.Begin();
             index != striped_voltage.End();
             ++index)
        {
            // assuming LR model has Ena = 54.4 and Ek = -77
            double Ena   =  54.4;   // mV
            double Ek    = -77.0;   // mV

            TS_ASSERT_LESS_THAN_EQUALS( voltage[index], Ena +  30);
            TS_ASSERT_LESS_THAN_EQUALS(-voltage[index] + (Ek-30), 0);

            std::vector<double> r_ode_vars = bidomain_problem.GetBidomainTissue()->GetCardiacCell(index.Global)->GetStdVecStateVariables();
            for (int j=0; j<8; j++)
            {
                // if not voltage or calcium ion conc, test whether between 0 and 1
                if ((j!=0) && (j!=7))
                {
                    TS_ASSERT_LESS_THAN_EQUALS( r_ode_vars[j], 1.0);
                    TS_ASSERT_LESS_THAN_EQUALS(-r_ode_vars[j], 0.0);
                }
            }

            // wave shouldn't have reached the second half of the mesh so
            // these should all be near the resting potential

            if (index.Global>50)
            {
                TS_ASSERT_DELTA(voltage[index], -83.85, 0.1);
            }

            // final voltages for nodes 0 to 5 produced with ksp_rtol=1e-9
            double test_values[6]={31.0335, 28.9214, 20.0279, -3.92649, -57.9395, -79.7754};

            for (unsigned node=0; node<=5; node++)
            {
                if (index.Global == node)
                {
                    // test against hardcoded value to check nothing has changed
                    TS_ASSERT_DELTA(voltage[index], test_values[node], 7e-3);
                    //With ksp_rtol set to 1e-6 the starting value may lead to changes of more that 1e-3 in final answer
                }
            }
        }
        DistributedVector::Stripe extracellular_potential(striped_voltage,1);
        if (striped_voltage.IsGlobalIndexLocal(100))
        {
            TS_ASSERT_DELTA(extracellular_potential[100], 0.0, 1e-6);
        }
    }

    // NOTE: This test uses NON-PHYSIOLOGICAL parameters values (conductivities,
    // surface-area-to-volume ratio, capacitance, stimulus amplitude). Essentially,
    // the equations have been divided through by the surface-area-to-volume ratio.
    // (Historical reasons...)
    void TestBidomainDg01DAveragePhiEOverDifferentRows()
    {
        HeartEventHandler::Disable();

        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1_100_elements");
        HeartConfig::Instance()->SetOutputDirectory("bidomainDg01d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        // Final values to test against have been produced with ksp_rtol=1e-9
        HeartConfig::Instance()->SetUseRelativeTolerance(1e-8);

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory;
        BidomainProblem<1> bidomain_problem( &bidomain_cell_factory );

        bidomain_problem.Initialise();

        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);

        // Apply the constraint 'Average phi = 0' to nodes 0, 50, 100.
        for (unsigned node=0; node<bidomain_problem.rGetMesh().GetNumNodes(); node+=50)
        {
            bidomain_problem.SetNodeForAverageOfPhiZeroed(node);

            try
            {
                bidomain_problem.Solve();
            }
            catch (Exception e)
            {
                TS_FAIL(e.GetMessage());
            }

            DistributedVector striped_voltage = bidomain_problem.GetSolutionDistributedVector();
            DistributedVector::Stripe voltage(striped_voltage,0);
            DistributedVector::Stripe phi_e(striped_voltage,1);

            for (DistributedVector::Iterator index = striped_voltage.Begin();
                 index != striped_voltage.End();
                 ++index)
            {
                // assuming LR model has Ena = 54.4 and Ek = -77
                double Ena   =  54.4;   // mV
                double Ek    = -77.0;   // mV

                TS_ASSERT_LESS_THAN_EQUALS( voltage[index], Ena +  30);
                TS_ASSERT_LESS_THAN_EQUALS(-voltage[index] + (Ek-30), 0);

                std::vector<double> r_ode_vars = bidomain_problem.GetBidomainTissue()->GetCardiacCell(index.Global)->GetStdVecStateVariables();
                for (int j=0; j<8; j++)
                {
                    // if not voltage or calcium ion conc, test whether between 0 and 1
                    if ((j!=0) && (j!=7))
                    {
                        TS_ASSERT_LESS_THAN_EQUALS( r_ode_vars[j], 1.0);
                        TS_ASSERT_LESS_THAN_EQUALS(-r_ode_vars[j], 0.0);
                    }
                }

                // wave shouldn't have reached the second half of the mesh so
                // these should all be near the resting potential

                if (index.Global>50)
                {
                    TS_ASSERT_DELTA(voltage[index], -83.85, 0.1);
                }

                // final voltages for nodes 0 to 5 produced with ksp_rtol=1e-9
                double voltage_test_values[6]={31.0335, 28.9214, 20.0279, -3.92649, -57.9395, -79.7754};

                if (index.Global<6)
                {
                    // test against hardcoded value to check nothing has changed
                    TS_ASSERT_DELTA(voltage[index], voltage_test_values[index.Global], 7e-3);
                }

                // final extracellular potencials for nodes 0 to 5 produced with ksp_rtol=1e-9
                double phi_e_test_values[6]={-55.2567, -54.2006, -49.7538, -37.7767, -10.7701, 0.148278};

                if (index.Global<6)
                {
                    // test against hardcoded value to check nothing has changed
                    TS_ASSERT_DELTA(phi_e[index], phi_e_test_values[index.Global], 7e-3);
                }
            }

            // check mean of extracellular potential is 0
            double local_phi_e=0.0;
            double total_phi_e=0.0;

            for (DistributedVector::Iterator index = striped_voltage.Begin();
                 index != striped_voltage.End();
                 ++index)
            {
                local_phi_e += phi_e[index];
            }

            int ierr = MPI_Allreduce(&local_phi_e, &total_phi_e, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
            TS_ASSERT_EQUALS(ierr, (int)MPI_SUCCESS);

            TS_ASSERT_DELTA(total_phi_e, 0, 1e-4);

            bidomain_problem.Initialise();
        }

        // Coverage of the exception in the solver itself
        BoundaryConditionsContainer<1,1,2> container;

        BidomainSolver<1,1> bidomain_solver(false,
                                            &bidomain_problem.rGetMesh(),
                                            bidomain_problem.GetBidomainTissue(),
                                            &container);

        TS_ASSERT_THROWS_THIS(bidomain_solver.SetRowForAverageOfPhiZeroed(0),
                "Row for applying the constraint \'Average of phi_e = zero\' should be odd in C++ like indexing");

        HeartEventHandler::Enable();
    }

    /*
     * The monodomain equations are obtained by taking the limit of the bidomain
     * equations as sigma_e tends to infinity (corresponding to the extracellular
     * space being grounded). Therefore, if we set sigma_e very large (relative to
     * sigma_i) in a bidomain simulation it should agree with a monodomain
     * simulation with the same parameters.
     *
     * NOTE: This test uses NON-PHYSIOLOGICAL parameters values (conductivities,
     * surface-area-to-volume ratio, capacitance, stimulus amplitude). Essentially,
     * the equations have been divided through by the surface-area-to-volume ratio.
     * (Historical reasons...)
     *
     */
    void TestCompareBidomainProblemWithMonodomain() throw(Exception)
    {
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1_100_elements");
        HeartConfig::Instance()->SetOutputDirectory("Monodomain1dVersusBidomain");
        HeartConfig::Instance()->SetOutputFilenamePrefix("monodomain1d");

        Vec monodomain_results;

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;

        // To avoid an issue with the Event handler only one simulation should be
        // in existance at a time: therefore monodomain simulation is defined in a block
        {
            ///////////////////////////////////////////////////////////////////
            // monodomain
            ///////////////////////////////////////////////////////////////////
            MonodomainProblem<1> monodomain_problem( &cell_factory );

            monodomain_problem.Initialise();

            HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
            HeartConfig::Instance()->SetCapacitance(1.0);

            // now solve
            monodomain_problem.Solve();

            VecDuplicate(monodomain_problem.GetSolution(), &monodomain_results);
            VecCopy(monodomain_problem.GetSolution(), monodomain_results);
        }

        ///////////////////////////////////////////////////////////////////
        // bidomain
        ///////////////////////////////////////////////////////////////////

        // keep the intra conductivity to be the same as monodomain
        // and the extra conductivity to be very large in comparison
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(1));
        HeartConfig::Instance()->SetOutputDirectory("Bidomain1d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain1d");

        BidomainProblem<1> bidomain_problem( &cell_factory );

        // Not really needed for coverage, but it's worth checking that logic works properly in both scenarios.
        TS_ASSERT(!bidomain_problem.GetHasBath());

        bidomain_problem.Initialise();

        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);

        // now solve
        bidomain_problem.Solve();

        ///////////////////////////////////////////////////////////////////
        // compare
        ///////////////////////////////////////////////////////////////////
        DistributedVector dist_bidomain_voltage = bidomain_problem.GetSolutionDistributedVector();
        DistributedVector monodomain_voltage = dist_bidomain_voltage.GetFactory()->CreateDistributedVector(monodomain_results);
        DistributedVector::Stripe bidomain_voltage(dist_bidomain_voltage, 0);
        DistributedVector::Stripe extracellular_potential(dist_bidomain_voltage, 1);

        for (DistributedVector::Iterator index = monodomain_voltage.Begin();
             index != monodomain_voltage.End();
             ++index)
        {
            if (index.Global==0)
            {
                TS_ASSERT_LESS_THAN(0, monodomain_voltage[index]);
            }
            // the mono and bidomains should agree closely
            TS_ASSERT_DELTA(monodomain_voltage[index], bidomain_voltage[index], 0.4);

            // the extracellular potential should be uniform
            TS_ASSERT_DELTA(extracellular_potential[index], 0, 0.06);
        }

        PetscTools::Destroy(monodomain_results);
    }

    ///////////////////////////////////////////////////////////////////
    // Solve a simple simulation and check the output was only
    // printed out at the correct times
    ///////////////////////////////////////////////////////////////////
    void TestBidomainProblemPrintsOnlyAtRequestedTimesAndOnlyRequestedNodes() throw (Exception)
    {
        HeartEventHandler::Disable();

        HeartConfig::Instance()->SetPrintingTimeStep(0.1);
        HeartConfig::Instance()->SetPdeTimeStep(0.01);
        HeartConfig::Instance()->SetSimulationDuration(0.3);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1mm_10_elements");
        HeartConfig::Instance()->SetOutputDirectory("Bidomain1d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_testPrintTimes");

        // run testing PrintingTimeSteps
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
        BidomainProblem<1>* p_bidomain_problem = new BidomainProblem<1>( &cell_factory );

        /* HOW_TO_TAG Cardiac/Output
         * Only output data for particular nodes
         */
        //Restrict the number of nodes
        std::vector<unsigned> nodes_to_be_output;
        nodes_to_be_output.push_back(0);
        nodes_to_be_output.push_back(5);
        nodes_to_be_output.push_back(10);
        p_bidomain_problem->SetOutputNodes(nodes_to_be_output);

        // Duplicate this with a single trace at node 5
        boost::shared_ptr<SingleTraceOutputModifier> trace_5(new SingleTraceOutputModifier("trace_5.txt", 5, 0.1));
        p_bidomain_problem->AddOutputModifier(trace_5);

        // for coverage:
        p_bidomain_problem->SetWriteInfo();

        p_bidomain_problem->Initialise();
        p_bidomain_problem->Solve();

        // Extra trace data
        OutputFileHandler handler("Bidomain1d", false);
        NumericFileComparison comp(handler.GetOutputDirectoryFullPath()+ "trace_5.txt", "heart/test/data/Bidomain1d/trace_5.txt");
        TS_ASSERT(comp.CompareFiles(5e-4));

        // read data entries for the time file and check correct
        Hdf5DataReader data_reader1=p_bidomain_problem->GetDataReader();
        std::vector<double> times = data_reader1.GetUnlimitedDimensionValues();

        TS_ASSERT_EQUALS( times.size(), (unsigned) 4);
        TS_ASSERT_DELTA( times[0], 0.00, 1e-12);
        TS_ASSERT_DELTA( times[1], 0.10, 1e-12);
        TS_ASSERT_DELTA( times[2], 0.20, 1e-12);
        TS_ASSERT_DELTA( times[3], 0.30, 1e-12);

        //Get back node over all times
        std::vector<double> node_0 = data_reader1.GetVariableOverTime("V", 0);
        TS_ASSERT_EQUALS( node_0.size(), 4U);
        TS_ASSERT_DELTA( node_0[0], -83.853, 1e-10);
        TS_ASSERT_DELTA( node_0[1], -83.8354, 5e-4);
        TS_ASSERT_DELTA( node_0[2], -83.8266, 3e-4);
        TS_ASSERT_DELTA( node_0[3], -83.8200, 3e-4);
        std::vector<double> node_5 = data_reader1.GetVariableOverTime("V", 5);
        TS_ASSERT_EQUALS( node_5.size(), 4U);
        std::vector<double> node_10 = data_reader1.GetVariableOverTime("V", 10);
        TS_ASSERT_EQUALS( node_10.size(), 4U);

        //Can't read back this node as it wasn't written
        TS_ASSERT_THROWS_THIS( data_reader1.GetVariableOverTime("V", 1),
                "The incomplete dataset 'Data' does not contain info of node 1");

        delete p_bidomain_problem;

        p_bidomain_problem = new BidomainProblem<1>( &cell_factory );

        // Now check that we can turn off output printing
        // Output should be the same as above: printing every 10th time step
        // because even though we set to print every time step...
        HeartConfig::Instance()->SetPrintingTimeStep(1);
        // ...we have output turned off
        p_bidomain_problem->PrintOutput(false);

        p_bidomain_problem->Initialise();
        p_bidomain_problem->Solve();

        Hdf5DataReader data_reader3=p_bidomain_problem->GetDataReader();
        times = data_reader3.GetUnlimitedDimensionValues();

        TS_ASSERT_EQUALS( times.size(), (unsigned) 4);
        TS_ASSERT_DELTA( times[0], 0.00,  1e-12);
        TS_ASSERT_DELTA( times[1], 0.10,  1e-12);
        TS_ASSERT_DELTA( times[2], 0.20,  1e-12);
        TS_ASSERT_DELTA( times[3], 0.30,  1e-12);

        delete p_bidomain_problem;
        HeartEventHandler::Enable();
    }

    void TestBidomainFallsOverProducesOutput() throw(Exception)
    {
#ifndef NDEBUG //Note that this test relies on the debug VerifyStateVariables() method throwing
        HeartConfig::Instance()->SetSimulationDuration(0.3);  //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1_100_elements");
        HeartConfig::Instance()->SetOutputDirectory("BidomainFallsOver");
        HeartConfig::Instance()->SetOutputFilenamePrefix("res");
        HeartConfig::Instance()->SetVisualizeWithMeshalyzer();

        // Something happens at 0.1ms

        // DelayedTotalStimCellFactory bidomain_cell_factory(-6e5); //Normal stimulus
        DelayedTotalStimCellFactory bidomain_cell_factory(-6e6); //Takes sodium out of range
        BidomainProblem<1> bidomain_problem( &bidomain_cell_factory );

        bidomain_problem.Initialise();

        TS_ASSERT_THROWS_CONTAINS(bidomain_problem.Solve(),
               "State variable fast_sodium_current_m_gate__m has gone out of range. Check numerical parameters, for example time and space stepsizes");

        // Make sure that there's time for the files to be written
        // (most files are only written by the master)
        PetscTools::Barrier();

        //Test for regular output
        Hdf5DataReader data_reader=bidomain_problem.GetDataReader();
        std::vector<double> times = data_reader.GetUnlimitedDimensionValues();
        //TS_ASSERT_EQUALS( times.size(), 31U);// For normal stimulation
        TS_ASSERT_EQUALS( times.size(), 31U);  // For fully allocated hdf5
        //TS_ASSERT_EQUALS( times.size(), 21U);// For over stimulation

        TS_ASSERT_DELTA( times[1], 0.01, 1e-12);
        //TS_ASSERT_DELTA( times.back(), 0.30, 1e-12);//For normal stimulation
        //TS_ASSERT_DELTA( times.back(), 0.20, 1e-12);//For over stimulation
        TS_ASSERT_DELTA( times.back(), 0, 1e-12);//For hdf5 writer allocated up to 30 slots
        TS_ASSERT_DELTA( times[20], 0.20, 1e-12);//For hdf5 writer allocated up to 30 slots
        TS_ASSERT_DELTA( times[21],    0, 1e-12);//For hdf5 writer allocated up to 30 slots

        // Test for post-processed output (and don't wipe the directory!)
        OutputFileHandler handler("BidomainFallsOver/output", false);

        std::string files[6] = {"res_mesh.pts","res_mesh.cnnx","ChasteParameters.xml",
                                "res_Phi_e.dat","res_V.dat","res_times.info"};

        for(unsigned i=0; i<6; i++)
        {
            TS_ASSERT(handler.FindFile(files[i]).Exists());
        }
#endif //NDEBUG Note that this test relies on the debug VerifyStateVariables() method throwing
    }

    void TestBidomainProblemExceptions() throw (Exception)
    {
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        // Throws because we've not called initialise
        TS_ASSERT_THROWS_THIS(bidomain_problem.Solve(), "Cardiac tissue is null, Initialise() probably hasn\'t been called");

        // Throws because mesh filename is unset
        TS_ASSERT_THROWS_CONTAINS(bidomain_problem.Initialise(),
                "No mesh given: define it in XML parameters file or call SetMesh()\n"
                "No XML element Simulation/Mesh found in parameters when calling");
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1mm_10_elements");
        TS_ASSERT_THROWS_NOTHING(bidomain_problem.Initialise());

        // Negative simulation duration
        HeartConfig::Instance()->SetSimulationDuration(-1.0);  //ms
        TS_ASSERT_THROWS_THIS(bidomain_problem.Solve(), "End time should be in the future");
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms

        // set output data to avoid their exceptions (which is covered in TestMonoDg0Assembler
        HeartConfig::Instance()->SetOutputDirectory("temp");
        HeartConfig::Instance()->SetOutputFilenamePrefix("temp");

        // Exception caused by relative tolerance and no clamping
        HeartConfig::Instance()->SetUseRelativeTolerance(2e-3);
        TS_ASSERT_THROWS_THIS(bidomain_problem.Solve(), "Bidomain external voltage is not bounded in this simulation - use KSP *absolute* tolerance");
        HeartConfig::Instance()->SetUseAbsoluteTolerance(2e-3);

        // Throws (in AbstractCardiacProblem) as dt does not divide end time
        HeartConfig::Instance()->SetPrintingTimeStep(0.15);
        HeartConfig::Instance()->SetPdeTimeStep(0.15);
        TS_ASSERT_THROWS_THIS( bidomain_problem.Solve(),"PDE timestep does not seem to divide end time - check parameters" );
        HeartConfig::Instance()->SetPdeTimeStep(0.01);
        HeartConfig::Instance()->SetPrintingTimeStep(0.01);

        //Throws because the node number is slightly bigger than the number of nodes in the mesh
        std::vector<unsigned> too_large;
        too_large.push_back(4358743);
        bidomain_problem.SetFixedExtracellularPotentialNodes(too_large);
        TS_ASSERT_THROWS_THIS(bidomain_problem.Solve(), "Fixed node number must be less than total number nodes");

        //Explicitly reset the counters so the next test in the test suite doesn't find on
        HeartEventHandler::Reset();
    }

    void TestCompareOrthotropicWithAxisymmetricBidomain() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetMeshFileName("heart/test/data/box_shaped_heart/box_heart", cp::media_type::Orthotropic);
        HeartConfig::Instance()->SetOutputDirectory("OrthotropicBidomain");
        HeartConfig::Instance()->SetOutputFilenamePrefix("ortho3d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 3> cell_factory;

        ///////////////////////////////////////////////////////////////////
        // orthotropic
        ///////////////////////////////////////////////////////////////////

        BidomainProblem<3> orthotropic_bido( &cell_factory );

        orthotropic_bido.Initialise();
        orthotropic_bido.Solve();

        ///////////////////////////////////////////////////////////////////
        // axisymmetric
        ///////////////////////////////////////////////////////////////////
        HeartConfig::Instance()->SetMeshFileName("heart/test/data/box_shaped_heart/box_heart", cp::media_type::Axisymmetric);
        HeartConfig::Instance()->SetOutputDirectory("AxisymmetricBidomain");
        HeartConfig::Instance()->SetOutputFilenamePrefix("axi3d");

        BidomainProblem<3> axisymmetric_bido( &cell_factory);
        HeartConfig::Instance()->SetVisualizeWithCmgui(true);
        HeartConfig::Instance()->SetVisualizeWithVtk(true);

        axisymmetric_bido.Initialise();
        axisymmetric_bido.Solve();

        ///////////////////////////////////////////////////////////////////
        // compare
        ///////////////////////////////////////////////////////////////////
        DistributedVector orthotropic_solution = orthotropic_bido.GetSolutionDistributedVector();
        DistributedVector axisymmetric_solution = axisymmetric_bido.GetSolutionDistributedVector();

        DistributedVector::Stripe ortho_voltage(orthotropic_solution, 0);
        DistributedVector::Stripe axi_voltage(axisymmetric_solution, 0);

        DistributedVector::Stripe ortho_ex_pot(orthotropic_solution, 1);
        DistributedVector::Stripe axi_ex_pot(axisymmetric_solution, 1);

        for (DistributedVector::Iterator index = orthotropic_solution.Begin();
             index != orthotropic_solution.End();
             ++index)
        {
            TS_ASSERT_DELTA(ortho_voltage[index], axi_voltage[index], 1e-7);
            TS_ASSERT_DELTA(ortho_ex_pot[index], axi_ex_pot[index], 1e-7);
        }

        // We check that the cmgui files generated by the convert command in the problem class are OK
        // We compare against mesh files and one data file that are known to be visualized correctly in Cmgui.

        // The files written in parallel are different from the ones written in sequential because of the different node
        // renumbering, therefore we test only the sequential case.
        // Note that the outputs of sequential and parallel simulations look the same when loaded with cmgui.
        // There are also minor rounding differences at the last decimal figure between sequential and parallel.
        EXIT_IF_PARALLEL

        std::string results_dir = OutputFileHandler::GetChasteTestOutputDirectory() + "AxisymmetricBidomain/cmgui_output/";

         //the mesh files...
        std::string node_file1 = results_dir + "/axi3d.exnode";
        std::string node_file2 = "heart/test/data/CmguiData/bidomain/bidomain3dValid.exnode";
        std::string elem_file1 = results_dir + "/axi3d.exelem";
        std::string elem_file2 = "heart/test/data/CmguiData/bidomain/bidomain3dValid.exelem";

        FileComparison comparer(node_file1,node_file2);
        TS_ASSERT(comparer.CompareFiles());

        FileComparison comparer2(elem_file1,elem_file2);
        TS_ASSERT(comparer2.CompareFiles());

        //...and a couple of data files as examples
        FileComparison comparer3(results_dir + "/axi3d_61.exnode",
                "heart/test/data/CmguiData/bidomain/bidomain3dValidData61.exnode");
        TS_ASSERT(comparer3.CompareFiles());

        FileComparison comparer4(results_dir + "/axi3d_100.exnode",
                "heart/test/data/CmguiData/bidomain/bidomain3dValidData100.exnode");
        TS_ASSERT(comparer4.CompareFiles());

        //info file
        FileComparison comparer5(results_dir + "/axi3d_times.info",
                "heart/test/data/CmguiData/bidomain/axi3d_times.info");
        TS_ASSERT(comparer5.CompareFiles());

        //HeartConfig XML
        TS_ASSERT(FileFinder(results_dir + "ChasteParameters.xml").Exists());

#ifdef CHASTE_VTK
// Requires  "sudo aptitude install libvtk5-dev" or similar
        results_dir = OutputFileHandler::GetChasteTestOutputDirectory() + "AxisymmetricBidomain/vtk_output/";

        VtkMeshReader<3,3> mesh_reader(results_dir + "axi3d.vtu");
        TS_ASSERT_EQUALS( mesh_reader.GetNumNodes(), 50U);
        TS_ASSERT_EQUALS( mesh_reader.GetNumElements(), 139U);

        std::vector<double> first_node = mesh_reader.GetNextNode();
        TS_ASSERT_DELTA( first_node[0] , 0.0 , 1e-6 );
        TS_ASSERT_DELTA( first_node[1] , 0.0, 1e-6 );
        TS_ASSERT_DELTA( first_node[2] , 0.0 , 1e-6 );

        std::vector<double> next_node = mesh_reader.GetNextNode();
        TS_ASSERT_DELTA( next_node[0] , 0.0 , 1e-6 );
        TS_ASSERT_DELTA( next_node[1] , 5.0 , 1e-6 );
        TS_ASSERT_DELTA( next_node[2] , 0.0 , 1e-6 );

        //V_m and phi_e samples
        std::vector<double> v_at_last, v_at_50, phi_at_50;
        mesh_reader.GetPointData( "V_000050", v_at_50);
        TS_ASSERT_DELTA( v_at_50[0],  -83.6547, 1e-3 );
        TS_ASSERT_DELTA( v_at_50[25], -83.8549, 1e-3 );
        TS_ASSERT_DELTA( v_at_50[49], -83.8557, 1e-3 );
        mesh_reader.GetPointData( "Phi_e_000050", phi_at_50);
        TS_ASSERT_DELTA( phi_at_50[0],  -0.0332, 1e-3 );
        TS_ASSERT_DELTA( phi_at_50[25],  0.0064, 1e-3 );
        TS_ASSERT_DELTA( phi_at_50[49],  0.0062, 1e-3 );
        mesh_reader.GetPointData( "V_000100", v_at_last);
        TS_ASSERT_DELTA( v_at_last[0],   -83.6863, 1e-3 );

        //HeartConfig XML
        TS_ASSERT(FileFinder(results_dir + "ChasteParameters.xml").Exists());
#else
        std::cout << "This test ran, but did not test VTK-dependent functions as VTK visualization is not enabled." << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_VTK
     }

    // Test the functionality for outputting the values of requested cell state variables
    void TestBidomainProblemPrintsMultipleVariables() throw (Exception)
    {
        // Get the singleton in a clean state
        HeartConfig::Instance()->Reset();

        // Set configuration file
        HeartConfig::Instance()->SetParametersFile("heart/test/data/xml/MultipleVariablesBidomain.xml");

        // Override the variables we are interested in writing.
        std::vector<std::string> output_variables;
        output_variables.push_back("calcium_dynamics__Ca_NSR");
        output_variables.push_back("ionic_concentrations__Nai");
        output_variables.push_back("fast_sodium_current_j_gate__j");
        output_variables.push_back("ionic_concentrations__Ki");

        HeartConfig::Instance()->SetOutputVariables( output_variables );

        // Set up problem
        PlaneStimulusCellFactory<CellFaberRudy2000FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        // Solve
        bidomain_problem.Initialise();
        bidomain_problem.Solve();

        // Get a reference to a reader object for the simulation results
        Hdf5DataReader data_reader1 = bidomain_problem.GetDataReader();
        std::vector<double> times = data_reader1.GetUnlimitedDimensionValues();

        // Check there is information about 11 timesteps (0, 0.01, 0.02, ...)
        unsigned num_steps = 11u;
        TS_ASSERT_EQUALS( times.size(), num_steps);
        TS_ASSERT_DELTA( times[0], 0.0, 1e-12);
        TS_ASSERT_DELTA( times[1], 0.01, 1e-12);
        TS_ASSERT_DELTA( times[2], 0.02, 1e-12);
        TS_ASSERT_DELTA( times[3], 0.03, 1e-12);

        // There should be 11 values per variable and node.
        std::vector<double> node_5_v = data_reader1.GetVariableOverTime("V", 5);
        TS_ASSERT_EQUALS( node_5_v.size(), num_steps);

        std::vector<double> node_5_phi = data_reader1.GetVariableOverTime("Phi_e", 5);
        TS_ASSERT_EQUALS( node_5_phi.size(), num_steps);

        for (unsigned i=0; i<output_variables.size(); i++)
        {
            unsigned global_index = 2+i*2;
            std::vector<double> values = data_reader1.GetVariableOverTime(output_variables[i], global_index);
            TS_ASSERT_EQUALS( values.size(), num_steps);

            // Check the last values match the cells' state
            if (bidomain_problem.rGetMesh().GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index))
            {
                AbstractCardiacCellInterface* p_cell = bidomain_problem.GetTissue()->GetCardiacCell(global_index);
                TS_ASSERT_DELTA(values.back(), p_cell->GetAnyVariable(output_variables[i],0), 1e-12);
            }
        }
    }

    /* HOW_TO_TAG Cardiac/Output
     * Output all cell model state variables for the cell model used in a particular simulation
     */
    void TestBidomainProblemPrintsAllStateVariables() throw (Exception)
    {
        // Get the singleton in a clean state
        HeartConfig::Instance()->Reset();

        // Set configuration file
        HeartConfig::Instance()->SetParametersFile("heart/test/data/xml/MultipleVariablesBidomain.xml");

        //Override the output directory so that this test goes into a fresh place
        HeartConfig::Instance()->SetOutputDirectory("AllStateVariablesBidomain");

        {
            //MultipleVariablesBidomain.xml uses FaberRudy.  Interrogate a model to get all the state variables
            boost::shared_ptr<ZeroStimulus> p_stimulus(new ZeroStimulus());
            boost::shared_ptr<EulerIvpOdeSolver> p_solver(new EulerIvpOdeSolver);
            CellFaberRudy2000FromCellML temporary_fr2000_ode_system(p_solver, p_stimulus);
            TS_ASSERT_EQUALS(temporary_fr2000_ode_system.rGetStateVariableNames().size(), 25u);

            // Copy all the state variable names from the temporary cell into HeartConfig so that they are all printed
            HeartConfig::Instance()->SetOutputVariables( temporary_fr2000_ode_system.rGetStateVariableNames() );
        }

        // Set up problem
        PlaneStimulusCellFactory<CellFaberRudy2000FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        // Solve
        bidomain_problem.Initialise();
        bidomain_problem.Solve();

        // Get a reference to a reader object for the simulation results
        Hdf5DataReader data_reader1 = bidomain_problem.GetDataReader();

        // Get an obscure state variable
        std::vector<double> obscure_at_5 = data_reader1.GetVariableOverTime("slow_delayed_rectifier_potassium_current_xs2_gate__xs2", 5);
        TS_ASSERT_EQUALS( obscure_at_5.size(), 11u);
        for (unsigned time_step=0; time_step<11; time_step++)
        {
            TS_ASSERT_DELTA( obscure_at_5[time_step], 0.0044, 1e-3); //Does not change over time at this node
        }

    }

    /*
     * Simple bidomain simulation to test against in the archiving tests below
     *
     * NOTE: This test uses NON-PHYSIOLOGICAL parameters values (conductivities,
     * surface-area-to-volume ratio, capacitance, stimulus amplitude). Essentially,
     * the equations have been divided through by the surface-area-to-volume ratio.
     * (Historical reasons...)
     */
    void TestSimpleBidomain1D() throw(Exception)
    {
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.1);
        HeartConfig::Instance()->SetSimulationDuration(2.0); //ms
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1mm_10_elements");
        HeartConfig::Instance()->SetOutputDirectory("BidomainSimple1dInOneGo");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        bidomain_problem.Initialise();
        bidomain_problem.Solve();

        // Check some voltages
        ReplicatableVector solution_replicated(bidomain_problem.GetSolution());

        double atol=5e-3;

        TS_ASSERT_DELTA(solution_replicated[1], -16.4861, atol);
        TS_ASSERT_DELTA(solution_replicated[2], 22.8117, atol);
        TS_ASSERT_DELTA(solution_replicated[3], -16.4893, atol);
        TS_ASSERT_DELTA(solution_replicated[5], -16.5617, atol);
        TS_ASSERT_DELTA(solution_replicated[7], -16.6761, atol);
        TS_ASSERT_DELTA(solution_replicated[9], -16.8344, atol);
        TS_ASSERT_DELTA(solution_replicated[10], 25.3148, atol);

        // Save solution for later comparison
        for (unsigned index=0; index<solution_replicated.GetSize(); index++)
        {
            mSolutionReplicated1d2ms.push_back(solution_replicated[index]);
        }

    }

    /*
     * Simple bidomain simulation with a straight permutation applied.
     * HOW_TO_TAG Cardiac/Output
     * Use the `SingleTraceOutputModifier` to output based on a global index
     * (index AFTER any permutation has been applied)
     *
     * NOTE: This test uses NON-PHYSIOLOGICAL parameters values.
     */
    void TestPermutedBidomain1D() throw(Exception)
    {

        TetrahedralMesh<1,1> mesh;
        TrianglesMeshReader<1,1> reader("mesh/test/data/1D_0_to_1mm_10_elements");
        mesh.ConstructFromMeshReader(reader);

        std::vector<unsigned> rotation_perm;

        unsigned number_nodes=11;
        for (unsigned index=0; index<(unsigned)number_nodes; index++)
        {
            rotation_perm.push_back( (index + 3) % number_nodes); // 3, 4, ... 0, 1, 2
        }

        //Rotate the permutation
        mesh.PermuteNodes(rotation_perm);

        HeartConfig::Instance()->SetOutputUsingOriginalNodeOrdering(true);

        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);

        HeartConfig::Instance()->SetSimulationDuration(0.5); //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainUnpermuted1d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");
        HeartConfig::Instance()->SetVisualizeWithMeshalyzer();

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        // Single trace at node 5
        unsigned new_index_for_5 = mesh.rGetNodePermutation()[5];
        TS_ASSERT_EQUALS(new_index_for_5, 8u);
        ChastePoint<1> point(0.05);
        unsigned index_for_location = mesh.GetNearestNodeIndex(point);
        TS_ASSERT_EQUALS(index_for_location, 8u);

        boost::shared_ptr<SingleTraceOutputModifier> trace_5(new SingleTraceOutputModifier("trace_5.txt", new_index_for_5, 0.1));
        bidomain_problem.AddOutputModifier(trace_5);

        bidomain_problem.SetMesh(&mesh);
        bidomain_problem.Initialise();
        bidomain_problem.Solve();

        // Extra trace data
        OutputFileHandler handler1("BidomainUnpermuted1d", false);
        NumericFileComparison comp1(handler1.GetOutputDirectoryFullPath()+ "trace_5.txt", "heart/test/data/BidomainUnpermuted1d/trace_5.txt");
        TS_ASSERT(comp1.CompareFiles(5e-4));

        //Can't read in the final mesh since it's a 1d example...
        OutputFileHandler handler("BidomainUnpermuted1d/output", false);

        //Mesh
        FileComparison comparer(handler.GetOutputDirectoryFullPath() + "/BidomainLR91_1d_mesh.pts",
                "heart/test/data/BidomainUnpermuted1d/BidomainLR91_1d_mesh.pts");
        TS_ASSERT(comparer.CompareFiles());

        //Transmembrane
        std::string file1=handler.GetOutputDirectoryFullPath()+ "/BidomainLR91_1d_V.dat";
        std::string file2="heart/test/data/BidomainUnpermuted1d/BidomainLR91_1d_V.dat";
        NumericFileComparison comp(file1, file2);
        TS_ASSERT(comp.CompareFiles(3e-3)); //This can be quite flexible since the permutation differences will be quite large
    }

    /*
     * This test is almost identical to TestSimpleBidomain1D
     * and relies on that test generating a h5 file to check against.
     *
     * NOTE: This test uses NON-PHYSIOLOGICAL parameters values (conductivities,
     * surface-area-to-volume ratio, capacitance, stimulus amplitude). Essentially,
     * the equations have been divided through by the surface-area-to-volume ratio.
     * (Historical reasons...)
     */
    void TestBidomainProblemInTwoHalves() throw (Exception)
    {
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.1);
        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1mm_10_elements");
        HeartConfig::Instance()->SetOutputDirectory("BidomainSimple1dInTwoHalves");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
        BidomainProblem<1> bidomain_problem( &cell_factory );

        bidomain_problem.Initialise();

        HeartConfig::Instance()->SetSimulationDuration(1.0);
        bidomain_problem.Solve();
        TS_ASSERT_DELTA(bidomain_problem.GetCurrentTime(), 1.0, 1e-12);

        HeartConfig::Instance()->SetSimulationDuration(2.0);
        bidomain_problem.Solve();
        TS_ASSERT_DELTA(bidomain_problem.GetCurrentTime(), 2.0, 1e-12);

        // Check some voltages
        ReplicatableVector solution_replicated(bidomain_problem.GetSolution());

        double atol=5e-3;

        TS_ASSERT_DELTA(solution_replicated[1], -16.4861, atol);
        TS_ASSERT_DELTA(solution_replicated[2], 22.8117, atol);
        TS_ASSERT_DELTA(solution_replicated[3], -16.4893, atol);
        TS_ASSERT_DELTA(solution_replicated[5], -16.5617, atol);
        TS_ASSERT_DELTA(solution_replicated[7], -16.6761, atol);
        TS_ASSERT_DELTA(solution_replicated[9], -16.8344, atol);
        TS_ASSERT_DELTA(solution_replicated[10], 25.3148, atol);

        // Compare against saved solution
        for (unsigned index=0; index<solution_replicated.GetSize(); index++)
        {
            TS_ASSERT_DELTA(solution_replicated[index], mSolutionReplicated1d2ms[index], 5e-11);
        }

        // Check output file contains results for the whole simulation and agree with normal test
        TS_ASSERT(CompareFilesViaHdf5DataReader("BidomainSimple1dInTwoHalves", "BidomainLR91_1d", true,
                                                "BidomainSimple1dInOneGo", "BidomainLR91_1d", true));

        // Test that we can keep solving even if the results have been deleted (i.e. by creating a new
        // .h5 file when we realize that there isn't one to extend)
        HeartConfig::Instance()->SetOutputDirectory("BidomainSimple1dInTwoHalves_2");
        OutputFileHandler file_handler("BidomainSimple1dInTwoHalves_2", true);
        FileFinder h5_file = file_handler.FindFile("BidomainLR91_1d.h5");
        TS_ASSERT(!h5_file.Exists());
        HeartConfig::Instance()->SetSimulationDuration(3.0);
        TS_ASSERT_THROWS_NOTHING( bidomain_problem.Solve() );
        TS_ASSERT(h5_file.Exists());
    }

    /**
     * Not a very thorough test yet - just checks we can load a problem, simulate it, and
     * get expected results.
     *
     * This test relies on the h5 file generated in TestSimpleBidomain1D. Always run after!
     */
    void TestArchiving() throw(Exception)
    {
        FileFinder archive_dir("bidomain_problem_archive", RelativeTo::ChasteTestOutput);
        std::string archive_file = "bidomain_problem.arch";

        // Values to test against after load
        unsigned num_cells;

        // Save
        {
            HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(0.0005));
            HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(0.0005));
            HeartConfig::Instance()->SetMeshFileName("mesh/test/data/1D_0_to_1mm_10_elements");
            HeartConfig::Instance()->SetOutputDirectory("BiProblemArchive");
            HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");
            HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
            HeartConfig::Instance()->SetCapacitance(1.0);
            HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.1);

            PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory;
            BidomainProblem<1> bidomain_problem( &cell_factory );

            bidomain_problem.Initialise();
            HeartConfig::Instance()->SetSimulationDuration(1.0); //ms
            bidomain_problem.Solve();

            num_cells = bidomain_problem.GetTissue()->rGetCellsDistributed().size();

            ArchiveOpener<boost::archive::text_oarchive, std::ofstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_oarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacProblem<1,1,2>* const p_bidomain_problem = &bidomain_problem;
            (*p_arch) & p_bidomain_problem;
        }

        // Load
        {
            ArchiveOpener<boost::archive::text_iarchive, std::ifstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_iarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacProblem<1,1,2> *p_bidomain_problem;
            (*p_arch) >> p_bidomain_problem;

            // Check values
            TS_ASSERT_EQUALS(p_bidomain_problem->GetTissue()->rGetCellsDistributed().size(),
                             num_cells);

            HeartConfig::Instance()->SetSimulationDuration(2.0); //ms
            p_bidomain_problem->Solve();

            // check some voltages
            ReplicatableVector solution_replicated(p_bidomain_problem->GetSolution());
            double atol=5e-3;
            TS_ASSERT_DELTA(solution_replicated[1], -16.4861, atol);
            TS_ASSERT_DELTA(solution_replicated[2], 22.8117, atol);
            TS_ASSERT_DELTA(solution_replicated[3], -16.4893, atol);
            TS_ASSERT_DELTA(solution_replicated[5], -16.5617, atol);
            TS_ASSERT_DELTA(solution_replicated[7], -16.6761, atol);
            TS_ASSERT_DELTA(solution_replicated[9], -16.8344, atol);
            TS_ASSERT_DELTA(solution_replicated[10], 25.3148, atol);

            for (unsigned index=0; index<solution_replicated.GetSize(); index++)
            {
                //Shouldn't differ from the original run at all
                TS_ASSERT_DELTA(solution_replicated[index], mSolutionReplicated1d2ms[index],  5e-11);
            }
            // check output file contains results for the whole simulation
            TS_ASSERT(CompareFilesViaHdf5DataReader("BiProblemArchive", "BidomainLR91_1d", true,
                                                    "BidomainSimple1dInOneGo", "BidomainLR91_1d", true));

            // Free memory
            delete p_bidomain_problem;
        }
    }

};

#endif /*TESTBIDOMAINPROBLEM_HPP_*/

max_stars_repo_name: ktunya/ChasteMod

max_issues_repo_name: ktunya/ChasteMod

max_forks_repo_name: ktunya/ChasteMod


text:
/*
this file should cover the following functions in hmm.h
- hmm_vec2 HMM_MultiplyVec2(hmm_vec2 Left, hmm_vec2 Right);
- hmm_vec2 HMM_MultiplyVec2f(hmm_vec2 Left, float Right);
- hmm_vec3 HMM_MultiplyVec3(hmm_vec3 Left, hmm_vec3 Right);
- hmm_vec3 HMM_MultiplyVec3f(hmm_vec3 Left, float Right);
- hmm_vec4 HMM_MultiplyVec4(hmm_vec4 Left, hmm_vec4 Right);
- hmm_vec4 HMM_MultiplyVec4f(hmm_vec4 Left, float Right);
- hmm_mat4 HMM_MultiplyMat4(hmm_mat4 Left, hmm_mat4 Right);
- hmm_mat4 HMM_MultiplyMat4f(hmm_mat4 Matrix, float Scalar);
- hmm_vec4 HMM_MultiplyMat4ByVec4(hmm_mat4 Matrix, hmm_vec4 Vector);

take care of operator overloading also for these functions:
- hmm_vec2 operator*(hmm_vec2 Left, hmm_vec2 Right);
- hmm_vec3 operator*(hmm_vec3 Left, hmm_vec3 Right);
- hmm_vec4 operator*(hmm_vec4 Left, hmm_vec4 Right);
- hmm_mat4 operator*(hmm_mat4 Left, hmm_mat4 Right);
- hmm_vec2 operator*(hmm_vec2 Left, float Right);
- hmm_vec3 operator*(hmm_vec3 Left, float Right);
- hmm_vec4 operator*(hmm_vec4 Left, float Right);
- hmm_mat4 operator*(hmm_mat4 Left, float Right);
- hmm_vec2 operator*(float Left, hmm_vec2 Right);
- hmm_vec3 operator*(float Left, hmm_vec3 Right);
- hmm_vec4 operator*(float Left, hmm_vec4 Right);
- hmm_mat4 operator*(float Left, hmm_mat4 Right);
- hmm_vec4 operator*(hmm_mat4 Matrix, hmm_vec4 Vector);
- hmm_vec2 &operator*=(hmm_vec2 &Left, hmm_vec2 Right);
- hmm_vec3 &operator*=(hmm_vec3 &Left, hmm_vec3 Right);
- hmm_vec4 &operator*=(hmm_vec4 &Left, hmm_vec4 Right);
- hmm_vec2 &operator*=(hmm_vec2 &Left, float Right);
- hmm_vec3 &operator*=(hmm_vec3 &Left, float Right);
- hmm_vec4 &operator*=(hmm_vec4 &Left, float Right);
- hmm_mat4 &operator*=(hmm_mat4 &Left, float Right);
*/
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include "HMM.h"
#include "test_helpers.h"
namespace utf = boost::unit_test;
namespace tt = boost::test_tools;
using namespace std;

BOOST_AUTO_TEST_SUITE(multiplication_test_suite)

BOOST_AUTO_TEST_SUITE(multiplication_vec2_by_int, *utf::label("trivial"))

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2_int_BothVecsZeros, *utf::tolerance(0.00001)) {
    //Arrange
    auto zeros = HMM_Vec2i(0, 0);
    //Act
    auto res = HMM_MultiplyVec2(zeros, zeros);
    //Assert
    float expectedRes[2] = {0, 0};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2_int_V1_zeros, *utf::tolerance(0.00001)) {
    //Arrange
    auto zeros = HMM_Vec2i(0, 0);
    auto right = HMM_Vec2i(5, -10);
    //Act
    auto res = HMM_MultiplyVec2(zeros, right);
    //Assert
    float expectedRes[2] = {0, 0};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2_int_postives, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2i(1, 5);
    auto right = HMM_Vec2i(3, 4);
    //Act
    auto res = HMM_MultiplyVec2(left, right);
    //Assert
    float expectedRes[2] = {left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2_int_pos_neg, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2i(1, 5);
    auto right = HMM_Vec2i(-3, -4);
    //Act
    auto res = HMM_MultiplyVec2(left, right);
    //Assert
    float expectedRes[2] = {left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]};
    BOOST_TEST(vector<float>(res.Elements, res.Elements + 2) == vector<float>({left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]}));
}

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2_float_pos_neg, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(.5, .2);
    auto right = HMM_Vec2(-2, -10.);
    //Act
    auto res = HMM_MultiplyVec2(left, right);
    //Assert
    float expectedRes[2] = {left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_SUITE_END()
//--------------------------------------------------------------

BOOST_AUTO_TEST_SUITE(multiplication_vec2_by_float)

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyVec2f_vec_x_posfloat, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(.5, .2);
    float right = 10.;
    //Act
    auto res = HMM_MultiplyVec2f(left, right);
    //Assert
    float expectedRes[2] = {left.Elements[0] * right, left.Elements[1] * right};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_SUITE_END()
//--------------------------------------------------------------
BOOST_AUTO_TEST_SUITE(mul_vec2_using_operators)

BOOST_AUTO_TEST_CASE(test_mul_vec2_by_vec2_using_mul_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(10., 1.);
    auto right = HMM_Vec2(-2, -5.);
    //Act
    auto res = left * right;
    //Assert
    float expectedRes[2] = {left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]};
    BOOST_TEST(res.Elements == expectedRes);
}

BOOST_AUTO_TEST_CASE(test_mul_vec2_by_float_using_mul_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(.5, .2);
    float right = 10.;
    //Act
    auto res = left * right;
    //Assert
    BOOST_TEST(vector<float>(res.Elements, res.Elements + 2) == vector<float>({left.Elements[0] * right, left.Elements[1] * right}));
}

BOOST_AUTO_TEST_CASE(test_mul_vec2_by_vec2_using_mul_equal_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(.5, .2);
    auto right = HMM_Vec2(-2, -5.);
    //Act
    auto res = left;
    res *= right;
    //Assert
    BOOST_TEST(vector<float>(res.Elements, res.Elements + 2) == vector<float>({left.Elements[0] * right.Elements[0], left.Elements[1] * right.Elements[1]}));
}

BOOST_AUTO_TEST_CASE(test_mul_vec2_by_float_using_mul_equal_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = HMM_Vec2(.5, .2);
    float right = 10.;
    //Act
    auto res = left;
    res *= right;
    //Assert
    BOOST_TEST(vector<float>(res.Elements, res.Elements + 2) == vector<float>({left.Elements[0] * right, left.Elements[1] * right}));
}

BOOST_AUTO_TEST_SUITE_END()

//--------------------------------------------------------------

BOOST_AUTO_TEST_SUITE(multiplication_mat4_by_mat4)

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyMat4, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    //Act
    hmm_mat4 res = HMM_MultiplyMat4(left, left);

    //Assert
    float expMat[4][4] = {{38, 80, 50, 56},
                          {38, 236, 50, 56},
                          {158, 248, 218, 248},
                          {75, -90, 105, 120}};

    BOOST_TEST(res.Elements == expMat);
}

BOOST_AUTO_TEST_SUITE_END()

//--------------------------------------------------------------

BOOST_AUTO_TEST_SUITE(multiplication_mat4_by_floatscalar)

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyMat4f, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    float right = .5;
    //Act
    hmm_mat4 res = HMM_MultiplyMat4f(left, right);

    //Assert
    float expMat[4][4] = {{arr[0][0] * right, arr[0][1] * right, arr[0][2] * right, arr[0][3] * right},
                          {arr[1][0] * right, arr[1][1] * right, arr[1][2] * right, arr[1][3] * right},
                          {arr[2][0] * right, arr[2][1] * right, arr[2][2] * right, arr[2][3] * right},
                          {arr[3][0] * right, arr[3][1] * right, arr[3][2] * right, arr[3][3] * right}};

    BOOST_TEST(res.Elements == expMat);
}

BOOST_AUTO_TEST_SUITE_END()

//--------------------------------------------------------------

BOOST_AUTO_TEST_SUITE(multiplication_mat4_by_vec4)

BOOST_AUTO_TEST_CASE(test_HMM_MultiplyMat4ByVec4, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    hmm_vec4 right = HMM_Vec4(.5, 1., 0., -3.);
    //Act
    hmm_vec4 res = HMM_MultiplyMat4ByVec4(left, right);
    //Assert
    float expRes[] = {5.5, -50, 8.5, 10};
    BOOST_TEST(res.Elements == expRes);
}

BOOST_AUTO_TEST_SUITE_END()
//--------------------------------------------------------------

BOOST_AUTO_TEST_SUITE(multiplication_using_operators)

BOOST_AUTO_TEST_CASE(test_mul_mat4_by_mat4_using_mul_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    //Act
    hmm_mat4 res = left * left;

    //Assert
    float expMat[4][4] = {{38, 80, 50, 56},
                          {38, 236, 50, 56},
                          {158, 248, 218, 248},
                          {75, -90, 105, 120}};

    BOOST_TEST(res.Elements == expMat);
}

BOOST_AUTO_TEST_CASE(test_mul_mat4_by_floatscalar_using_mul_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    float right = .5;
    //Act
    hmm_mat4 res = left * right;
    //Assert
    float expMat[4][4] = {{arr[0][0] * right, arr[0][1] * right, arr[0][2] * right, arr[0][3] * right},
                          {arr[1][0] * right, arr[1][1] * right, arr[1][2] * right, arr[1][3] * right},
                          {arr[2][0] * right, arr[2][1] * right, arr[2][2] * right, arr[2][3] * right},
                          {arr[3][0] * right, arr[3][1] * right, arr[3][2] * right, arr[3][3] * right}};

    BOOST_TEST(res.Elements == expMat);
}

BOOST_AUTO_TEST_CASE(test_mul_mat4_by_vec4_using_mul_operator, *utf::tolerance(0.00001)) {
    //Arrange
    auto left = hmm_mat4();
    float arr[4][4] = {{1, 2, 3, 4},
                       {5, -6, 7, 8},
                       {9, 10, 11, 12},
                       {0, 15, 0, 0}};
    copy(&arr[0][0], &arr[0][0] + 4 * 4, &left.Elements[0][0]);

    hmm_vec4 right = HMM_Vec4(.5, 1., 0., -3.);
    //Act
    hmm_vec4 res = left * right;
    //Assert
    float epxRes[] = {5.5, -50, 8.5, 10};
    BOOST_TEST(res.Elements == epxRes);
}

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE_END()

ext: cpp

max_stars_repo_name: AmmarRabie/Test-HMM

max_issues_repo_name: AmmarRabie/Test-HMM

max_forks_repo_name: AmmarRabie/Test-HMM


text:
// Copyright Abel Sinkovics (abel@sinkovics.hu) 2013.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#include <mpllibs/metamonad/instantiate.hpp>
#include <mpllibs/metamonad/tmp_value.hpp>
#include <mpllibs/metamonad/apply.hpp>

#include <boost/test/unit_test.hpp>

#include <boost/mpl/vector.hpp>
#include <boost/mpl/assert.hpp>

#include <boost/type_traits.hpp>

using mpllibs::metamonad::tmp_value;

namespace
{
  template <class A, class B>
  struct some_template : tmp_value<some_template<A, B> > {};
}

BOOST_AUTO_TEST_CASE(test_instantiate)
{
  using mpllibs::metamonad::instantiate2;
  using mpllibs::metamonad::apply;

  using boost::mpl::vector;
  using boost::mpl::na;

  using boost::is_same;

  // test_instantiation
  BOOST_MPL_ASSERT((
    is_same<
      some_template<int, double>,
      apply<instantiate2<some_template>, vector<int, double> >::type
    >
  ));

  // test_instantiation_with_na
  BOOST_MPL_ASSERT((
    is_same<
      some_template<int, na>,
      apply<instantiate2<some_template>, vector<int> >::type
    >
  ));
}


text:
/* ---------------------------------------------------------------------
 *
 * Copyright (C) 1999 - 2019 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level of the deal.II distribution.
 *
 * ---------------------------------------------------------------------
 *
 * based on deal.II step-1
 */

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>

#include <cmath>
#include <fstream>
#include <iostream>

using namespace dealii;

//! Generate a hypercube, and output it as an svg file.
void first_grid(Triangulation<2> &triangulation)
{
  GridGenerator::hyper_cube(triangulation);
  triangulation.refine_global(4);

  std::ofstream out("grid-1.svg");
  GridOut       grid_out;
  grid_out.write_svg(triangulation, out);
  std::cout << "Grid written to grid-1.svg" << std::endl;
}

//! Generate a locally refined hyper_shell, and output it as an svg file.
void second_grid(Triangulation<2> &triangulation)
{
  const Point<2> center(1, 0);
  const double   inner_radius = 0.5, outer_radius = 1.0;
  GridGenerator::hyper_shell(
    triangulation, center, inner_radius, outer_radius, 10);

  // triangulation.reset_manifold(0);

  for (unsigned int step = 0; step < 5; ++step)
    {
      for (auto &cell : triangulation.active_cell_iterators())
        {
          for (const auto v : cell->vertex_indices())
            {
              const double distance_from_center =
                center.distance(cell->vertex(v));

              if (std::fabs(distance_from_center - inner_radius) <=
                  1e-6 * inner_radius)
                {
                  cell->set_refine_flag();
                  break;
                }
            }
        }

      triangulation.execute_coarsening_and_refinement();
    }

  std::ofstream out("grid-2.svg");
  GridOut       grid_out;
  grid_out.write_svg(triangulation, out);

  std::cout << "Grid written to grid-2.svg" << std::endl;
}

//! Create an L-shaped domain with one global refinement, and write it on
// `third_grid.vtk`.  Refine the L-shaped mesh adaptively around the re-entrant
// corner three times (after the global refinement you already did), but with a
// twist: refine all cells with the distance between the center of the cell and
// re-entrant corner is smaller than 1/3.
void third_grid(Triangulation<2> &tria)
{
  // Insert code here
  // last year code in comments
  /* Triangulation<2> tr1,tr2,tr3,tr_final;
   const Point<2> p1(0, 0);
   const Point<2> p2(1,2);
   const Point<2> p3(2,1);
   const Point<2> p4(1,0);
   const Point<2> p5(1,1);
   const Point<2> p6(0,1);
   GridGenerator::hyper_rectangle(tr1,p1,p5);
   GridGenerator::hyper_rectangle(tr2,p6,p2);
   GridGenerator::hyper_rectangle(tr3,p4,p3);
   //they cannot be refined!
   GridGenerator::merge_triangulations ({&tr1,&tr2,&tr3},tr_final);*/
  GridGenerator::hyper_L(tria);
  tria.refine_global(1);
  const Point<2> p5(0, 0);
  for (unsigned int step = 0; step < 3; ++step)
    {
      // Active cells are those that are not further refined
      // we need to mark cells for refinement
      for (auto &cell : tria.active_cell_iterators())
        {
          const double distance_from_corner = cell->center().distance({0, 0});
          // choose whatever refinement condition
          if (distance_from_corner < 2.0 / 3.0)
            {
              cell->set_refine_flag();
            }
        }
      // refine global calls this function too
      tria.execute_coarsening_and_refinement();
    } // for steps loop

  // tr_final.refine_global(2); //we want a nice picture :)

  std::ofstream file_var("grid-3.vtk");
  GridOut       grid_out;
  grid_out.write_vtk(tria, file_var);
  std::cout << "Grid written to grid-3.vtk" << std::endl;
}

//! Returns a tuple with number of levels, number of cells, number of active
// cells. Test this with all of  your meshes.
std::tuple<unsigned int, unsigned int, unsigned int>
get_info(const Triangulation<2> &tria)
{
  // Insert code here
  return std::make_tuple(tria.n_levels(),
                         tria.n_cells(),
                         tria.n_active_cells());
}

void
torus_grid()
{
  Triangulation<2, 3> tria;
  GridGenerator::torus(tria, 2, 1);
  tria.refine_global(2);

  std::ofstream out("grid-torus.vtk");
  GridOut       grid_out;
  grid_out.write_vtk(tria, out);
  std::cout << "Grid written to grid-torus.vtk" << std::endl;
}

int
main()
{
  Triangulation<2> triangulation;
  first_grid(triangulation);
  triangulation.clear();
  second_grid(triangulation);
  triangulation.clear();
  third_grid(triangulation);

  torus_grid();
}

ext: cc


text:
//==============================================================================
//         Copyright 2003 - 2011   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_DISPATCH_META_PRINT_HPP_INCLUDED
#define BOOST_DISPATCH_META_PRINT_HPP_INCLUDED

#include <boost/config.hpp>

/*!
 * \file
 * \brief Defines and implements the \ref boost::dispatch::meta::print \metafunction
 */

#if defined(__GNUC__) && !defined(BOOST_CLANG) && !defined(BOOST_INTEL) && !defined(__CUDACC__)

namespace boost { namespace dispatch { namespace meta
{
  template<class T>
  struct print
  {
    typedef T type;
    unsigned : 80;
  };
} } }

#else

#include <boost/mpl/print.hpp>
namespace boost { namespace dispatch { namespace meta
{
  using boost::mpl::print;
} } }

#endif

#endif

ext: hpp


text:
//
// Copyright 2005-2007 Adobe Systems Incorporated
//
// Distributed under the Boost Software License, Version 1.0
// See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt
//
#ifndef BOOST_GIL_EXTENSION_NUMERIC_PIXEL_NUMERIC_OPERATIONS_HPP
#define BOOST_GIL_EXTENSION_NUMERIC_PIXEL_NUMERIC_OPERATIONS_HPP

#include <boost/gil/extension/numeric/channel_numeric_operations.hpp>

#include <boost/gil/color_base_algorithm.hpp>
#include <boost/gil/pixel.hpp>

namespace boost { namespace gil {

// Function objects and utilities for pixel-wise numeric operations.
//
// List of currently defined functors:
//   pixel_plus_t (+)
//   pixel_minus_t (-)
//   pixel_multiplies_scalar_t (*)
//   pixel_divides_scalar_t (/)
//   pixel_halves_t (/=2),
//   pixel_zeros_t (=0)
//   pixel_assigns_t (=)

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise addition of two pixels.
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelRef2 - models PixelConcept
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef1, typename PixelRef2, typename PixelResult>
struct pixel_plus_t
{
    auto operator()(PixelRef1 const& p1, PixelRef2 const& p2) const -> PixelResult
    {
        PixelResult result;
        static_transform(p1, p2, result,
            channel_plus_t
            <
                typename channel_type<PixelRef1>::type,
                typename channel_type<PixelRef2>::type,
                typename channel_type<PixelResult>::type
            >());
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise subtraction of two pixels.
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelRef2 - models PixelConcept
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef1, typename PixelRef2, typename PixelResult>
struct pixel_minus_t
{
    auto operator()(PixelRef1 const& p1, PixelRef2 const& p2) const -> PixelResult
    {
        PixelResult result;
        static_transform(p1, p2, result,
            channel_minus_t
            <
                typename channel_type<PixelRef1>::type,
                typename channel_type<PixelRef2>::type,
                typename channel_type<PixelResult>::type
            >());
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise multiplication of pixel elements by scalar.
/// \tparam PixelRef - models PixelConcept
/// \tparam Scalar - models a scalar type
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef, typename Scalar, typename PixelResult>
struct pixel_multiplies_scalar_t
{
    auto operator()(PixelRef const& p, Scalar const& s) const -> PixelResult
    {
        PixelResult result;
        static_transform(p, result,
            std::bind(
                channel_multiplies_scalar_t<typename channel_type<PixelRef>::type,
                Scalar,
                typename channel_type<PixelResult>::type>(),
                std::placeholders::_1, s));
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise multiplication of two pixels.
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef1, typename PixelRef2, typename PixelResult>
struct pixel_multiply_t
{
    auto operator()(PixelRef1 const& p1, PixelRef2 const& p2) const -> PixelResult
    {
        PixelResult result;
        static_transform(p1, p2, result,
            channel_multiplies_t
            <
                typename channel_type<PixelRef1>::type,
                typename channel_type<PixelRef2>::type,
                typename channel_type<PixelResult>::type
            >());
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise division of pixel elements by scalar.
/// \tparam PixelRef - models PixelConcept
/// \tparam Scalar - models a scalar type
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef, typename Scalar, typename PixelResult>
struct pixel_divides_scalar_t
{
    auto operator()(PixelRef const& p, Scalar const& s) const -> PixelResult
    {
        PixelResult result;
        static_transform(p, result,
            std::bind(channel_divides_scalar_t<typename channel_type<PixelRef>::type,
                Scalar,
                typename channel_type<PixelResult>::type>(),
                std::placeholders::_1, s));
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise division of two pixels.
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelRef1 - models PixelConcept
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef1, typename PixelRef2, typename PixelResult>
struct pixel_divide_t
{
    auto operator()(PixelRef1 const& p1, PixelRef2 const& p2) const -> PixelResult
    {
        PixelResult result;
        static_transform(p1, p2, result,
            channel_divides_t
            <
                typename channel_type<PixelRef1>::type,
                typename channel_type<PixelRef2>::type,
                typename channel_type<PixelResult>::type
            >());
        return result;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Performs channel-wise division by 2
/// \tparam PixelRef - models PixelConcept
template <typename PixelRef>
struct pixel_halves_t
{
    auto operator()(PixelRef& p) const -> PixelRef&
    {
        static_for_each(p, channel_halves_t<typename channel_type<PixelRef>::type>());
        return p;
    }
};

/// \ingroup PixelNumericOperations
/// \brief Sets pixel elements to zero (for whatever zero means)
/// \tparam PixelRef - models PixelConcept
template <typename PixelRef>
struct pixel_zeros_t
{
    auto operator()(PixelRef& p) const -> PixelRef&
    {
        static_for_each(p, channel_zeros_t<typename channel_type<PixelRef>::type>());
        return p;
    }
};

/// \brief Sets pixel elements to zero (for whatever zero means)
/// \tparam Pixel - models PixelConcept
template <typename Pixel>
void zero_channels(Pixel& p)
{
    static_for_each(p, channel_zeros_t<typename channel_type<Pixel>::type>());
}

/// \ingroup PixelNumericOperations
/// \brief Casts and assigns a pixel to another
///
/// A generic implementation for casting and assigning a pixel to another.
/// User should specialize it for better performance.
///
/// \tparam PixelRef - models PixelConcept
/// \tparam PixelResult - models PixelValueConcept
template <typename PixelRef, typename PixelResult>
struct pixel_assigns_t
{
    auto operator()(PixelRef const& src, PixelResult& dst) const -> PixelResult
    {
        static_for_each(src, dst,
            channel_assigns_t
            <
                typename channel_type<PixelRef>::type,
                typename channel_type<PixelResult>::type
            >());
        return dst;
    }
};

}} // namespace boost::gil

#endif

max_stars_repo_name: kindlychung/mediasoup-sfu-cpp

max_issues_repo_name: kindlychung/mediasoup-sfu-cpp

max_forks_repo_name: kindlychung/mediasoup-sfu-cpp


text:
// System
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <fstream>
#include <vector>
#include <chrono>
#include <Eigen/Dense>
#include <pcl/point_cloud.h>
#include <pcl/common/transforms.h>
#include <pcl/filters/crop_box.h>

// Custom
#include <gpg/candidates_generator.h>
#include <gpg/hand_search.h>
#include <gpg/finger_hand.h>
#include <gpg/config_file.h>
#include <gpg/plot.h>

const double PI = M_PI;
const double PI_2 = M_PI*2;

// function to read in a double array from a single line of a configuration file
std::vector<double> stringToDouble(const std::string& str)
{
  std::vector<double> values;
  std::stringstream ss(str);
  double v;

  while (ss >> v)
  {
    values.push_back(v);
    if (ss.peek() == ' ')
    {
      ss.ignore();
    }
  }

  return values;
}

// function to load registrated object pose between mesh and point cloud
Eigen::Affine3d loadRegistration(const std::string& filepath)
{
  std::ifstream ifs(filepath, std::ifstream::in);
  std::string s;
  ifs >> s;
  Eigen::Matrix4d mesh2cloud;
  mesh2cloud.row(3) << 0,0,0,1;
  std::stringstream ss;
  ss << s;
  std::string substr;
  for (int i=0; i<12; ++i) {
      std::getline(ss, substr, ',');
      mesh2cloud.col(i/3).row(i%3) << std::stod(substr);
  }
  return Eigen::Affine3d(mesh2cloud);
}

Grasp transformGrasp(const Grasp& grasp, const FingerHand &finger_hand, Eigen::Affine3d transform)
{
  Eigen::Matrix3d I;
  I.setIdentity();
  Grasp grasp_I(transform*grasp.getSample(), I, finger_hand);
  grasp_I.setGraspWidth(grasp.getGraspWidth());
  grasp_I.setGraspBottom(transform*grasp.getGraspBottom());
  grasp_I.setGraspTop(transform*grasp.getGraspTop());
  grasp_I.setGraspSurface(transform*grasp.getGraspSurface());
  return grasp_I;
}

/*
inline double random_uniform(void)
{
  return std::max(1e-8, std::min(1.0-1e-8,(double)rand() / (double)RAND_MAX));
}

inline Eigen::Quaterniond randomQuaternion_sample(double u1, double u2, double u3)
{
  Eigen::Quaterniond random_quat(std::sqrt(1.-u1)*std::sin(PI_2*u2), std::sqrt(1.-u1)*std::cos(PI_2*u2), \
                                 std::sqrt(u1   )*std::sin(PI_2*u3), std::sqrt(u1   )*std::cos(PI_2*u3));
  return random_quat;
}

Eigen::Quaterniond randomQuaternion(void)
{
  // Random Rotation
  double u1 = random_uniform();
  double u2 = random_uniform();
  double u3 = random_uniform();
  return randomQuaternion_sample(u1, u2, u3);
}
*/

int main(int argc, char* argv[])
{
  srand(time(NULL));
  // Read parameters from configuration file.
  ConfigFile config_file(argv[1]);

  double finger_width = config_file.getValueOfKey<double>("finger_width", 0.01);
  double hand_outer_diameter  = config_file.getValueOfKey<double>("hand_outer_diameter", 0.12);
  double hand_depth = config_file.getValueOfKey<double>("hand_depth", 0.06);
  double hand_height  = config_file.getValueOfKey<double>("hand_height", 0.02);
  double init_bite  = config_file.getValueOfKey<double>("init_bite", 0.01);
  double friction_coeff = config_file.getValueOfKey<double>("friction_coeff", 20.0);
  int viable_thresh = config_file.getValueOfKey<int>("viable_thresh", 6);
  bool negative_sample = config_file.getValueOfKey<bool>("negative_sample", false);

  std::cout << "finger_width: " << finger_width << "\n";
  std::cout << "hand_outer_diameter: " << hand_outer_diameter << "\n";
  std::cout << "hand_depth: " << hand_depth << "\n";
  std::cout << "hand_height: " << hand_height << "\n";
  std::cout << "init_bite: " << init_bite << "\n";
  std::cout << "friction_coeff: " << friction_coeff << "\n";
  std::cout << "viable_thresh: " << viable_thresh << "\n";
  std::cout << "negative_sample: " << negative_sample << "\n";

  bool voxelize = config_file.getValueOfKey<bool>("voxelize", true);
  bool remove_outliers = config_file.getValueOfKey<bool>("remove_outliers", false);
  std::string workspace_str = config_file.getValueOfKeyAsString("workspace", "");
  std::string camera_pose_str = config_file.getValueOfKeyAsString("camera_pose", "");
  std::vector<double> workspace = stringToDouble(workspace_str);
  std::vector<double> camera_pose = stringToDouble(camera_pose_str);
  std::cout << "voxelize: " << voxelize << "\n";
  std::cout << "remove_outliers: " << remove_outliers << "\n";
  std::cout << "workspace: " << workspace_str << "\n";
  std::cout << "camera_pose: " << camera_pose_str << "\n";

  int num_samples = config_file.getValueOfKey<int>("num_samples", 1000);
  int num_threads = config_file.getValueOfKey<int>("num_threads", 1);
  double nn_radius = config_file.getValueOfKey<double>("nn_radius", 0.01);
  int num_orientations = config_file.getValueOfKey<int>("num_orientations", 8);
  int rotation_axis = config_file.getValueOfKey<int>("rotation_axis", 2);
  std::cout << "num_samples: " << num_samples << "\n";
  std::cout << "num_threads: " << num_threads << "\n";
  std::cout << "nn_radius: " << nn_radius << "\n";
  std::cout << "num_orientations: " << num_orientations << "\n";
  std::cout << "rotation_axis: " << rotation_axis << "\n";

  bool plot_grasps = config_file.getValueOfKey<bool>("plot_grasps", true);
  bool plot_candidates = config_file.getValueOfKey<bool>("plot_candidates", true);
  bool plot_normals = config_file.getValueOfKey<bool>("plot_normals", false);
  int max_samples = config_file.getValueOfKey<int>("max_samples", 10);
  std::cout << "plot_grasps: " << plot_grasps << "\n";
  std::cout << "plot_candidates: " << plot_candidates << "\n";
  std::cout << "plot_normals: " << plot_normals << "\n";
  std::cout << "max_samples: " << max_samples << "\n";

  // Create object to generate grasp candidates.
  CandidatesGenerator::Parameters generator_params;
  generator_params.num_samples_ = num_samples;
  generator_params.num_threads_ = num_threads;
  generator_params.plot_normals_ = plot_normals;
  generator_params.plot_grasps_ = plot_candidates;
  generator_params.remove_statistical_outliers_ = remove_outliers;
  generator_params.voxelize_ = voxelize;
  generator_params.workspace_ = workspace;
  HandSearch::Parameters hand_search_params;
  hand_search_params.finger_width_ = finger_width;
  hand_search_params.hand_outer_diameter_ = hand_outer_diameter;
  hand_search_params.hand_depth_ = hand_depth;
  hand_search_params.hand_height_ = hand_height;
  hand_search_params.init_bite_ = init_bite;
  hand_search_params.nn_radius_frames_ = nn_radius;
  hand_search_params.num_orientations_ = num_orientations;
  hand_search_params.num_samples_ = num_samples;
  hand_search_params.num_threads_ = num_threads;
  hand_search_params.rotation_axis_ = rotation_axis;
  hand_search_params.friction_coeff_ = friction_coeff;
  hand_search_params.viable_thresh_ = viable_thresh;
  CandidatesGenerator candidates_generator(generator_params, hand_search_params);
  HandSearch handsearch(hand_search_params);

  // Set the camera pose.
  Eigen::Matrix3Xd view_points(3,1);
  view_points << camera_pose[3], camera_pose[6], camera_pose[9];

  // Random Rotation
  /*
  for (double u1=1e-8; u1<1-1e-8; u1+=0.05)
    for (double u2=1e-8; u2<1-1e-8; u2+=0.05)
      for (double u3=1e-8; u3<1-1e-8; u3+=0.05)
      {
        Eigen::Quaterniond random_quat(randomQuaternion());
        std::cout << "Random Euler Angle: " << random_quat.toRotationMatrix().eulerAngles(0, 1, 2).transpose() << std::endl;
      }
  */

  Eigen::Affine3d mesh2cloud(loadRegistration(argv[3]));
  // pcl::transformPointCloud(*mesh_cam.getCloudOriginal(), *new_pc, InvTrans);

  // Create object to load point cloud from file.
  // CloudCamera cloud_cam(argv[3], view_points);
  PointCloudRGB::Ptr mesh_cam_pts = CloudCamera::loadPointCloudFromFile(argv[2]);
  pcl::transformPointCloud(*mesh_cam_pts, *mesh_cam_pts, mesh2cloud);
  CloudCamera mesh_cam(mesh_cam_pts, 0, view_points);
  if (mesh_cam.getCloudOriginal()->size() == 0)
  {
    std::cout << "mesh_cam: Input point cloud is empty or does not exist!\n";
    return (-1);
  }
  /*
  if (cloud_cam.getCloudOriginal()->size() == 0)
  {
    std::cout << "cloud_cam: Input point cloud is empty or does not exist!\n";
    return (-1);
  }
  */
  std::chrono::steady_clock::time_point total_begin = std::chrono::steady_clock::now();

  // Point cloud preprocessing: voxelize, remove statistical outliers, workspace filter, compute normals, subsample.
  candidates_generator.preprocessPointCloud(mesh_cam);

  // Generate a list of grasp candidates.
  std::vector<Grasp> candidates = candidates_generator.generateGraspCandidates(mesh_cam);
  std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
  std::vector<int> labels = handsearch.reevaluateHypotheses(mesh_cam, candidates);
  std::vector<int> good_index;
  for (int i=0; i<labels.size(); ++i) {
      bool b_good = (labels[i]==2);
      if (negative_sample) b_good = !b_good;
      if (b_good) good_index.push_back(i); // good_grasps.push_back(candidates[i]);
  }
  std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
  std::cout << "Evaluation => " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count()*1e-6 << " seconds" << std::endl;
  std::cout << "Total => " << std::chrono::duration_cast<std::chrono::microseconds>(end - total_begin).count()*1e-6 << " seconds" << std::endl;
  std::cout << "Generated " << good_index.size() << (negative_sample?" bad":" good") << " grasps." << std::endl;

  if (plot_grasps) {
      std::vector<Grasp> good_grasps;
      std::random_shuffle(good_index.begin(), good_index.end());

      for (int i=0; i<std::min(max_samples, (int)good_index.size()); ++i) {
          good_grasps.push_back(candidates[good_index[i]]);
      }

      std::cout << "Selected " << good_grasps.size() << " grasps." << std::endl;
          Plot plotter;
          plotter.plotFingers3D(good_grasps, mesh_cam.getCloudOriginal(), "Selected Grasps", hand_search_params.hand_outer_diameter_,
                  hand_search_params.finger_width_, hand_search_params.hand_depth_, hand_search_params.hand_height_);
  }
  // Eigen::Vector3d euler_mean, euler2_mean;
  std::vector<int> roll_h(32,0), pitch_h(32,0), yaw_h(32,0);
  std::ofstream output_fs(argv[4]);
  // euler_mean.setZero();
  // euler2_mean.setZero();
  /*
  for (int i=0; i<good_grasps.size(); ++i) {
      Grasp single_grasp = good_grasps[i];
  */
  for (auto idx : good_index) {
      Grasp single_grasp = candidates[idx];
      std::vector<Grasp> grasp_vec, grasp_vec_centered;
      Eigen::Vector3d trans = single_grasp.getGraspBottom().cast<double>();
      Eigen::Vector3d top_vec = single_grasp.getGraspTop();
      Eigen::Matrix3d rot   = single_grasp.getFrame().cast<double>();
      double top = single_grasp.getTop();
      double bottom = single_grasp.getBottom();
      grasp_vec.push_back(single_grasp);
      /*
      if (plot_grasps) {
          Plot plotter;
          plotter.plotFingers3D(grasp_vec, mesh_cam.getCloudOriginal(), "Good Grasps (mesh)", hand_search_params.hand_outer_diameter_,
                  hand_search_params.finger_width_, hand_search_params.hand_depth_, hand_search_params.hand_height_);
      }
      if (plot_grasps) {
          Plot plotter;
          plotter.plotFingers3D(grasp_vec, cloud_cam.getCloudOriginal(), "Good Grasps (cloud)", hand_search_params.hand_outer_diameter_,
                  hand_search_params.finger_width_, hand_search_params.hand_depth_, hand_search_params.hand_height_);
      }
      */
      Eigen::Matrix4d Trans_m, InvTrans_m;
      Eigen::Affine3d Trans, InvTrans;
      Eigen::Vector3d euler;
      if (rot(2,2)<0) // up-side-down
      {
        rot = rot * Eigen::AngleAxisd(PI, Eigen::Vector3d::UnitX());
      }
      euler = rot.eulerAngles(0, 1, 2); // roll, pitch, yaw
      Trans_m.setIdentity();
      Trans_m.block<3,3>(0,0) = rot;
      Trans_m.block<3,1>(0,3) = trans;
      Trans.matrix() = Trans_m;
      InvTrans = Trans.inverse();
      InvTrans_m = InvTrans.matrix();
      double roll = euler(0);
      double pitch = euler(1);
      double yaw = euler(2);

      // std::cout << "Euler Angles (x-roll, y-pitch, z-yaw): " << euler.transpose() << std::endl;
      // std::cout << "Translation: " << trans.transpose() << std::endl;
      // euler_mean += euler;
      // euler2_mean += (Eigen::Vector3d)(euler.array()*euler.array());

      if (roll<0) roll = std::min(PI_2, std::max(0.0, roll+PI_2));
      if (pitch<0) pitch = std::min(PI_2, std::max(0.0, pitch+PI_2));
      if (yaw<0) yaw = std::min(PI_2, std::max(0.0, yaw+PI_2));
      ++roll_h[std::min((int)roll_h.size()-1,(int)(roll/PI_2*roll_h.size()))];
      ++pitch_h[std::min((int)pitch_h.size()-1,(int)(pitch/PI_2*pitch_h.size()))];
      ++yaw_h[std::min((int)yaw_h.size()-1,(int)(yaw/PI_2*yaw_h.size()))];

      output_fs << std::setprecision(8) \
                << Trans_m.row(0)         << ' ' \
                << Trans_m.row(1)         << ' ' \
                << Trans_m.row(2)         << ' ' \
                << InvTrans_m.row(0)      << ' ' \
                << InvTrans_m.row(1)      << ' ' \
                << InvTrans_m.row(2)      << std::endl;

      /*
      std::cout << Trans.matrix()    << std::endl \
                << InvTrans.matrix() << std::endl;

      std::cout << std::endl;

      std::cout << top     << std::endl \
                << top_vec << std::endl \
                << std::endl \
                << bottom  << std::endl \
                << trans   << std::endl;

      std::cout << std::endl;

      std::vector<int> enclosed_indecies;
      FingerHand finger_hand(finger_width, hand_outer_diameter, hand_depth);
      Grasp single_grasp_I = transformGrasp(single_grasp, finger_hand, InvTrans);
      grasp_vec_centered.push_back(single_grasp_I);

      pcl::CropBox<pcl::PointXYZRGBA> boxFilter(true);
      boxFilter.setMin(Eigen::Vector4f(0, -(hand_outer_diameter-finger_width*2)/2.0, -hand_height/2.0, 1.0));
      boxFilter.setMax(Eigen::Vector4f(top-bottom, (hand_outer_diameter-finger_width*2)/2.0, hand_height/2.0, 1.0));
      boxFilter.setTransform(InvTrans.cast<float>());

      PointCloudRGB::Ptr new_pc(new PointCloudRGB), new_pc_crop(new PointCloudRGB);
      boxFilter.setInputCloud(cloud_cam.getCloudOriginal());
      boxFilter.filter(*new_pc_crop);
      boxFilter.filter(enclosed_indecies);
      // for (auto v : enclosed_indecies) std::cout << v << ' ';
      // std::cout << std::endl;
      pcl::transformPointCloud(*new_pc_crop, *new_pc_crop, InvTrans);
      pcl::transformPointCloud(*cloud_cam.getCloudOriginal(), *new_pc, InvTrans);
      if (plot_grasps) {
          Plot plotter;
          plotter.plotFingers3D(grasp_vec_centered, new_pc, "Good Grasps (centered)", hand_search_params.hand_outer_diameter_,
                  hand_search_params.finger_width_, hand_search_params.hand_depth_, hand_search_params.hand_height_);
      }

      if (plot_grasps) {
          Plot plotter;
          plotter.plotFingers3D(grasp_vec_centered, new_pc_crop, "Good Grasps (centered-cropped)", hand_search_params.hand_outer_diameter_,
                  hand_search_params.finger_width_, hand_search_params.hand_depth_, hand_search_params.hand_height_);
      }
      new_pc->clear();
      new_pc_crop->clear();
      */
  }

  // euler_mean /= (double)good_grasps.size();
  // euler2_mean /= (double)good_grasps.size();

  // std::cout << "Euler Means: " << euler_mean.transpose() << std::endl;
  // std::cout << "Euler Vars : " << (euler2_mean-(Eigen::Vector3d)(euler_mean.array()*euler_mean.array())).transpose() << std::endl;

  std::cout << "Roll histogram (0~2*pi):" << std::endl;
  for (auto v : roll_h) std::cout << v << ' ';
  std::cout << std::endl;

  std::cout << "Pitch histogram (0~2*pi):" << std::endl;
  for (auto v : pitch_h) std::cout << v << ' ';
  std::cout << std::endl;

  std::cout << "Yaw histogram (0~2*pi):" << std::endl;
  for (auto v : yaw_h) std::cout << v << ' ';
  std::cout << std::endl;

  output_fs.close();

  return 0;
}