text:
//---------------------------------------------------------------------------//
//!
//! \file   tstHistogramPartiallyTabularBivariateDistribution.cpp
//! \author Alex Robinson
//! \brief  The histogram partially tabular two-dimensional dist. unit tests
//!
//---------------------------------------------------------------------------//

// Std Lib Includes
#include <iostream>
#include <memory>

// Boost Includes
#include <boost/units/systems/cgs.hpp>
#include <boost/units/io.hpp>

// FRENSIE Includes
#include "Utility_HistogramPartiallyTabularBasicBivariateDistribution.hpp"
#include "Utility_DeltaDistribution.hpp"
#include "Utility_UniformDistribution.hpp"
#include "Utility_ExponentialDistribution.hpp"
#include "Utility_ElectronVoltUnit.hpp"
#include "Utility_BarnUnit.hpp"
#include "Utility_UnitTestHarnessWithMain.hpp"
#include "ArchiveTestHelpers.hpp"

//---------------------------------------------------------------------------//
// Testing Types
//---------------------------------------------------------------------------//

using boost::units::quantity;
using Utility::Units::MegaElectronVolt;
using Utility::Units::MeV;
using Utility::Units::Barn;
using Utility::Units::barn;
using Utility::Units::barns;
namespace cgs = boost::units::cgs;

typedef TestArchiveHelper::TestArchives TestArchives;

//---------------------------------------------------------------------------//
// Testing Variables
//---------------------------------------------------------------------------//
std::shared_ptr<Utility::UnitAwarePartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> >
  unit_aware_distribution;

std::shared_ptr<Utility::PartiallyTabularBasicBivariateDistribution> distribution;

//---------------------------------------------------------------------------//
// Tests.
//---------------------------------------------------------------------------//
// Check that the upper bound of the distribution primary independent variable
// can be returned
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   getUpperBoundOfPrimaryIndepVar )
{
  FRENSIE_CHECK_EQUAL( distribution->getUpperBoundOfPrimaryIndepVar(), 2.0 );
}

//---------------------------------------------------------------------------//
// Check that the upper bound of the unit-aware distribution primary
// independent variable can be returned
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   getUpperBoundOfPrimaryIndepVar )
{
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->getUpperBoundOfPrimaryIndepVar(), 2.0*MeV );
}

//---------------------------------------------------------------------------//
// Check that the lower bound of the distribution primary independent variable
// can be returned
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   getLowerBoundOfPrimaryIndepVar )
{
  FRENSIE_CHECK_EQUAL( distribution->getLowerBoundOfPrimaryIndepVar(), 0.0 );
}

//---------------------------------------------------------------------------//
// Check that the lower bound of the unit-aware distribution primary
// independent variable can be returned
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   getLowerBoundOfPrimaryIndepVar )
{
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->getLowerBoundOfPrimaryIndepVar(), 0.0*MeV );
}

//---------------------------------------------------------------------------//
// Check that the upper bound of the conditional distribution independent
// variable can be returned
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   getUpperBoundOfSecondaryConditionalIndepVar )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL(
                 distribution->getUpperBoundOfSecondaryConditionalIndepVar(-1.0), 0.0 );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( distribution->getLowerBoundOfSecondaryConditionalIndepVar(-1.0),
                       0.0 );

  distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL( distribution->getUpperBoundOfSecondaryConditionalIndepVar(0.0),
                       std::numeric_limits<double>::infinity() );

  // In the second bin
  FRENSIE_CHECK_EQUAL( distribution->getUpperBoundOfSecondaryConditionalIndepVar(0.5),
                       std::numeric_limits<double>::infinity() );

  // On the third bin
  FRENSIE_CHECK_EQUAL(
                 distribution->getUpperBoundOfSecondaryConditionalIndepVar(1.0), 10.0 );

  // In the third bin
  FRENSIE_CHECK_EQUAL(
                 distribution->getUpperBoundOfSecondaryConditionalIndepVar(1.5), 10.0 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL(
                 distribution->getUpperBoundOfSecondaryConditionalIndepVar(2.0), 10.0 );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL(
                  distribution->getUpperBoundOfSecondaryConditionalIndepVar(3.0), 0.0 );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
                  distribution->getUpperBoundOfSecondaryConditionalIndepVar(3.0), 0.0 );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the upper bound of the unit-aware conditional distribution
// independent variable can be returned
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   getUpperBoundOfSecondaryConditionalIndepVar )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(-1.0*MeV),
         0.0*cgs::centimeter );

  // Before the first bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(-1.0*MeV),
         0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(0.0*MeV),
         std::numeric_limits<double>::infinity()*cgs::centimeter );

  // In the second bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(0.5*MeV),
         std::numeric_limits<double>::infinity()*cgs::centimeter );

  // On the third bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(1.0*MeV),
         10.0*cgs::centimeter );

  // In the third bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(1.5*MeV),
         10.0*cgs::centimeter );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(2.0*MeV),
         10.0*cgs::centimeter );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(3.0*MeV),
         0.0*cgs::centimeter );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getUpperBoundOfSecondaryConditionalIndepVar(3.0*MeV),
         0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the lower bound of the conditional distribution independent
// variable can be returned
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   getLowerBoundOfSecondaryConditionalIndepVar )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL(
                 distribution->getLowerBoundOfSecondaryConditionalIndepVar(-1.0), 0.0 );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
                 distribution->getLowerBoundOfSecondaryConditionalIndepVar(-1.0), 0.0 );

  distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(0.0), 0.0 );

  // In the second bin
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(0.5), 0.0 );

  // On the third bin
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(1.0), 0.0 );

  // In the third bin
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(1.5), 0.0 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(2.0), 0.0 );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(3.0), 0.0 );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
                  distribution->getLowerBoundOfSecondaryConditionalIndepVar(3.0), 0.0 );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the lower bound of the unit-aware conditional distribution
// independent variable can be returned
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   getLowerBoundOfSecondaryConditionalIndepVar )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(-1.0*MeV),
         0.0*cgs::centimeter );

  // Before the first bin - with extensino
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(-1.0*MeV),
         0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(0.0*MeV),
         0.0*cgs::centimeter );

  // In the second bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(0.5*MeV),
         0.0*cgs::centimeter );

  // On the third bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(1.0*MeV),
         0.0*cgs::centimeter );

  // In the third bin
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(1.5*MeV),
         0.0*cgs::centimeter );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(2.0*MeV),
         0.0*cgs::centimeter );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(3.0*MeV),
         0.0*cgs::centimeter );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL(
         unit_aware_distribution->getLowerBoundOfSecondaryConditionalIndepVar(2.0*MeV),
         0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check if the distribution is tabular in the primary dimension
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   isPrimaryDimensionTabular )
{
  FRENSIE_CHECK( distribution->isPrimaryDimensionTabular() );
}

//---------------------------------------------------------------------------//
// Check if the unit-aware distribution is tabular in the primary dimension
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   isPrimaryDimensionTabular )
{
  FRENSIE_CHECK( unit_aware_distribution->isPrimaryDimensionTabular() );
}

//---------------------------------------------------------------------------//
// Check if the distribution is continuous in the primary dimension
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   isPrimaryDimensionContinuous )
{
  FRENSIE_CHECK( distribution->isPrimaryDimensionContinuous() );
}

//---------------------------------------------------------------------------//
// Check if the unit-aware distribution is continuous in the primary dimension
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   isPrimaryDimensionContinuous )
{
  FRENSIE_CHECK( unit_aware_distribution->isPrimaryDimensionContinuous() );
}

