tstGammaResourceIO.cpp

./Core/io/tests/tstGammaResourceIO.cpp

#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <istream>
#include <sstream>
#include <string>
#include <vector>
#include "Origen/Core/fn/print.h"
#include "Origen/Core/io/GammaResourceIO.h"
#include "Origen/Core/io/GammaResourceIO_cimg.h"
#include "Origen/Core/re/GammaResource.h"

#include "Nemesis/gtest/nemesis_gtest.hh"
#include "ScaleData/Core/Utils.h"
#include "ScaleUtils/IO/Utils.h"

#include "Origen/Core/config.h"

#include "Origen/Core/TestPaths.h"

#include <string>

using namespace Origen;
using namespace nemesis;

using ScaleUtils::IO::nprintf;

const std::string snippet =
    R"(# new comment character!
  40070           1.        0.        1.        0.        0.4.9862E-02     be     7 endf half=4.5982E+06
          4.7760E-01 1.0440E-01
 110250          11.        1.       10.        0.        0.4.3596E-01     na    25 endf half=5.9100E+01
          1.2540E-03 5.1562E-07 3.8971E-01 1.2675E-01 5.8503E-01 1.3000E-01
          8.3684E-01 1.0400E-03 9.7474E-01 1.4950E-01 9.8987E-01 1.6640E-03
          1.3795E+00 2.3140E-03 1.6117E+00 9.4770E-02 1.9645E+00 1.4664E-03
          2.2163E+00 9.3470E-04 2.8013E+00 4.9400E-04
 922380          17.       15.        2.        0.        0.1.2817E-03      u   238 endf half=1.4100E+17
          3.0613E-03 1.0808E-02 1.2960E-02 3.0384E-02 1.3440E-02 3.9535E-03
          1.6150E-02 3.7072E-02 1.9149E-02 8.2288E-03 9.0330E-02 6.7129E-06
          9.3795E-02 1.1074E-05 1.0528E-01 1.2914E-06 1.0607E-01 2.5224E-06
          1.0661E-01 4.3217E-08 1.0677E-01 4.9923E-08 1.0895E-01 3.2125E-07
          1.0915E-01 6.4549E-07 1.0939E-01 1.2493E-08 1.0943E-01 1.4463E-08
          4.9550E-02 6.4000E-04 1.1350E-01 1.0200E-04
 581510         138.        0.        0.      138.        0.1.0269E+00     ce   151 endf half=1.7600E+00
          2.2500E-02 1.3548E+00 6.2500E-02 3.3108E-01 1.0250E-01 1.9225E-01
          1.4250E-01 1.3092E-01 1.8250E-01 9.8551E-02 2.2250E-01 7.9834E-02
          2.6250E-01 6.9395E-02 3.0250E-01 6.5063E-02 3.4250E-01 6.5528E-02
          3.8250E-01 7.0008E-02 4.2250E-01 7.7503E-02 4.6250E-01 8.6724E-02
          5.0250E-01 9.4915E-02 5.4250E-01 9.7980E-02 5.8250E-01 9.1299E-02
          6.2250E-01 7.3918E-02 6.6250E-01 5.1067E-02 7.0250E-01 3.0740E-02
          7.4250E-01 1.7830E-02 7.8250E-01 1.1763E-02 8.2250E-01 9.5956E-03
          8.6250E-01 8.9293E-03 9.0250E-01 8.6778E-03 9.4250E-01 8.5047E-03
          9.8250E-01 8.3443E-03 1.0225E+00 8.1860E-03 1.0625E+00 8.0259E-03
          1.1025E+00 7.8665E-03 1.1425E+00 7.7099E-03 1.1825E+00 7.5550E-03
          1.2225E+00 7.3983E-03 1.2625E+00 7.2461E-03 1.3025E+00 7.0921E-03
          1.3425E+00 6.9407E-03 1.3825E+00 6.7919E-03 1.4225E+00 6.6462E-03
          1.4625E+00 6.5032E-03 1.5025E+00 6.3658E-03 1.5425E+00 6.2319E-03
          1.5825E+00 6.1025E-03 1.6225E+00 5.9774E-03 1.6625E+00 5.8535E-03
          1.7025E+00 5.7240E-03 1.7425E+00 5.5928E-03 1.7825E+00 5.4601E-03
          1.8225E+00 5.3296E-03 1.8625E+00 5.2010E-03 1.9025E+00 5.0761E-03
          1.9425E+00 4.9522E-03 1.9825E+00 4.8311E-03 2.0225E+00 4.7079E-03
          2.0625E+00 4.5800E-03 2.1025E+00 4.4512E-03 2.1425E+00 4.3273E-03
          2.1825E+00 4.2108E-03 2.2225E+00 4.1015E-03 2.2625E+00 3.9939E-03
          2.3025E+00 3.8715E-03 2.3425E+00 3.7174E-03 2.3825E+00 3.5146E-03
          2.4225E+00 3.2461E-03 2.4625E+00 2.9200E-03 2.5025E+00 2.5694E-03
          2.5425E+00 2.2613E-03 2.5825E+00 2.0431E-03 2.6225E+00 1.9106E-03
          2.6625E+00 1.8280E-03 2.7025E+00 1.7550E-03 2.7425E+00 1.6689E-03
          2.7825E+00 1.5694E-03 2.8225E+00 1.4721E-03 2.8625E+00 1.3903E-03
          2.9025E+00 1.3311E-03 2.9425E+00 1.2878E-03 2.9825E+00 1.2553E-03
          3.0225E+00 1.2358E-03 3.0625E+00 1.2339E-03 3.1025E+00 1.2310E-03
          3.1425E+00 1.1674E-03 3.1825E+00 9.8368E-04 3.2225E+00 6.9676E-04
          3.2625E+00 4.0459E-04 3.3025E+00 1.9616E-04 3.3425E+00 8.5655E-05
          3.3825E+00 3.8926E-05 3.4225E+00 2.1949E-05 3.4625E+00 1.6267E-05
          3.5025E+00 1.4034E-05 3.5425E+00 1.2536E-05 3.5825E+00 1.1314E-05
          3.6225E+00 1.0337E-05 3.6625E+00 9.2866E-06 3.7025E+00 7.7799E-06
          3.7425E+00 5.9855E-06 3.7825E+00 4.4653E-06 3.8225E+00 3.5274E-06
          3.8625E+00 3.0718E-06 3.9025E+00 2.8621E-06 3.9425E+00 2.7342E-06
          3.9825E+00 2.6289E-06 4.0225E+00 2.5313E-06 4.0625E+00 2.4366E-06
          4.1025E+00 2.3476E-06 4.1425E+00 2.2628E-06 4.1825E+00 2.1817E-06
          4.2225E+00 2.1057E-06 4.2625E+00 2.0365E-06 4.3025E+00 1.9826E-06
          4.3425E+00 1.9460E-06 4.3825E+00 1.9259E-06 4.4225E+00 1.9050E-06
          4.4625E+00 1.8630E-06 4.5025E+00 1.7765E-06 4.5425E+00 1.6420E-06
          4.5825E+00 1.5019E-06 4.6225E+00 1.4288E-06 4.6625E+00 1.4325E-06
          4.7025E+00 1.4139E-06 4.7425E+00 1.2703E-06 4.7825E+00 1.0185E-06
          4.8225E+00 7.6317E-07 4.8625E+00 5.7343E-07 4.9025E+00 4.5278E-07
          4.9425E+00 3.8173E-07 4.9825E+00 3.3967E-07 5.0225E+00 2.9425E-07
          5.0625E+00 2.2902E-07 5.1025E+00 1.5401E-07 5.1425E+00 9.1661E-08
          5.1825E+00 4.9117E-08 5.2225E+00 2.3192E-08 5.2625E+00 9.0785E-09
          5.3025E+00 2.8004E-09 5.3425E+00 6.5733E-10 5.3825E+00 1.1557E-10
          5.4225E+00 1.5076E-11 5.4625E+00 1.4490E-12 5.4862E+00 4.6768E-14
 711711          15.       14.        1.        0.        0.2.0209E-03     lu   171 endf half=7.9000E+01
          7.1715E-03 7.8516E-03 7.6402E-03 8.2277E-02 8.7704E-03 9.7394E-02
          1.0158E-02 1.5493E-02 5.3069E-02 8.2776E-04 5.4199E-02 1.6239E-03
          6.1184E-02 1.8795E-04 6.1426E-02 3.6538E-04 6.1791E-02 3.7318E-06
          6.1844E-02 4.7723E-06 6.3032E-02 4.2509E-05 6.3087E-02 8.2187E-05
          6.3233E-02 8.6910E-07 6.3243E-02 1.1052E-06 7.1100E-02 2.1100E-03
)";