//---------------------------------------------------------------------------//
// Check if the distribution has the same bounds as another distribution
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   hasSamePrimaryBounds )
{
  std::unique_ptr<Utility::BasicBivariateDistribution> other_distribution;

  {
    std::vector<double> primary_grid( {0.0, 2.0} );

    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::DeltaDistribution( 0.0 ) );
    secondary_dists[1].reset( new Utility::DeltaDistribution( 1.0 ) );

    other_distribution.reset(
             new Utility::HistogramPartiallyTabularBasicBivariateDistribution(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<double> primary_grid( {-1.0, 2.0} );

    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::DeltaDistribution( 0.0 ) );
    secondary_dists[1].reset( new Utility::DeltaDistribution( 1.0 ) );

    other_distribution.reset(
              new Utility::HistogramPartiallyTabularBasicBivariateDistribution(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<double> primary_grid( {0.0, 3.0} );

    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::DeltaDistribution( 0.0 ) );
    secondary_dists[1].reset( new Utility::DeltaDistribution( 1.0 ) );

    other_distribution.reset(
             new Utility::HistogramPartiallyTabularBasicBivariateDistribution(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<double> primary_grid( {-1.0, 3.0} );

    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::DeltaDistribution( 0.0 ) );
    secondary_dists[1].reset( new Utility::DeltaDistribution( 1.0 ) );

    other_distribution.reset(
              new Utility::HistogramPartiallyTabularBasicBivariateDistribution(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !distribution->hasSamePrimaryBounds( *other_distribution ) );
}

//---------------------------------------------------------------------------//
// Check if the unit-aware distribution has the same bounds as another
// distribution
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   hasSamePrimaryBounds )
{
  std::unique_ptr<Utility::UnitAwareBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> > other_distribution;

  {
    std::vector<quantity<MegaElectronVolt> >
      primary_grid( {0.0*MeV, 2.0*MeV} );

    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > > secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 0.0*cgs::centimeter ) );
    secondary_dists[1].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 1.0*cgs::centimeter ) );

    other_distribution.reset(
          new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( unit_aware_distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<quantity<MegaElectronVolt> >
      primary_grid( {-1.0*MeV, 2.0*MeV} );

    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > > secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 0.0*cgs::centimeter ) );
    secondary_dists[1].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 1.0*cgs::centimeter ) );

    other_distribution.reset(
          new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !unit_aware_distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<quantity<MegaElectronVolt> >
      primary_grid( {0.0*MeV, 3.0*MeV} );

    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > > secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 0.0*cgs::centimeter ) );
    secondary_dists[1].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 1.0*cgs::centimeter ) );

    other_distribution.reset(
          new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !unit_aware_distribution->hasSamePrimaryBounds( *other_distribution ) );

  {
    std::vector<quantity<MegaElectronVolt> >
      primary_grid( {-1.0*MeV, 3.0*MeV} );

    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > > secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 0.0*cgs::centimeter ) );
    secondary_dists[1].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 1.0*cgs::centimeter ) );

    other_distribution.reset(
          new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>(
                                             primary_grid, secondary_dists ) );
  }

  FRENSIE_CHECK( !unit_aware_distribution->hasSamePrimaryBounds( *other_distribution ) );
}

//---------------------------------------------------------------------------//
// Check that the distribution can be evaluated
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution, evaluate )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, 0.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, 1.0 ), 0.0 );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( -1.0, 1.0 ), 0.0 );

  distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.0, 1.0 ), exp( -1.0 ) );

  // In the second bin
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.5, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.5, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 0.5, 1.0 ), exp( -1.0 ) );

  // On the third bin
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.0, 1.0 ), 1.0 );

  // In the third bin
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.5, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.5, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 1.5, 1.0 ), 1.0 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 2.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 2.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 2.0, 1.0 ), 1.0 );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, 0.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, 1.0 ), 0.0 );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluate( 3.0, 1.0 ), 0.0 );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the unit-aware distribution can be evaluated
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   evaluate )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, 0.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, 1.0*cgs::centimeter ), 0.0*barn );

  // Before the first bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( -1.0*MeV, 1.0*cgs::centimeter ), 0.0*barn );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.0*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.0*MeV, 1.0*cgs::centimeter ), exp( -1.0 )*barns );

  // In the second bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.5*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.5*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 0.5*MeV, 1.0*cgs::centimeter ), exp( -1.0 )*barns );

  // On the third bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.0*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.0*MeV, 1.0*cgs::centimeter ), 1.0*barn );

  // In the third bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.5*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.5*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 1.5*MeV, 1.0*cgs::centimeter ), 1.0*barn );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 2.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 2.0*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 2.0*MeV, 1.0*cgs::centimeter ), 1.0*barn );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, 0.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, 1.0*cgs::centimeter ), 0.0*barn );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, -1.0*cgs::centimeter ), 0.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, 0.0*cgs::centimeter ), 1.0*barn );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluate( 3.0*MeV, 1.0*cgs::centimeter ), 0.0*barn );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the secondary conditional PDF can be evaluated
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   evaluateSecondaryConditionalPDF )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, 0.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, 1.0 ), 0.0 );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( -1.0, 1.0 ), 0.0 );

  distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.0, 1.0 ), exp( -1.0 ) );

  // In the second bin
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.5, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.5, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 0.5, 1.0 ), exp( -1.0 ) );

  // On the third bin
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.0, 0.0 ), 0.1 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.0, 1.0 ), 0.1 );

  // In the third bin
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.5, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.5, 0.0 ), 0.1 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 1.5, 1.0 ), 0.1 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 2.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 2.0, 0.0 ), 0.1 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 2.0, 1.0 ), 0.1 );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, 0.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, 1.0 ), 0.0 );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, -1.0 ), 0.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, 0.0 ), 1.0 );
  FRENSIE_CHECK_EQUAL( distribution->evaluateSecondaryConditionalPDF( 3.0, 1.0 ), 0.0 );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the unit-aware secondary conditional PDF can be evaluated
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   evaluateSecondaryConditionalPDF )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, 0.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, 1.0*cgs::centimeter ), 0.0/cgs::centimeter );

  // Before the first bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, 0.0*cgs::centimeter ), 1.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( -1.0*MeV, 1.0*cgs::centimeter ), 0.0/cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.0*MeV, 0.0*cgs::centimeter ), 1.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.0*MeV, 1.0*cgs::centimeter ), exp( -1.0 )/cgs::centimeter );

  // In the second bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.5*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.5*MeV, 0.0*cgs::centimeter ), 1.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 0.5*MeV, 1.0*cgs::centimeter ), exp( -1.0 )/cgs::centimeter );

  // On the third bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.0*MeV, 0.0*cgs::centimeter ), 0.1/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.0*MeV, 1.0*cgs::centimeter ), 0.1/cgs::centimeter );

  // In the third bin
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.5*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.5*MeV, 0.0*cgs::centimeter ), 0.1/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 1.5*MeV, 1.0*cgs::centimeter ), 0.1/cgs::centimeter );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 2.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 2.0*MeV, 0.0*cgs::centimeter ), 0.1/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 2.0*MeV, 1.0*cgs::centimeter ), 0.1/cgs::centimeter );

  // After the third bin - no extension
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, 0.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, 1.0*cgs::centimeter ), 0.0/cgs::centimeter );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, -1.0*cgs::centimeter ), 0.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, 0.0*cgs::centimeter ), 1.0/cgs::centimeter );
  FRENSIE_CHECK_EQUAL( unit_aware_distribution->evaluateSecondaryConditionalPDF( 3.0*MeV, 1.0*cgs::centimeter ), 0.0/cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that a secondary conditional PDF can be sampled
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   sampleSecondaryConditional )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_THROW( distribution->sampleSecondaryConditional( -1.0 ),
              std::logic_error );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  double sample = distribution->sampleSecondaryConditional( -1.0 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  std::vector<double> fake_stream( 2 );
  fake_stream[0] = 0.0;
  fake_stream[1] = 0.5;

  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditional( 0.0 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  sample = distribution->sampleSecondaryConditional( 0.0 );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453, 1e-9 );

  // In the second bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditional( 0.5 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  sample = distribution->sampleSecondaryConditional( 0.5 );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453, 1e-9 );

  // On the third bin
  fake_stream.push_back( 1.0 - 1e-15 );
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditional( 1.0 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  sample = distribution->sampleSecondaryConditional( 1.0 );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );

  sample = distribution->sampleSecondaryConditional( 1.0 );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );

  // In the third bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditional( 1.5 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  sample = distribution->sampleSecondaryConditional( 1.5 );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );

  sample = distribution->sampleSecondaryConditional( 1.5 );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditional( 2.0 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  sample = distribution->sampleSecondaryConditional( 2.0 );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );

  sample = distribution->sampleSecondaryConditional( 2.0 );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );

  // After the third bin - no extension
  Utility::RandomNumberGenerator::unsetFakeStream();

  FRENSIE_CHECK_THROW( distribution->sampleSecondaryConditional( 3.0 ),
              std::logic_error );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  sample = distribution->sampleSecondaryConditional( 3.0 );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that a unit-aware secondary conditional PDF can be sampled
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   sampleSecondaryConditional )
{
  // Before the first bin - no extension
  FRENSIE_CHECK_THROW( unit_aware_distribution->sampleSecondaryConditional( -1.0*MeV ),
              std::logic_error );

  // Before the first bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  quantity<cgs::length> sample = unit_aware_distribution->sampleSecondaryConditional( -1.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  // On the second bin (first bin boundary = second bin boundary)
  std::vector<double> fake_stream( 2 );
  fake_stream[0] = 0.0;
  fake_stream[1] = 0.5;

  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditional( 0.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 0.0*MeV );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453*cgs::centimeter, 1e-9 );

  // In the second bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditional( 0.5*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 0.5*MeV );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453*cgs::centimeter, 1e-9 );

  // On the third bin
  fake_stream.push_back( 1.0 - 1e-15 );
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.0*MeV );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );

  // In the third bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.5*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.5*MeV );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 1.5*MeV );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditional( 2.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 2.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );

  sample = unit_aware_distribution->sampleSecondaryConditional( 2.0*MeV );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );

  // After the third bin - no extension
  Utility::RandomNumberGenerator::unsetFakeStream();

  FRENSIE_CHECK_THROW( unit_aware_distribution->sampleSecondaryConditional( 3.0*MeV ),
              std::logic_error );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  sample = unit_aware_distribution->sampleSecondaryConditional( 3.0*MeV );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that a secondary conditional PDF can be sampled
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   sampleSecondaryConditionalAndRecordTrials )
{
  Utility::DistributionTraits::Counter trials = 0u;

  // Before the first bin - no extension
  FRENSIE_CHECK_THROW( distribution->sampleSecondaryConditionalAndRecordTrials( -1.0, trials ),
              std::logic_error );
  FRENSIE_CHECK_EQUAL( trials, 0u );

  // Before the first bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  double sample = distribution->sampleSecondaryConditionalAndRecordTrials( -1.0 , trials );

  distribution->limitToPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 1u );

  // On the second bin (first bin boundary = second bin boundary)
  std::vector<double> fake_stream( 2 );
  fake_stream[0] = 0.0;
  fake_stream[1] = 0.5;

  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 0.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 2u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 0.0, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 3u );

  // In the second bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 0.5, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 4u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 0.5, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 5u );

  // On the third bin
  fake_stream.push_back( 1.0 - 1e-15 );
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 6u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );
  FRENSIE_CHECK_EQUAL( trials, 7u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.0, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 8u );

  // In the third bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.5, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 9u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.5, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );
  FRENSIE_CHECK_EQUAL( trials, 10u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 1.5, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 11u );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 2.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 12u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 2.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0 );
  FRENSIE_CHECK_EQUAL( trials, 13u );

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 2.0, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 14u );

  // After the third bin - no extension
  Utility::RandomNumberGenerator::unsetFakeStream();

  FRENSIE_CHECK_THROW( distribution->sampleSecondaryConditionalAndRecordTrials( 3.0, trials ),
              std::logic_error );

  // After the third bin - with extension
  distribution->extendBeyondPrimaryIndepLimits();

  sample = distribution->sampleSecondaryConditionalAndRecordTrials( 3.0, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0 );
  FRENSIE_CHECK_EQUAL( trials, 15u );

  distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that a unit-aware secondary conditional PDF can be sampled
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   sampleSecondaryConditionalAndRecordTrials )
{
  Utility::DistributionTraits::Counter trials = 0u;

  // Before the first bin - no extension
  FRENSIE_CHECK_THROW( unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( -1.0*MeV, trials ),
              std::logic_error );
  FRENSIE_CHECK_EQUAL( trials, 0u );

  // Before the first bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  quantity<cgs::length> sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( -1.0*MeV, trials );

  unit_aware_distribution->limitToPrimaryIndepLimits();

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 1u );

  // On the second bin (first bin boundary = second bin boundary)
  std::vector<double> fake_stream( 2 );
  fake_stream[0] = 0.0;
  fake_stream[1] = 0.5;

  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 0.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 2u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 0.0*MeV, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453*cgs::centimeter, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 3u );

  // In the second bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 0.5*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 4u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 0.5*MeV, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 0.6931471805599453*cgs::centimeter, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 5u );

  // On the third bin
  fake_stream.push_back( 1.0 - 1e-15 );
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 6u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 7u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.0*MeV, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 8u );

  // In the third bin
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.5*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 9u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.5*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 10u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 1.5*MeV, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 11u );

  // On the upper bin boundary - this should be treated as part of the third
  // bin (special case)
  Utility::RandomNumberGenerator::setFakeStream( fake_stream );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 2.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 12u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 2.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 5.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 13u );

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 2.0*MeV, trials );

  FRENSIE_CHECK_FLOATING_EQUALITY( sample, 10.0*cgs::centimeter, 1e-9 );
  FRENSIE_CHECK_EQUAL( trials, 14u );

  // After the third bin - no extension
  Utility::RandomNumberGenerator::unsetFakeStream();

  FRENSIE_CHECK_THROW( unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 3.0*MeV, trials ),
              std::logic_error );

  // After the third bin - with extension
  unit_aware_distribution->extendBeyondPrimaryIndepLimits();

  sample = unit_aware_distribution->sampleSecondaryConditionalAndRecordTrials( 3.0*MeV, trials );

  FRENSIE_CHECK_EQUAL( sample, 0.0*cgs::centimeter );
  FRENSIE_CHECK_EQUAL( trials, 15u );

  unit_aware_distribution->limitToPrimaryIndepLimits();
}

//---------------------------------------------------------------------------//
// Check that the distribution can be placed in a stream
FRENSIE_UNIT_TEST( HistogramPartiallyTabularBasicBivariateDistribution,
                   ostream_operator )
{
  std::ostringstream oss;

  oss << *distribution;

  Utility::VariantMap dist_data =
    Utility::fromString<Utility::VariantMap>( oss.str() );

  FRENSIE_CHECK_EQUAL( dist_data["type"].toString(),
                       "HistogramPartiallyTabularBasicBivariateDistribution",
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["primary independent unit"].toString(),
                       "void",
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["secondary independent unit"].toString(),
                       "void",
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["dependent unit"].toString(),
                       "void",
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["primary grid"].toType<std::vector<double> >(),
                       std::vector<double>( {0.0, 0.0, 1.0, 2.0} ),
                       SHOW_LHS );

  Utility::VariantVector secondary_dists =
    dist_data["secondary dists"].toVector();

  FRENSIE_REQUIRE_EQUAL( secondary_dists.size(), 4, SHOW_LHS );

  Utility::VariantMap secondary_dist_data = secondary_dists[0].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Delta Distribution",
                       SHOW_LHS );

  secondary_dist_data = secondary_dists[1].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Exponential Distribution",
                       SHOW_LHS );

  secondary_dist_data = secondary_dists[2].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Uniform Distribution",
                       SHOW_LHS );

  FRENSIE_CHECK_EQUAL( secondary_dists[0].toMap(),
                       secondary_dists[3].toMap(),
                       SHOW_BOTH );
}

//---------------------------------------------------------------------------//
// Check that the unit-aware distribution can be placed in a stream
FRENSIE_UNIT_TEST( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                   ostream_operator )
{
  std::ostringstream oss;

  oss << *unit_aware_distribution;

  Utility::VariantMap dist_data =
    Utility::fromString<Utility::VariantMap>( oss.str() );

  FRENSIE_CHECK_EQUAL( dist_data["type"].toString(),
                       "HistogramPartiallyTabularBasicBivariateDistribution",
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["primary independent unit"].toString(),
                       Utility::UnitTraits<MegaElectronVolt>::name(),
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["secondary independent unit"].toString(),
                       Utility::UnitTraits<cgs::length>::name(),
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["dependent unit"].toString(),
                       Utility::UnitTraits<Barn>::name(),
                       SHOW_LHS );
  FRENSIE_CHECK_EQUAL( dist_data["primary grid"].toType<std::vector<quantity<MegaElectronVolt> > >(),
                       std::vector<quantity<MegaElectronVolt> >( {0.0*MeV, 0.0*MeV, 1.0*MeV, 2.0*MeV} ),
                       SHOW_LHS );

  Utility::VariantVector secondary_dists =
    dist_data["secondary dists"].toVector();

  FRENSIE_REQUIRE_EQUAL( secondary_dists.size(), 4, SHOW_LHS );

  Utility::VariantMap secondary_dist_data = secondary_dists[0].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Delta Distribution",
                       SHOW_LHS );

  secondary_dist_data = secondary_dists[1].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Exponential Distribution",
                       SHOW_LHS );

  secondary_dist_data = secondary_dists[2].toMap();

  FRENSIE_CHECK_EQUAL( secondary_dist_data["type"].toString(),
                       "Uniform Distribution",
                       SHOW_LHS );

  FRENSIE_CHECK_EQUAL( secondary_dists[0].toMap(),
                       secondary_dists[3].toMap(),
                       SHOW_BOTH );
}