void gnuplot_binning_output( std::string prefix,
                             double emax,
                             double emin,
                             const EmissionSpectrum& em )
{
    std::string dtxt( prefix + "_discrete.txt" );
    std::string ctxt( prefix + "_continuous.txt" );

    // output discrete line spectrum
    {
        std::ofstream ofs( dtxt );
        ofs << em.print_discrete();
    }
    // output continuous spectrum
    {
        std::ofstream ofs( ctxt );
        ofs << em.print_continuous();
    }

    // output driver
    {
        std::ofstream ofs( prefix + "_spectrum.gpt" );
        ofs << "set logscale xy\n";
        ofs << "set xrange [" << emin << ":" << emax << "] reverse\n";
        ofs << "plot '" << dtxt << "' u 1:2:(0):(-$2) w vectors nohead,";
        ofs << "     '" << ctxt << "' u 1:2 w l\n";
    }
}
class GammaResourceData : public ::testing::TestWithParam<int>
{
  protected:
    void SetUp()
    {
        SCOPED_TIMER( "load" );

        Vec_Str paths( {"origen.rev00.mpbrh2om.data",
                        "origen.rev00.mpbrh2op.data",
                        "origen.rev00.mpbruo2m.data",
                        "origen.rev00.mpbruo2p.data",
                        "origen.rev00.mpsfangm.data",
                        "origen.rev04.mpdkxgam.data"} );

        std::string scale_data = "/Users/ww5/lib/scale_dev_data/origen_data/";
        std::cout << nprintf( "%30s %7s %3s %3s %1s %7s %7s %7s\n",
                              "file",
                              "id",
                              "nc",
                              "nd",
                              "?",
                              "totI",
                              "totEI",
                              "meanE" );

        for( auto path : paths )
        {
            std::string full_path( scale_data + path );
            GammaResourceIO io;
            DataContainer opts;

            SP_GammaResource resource( new GammaResource() );
            bool read = io.load( *resource, full_path, opts );
            EXPECT_TRUE( read );
            resource_map[path] = resource;

            for( size_t i = 0; i < resource->emission_size(); ++i )
            {
                EmissionSpectrum em = resource->emission_at( i );
                std::cout << nprintf(
                    "%30s %7d %7s %3d %3d %1s %7.1e %7.1e %7.1e\n",
                    path.c_str(),
                    em.id(),
                    ScaleData::Utils::pizzzaaa_to_symbol( em.id() ).c_str(),
                    em.num_continuous(),
                    em.num_discrete(),
                    ( em.num_continuous() == 0 || em.num_discrete() == 0 )
                        ? "N"
                        : "Y",
                    em.total_intensity(),
                    em.total_energy_intensity(),
                    em.mean_energy() );
            }
        }
    }
    std::map<std::string, SP_GammaResource> resource_map;
};