//---------------------------------------------------------------------------//
// Check that the distribution can be archived
FRENSIE_UNIT_TEST_TEMPLATE_EXPAND( HistogramPartiallyTabularBasicBivariateDistribution,
                                   archive,
                                   TestArchives )
{
  FETCH_TEMPLATE_PARAM( 0, RawOArchive );
  FETCH_TEMPLATE_PARAM( 1, RawIArchive );

  typedef typename std::remove_pointer<RawOArchive>::type OArchive;
  typedef typename std::remove_pointer<RawIArchive>::type IArchive;

  std::string archive_base_name( "test_histogram_partially_tabular_basic_bivariate_dist" );
  std::ostringstream archive_ostream;

  std::string test_dist_string;

  // Create and archive some distributions
  {
    std::unique_ptr<OArchive> oarchive;

    createOArchive( archive_base_name, archive_ostream, oarchive );

    std::vector<double> primary_grid( {-1.0, 1.0} );
    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::DeltaDistribution( -1.0 ) );
    secondary_dists[1].reset( new Utility::DeltaDistribution( 1.0 ) );

    std::shared_ptr<Utility::HistogramPartiallyTabularBasicBivariateDistribution>
      concrete_dist( new Utility::HistogramPartiallyTabularBasicBivariateDistribution( primary_grid, secondary_dists ) );

    test_dist_string = Utility::toString( *concrete_dist );

    std::shared_ptr<Utility::BasicBivariateDistribution>
      basic_dist( concrete_dist );

    std::shared_ptr<Utility::PartiallyTabularBasicBivariateDistribution>
      partially_tabular_dist( concrete_dist );

    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( concrete_dist ) );
    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( basic_dist ) );
    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( partially_tabular_dist ) );
  }

  // Copy the archive ostream to an istream
  std::istringstream archive_istream( archive_ostream.str() );

  // Load the archived distributions
  std::unique_ptr<IArchive> iarchive;

  createIArchive( archive_istream, iarchive );

  std::shared_ptr<Utility::HistogramPartiallyTabularBasicBivariateDistribution> concrete_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( concrete_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *concrete_dist ),
                       test_dist_string,
                       SHOW_LHS );

  std::shared_ptr<Utility::BasicBivariateDistribution> basic_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( basic_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *basic_dist ),
                       test_dist_string,
                       SHOW_LHS );

  std::shared_ptr<Utility::PartiallyTabularBasicBivariateDistribution> partially_tabular_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( partially_tabular_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *partially_tabular_dist ),
                       test_dist_string,
                       SHOW_LHS );
}

//---------------------------------------------------------------------------//
// Check that the distribution can be archived
FRENSIE_UNIT_TEST_TEMPLATE_EXPAND( UnitAwareHistogramPartiallyTabularBasicBivariateDistribution,
                                   archive,
                                   TestArchives )
{
  FETCH_TEMPLATE_PARAM( 0, RawOArchive );
  FETCH_TEMPLATE_PARAM( 1, RawIArchive );

  typedef typename std::remove_pointer<RawOArchive>::type OArchive;
  typedef typename std::remove_pointer<RawIArchive>::type IArchive;

  std::string archive_base_name( "test_unit_aware_histogram_partially_tabular_basic_bivariate_dist" );
  std::ostringstream archive_ostream;

  std::string test_dist_string;

  // Create and archive some distributions
  {
    std::unique_ptr<OArchive> oarchive;

    createOArchive( archive_base_name, archive_ostream, oarchive );

    std::vector<quantity<MegaElectronVolt> > primary_grid( {0.0*MeV, 1.0*MeV} );
    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > >
      secondary_dists( 2 );

    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( -1.0*cgs::centimeter ) );
    secondary_dists[1].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 1.0*cgs::centimeter ) );

    std::shared_ptr<Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> >
      concrete_dist( new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>( primary_grid, secondary_dists ) );

    test_dist_string = Utility::toString( *concrete_dist );

    std::shared_ptr<Utility::UnitAwareBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> >
      basic_dist( concrete_dist );

    std::shared_ptr<Utility::UnitAwarePartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> >
      partially_tabular_dist( concrete_dist );

    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( concrete_dist ) );
    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( basic_dist ) );
    FRENSIE_REQUIRE_NO_THROW( (*oarchive) << BOOST_SERIALIZATION_NVP( partially_tabular_dist ) );
  }

  // Copy the archive ostream to an istream
  std::istringstream archive_istream( archive_ostream.str() );

  // Load the archived distributions
  std::unique_ptr<IArchive> iarchive;

  createIArchive( archive_istream, iarchive );

  std::shared_ptr<Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> > concrete_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( concrete_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *concrete_dist ),
                       test_dist_string,
                       SHOW_LHS );

  std::shared_ptr<Utility::UnitAwareBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> > basic_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( basic_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *basic_dist ),
                       test_dist_string,
                       SHOW_LHS );

  std::shared_ptr<Utility::UnitAwarePartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn> > partially_tabular_dist;

  FRENSIE_REQUIRE_NO_THROW( (*iarchive) >> BOOST_SERIALIZATION_NVP( partially_tabular_dist ) );
  FRENSIE_CHECK_EQUAL( Utility::toString( *partially_tabular_dist ),
                       test_dist_string,
                       SHOW_LHS );
}

//---------------------------------------------------------------------------//
// Custom setup
//---------------------------------------------------------------------------//
FRENSIE_CUSTOM_UNIT_TEST_SETUP_BEGIN();

FRENSIE_CUSTOM_UNIT_TEST_INIT()
{
  // Create the two-dimensional distribution
  {
    std::vector<double> primary_grid( 4 );
    std::vector<std::shared_ptr<const Utility::UnivariateDistribution> >
      secondary_dists( 4 );

    // Create the secondary distribution in the first bin
    primary_grid[0] = 0.0;
    secondary_dists[0].reset( new Utility::DeltaDistribution( 0.0 ) );

    // Create the secondary distribution in the second bin
    primary_grid[1] = 0.0;
    secondary_dists[1].reset( new Utility::ExponentialDistribution( 1.0, 1.0 ) );

    // Create the secondary distribution in the third bin
    primary_grid[2] = 1.0;
    secondary_dists[2].reset( new Utility::UniformDistribution( 0.0, 10.0, 1.0 ) );
    // Create the secondary distribution beyond the third bin
    primary_grid[3] = 2.0;
    secondary_dists[3] = secondary_dists[0];

    distribution.reset(
              new Utility::HistogramPartiallyTabularBasicBivariateDistribution(
                                                           primary_grid,
                                                           secondary_dists ) );
  }

  // Create the unit-aware two-dimensional distribution
  {
    std::vector<quantity<MegaElectronVolt> > primary_bins( 4 );

    std::vector<std::shared_ptr<const Utility::UnitAwareUnivariateDistribution<cgs::length,Barn> > > secondary_dists( 4 );

    // Create the secondary distribution in the first bin
    primary_bins[0] = 0.0*MeV;
    secondary_dists[0].reset( new Utility::UnitAwareDeltaDistribution<cgs::length,Barn>( 0.0*cgs::centimeter ) );

    // Create the secondary distribution in the second bin
    primary_bins[1] = 0.0*MeV;
    secondary_dists[1].reset( new Utility::UnitAwareExponentialDistribution<cgs::length,Barn>( 1.0*barn, 1.0/cgs::centimeter ) );

    // Create the secondary distribution in the third bin
    primary_bins[2] = 1.0*MeV;
    secondary_dists[2].reset( new Utility::UnitAwareUniformDistribution<cgs::length,Barn>( 0.0*cgs::centimeter, 10.0*cgs::centimeter, 1.0*barn ) );

    // Create the secondary distribution beyond the third bin
    primary_bins[3] = 2.0*MeV;
    secondary_dists[3] = secondary_dists[0];

    unit_aware_distribution.reset( new Utility::UnitAwareHistogramPartiallyTabularBasicBivariateDistribution<MegaElectronVolt,cgs::length,Barn>( primary_bins, secondary_dists ) );
  }

  // Initialize the random number generator
  Utility::RandomNumberGenerator::createStreams();
}

FRENSIE_CUSTOM_UNIT_TEST_SETUP_END();

//---------------------------------------------------------------------------//
// end tstHistogramPartiallyTabularBasicBivariateDistribution.cpp
//---------------------------------------------------------------------------//


text:
#include <Python.h>

#pragma once

#include <boost/asio.hpp>
#include <fc/exception/exception.hpp>
#include <fc/io/json.hpp>
#include <fc/variant.hpp>
#include <iostream>
#include <string>
#include <vector>

//#include <eosio/account_history_api_plugin/account_history_api_plugin.hpp>
#include <eosio/chain/config.hpp>
#include <eosio/chain_plugin/chain_plugin.hpp>

#include <boost/algorithm/string/split.hpp>
#include <boost/range/adaptor/transformed.hpp>
#include <boost/range/algorithm/sort.hpp>
#include <eosio/utilities/key_conversion.hpp>

#include <IR/Module.h>
#include <IR/Validate.h>
#include <Inline/BasicTypes.h>
#include <Runtime/Runtime.h>
#include <WASM/WASM.h>
#include <WAST/WAST.h>

#include <fc/io/fstream.hpp>

#include <Python.h>
#include <boost/thread.hpp>
#include <fc/log/logger_config.hpp>

#include <eosio/chain/exceptions.hpp>
#include <eosio/chain_api_plugin/chain_api_plugin.hpp>
#include <eosio/producer_plugin/producer_plugin.hpp>
#include <fc/io/json.hpp>

#include <eosio/chain/wast_to_wasm.hpp>
#include <eosio/chain_plugin/chain_plugin.hpp>
#include <eosio/wallet_plugin/wallet_manager.hpp>
#include <eosio/wallet_plugin/wallet_plugin.hpp>

#include <eosio/chain/db_api.hpp>

using namespace std;
using namespace eosio;
using namespace eosio::chain;
using namespace eosio::utilities;
using namespace eosio::wallet;
using namespace eosio::chain_apis;

bool app_isshutdown_();
uint32_t now2_();

controller& db();
string to_bitset(uint64_t val);
wallet_manager& get_wm();

#define INT_65 65

int produce_block_();
int produce_block_start_();
int produce_block_end_();

PyObject* get_info_();
PyObject* get_block_(char* num_or_id);
PyObject* get_account_(const char* name);
bool is_account_(const char* _name);

PyObject* get_currency_balance_(string& _code, string& _account, string& _symbol);

PyObject* create_account_(string creator, string newaccount, string owner,
                          string active, int sign);

PyObject* create_key_();
PyObject* get_public_key_(string& wif_key);

int get_transaction_(string& id, string& result);

int get_code_(string& name, string& wast, string& abi, string& code_hash, int& vm_type);
int get_table_(string& scope, string& code, string& table, string& result);

uint64_t string_to_uint64_(string str);
string uint64_to_string_(uint64_t n);
string convert_to_eth_address(string& name);
string convert_from_eth_address(string& eth_address);

PyObject* push_transactions_(vector<vector<chain::action>>& vv, bool sign, uint64_t skip_flag = 0, bool async = false, bool compress = false, int32_t max_ram_usage=std::numeric_limits<int32_t>::max());

void wast2wasm_( string& wast ,string& result);
void wasm2wast_(string& wasm, string& result);

bool is_replay_();

void fc_pack_bytes_(string& in, string& out);
void fc_unpack_bytes_(string& in, string& out);

void fc_pack_setabi_(string& abiPath, uint64_t account, string& out);
void fc_pack_setconfig_(string& abiPath, uint64_t account, string& out);

void fc_pack_updateauth(string& _account, string& _permission, string& _parent, string& _auth, uint32_t _delay, string& result);
void fc_pack_args(uint64_t code, uint64_t action, string& json, string& bin);
PyObject* fc_unpack_args(uint64_t code, uint64_t action, string& bin);

PyObject* gen_transaction_(vector<chain::action>& v, int expiration);
PyObject* sign_transaction_(string& trx_json_to_sign, string& str_private_key);
PyObject* push_raw_transaction_(string& signed_trx);

PyObject* get_actions_(uint64_t account, int pos, int offset);
PyObject* get_transaction_(string& id);

void n_to_symbol_(uint64_t n, string& out);
uint64_t symbol_to_n_(string& n);

bool debug_mode_();

void get_code_hash_(string& name, string& code_hash); // @suppress("Ambiguous problem")
void sha256_(string& data, string& hash); // @suppress("Ambiguous problem")

uint64_t ethaddr2n_(string& addr);
void n2ethaddr_(uint64_t n, string& addr);

uint64_t get_code_update_time_ms_(string& name);

void transaction_listen_();
void transaction_disconnect_();

void get_active_producers_(vector<string>& producers);

extern "C" {
	int has_option(const char* _option);
	int get_option(const char* option, char *result, int size);
}

void producer_pause_();
void producer_resume_();

ext: hpp

max_stars_repo_name: learnforpractice/pyeos

max_issues_repo_name: learnforpractice/pyeos

max_forks_repo_name: learnforpractice/pyeos


text:
/*
 * server.cpp
 *
 *  Created on: Oct 21, 2011
 *      Author: Vincent Zhang, ivincent.zhang@gmail.com
 * */

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

#include "config.h"

#include <signal.h>

#include <cstdlib>
#include <memory>
#include <thread>

#include <boost/filesystem.hpp>
#include <boost/program_options.hpp>

#include <muduo/net/EventLoop.h>

#include <pioneer/net/net.h>
#include <pioneer/net/net_pools.h>
#include <pioneer/net/net_handlers.h>
#include <pioneer/net/multicast.h>
#include <pioneer/net/rpc_clients.h>

#include "service/rfc_func.h"
#include "service/rfc_func.server.ipp"

namespace bf = boost::filesystem;

using muduo::net::EventLoop;
using muduo::net::InetAddress;
using pioneer::net::connection_handler;
using pioneer::net::message_handler;
using pioneer::net::inward_tag;
using pioneer::net::outward_tag;

using namespace pioneer;

std::shared_ptr<EventLoop> g_report_server_base_loop;
std::shared_ptr<EventLoop> g_inward_server_base_loop;
std::shared_ptr<EventLoop> g_outward_server_base_loop;
std::shared_ptr<net::mcast_server> g_mcast_server;

void at_signal() {
  if (system::context::system_quitting) {
    LOG(INFO) << "the system is already quitting";

    return;
  }

  system::context::system_quitting = true;

  net::inward_client_pool::ref().stop();
  net::mcast_client::ref().stop();

  if (g_mcast_server) g_mcast_server->stop();
  if (g_report_server_base_loop) g_report_server_base_loop->quit();
  if (g_inward_server_base_loop) g_inward_server_base_loop->quit();
  if (g_outward_server_base_loop) g_outward_server_base_loop->quit();
}

void signal_handler(int signal_no) {
  switch (signal_no) {
  case SIGHUP:
  case SIGTERM:
  case SIGQUIT:
  case SIGINT:
    at_signal();
    break;
  default:
    break;
  }
}

class pioneer_server {
public:

  pioneer_server(int outward_port, int inward_port, int reporter_port,
      int outward_server_threads, int inward_server_threads, int icp_threads, bool logtostderr) :
    _outward_server_address(outward_port), _inward_server_address(inward_port), _report_server_address(reporter_port),
    _outward_server_threads(outward_server_threads), _inward_server_threads(inward_server_threads), _icp_threads(icp_threads),
    _logtostderr(logtostderr)
  {
  }

  ~pioneer_server() {
    try {
      at_exit();
    }
    catch (...) {
      LOG(ERROR) << "exception at exit";
    }
  }

  pioneer_server(const pioneer_server&) = delete;
  const pioneer_server& operator=(const pioneer_server&) = delete;

public:

  void start() {
    // ****************************** local system ****************************
    init_glog();

    install_signal_handlers();

    // ****************************** report server ********************************
    start_report_server();

    // ****************************** main UDP service *****************************
    // start the mcast server so that we can receive UDP messages from the cluster
    start_mcast_server();
    // init mcast client so that we can send UDP message into the cluster
    init_mcast_client();

    // ****************************** main TCP server ******************************
    // start a TCP server in a standalone thread for TCP requests from outward the cluster
    start_outward_server();
    // start a TCP server in a standalone thread for TCP requests from inward the cluster
    start_inward_server();

    // ****************************** main TCP client service ***********************
    // init inner client pool so that we can establish connections to other inner nodes
    init_inward_client_pool();

    sleep(5); // TODO : avoid hard coding

    LOG(INFO) << "\n\n====================let's go====================\n\n";

    LOG(INFO) << "press Ctrl+c to exit";

    for (auto t : _main_threads) {
      t.second->join();
    }

    LOG(INFO) << "\n\n=========================exit====================\n\n";

    _main_threads.clear();

    if (g_outward_server_base_loop) g_outward_server_base_loop.reset();
    if (g_inward_server_base_loop) g_inward_server_base_loop.reset();
    if (g_report_server_base_loop) g_report_server_base_loop.reset();
  }

protected:

  void init_glog() {
    if (!_logtostderr) google::InitGoogleLogging("pioneer");

    LOG(INFO) << "glog has been initialized";
  }

  void install_signal_handlers() {
    ::signal(SIGHUP, &signal_handler); // terminal quit
    ::signal(SIGTERM, &signal_handler); // kill
    ::signal(SIGQUIT, &signal_handler); /* ctrl-\ */
    ::signal(SIGINT, &signal_handler); // ctrl-c
  }

  void start_report_server() {
    auto f = [this]() {
      if (g_report_server_base_loop) return;

      LOG(INFO) << "staring report server, listening at " << _report_server_address.toIpPort().c_str() << "...";

      g_report_server_base_loop.reset(new EventLoop);
      net::report_server server(g_report_server_base_loop.get(), _report_server_address, "report server");

      server.setHttpCallback(boost::bind(message_handler::on_report_server_message, _1, _2));

      server.start();
      g_report_server_base_loop->loop();

      LOG(INFO) << "quit report server";
    };

    // run in a new thread
    _main_threads["report_server"] = std::make_shared<std::thread>(f);
  }

  void start_mcast_server() {
    auto f = [this]() {
      if (g_mcast_server) return;

      LOG(INFO) << "starting mcast server...";

      g_mcast_server.reset(new net::mcast_server(PIONEER_MULTIGROUP));

      g_mcast_server->set_message_callback(net::message_handler::on_mcast_message);
      g_mcast_server->start();

      LOG(INFO) << "quit mcast server";
    };

    // run in a new thread
    _main_threads["mcast_server"] = std::make_shared<std::thread>(f);
  }

  void init_mcast_client() {
    LOG(INFO) << "initializing mcast client...";

    net::mcast_client::ref().init(PIONEER_MULTIGROUP);
  }

  void start_outward_server() {
    auto f = [this]() {
      if (g_outward_server_base_loop) return;

      LOG(INFO) << "starting outward server, listening at " << _outward_server_address.toIpPort().c_str() << "...";

      g_outward_server_base_loop.reset(new EventLoop);
      net::outward_server server(g_outward_server_base_loop.get(), _outward_server_address, "outward server");
      server.setThreadNum(_outward_server_threads);

      server.setConnectionCallback(boost::bind(connection_handler::on_outward_server_connection, _1));
      server.setMessageCallback(boost::bind(message_handler::on_outward_server_message, _1, _2, _3));
      server.setWriteCompleteCallback(boost::bind(connection_handler::on_write_complete<outward_tag>, _1));

      server.start();
      g_outward_server_base_loop->loop();

      LOG(INFO) << "quit outward server";
    };

    // run in a new thread
    _main_threads["outward_server"] = std::make_shared<std::thread>(f);
  }

  void start_inward_server() {
    auto f = [this]() {
      if (g_inward_server_base_loop) return;

      LOG(INFO) << "starting inner server, listening at " << _inward_server_address.toIpPort().c_str() << "...";

      g_inward_server_base_loop.reset(new EventLoop);
      net::inward_server server(g_inward_server_base_loop.get(), _inward_server_address, "inward server");
      server.setThreadNum(_inward_server_threads);

      server.setConnectionCallback(boost::bind(connection_handler::on_inward_server_connection, _1));
      server.setMessageCallback(boost::bind(message_handler::on_inward_server_message, _1, _2, _3));
      server.setWriteCompleteCallback(boost::bind(connection_handler::on_write_complete<inward_tag>, _1));

      server.start();
      g_inward_server_base_loop->loop();

      LOG(INFO) << "quit inward server";
    };

    // run in a new thread
    _main_threads["inward_server"] = std::make_shared<std::thread>(f);
  }

  void init_inward_client_pool() {
    auto f = [this]() {
      LOG(INFO) << "starting inner node client pool service...";

      auto& tcp_client_pool = net::inward_client_pool::ref();

      tcp_client_pool.set_server_port(PIONEER_INWARD_SERVER_PORT); // TODO : parameterize this
      tcp_client_pool.set_thread_num(_icp_threads);

      tcp_client_pool.set_connection_callback(boost::bind(net::connection_handler::on_inward_client_connection, _1));
      tcp_client_pool.set_message_callback(boost::bind(net::message_handler::on_inward_client_message, _1, _2, _3));
      tcp_client_pool.set_write_complete_callback(boost::bind(net::connection_handler::on_write_complete<net::inward_tag>, _1));

      tcp_client_pool.init();
      tcp_client_pool.start();

      LOG(INFO)<<"quit inner client pool service";
    };

    // run in a new thread
    _main_threads["inward_client_pool"] = std::make_shared<std::thread>(f);
  }

  void at_exit() {
    LOG(INFO) << "all services are stopped, do the cleaning";
  }

private:

  InetAddress _outward_server_address; // outward server address
  InetAddress _inward_server_address; // inner server address
  InetAddress _report_server_address; // report server address

  int _outward_server_threads; // outward server thread number
  int _inward_server_threads; // inner server thread number
  int _icp_threads;  // inner client pool thread number

  bool _logtostderr;

  std::map<std::string, std::shared_ptr<std::thread>> _main_threads;
};

int main(int argc, char* argv[]) {
  namespace po = boost::program_options;

  std::string help = std::string("usage : ") + argv[0] + " [options]...";

  po::options_description desc("allowed options:");
  desc.add_options()
      ("help", "Usage : [options]...")
      ("outward_port", po::value<int>()->default_value(PIONEER_OUTWARD_SERVER_PORT), "outward server port")
      ("inward_port", po::value<int>()->default_value(PIONEER_INWARD_SERVER_PORT), "inward server port")
      ("reporter_port", po::value<int>()->default_value(PIONEER_REPORT_SERVER_PORT), "report server port")
      ("outward_server_threads", po::value<int>()->default_value(OUTWARD_SERVER_THREADS), "outward server thread number")
      ("inward_server_threads", po::value<int>()->default_value(INWARD_SERVER_THREADS), "inward server thread number")
      ("icp_threads", po::value<int>()->default_value(INWARD_CLIENT_POOL_THREADS), "inward client pool thread number")
      ("logtostderr", po::value<bool>()->default_value(true), "all logs are written to stderr instead of file")
      ;

  po::variables_map vm;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify(vm);

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

  // make it a local variable to watch the destruction
  {
    pioneer_server server(
        vm["outward_port"].as<int>(),
        vm["inward_port"].as<int>(),
        vm["reporter_port"].as<int>(),
        vm["outward_server_threads"].as<int>(),
        vm["inward_server_threads"].as<int>(),
        vm["icp_threads"].as<int>(),
        vm["logtostderr"].as<bool>());

    server.start();
  }

  LOG(INFO) << "bye!";
}

ext: cpp

max_stars_repo_name: galaxyeye/pioneer

max_issues_repo_name: galaxyeye/pioneer

max_forks_repo_name: galaxyeye/pioneer


text:
#include <pybind11/pybind11.h>

#include <boost/lexical_cast.hpp>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/Triangulation_conformer_2.h>
#include <CGAL/lloyd_optimize_mesh_2.h>

#include <type_traits>
#include <typeinfo>
#ifndef _MSC_VER
#   include <cxxabi.h>
#endif
#include <memory>
#include <string>
#include <cstdlib>

#include <functional>

template <typename...> struct hash;

template<typename T>
struct hash<T>
        : public std::hash<T>
{
    using std::hash<T>::hash;
};