TEST_P( GammaResourceData, DISABLED_BinningPerformance )
{
    // get an energy grid
    double emax = 1e+1;
    double emin = 1e-4;

    size_t ndiv = GetParam();
    Vec_Dbl e2 = EmissionSpectrum::create_bounds_log( emax, emin, ndiv );
    Vec_Dbl b2( ndiv + 1, 0.0 );

    // request a special print of an id
    int print_id = 27060;  // 52137;

    // weighting map
    std::map<int, double> map;
    map[37094] = 1.0;
    map[92235] = 1.0;

    // running sum
    double running_sum = 0;

    for( auto entry : resource_map )
    {
        const GammaResource& resource = *entry.second;
        const std::string& path = entry.first;

        for( size_t i = 0; i < resource.emission_size(); ++i )
        {
            EmissionSpectrum em = resource.emission_at( i );

            if( em.id() == print_id )
            {
                std::string out( path + "_" + std::to_string( em.id() ) );
                std::cerr << "printing " << out << std::endl;
                gnuplot_binning_output( out, emax, emin, em );
            }

            {
                // bin with weight in map or something small
                // if not found so is not a null op
                double id_wt = map.count( em.id() ) > 0 ? map[em.id()] : 1e-20;
                em.bin( e2, &b2, id_wt );

                // check that we are always increasing the total intensity
                double this_sum = ScaleSTL::sum( b2 );
                jDebugLine( "  intensity sum=" << this_sum );
                EXPECT_GT( this_sum, running_sum )
                    << " failed to increase intensity with id=" << em.id();
                if( this_sum <= running_sum )
                {
                    std::string out( path + "_" + std::to_string( em.id() ) );
                    jDebugLine( "   printing " << out );
                    gnuplot_binning_output( out, emax, emin, em );
                }
                running_sum = this_sum;
            }
        }
    }

    printTimingReport();

    // print out complete bins
    {
        std::ofstream ofs( "bins" + std::to_string( ndiv ) + ".txt" );
        ofs << Spectrum::print_bins( e2, b2 );
    }
    {
        std::ofstream ofs( "bins" + std::to_string( ndiv ) + ".gpt" );
        ofs << "set logscale xy\n";
        ofs << "set xrange [" << emin << ":" << emax << "] reverse\n";
        ofs << "plot 'bins" + std::to_string( ndiv ) + ".txt' u 1:2 w l\n";
    }
    return;
}