template <typename T, typename... Rest>
struct hash<T, Rest...>
{
    inline std::size_t operator()(const T& v, const Rest&... rest) {
        std::size_t seed = hash<Rest...>{}(rest...);
        seed ^= hash<T>{}(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
        return seed;
    }
};

template <class T>
std::string
type_name()
{
    typedef typename std::remove_reference<T>::type TR;
    std::unique_ptr<char, void(*)(void*)> own
            (
            nullptr,
            std::free
    );
    std::string r = own != nullptr ? own.get() : typeid(TR).name();
    if (std::is_const<TR>::value)
        r += " const";
    if (std::is_volatile<TR>::value)
        r += " volatile";
    if (std::is_lvalue_reference<T>::value)
        r += "&";
    else if (std::is_rvalue_reference<T>::value)
        r += "&&";
    return r;
}

namespace py = pybind11;

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Vb = CGAL::Delaunay_mesh_vertex_base_2<K>;
using Fb = CGAL::Delaunay_mesh_face_base_2<K>;
using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
using CDT = CGAL::Constrained_Delaunay_triangulation_2<K, Tds>;
using Point = CDT::Point;
using Vertex_handle = CDT::Vertex_handle;
using Criteria = CGAL::Delaunay_mesh_size_criteria_2<CDT>;
using Mesher = CGAL::Delaunay_mesher_2<CDT, Criteria>;

template <typename T>
class TypedInputIterator
{
public:
    using iterator_category = std::input_iterator_tag;
    using difference_type = std::ptrdiff_t;
    using value_type = T;
    using pointer = T*;
    using reference = T&;

    explicit TypedInputIterator(py::iterator& py_iter) :
            py_iter_(py_iter)
    {
    }

    explicit TypedInputIterator(py::iterator&& py_iter) :
            py_iter_(py_iter)
    {
    }

    value_type operator*()
    {
        return (*py_iter_).template cast<value_type>();
    }

    TypedInputIterator operator++(int)
    {
        auto copy = *this;
        ++py_iter_;
        return copy;
    }

    TypedInputIterator& operator++()
    {
        ++py_iter_;
        return *this;
    }

    bool operator!=(TypedInputIterator &rhs)
    {
        return py_iter_ != rhs.py_iter_;
    }

    bool operator==(TypedInputIterator &rhs)
    {
        return py_iter_ == rhs.py_iter_;
    }

private:
    py::iterator py_iter_;
};

PYBIND11_MODULE(daidalus, m)
{

    py::class_<Point>(m, "Point")
        .def(py::init<int, int>(),  py::arg("x"), py::arg("y"))
        .def(py::init<double, double>(), py::arg("x"), py::arg("y"))
        .def_property_readonly("x", &Point::x)
        .def_property_readonly("y", &Point::y)
        .def("__repr__",
             [](const Point &p) {
                 std::string r("Point(");
                 r += boost::lexical_cast<std::string>(p.x());
                 r += ", ";
                 r += boost::lexical_cast<std::string>(p.y());
                 r += ")";
                 return r;
             })
        .def("__hash__",
             [](const Point &p) {
                 std::hash<double> double_hash;
                 auto x_hash = double_hash(p.x());
                 auto y_hash = double_hash(p.y());
                 return y_hash ^ x_hash + 0x9e3779b9 + (y_hash << 6) + (y_hash >> 2);
             }
        )
        .def("__eq__",
            [](const Point &p, const Point & q) {
                return p == q;
            }
        )
        ;

    py::class_<Vertex_handle>(m, "VertexHandle")
        .def_property_readonly(
            "point",
            [](const Vertex_handle& vertex_handle)
            {
                return vertex_handle->point();
            })
            ;

    py::class_<CDT::Finite_vertices_iterator::value_type>(m, "Vertex")
            .def_property_readonly(
                "point", [](CDT::Finite_vertices_iterator::value_type& vertex)
                {
                    return vertex.point();
                }
            )
            ;

    py::class_<CDT::Finite_faces_iterator::value_type>(m, "Face")
            .def("vertex_handle",
                [](CDT::Finite_faces_iterator::value_type& face, int index)
                {
                    return face.vertex(index);
                },
                py::arg("index")
            )
            ;

    m.def("print_faces_iterator_value_type", [](){
        std::cout << type_name<CDT::Finite_faces_iterator::value_type>();
    });

    py::class_<CDT>(m, "ConstrainedDelaunayTriangulation")
        .def(py::init())
        .def("insert", [](CDT & cdt, const Point & p) { return cdt.insert(p); })
        .def("insert_constraint",
             [](CDT & cdt, Vertex_handle a, Vertex_handle b)
             {
                 cdt.insert_constraint(a, b);
             })
        .def("number_of_vertices", &CDT::number_of_vertices)
        .def("number_of_faces", &CDT::number_of_faces)
        .def("finite_vertices", [](CDT & cdt) -> py::iterator
        {
            return py::make_iterator(cdt.finite_vertices_begin(), cdt.finite_vertices_end());
        })
        .def("finite_faces", [](CDT & cdt) -> py::iterator
        {
            return py::make_iterator(cdt.finite_faces_begin(), cdt.finite_faces_end());
        })
        ;

    py::class_<Criteria>(m, "Criteria")
            .def(py::init<double, double>(),
                 py::arg("aspect_bound") = 0.125,
                 py::arg("size_bound") = 0.0)
            .def_property("size_bound", &Criteria::size_bound, &Criteria::set_size_bound)
            .def_property("aspect_bound",
                          [](const Criteria & c) { c.bound(); },
                          [](Criteria & c, double bound) { c.set_bound(bound); })
            ;

    py::class_<Mesher>(m, "Mesher")
        .def(py::init<CDT&>())
        .def("seeds_from", [](Mesher & mesher, py::iterable iterable)
        {
            py::iterator iterator = py::iter(iterable);
            TypedInputIterator<Point> points_begin(iterator);
            TypedInputIterator<Point> points_end(py::iterator::sentinel());
            mesher.set_seeds(points_begin, points_end);
        })
        .def("refine_mesh", &Mesher::refine_mesh)
        .def_property(
                "criteria",
                &Mesher::get_criteria,
                [](Mesher& mesher, const Criteria & criteria)
                {
                    mesher.set_criteria(criteria);
                }
        )
        ;

    m.def("make_conforming_delaunay",
          &CGAL::make_conforming_Delaunay_2<CDT>,
          py::arg("cdt")
    );

    m.def("make_conforming_gabriel",
          &CGAL::make_conforming_Gabriel_2<CDT>,
          py::arg("cdt")
    );

    m.def("lloyd_optimize", [](
                  CDT& cdt,
                  int max_iteration_number,
                  double time_limit,
                  double convergence,
                  double freeze_bound)
          {
              CGAL::lloyd_optimize_mesh_2(cdt,
                                          CGAL::parameters::max_iteration_number = max_iteration_number,
                                          CGAL::parameters::time_limit = time_limit,
                                          CGAL::parameters::convergence = convergence,
                                          CGAL::parameters::freeze_bound = freeze_bound
              );
          },
          py::arg("cdt"),
          py::arg("max_iteration_number") = 0,
          py::arg("time_limit") = 0.0,
          py::arg("convergence") = 0.001,
          py::arg("freeze_bound") = 0.001
    );

}

max_stars_repo_name: josuehfa/DAASystem

max_issues_repo_name: josuehfa/DAASystem

max_forks_repo_name: josuehfa/DAASystem


text:
#include "relativecm_rme.h"

#include <Eigen/Dense>
#include <cmath>
#include <vector>

#include "basis_func/ho.h"
#include "chime.h"
#include "constants.h"
#include "quadpp/quadpp.h"
#include "quadpp/spline.h"
#include "tprme.h"

namespace chime {
namespace relcm {

constexpr double mPi = constants::pion_mass_fm;
constexpr double mN = constants::nucleon_mass_fm;
constexpr double FPi = constants::pion_decay_constant_fm;
constexpr double gA = constants::gA;

void ConstructMu2nNLOOperator(
    const basis::RelativeCMOperatorParametersLSJT& op_params,
    const basis::RelativeCMSpaceLSJT& relcm_space,
    std::array<basis::RelativeCMSectorsLSJT, 3>& relcm_sectors,
    std::array<basis::OperatorBlocks<double>, 3>& relcm_matrices,
    const double& oscillator_energy, const double& R)
{
  std::cout << " Constructing M1 operator...\n";
  assert(op_params.J0 == 1);
  assert(op_params.g0 == 0);
  assert((op_params.T0_min == 1) && (op_params.T0_max == 1));

  // Alias isospin rank.
  int T0 = op_params.T0_min;

  // Generate required harmonic oscillator basis functions, and radial
  // integral weights.
  int Nmax = op_params.Nmax;

  const int npts = 3001;
  const double low = 0, high = 1;
  Eigen::ArrayXd x, r, jac, wt;
  quadpp::SemiInfiniteIntegralMesh(npts, low, high, x, r, jac);
  wt = r * r * jac;  // weights for radial integral with transformed variable

  std::cout << "  Generating basis functions...\n";
  std::vector<Eigen::ArrayXXd> ho_wfs;
  double brel = chime::RelativeOscillatorLength(oscillator_energy);
  double bcm = chime::CMOscillatorLength(oscillator_energy);
  basis_func::ho::WaveFunctionsUptoMaxL(ho_wfs, r, Nmax, Nmax, brel,
                                        basis_func::Space::coordinate);

  // Relative radial integral kernels.
  std::cout << "  Generating integration kernels...\n";
  Eigen::ArrayXd mpir = mPi * r;
  Eigen::ArrayXd expmpir = Eigen::exp(-mpir);
  Eigen::ArrayXd ypir = expmpir / mpir;
  Eigen::ArrayXd zpir = (1. + mpir);
  Eigen::ArrayXd tpir = (-1. + 2 * mpir);
  Eigen::ArrayXd wpir = (1. + (3 * zpir / mpir.square()));

  // Semilocal coordinate space regulator.
  Eigen::ArrayXd scs_reg = Eigen::ArrayXd::Ones(npts);
  if (R != 0) {
    scs_reg *= Eigen::pow(1. - Eigen::exp(-(r * r) / (R * R)), 6);
  }

  // Zero initialize operator.
  std::cout << "  Zero initializing operator...\n";
  for (int T = op_params.T0_min; T <= op_params.T0_max; ++T) {
    relcm_sectors[T] = basis::RelativeCMSectorsLSJT(relcm_space, op_params.J0,
                                                    T, op_params.g0);
    basis::SetOperatorToZero(relcm_sectors[T], relcm_matrices[T]);
  }

  // Select T0 component.
  const basis::RelativeCMSectorsLSJT& sectors = relcm_sectors[T0];
  basis::OperatorBlocks<double>& matrices = relcm_matrices[T0];

  // Reduced matrix element calculation.
  std::cout << "  Starting matrix element calculation...\n";
  for (std::size_t sector_index = 0; sector_index < sectors.size();
       ++sector_index) {
    const basis::RelativeCMSectorsLSJT::SectorType& sector =
        sectors.GetSector(sector_index);
    const basis::RelativeCMSubspaceLSJT& bra_subspace = sector.bra_subspace();
    const basis::RelativeCMSubspaceLSJT& ket_subspace = sector.ket_subspace();

    // Alias for matrix.
    basis::OperatorBlock<double>& matrix = matrices[sector_index];

    // Extract subspace labels.
    int bra_L = bra_subspace.L();
    int bra_S = bra_subspace.S();
    int bra_J = bra_subspace.J();
    int bra_T = bra_subspace.T();
    int ket_L = ket_subspace.L();
    int ket_S = ket_subspace.S();
    int ket_J = ket_subspace.J();
    int ket_T = ket_subspace.T();

    if (bra_T == ket_T) {
      continue;
    }

    // Loop over bra and ket states.
    const std::size_t bra_subspace_size = bra_subspace.size();
    const std::size_t ket_subspace_size = ket_subspace.size();
#pragma omp parallel for collapse(2)
    for (std::size_t bra_index = 0; bra_index < bra_subspace_size;
         ++bra_index) {
      for (std::size_t ket_index = 0; ket_index < ket_subspace_size;
           ++ket_index) {
        const basis::RelativeCMStateLSJT bra_state(bra_subspace, bra_index);
        const basis::RelativeCMStateLSJT ket_state(ket_subspace, ket_index);

        // Extract state labels.
        int bra_nr = bra_state.Nr();
        int bra_lr = bra_state.lr();
        int bra_nc = bra_state.Nc();
        int bra_lc = bra_state.lc();
        int ket_nr = ket_state.Nr();
        int ket_lr = ket_state.lr();
        int ket_nc = ket_state.Nc();
        int ket_lc = ket_state.lc();

        // Common part of all radial integrals.
        Eigen::ArrayXd common_integrand = wt * scs_reg;
        common_integrand *=
            ho_wfs.at(bra_lr).row(bra_nr) * ho_wfs.at(ket_lr).row(ket_nr);

        // Reduced matrix element calculation.
        // Pauli matrix tensor product in spin space enforces the bra and
        // ket spins to be the same for the relative-cm part, and to be
        // different for the purely relative part.
        double rme = 0;

        if (bra_S == ket_S) {
          // Relative-cm part.
          double tp_a =
              tp::CCSpinTensorProductRME(bra_state, ket_state, 1, 1, 1, 0, 1);
          tp_a *= -std::sqrt(3.);

          Eigen::ArrayXd y = common_integrand * expmpir;
          y.head(1) = 0;  // Required to avoid divide by 0.
          y.tail(1) = 0;  // Required to avoid divide by 0.
          double integ_expmpir = quadpp::spline::Integrate(x, y);

          rme = tp_a * integ_expmpir;

          if (bra_S == 1) {
            // Rank 2 Pauli Matrix tensor product.

            double tp_b =
                tp::CCSpinTensorProductRME(bra_state, ket_state, 1, 1, 1, 2, 1);
            tp_b *= std::sqrt(3. / 5.);

            double tp_c =
                tp::CCSpinTensorProductRME(bra_state, ket_state, 1, 1, 2, 2, 1);
            tp_c *= std::sqrt(9. / 5.);

            double tp_d =
                tp::CCSpinTensorProductRME(bra_state, ket_state, 3, 1, 2, 2, 1);
            tp_d *= std::sqrt(14. / 5.);

            double tp_e =
                tp::CCSpinTensorProductRME(bra_state, ket_state, 3, 1, 3, 2, 1);
            tp_e *= std::sqrt(28. / 5.);

            y *= wpir;
            y.head(1) = 0;  // Required to avoid divide by 0.
            y.tail(1) = 0;  // Required to avoid divide by 0.
            double integ_expmpir_wpir = quadpp::spline::Integrate(x, y);

            rme += (tp_b + tp_c + tp_d + tp_e) * integ_expmpir_wpir;
          }

          double integ_cm = 0;  // CM coordinate integral; analytical result.
          if (bra_lc == ket_lc + 1) {
            integ_cm = ((std::sqrt(ket_nc + ket_lc + 1.5) * (bra_nc == ket_nc))
                        + (std::sqrt(ket_nc) * (bra_nc + 1 == ket_nc)));
          }
          else if (bra_lc + 1 == ket_lc) {
            integ_cm = ((std::sqrt(bra_nc + ket_nc + 1.5)) * (bra_nc == ket_nc)
                        + (std::sqrt(bra_nc) * (bra_nc == ket_nc + 1)));
          }
          integ_cm *= mPi * bcm;

          rme *= integ_cm;
        }
        else {
          // Purely relative part. The cm labels for the bra and ket
          // must be the same.
          if ((bra_nc == ket_nc) && (bra_lc == ket_lc)) {
            double tp_f =
                tp::CCSpinTensorProductRME(bra_state, ket_state, 2, 0, 2, 1, 1);
            tp_f *= std::sqrt(10.);

            Eigen::ArrayXd y = common_integrand * zpir * ypir;
            y.head(1) = 0;  // Required to avoid divide by 0.
            y.tail(1) = 0;  // Required to avoid divide by 0.
            double integ_zpir_ypir = quadpp::spline::Integrate(x, y);

            rme = tp_f * integ_zpir_ypir;

            if (bra_lr == ket_lr) {
              // Rank 0 spherical harmonic.
              double tp_g = tp::CCSpinTensorProductRME(bra_state, ket_state, 0,
                                                       0, 0, 1, 1);

              y = common_integrand * tpir * ypir;
              y.head(1) = 0;  // Required to avoid divide by 0.
              y.tail(1) = 0;  // Required to avoid divide by 0.
              double integ_tpir_ypir = quadpp::spline::Integrate(x, y);

              rme += tp_g * integ_tpir_ypir;
            }
          }
        }
        rme *= tp::SpinTensorProductRME(bra_T, ket_T, 1);  // Isospin.
        rme *= -(mN * mPi * gA * gA) / (24 * constants::pi * FPi * FPi);

        matrix(bra_index, ket_index) = rme;
      }
    }
  }
}

}  // namespace relcm
}  // namespace chime

max_stars_repo_name: e-eight/chime

max_issues_repo_name: e-eight/chime

max_forks_repo_name: e-eight/chime