INSTANTIATE_TEST_CASE_P( Ndiv,
                         GammaResourceData,
                         ::testing::Values( 1, 1e1, 1e2, 1e3, 1e4, 1e5 ) );

TEST( GammaResourceIO, cimg )
{
    GammaResource resource;

    bool test_local_full = false;
    if( test_local_full )
    {
        std::string path(
            "/Users/ww5/lib/scale_dev_data/origen_data/"
            "origen.rev04.mpdkxgam.data" );
        GammaResourceIO io;
        DataContainer opts;
        bool read = io.load( resource, path, opts );
        EXPECT_TRUE( read );
    }
    else
    {
        GammaResourceIO_cimg io;
        DataContainer opts;
        std::istringstream stream( snippet );
        bool read = io.read( resource, stream, opts );
        EXPECT_TRUE( read );
    }

    EXPECT_TRUE( resource.has_emission( 4007 ) );
    EXPECT_TRUE( resource.has_emission( 11025 ) );
    EXPECT_TRUE( resource.has_emission( 92238 ) );
    EXPECT_TRUE( resource.has_emission( 58151 ) );
    EXPECT_TRUE( resource.has_emission( 1071171 ) );
    EXPECT_FALSE( resource.has_emission( 1001 ) );

    // Lines test
    {
        EmissionSpectrum emission = resource.emission( 92238 );
        EXPECT_EQ( 0, emission.num_continuous() );
        EXPECT_EQ( 17, emission.num_discrete() );

        DiscreteSpectrum lines = emission.discrete();
        EXPECT_DOUBLE_EQ( 1.0808E-02, lines.search( 3.0613E-03 ) );
        EXPECT_DOUBLE_EQ( BAD_DOUBLE, lines.search( 3.0615E-03 ) );
    }

    // Continuum test
    {
        EmissionSpectrum emission = resource.emission( 58151 );
        EXPECT_EQ( 138, emission.num_continuous() );
        EXPECT_EQ( 0, emission.num_discrete() );
        Spectrum continuous = emission.continuous();
        EXPECT_EQ( 139, continuous.bounds_size() );
        double e1 = continuous.bounds_at( 1 );
        EXPECT_FLOAT_EQ( 5.4862E+00, e1 );
        double e2 = continuous.bounds_at( 2 );
        EXPECT_FLOAT_EQ( 5.4625E+00, e2 );
        double e0 = continuous.bounds_at( 0 );
        double e0_ref = e1 + ( e1 - e2 );  // use delta to get last bound
        EXPECT_FLOAT_EQ( e0_ref, e0 );
    }

    // accumulate source between 1.0 and 0.1 MeV
    Vec_Dbl bounds( {1.0, 0.1} );
    Vec_Dbl intensity( 1, 0.0 );
    for( size_t i = 0; i < resource.emission_size(); ++i )
    {
        EmissionSpectrum emission = resource.emission_at( i );
        emission.bin( bounds, &intensity, 1.2 );
    }
}