text:
#include "hoCuNDArray_utils.h"
#include "radial_utilities.h"
#include "cuNDArray_fileio.h"
#include "cuNDArray_math.h"
#include "imageOperator.h"
#include "identityOperator.h"
#include "hoPartialDerivativeOperator.h"
#include "hoCuConebeamProjectionOperator.h"
#include "cuConvolutionOperator.h"
#include "hoCuNDArray_blas.h"
#include "hoCuNDArray_elemwise.h"
#include "hoCuNDArray_blas.h"
#include "cgSolver.h"
#include "CBCT_acquisition.h"
#include "complext.h"
#include "encodingOperatorContainer.h"
#include "vector_td_io.h"
#include "hoPartialDerivativeOperator.h"
#include "hoCuGPBBSolver.h"
#include "hoCuTvOperator.h"
#include "hoCuTvPicsOperator.h"
#include "hoCuNCGSolver.h"
#include "hoCuCgDescentSolver.h"
#include "hoNDArray_utils.h"
#include "cuPartialDifferenceOperator.h"
#include "cuTvOperator.h"
#include "cuTv1dOperator.h"
#include "cuTvPicsOperator.h"
#include "CBSubsetOperator.h"
#include "osSPSSolver.h"
#include "osMOMSolver.h"
#include "cuOSMOMSolverD.h"
#include "osMOMSolverD2.h"
#include "hoCuNCGSolver.h"
#include "osMOMSolverD3.h"
#include "osMOMSolverL1.h"
#include "osMOMSolverF.h"
#include "osAHZCSolver.h"
#include "ADMMSolver.h"
#include <iostream>
#include <algorithm>
#include <sstream>
#include <math_constants.h>
#include <boost/program_options.hpp>
#include <boost/make_shared.hpp>
#include <GPUTimer.h>
#include <operators/cuGaussianFilterOperator.h>
#include <multiplicationOperatorContainer.h>
#include <cuDownsampleOperator.h>
#include <operators/cuSmallConvOperator.h>
#include <denoise/nonlocalMeans.h>
#include <solvers/osMOMSolverW.h>
#include <operators/cuTFFT.h>
#include "cuSolverUtils.h"
#include "osPDsolver.h"
#include "osLALMSolver.h"
#include "osLALMSolver2.h"
#include "cuATrousOperator.h"
#include "hdf5_utils.h"
#include "cuEdgeATrousOperator.h"
#include "cuDCTOperator.h"
#include "cuDCTDerivativeOperator.h"
#include "dicomWriter.h"
#include "conebeam_projection.h"
#include "weightingOperator.h"
#include "hoNDArray_math.h"
#include "cuNCGSolver.h"
#include "CT_acquisition.h"
#include "cuScaleOperator.h"
#include "cuTVPrimalDualOperator.h"
#include "cuWTVPrimalDualOperator.h"
#include "cuATVPrimalDualOperator.h"
#include "cuSSTVPrimalDualOperator.h"
#include "cuBilateralPriorPrimalDualOperator.h"
#include "cuTV4DPrimalDualOperator.h"
#include "CBSubsetWeightOperator.h"
#include "BilateralPriorOperator.h"
#include "cuPICSDualOperator.h"
#include "cuTVTFFT.h"

using namespace std;
using namespace Gadgetron;

namespace po = boost::program_options;

boost::shared_ptr<hoCuNDArray<float>> downsample_projections(hoCuNDArray<float>* projections, unsigned int num_downsamples )
{

    if (num_downsamples == 0) return boost::make_shared<hoCuNDArray<float>>(*projections);

    auto tmp = Gadgetron::downsample<float,2>(projections);

    for (int k = 1; k < num_downsamples; k++)
        tmp = Gadgetron::downsample<float,2>(tmp.get());

    return boost::make_shared<hoCuNDArray<float>>(*tmp);
}

boost::shared_ptr<cuNDArray<float> > calculate_prior(boost::shared_ptr<CBCT_binning>  binning,boost::shared_ptr<CBCT_acquisition> ps, hoCuNDArray<float>& projections, std::vector<size_t> is_dims, floatd3 imageDimensions){
    std::cout << "Calculating FDK prior" << std::endl;
    boost::shared_ptr<CBCT_binning> binning_pics=binning->get_3d_binning();
    std::vector<size_t> is_dims3d = is_dims;
    is_dims3d.pop_back();
    boost::shared_ptr< hoCuConebeamProjectionOperator >
            Ep( new hoCuConebeamProjectionOperator() );
    Ep->setup(ps,binning_pics,imageDimensions);
    Ep->set_codomain_dimensions(ps->get_projections()->get_dimensions().get());
    Ep->set_domain_dimensions(&is_dims3d);
    Ep->set_use_filtered_backprojection(true);
    boost::shared_ptr<hoCuNDArray<float> > prior3d(new hoCuNDArray<float>(&is_dims3d));
    Ep->mult_MH(&projections,prior3d.get());

    hoCuNDArray<float> tmp_proj(*ps->get_projections());
    Ep->mult_M(prior3d.get(),&tmp_proj);
    float s = dot(ps->get_projections().get(),&tmp_proj)/dot(&tmp_proj,&tmp_proj);
    *prior3d *= s;
    auto prior = boost::make_shared<cuNDArray<float>>(*prior3d);
    return prior;
}

int main(int argc, char** argv)
{

    string acquisition_filename;
    string outputFile;

    uintd3 imageSize;
    floatd3 voxelSize;
    int device;
    floatd2 scale_factor;
    unsigned int iterations;
    unsigned int subsets;
    float rho,tau,nlm_noise,bil_weight;
    float tv_weight,pics_weight, wavelet_weight,huber,sigma,dct_weight,sfr_weight,framelet_weight,atv_weight, sstv_weight;
    float framelet_weight4d;
    float tv_4d,atv_4d;
    bool use_non_negativity;
    int reg_iter;

    po::options_description desc("Allowed options");

    desc.add_options()
            ("help", "produce help message")
            ("acquisition,a", po::value<string>(&acquisition_filename)->default_value("acquisition.hdf5"), "Acquisition data")
            ("samples,n",po::value<unsigned int>(),"Number of samples per ray")
    ("output,f", po::value<string>(&outputFile)->default_value("reconstruction.hdf5"), "Output filename")
            ("size,s",po::value<uintd3>(&imageSize)->default_value(uintd3(512,512,1)),"Image size in pixels")
            ("binning,b",po::value<string>(),"Binning file for 4d reconstruction")
            ("SAG","Use exact SAG correction if present")
            ("voxelSize,v",po::value<floatd3>(&voxelSize)->default_value(floatd3(0.488f,0.488f,1.0f)),"Voxel size in mm")
            ("dimensions,d",po::value<floatd3>(),"Image dimensions in mm. Overwrites voxelSize.")
            ("iterations,i",po::value<unsigned int>(&iterations)->default_value(10),"Number of iterations")
    ("device",po::value<int>(&device)->default_value(0),"Number of the device to use (0 indexed)")
            ("downsample,D",po::value<floatd2>(&scale_factor)->default_value(floatd2(1,1)),"Downsample projections this factor")
            ("subsets,u",po::value<unsigned int>(&subsets)->default_value(10),"Number of subsets to use")
    ("TV",po::value<float>(&tv_weight)->default_value(0),"Total variation weight in spatial dimensions")
            ("TV4D",po::value<float>(&tv_4d)->default_value(0),"Total variation weight in temporal dimensions")
            ("ATV4D",po::value<float>(&atv_4d)->default_value(0),"Advanced Total variation weight in temporal dimensions")
            ("ATV",po::value<float>(&atv_weight)->default_value(0),"Advanced Total variation weight ")
            ("SSTV",po::value<float>(&sstv_weight)->default_value(0),"Scale space total variation weight ")
            ("PICS",po::value<float>(&pics_weight)->default_value(0),"PICS weight")
            ("Wavelet,W",po::value<float>(&wavelet_weight)->default_value(0),"Weight of the wavelet operator")
            ("Framelet",po::value<float>(&framelet_weight)->default_value(0),"Weight of the framelet operator")
            ("Framelet4D",po::value<float>(&framelet_weight4d)->default_value(0),"Weight of the framelet operator")
            ("Huber",po::value<float>(&huber)->default_value(0),"Huber weight")
            ("use_prior","Use an FDK prior")
            ("use_non_negativity",po::value<bool>(&use_non_negativity)->default_value(true),"Prevent image from having negative attenuation")
            ("sigma",po::value<float>(&sigma)->default_value(0.1),"Sigma for billateral filter")
            ("DCT",po::value<float>(&dct_weight)->default_value(0),"DCT regularization")
            ("SFR",po::value<float>(&sfr_weight)->default_value(0),"SFR regularization")
            ("3D","Only use binning for selecting valid projections")
            ("tau",po::value<float>(&tau)->default_value(1e-5),"Tau value for solver")
            ("reg_iter",po::value<int>(&reg_iter)->default_value(2))
            ("bilateral-weight",po::value<float>(&bil_weight)->default_value(0),"Bilateral weight")
            ("prior",po::value<string>(),"prior image")
            ("projection_weights",po::value<string>(),"Array containing weights to be applied to the projections.")

            ("NLM",po::value<float>(&nlm_noise)->default_value(0),"Use non-local means based on the 3d image")
            ;

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

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

    std::stringstream command_line_string;
    std::cout << "Command line options:" << std::endl;
    for (po::variables_map::iterator it = vm.begin(); it != vm.end(); ++it){
        boost::any a = it->second.value();
        command_line_string << it->first << ": ";
        if (a.type() == typeid(std::string)) command_line_string << it->second.as<std::string>();
        else if (a.type() == typeid(int)) command_line_string << it->second.as<int>();
        else if (a.type() == typeid(unsigned int)) command_line_string << it->second.as<unsigned int>();
        else if (a.type() == typeid(float)) command_line_string << it->second.as<float>();
        else if (a.type() == typeid(vector_td<float,3>)) command_line_string << it->second.as<vector_td<float,3> >();
        else if (a.type() == typeid(vector_td<float,2>)) command_line_string << it->second.as<vector_td<float,2> >();
        else if (a.type() == typeid(vector_td<int,3>)) command_line_string << it->second.as<vector_td<int,3> >();
        else if (a.type() == typeid(vector_td<unsigned int,3>)) command_line_string << it->second.as<vector_td<unsigned int,3> >();
        else if (a.type() == typeid(bool)) command_line_string << it->second.as<bool>();
        else command_line_string << "Unknown type" << std::endl;
        command_line_string << std::endl;
    }
    std::cout << command_line_string.str();

    cudaSetDevice(device);
    cudaDeviceReset();

    //Really weird stuff. Needed to initialize the device?? Should find real bug.
    cudaDeviceManager::Instance()->lockHandle();
    cudaDeviceManager::Instance()->unlockHandle();

    boost::shared_ptr<CBCT_acquisition> ps(new CBCT_acquisition());
    ps->load(acquisition_filename);
    ps->get_geometry()->print(std::cout);

    float SDD = ps->get_geometry()->get_SDD();
    float SAD = ps->get_geometry()->get_SAD();

    boost::shared_ptr<CBCT_binning> binning(new CBCT_binning());
    if (vm.count("binning")){
        std::cout << "Loading binning data" << std::endl;
        binning->load(vm["binning"].as<string>());
        if (vm.count("3D"))
            binning = binning->get_3d_binning();
    } else binning->set_as_default_3d_bin(ps->get_projections()->get_size(2));
    binning->print(std::cout);

    floatd3 imageDimensions;
    if (vm.count("dimensions")){
        imageDimensions = vm["dimensions"].as<floatd3>();
        voxelSize = imageDimensions/imageSize;
    }
    else imageDimensions = voxelSize*imageSize;

    float lengthOfRay_in_mm = norm(imageDimensions);
    unsigned int numSamplesPerPixel = 3;
    float minSpacing = min(voxelSize)/numSamplesPerPixel;

    unsigned int numSamplesPerRay;
    if (vm.count("samples")) numSamplesPerRay = vm["samples"].as<unsigned int>();
    else numSamplesPerRay = ceil( lengthOfRay_in_mm / minSpacing );

    float step_size_in_mm = lengthOfRay_in_mm / numSamplesPerRay;
    size_t numProjs = ps->get_projections()->get_size(2);
    size_t needed_bytes = 2 * prod(imageSize) * sizeof(float);
    std::vector<size_t> is_dims = to_std_vector((uint64d3)imageSize);

    std::cout << "IS dimensions " << is_dims[0] << " " << is_dims[1] << " " << is_dims[2] << std::endl;
    std::cout << "Image size " << imageDimensions << std::endl;
    std::cout << "FRAMELET WEIGHT " << framelet_weight << std::endl;
    is_dims.push_back(binning->get_number_of_bins());

    //scatter_correct(binning,ps,ps->get_projections().get(),is_dims,imageDimensions);

    if (scale_factor[0] != 1 || scale_factor[1] != 1)
        ps->downsample(scale_factor[0],scale_factor[1]);

    //osLALMSolver<cuNDArray<float>> solver;
    osMOMSolverD<cuNDArray<float>> solver;
    //osMOMSolverW<cuNDArray<float>> solver;
    //osMOMSolverL1<cuNDArray<float>> solver;
    //osAHZCSolver<cuNDArray<float>> solver;
    //osMOMSolverF<cuNDArray<float>> solver;
    //ADMMSolver<cuNDArray<float>> solver;
    solver.set_dump(false);

    boost::shared_ptr<cuNDArray<float>> prior;
    if (vm.count("use_prior") ) {
        auto projections = *ps->get_projections();
        prior = calculate_prior(binning,ps,projections,is_dims,imageDimensions);
        solver.set_x0(expand(prior.get(),binning->get_number_of_bins()));
        //prior = calculate_weightImage(binning,ps,projections,is_dims,imageDimensions);
//        solver.set_x0(prior);
    }
    if (vm.count("prior")){
        prior = boost::make_shared<cuNDArray<float>>(*read_nd_array<float>(vm["prior"].as<string>().c_str()));
    }

    solver.set_max_iterations(iterations);
    solver.set_output_mode(osSPSSolver<cuNDArray<float>>::OUTPUT_VERBOSE);
    solver.set_tau(tau);
    solver.set_non_negativity_constraint(use_non_negativity);
    solver.set_huber(huber);
    solver.set_reg_steps(reg_iter);
    //solver.set_rho(rho);

    if (atv_weight > 0) {

        auto ATV = boost::make_shared<cuATVPrimalDualOperator<float>>(0.0);
        ATV->set_weight(atv_weight);
        solver.add_regularization_operator(ATV);

    }

    if (sstv_weight > 0){
        auto SSTV = boost::make_shared<cuSSTVPrimalDualOperator<float>>();
        SSTV->set_weight(sstv_weight);
        solver.add_regularization_operator(SSTV);
    }

    if (tv_weight > 0) {

        //auto TV = boost::make_shared<cuWTVPrimalDualOperator<float>>(0.0,1e-3);
//

        auto TV = boost::make_shared<cuTVPrimalDualOperator<float>>(0.0);
        TV->set_weight(tv_weight);

        solver.add_regularization_operator(TV);

/*
      auto Dx = boost::make_shared<cuPartialDifferenceOperator<float, 3>>(0);
      Dx->set_weight(tv_weight);
      Dx->set_domain_dimensions(&is_dims);
      Dx->set_codomain_dimensions(&is_dims);

      auto Dy = boost::make_shared<cuPartialDifferenceOperator<float, 3>>(1);
      Dy->set_weight(tv_weight);
      Dy->set_domain_dimensions(&is_dims);
      Dy->set_codomain_dimensions(&is_dims);

      auto Dz = boost::make_shared<cuPartialDifferenceOperator<float, 3>>(2);
      Dz->set_weight(tv_weight);
      Dz->set_domain_dimensions(&is_dims);
      Dz->set_codomain_dimensions(&is_dims);
      solver.add_regularization_group({Dx, Dy, Dz});
*/

    }
    if (tv_4d > 0) {

        auto Dt = boost::make_shared<cuPartialDifferenceOperator<float, 4>>(3);
        Dt->set_weight(tv_4d);
        Dt->set_domain_dimensions(&is_dims);
        Dt->set_codomain_dimensions(&is_dims);
        solver.add_regularization_operator(Dt);
/**/

    }
    if (atv_4d > 0) {

        auto TVF =  boost::make_shared<cuTVTFFT>();
        TVF->set_weight(atv_4d);
        TVF->set_domain_dimensions(&is_dims);
        solver.add_regularization_operator(TVF);

/*
		  auto Dt = boost::make_shared<cuTV4DPrimalDualOperator<float>>(0);
		  Dt->set_weight(atv_4d);
		  solver.add_regularization_operator(Dt);
*/

    }

    if (framelet_weight > 0){
        auto stencils = std::vector<vector_td<float,3>>({
                                                                vector_td<float,3>(-std::sqrt(2.0f),0,std::sqrt(2.0f)),vector_td<float,3>(-1,2,-1),vector_td<float,3>(1,2,1) });
        std::cout << "Framelet weight " << framelet_weight << std::endl;
        for (int step_size = 0; step_size < 3; step_size++) {
            for (auto stencil : stencils) {
                std::cout << "Adding operators " << std::endl;
                auto Rx = boost::make_shared<cuSmallConvOperator<float, 4, 3>>(stencil, 0,std::pow(2,step_size));
                Rx->set_weight(framelet_weight);
                Rx->set_domain_dimensions(&is_dims);
                Rx->set_codomain_dimensions(&is_dims);

                auto Ry = boost::make_shared<cuSmallConvOperator<float, 4, 3>>(stencil, 1,std::pow(2,step_size));
                Ry->set_weight(framelet_weight);
                Ry->set_domain_dimensions(&is_dims);
                Ry->set_codomain_dimensions(&is_dims);

                auto Rz = boost::make_shared<cuSmallConvOperator<float, 4, 3>>(stencil, 2,std::pow(2,step_size));
                Rz->set_weight(framelet_weight);
                Rz->set_domain_dimensions(&is_dims);
                Rz->set_codomain_dimensions(&is_dims);
                solver.add_regularization_group({Rx, Ry, Rz});
//                solver.add_regularization_operator(Rx);
//                solver.add_regularization_operator(Ry);
//                solver.add_regularization_operator(Rz);

            }
        }
    }
    if (framelet_weight4d > 0){
        auto stencils = std::vector<vector_td<float,3>>({
                                                                vector_td<float,3>(-std::sqrt(2.0f),0,std::sqrt(2.0f)),vector_td<float,3>(-1,2,-1),vector_td<float,3>(1,2,1) });
        for (int step_size = 0; step_size < 3; step_size++) {
            for (auto stencil : stencils) {
                auto Rt = boost::make_shared<cuSmallConvOperator<float, 4, 3>>(stencil, 3,std::pow(2,step_size));
                Rt->set_weight(framelet_weight4d);
                Rt->set_domain_dimensions(&is_dims);
                Rt->set_codomain_dimensions(&is_dims);
                solver.add_regularization_operator(Rt);
            }
        }
    }

    if (dct_weight > 0){
        auto dctOp = boost::make_shared<cuDCTOperator<float>>();
//        auto dctOp = boost::make_shared<cuTFFT>();
        dctOp->set_domain_dimensions(&is_dims);
        dctOp->set_weight(dct_weight);
        solver.add_regularization_operator(dctOp);
    }

    if (sfr_weight > 0){
        auto sfrOp = boost::make_shared<cuTFFT>();
        sfrOp->set_domain_dimensions(&is_dims);
        sfrOp->set_weight(sfr_weight);
        solver.add_regularization_operator(sfrOp);
    }

    if (bil_weight > 0){
        auto bilOp = boost::make_shared<cuBilateralPriorPrimalDualOperator>(0.002,3.0,prior);
        bilOp->set_weight(bil_weight);
        solver.add_regularization_operator(bilOp);

        /*
        auto bilOp = boost::make_shared<BilateralPriorOperator>();
        bilOp->set_domain_dimensions(&is_dims);
        bilOp->set_codomain_dimensions(&is_dims);
        bilOp->set_weight(bil_weight);
        bilOp->set_sigma_spatial(3.0);
        bilOp->set_sigma_int(0.01);

        bilOp->set_prior(cuBilPrior);
        solver.add_regularization_operator(bilOp,0.01);
         */
    }

    if (pics_weight > 0){
        auto pics = boost::make_shared<cuPICSPrimalDualOperator<float>>();
        pics->set_weight(pics_weight);
        pics->set_prior(prior);
        solver.add_regularization_operator(pics);
    }

    boost::shared_ptr<CBSubsetOperator<cuNDArray>> E;
    if ( vm.count("projection_weights")) {
        auto weights = read_nd_array<float>(vm["projection_weights"].as<string>().c_str());
        if (scale_factor[0] != 1 || scale_factor[1] != 1) {
            cuNDArray<float> tmp_weights(*weights);
            auto dims = *tmp_weights.get_dimensions();

            tmp_weights = downsample_projections(&tmp_weights,scale_factor[0],scale_factor[1]);

            weights = tmp_weights.to_host();
        }

        if (!ps->get_projections()->dimensions_equal(weights.get()))
            throw std::runtime_error("Weight dimensions must match that of the projection data");
        auto EW  = boost::make_shared<CBSubsetWeightOperator<cuNDArray>>(subsets);
        EW->setup(ps,binning,imageDimensions,weights);
        E = EW;
    } else {
        std::cout <<"Normal projections" << std::endl;
        E = boost::make_shared<CBSubsetOperator<cuNDArray> >(subsets);

        E->setup(ps, binning, imageDimensions);
    }
    E->set_domain_dimensions(&is_dims);
    E->set_codomain_dimensions(ps->get_projections()->get_dimensions().get());

    solver.set_encoding_operator(E);

    auto projections = boost::make_shared<cuNDArray<float>>(*ps->get_projections());
    std::cout << "Projection norm:" << nrm2(projections.get()) << std::endl;

    boost::shared_ptr<cuNDArray<float>> result;
    {
        GPUTimer tim("Solver");
        result = solver.solve(projections.get());
    }
//	global_timer.reset();
    std::cout << "Penguin" << nrm2(result.get()) << std::endl;

    std::cout << "Result sum " << asum(result.get()) << std::endl;

    //apply_mask(result.get(),mask.get());

    std::cout << "Result sum " << asum(result.get()) << std::endl;
    //saveNDArray2HDF5(result.get(),outputFile,imageDimensions,vector_td<float,3>(0),command_line_string.str(),iterations);

    write_nd_array(result.get(),"reconstruction.real");

    if (wavelet_weight > 0){
        osMOMSolverF<cuNDArray<float>> solverF;
        solverF.set_max_iterations(10);
        solverF.set_x0(result);
        solverF.set_encoding_operator(E);
        solverF.set_non_negativity_constraint(true);

        auto wave = boost::make_shared<cuEdgeATrousOperator<float>>();

        wave->set_sigma(sigma);
        wave->set_domain_dimensions(&is_dims);
        if (binning->get_number_of_bins() == 1)
            wave->set_levels({2,2,2});
        else
            wave->set_levels({2,2,2,2});
        wave->set_weight(wavelet_weight);
        solverF.add_regularization_operator(wave);

        result = solverF.solve(projections.get());

    }
/*
	if (nlm_noise > 0){

		float* result_data = result->get_data_ptr();

		for (int i = 0; i < result->get_size(3); i++){
			cuNDArray<float> result_view(prior->get_dimensions(),result_data);

			nonlocal_means_ref(&result_view,prior.get(),nlm_noise);
			result_data += result_view.get_number_of_elements();
		}
	}
*/
    saveNDArray2HDF5(result.get(),outputFile,imageDimensions,floatd3(0,0,0),command_line_string.str(),iterations);
//	write_nd_array(result.get(),"reconstruction.real");
    //write_dicom(result.get(),command_line_string.str(),imageDimensions);

}

max_stars_repo_name: ahsanjav/gt-tomography

max_issues_repo_name: ahsanjav/gt-tomography

max_forks_repo_name: ahsanjav/gt-tomography


text:
/*
// Copyright (c) 2000-2009, Texas Engineering Experiment Station (TEES), a
// component of the Texas A&M University System.

// All rights reserved.

// The information and source code contained herein is the exclusive
// property of TEES and may not be disclosed, examined or reproduced
// in whole or in part without explicit written authorization from TEES.
*/

#ifndef STAPL_VIEWS_SYSTEM_VIEW_HPP
#define STAPL_VIEWS_SYSTEM_VIEW_HPP

#include <stapl/domains/interval.hpp>
#include <stapl/views/distribution_spec_view.hpp>
#include <stapl/containers/distribution/specification_functors.hpp>
#include <stapl/views/proxy/trivial_accessor.hpp>
#include <stapl/views/metadata/extraction/generator.hpp>
#include <boost/optional.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/utility/typed_in_place_factory.hpp>
#include <vector>

namespace stapl {

namespace dist_view_impl {

//////////////////////////////////////////////////////////////////////
/// @brief Mapping from location id in the creator gang to a location id
/// in the created gang based on an explicit set of locations that will
/// be part of the created gang.
///
/// @tparam LocDom Type of the domain used for the explicit location id set.
//////////////////////////////////////////////////////////////////////
template <typename LocDom>
struct explicit_loc_mapping_function
{
private:
  LocDom m_locs;

public:
  explicit_loc_mapping_function(LocDom const& locs)
    : m_locs(locs)
  { }

  stapl::location_type operator()(stapl::location_type id) const
  { return m_locs.distance(m_locs.first(), id); }

  bool operator==(explicit_loc_mapping_function const& other) const
  { return m_locs == other.m_locs; }

  void define_type(typer& t)
  { t.member(m_locs); }
};

//////////////////////////////////////////////////////////////////////
/// @brief Mapping from location id in the created gang to a location id
/// in the creator gang based on an explicit set of locations that will
/// be part of the created gang.
///
/// This mapping is the inverse of the mapping provided by
/// @ref explicit_loc_mapping_function.
///
/// @tparam LocDom Type of the domain used for the explicit location id set.
//////////////////////////////////////////////////////////////////////
template <typename LocDom>
struct explicit_loc_resolution_function
{
private:
  LocDom m_locs;

public:
  explicit_loc_resolution_function(LocDom const& locs)
    : m_locs(locs)
  { }

  stapl::location_type operator()(stapl::location_type id) const
  { return m_locs.advance(m_locs.first(), id); }

  bool operator==(explicit_loc_resolution_function const& other) const
  { return m_locs == other.m_locs; }

  void define_type(typer& t)
  { t.member(m_locs); }
};

//////////////////////////////////////////////////////////////////////
/// @brief Basic @ref p_object used by @ref system_container to initialize
///  a gang over its set of locations.
//////////////////////////////////////////////////////////////////////
struct trivial_p_object
  : public p_object
{ };

struct system_container;

//////////////////////////////////////////////////////////////////////
/// @todo either remove metadata from system view (does it make sense?)
///  or properly compute affinity information.
//////////////////////////////////////////////////////////////////////
struct system_container_distribution
  : public p_object
{
  typedef indexed_domain<location_type>            domain_type;
  typedef system_container*                        component_type;
  typedef metadata_entry<
    domain_type, component_type, unsigned int
  >                                                dom_info_type;

  future<dom_info_type> metadata_at(size_t gid)
  {
    const unsigned int id = 0;

    return make_ready_future(dom_info_type(
             id, domain_type(location_type(gid), location_type(gid)), 0,
             LQ_CERTAIN, affinity_tag(), this->get_rmi_handle(), gid));
  }
};

//////////////////////////////////////////////////////////////////////
/// @brief Container that represents the set of locations in the system.
///
/// The @ref identity_container was considered as an alternative to this
/// implementation, but lacked the finite size required to represent the
/// finite set of locations being used to execute an application.
/// @ingroup system_view
/// @internal
//////////////////////////////////////////////////////////////////////
struct system_container
{
  typedef location_type                                     value_type;
  typedef proxy<location_type,
            trivial_accessor<location_type> >               reference;
  // trivial_accessor provides read only.
  // reference and const_reference are the same type.
  typedef reference                                         const_reference;
  typedef domset1D<location_type>                           domain_type;
  typedef location_type                                     gid_type;
  typedef unsigned int                                      cid_type;
  typedef system_container_distribution                     distribution_type;
  typedef metadata::generator_extractor<system_container>   loc_dist_metadata;

private:
  location_type                                  m_size;
  distribution_type                              m_dist;
  level                                          m_level;
  boost::shared_ptr<std::vector<location_type>>  m_locs;
  boost::optional<trivial_p_object>              m_gang;

public:
  system_container(location_type nlocs)
    : m_size(nlocs), m_level(level::invalid), m_locs()
  { }

  system_container(level const& lvl)
    : m_size(0), m_level(lvl), m_locs()
  {
    stapl_assert(lvl == current_level || lvl == lowest_level,
     "unsupported level specifier in system container");
  }

  system_container(std::vector<location_type> const& locs)
    : m_size(locs.size()), m_level(level::invalid),
      m_locs(new std::vector<location_type>(locs))
  { }

  distribution_type& distribution(void)
  { return m_dist; }

  value_type get_element(gid_type const& index)
  {
    if (!m_size)
      m_size = stapl::get_num_locations();

    stapl_assert(index < m_size,
      "system_container::get_element received out of range index");
    if (!m_locs)
      return index;
    else
      return (*m_locs)[index];
  }

  const_reference operator[](gid_type const& index)
  {
    if (!m_size)
      m_size = stapl::get_num_locations();

    stapl_assert(index < m_size,
      "system_container::operator[] received out of range index");
    if (!m_locs)
      return reference(trivial_accessor<value_type>(index));
    else
      return reference(trivial_accessor<value_type>((*m_locs)[index]));
  }

  //////////////////////////////////////////////////////////////////////
  /// @brief Returns a reference to the element indexed.
  ///
  /// @param index of element to return
  /// @return element reference
  //////////////////////////////////////////////////////////////////////
  const_reference make_reference(gid_type const& index)
  {
    return this->operator[](index);
  }

  //////////////////////////////////////////////////////////////////////
  /// @internal
  //////////////////////////////////////////////////////////////////////
  size_t version(void) const
  {
    return 0;
  }

  size_t size(void) const
  {
    stapl_assert(m_size, "system_container uninitialized when size called");
    return m_size;
  }

  bool explicit_locs(void) const
  { return m_locs != nullptr; }

  std::vector<location_type> const& get_locs(void) const
  {
    if (!this->explicit_locs())
      stapl::abort("set of locations in system view is not set");

    return *m_locs;
  }

  bool level_specified(void) const
  { return m_level != level::invalid; }

  level level_spec(void) const
  { return m_level; }

  domain_type domain(void) const
  {
    if (!m_locs)
      return domain_type(0, m_size-1);
    else
    {
      if (m_locs->back() - m_locs->front() + 1 == m_size)
        return domain_type(m_locs->front(), m_locs->back());
      else
      {
        domain_type d;
        for (auto loc : *m_locs)
          d += loc;
        return d;
      }
    }
  }

  //////////////////////////////////////////////////////////////////////
  /// @brief Create a @ref gang representing the locations this system
  ///   container instance refers to, on demand and return a reference
  ///   to an @ref p_object in that communication group, so that other
  ///   p_objects can be initialized within it by utilizing gang switching.
  //////////////////////////////////////////////////////////////////////
  p_object const& group(void)
  {
    if (level_specified())
    {
       if (m_level == current_level)
       {
         if (!m_gang)
           m_gang = boost::in_place<trivial_p_object>();

         return *m_gang;
       }

       // else
       if (m_level == lowest_level)
       {
         if (!m_gang)
         {
           gang g;
           m_gang = boost::in_place<trivial_p_object>();
         }

        return *m_gang;
       }
    }

    stapl_assert(this->domain().contains(get_location_id()),
      "location shouldn't be participating in sub gang creation");

    if (!m_gang)
    {
      domain_type loc_domain = this->domain();

      gang g(loc_domain.size(),
             explicit_loc_mapping_function<domain_type>(loc_domain),
             explicit_loc_resolution_function<domain_type>(loc_domain));

      m_gang = boost::in_place<trivial_p_object>();
    }

    return *m_gang;
  }

  template <typename F>
  typename F::result_type apply_get(gid_type const& gid, F const& f) const
  {
    if (!m_locs)
      return f(gid);
    else
      return f((*m_locs)[gid]);
  }

  //////////////////////////////////////////////////////////////////////
  /// @internal
  //////////////////////////////////////////////////////////////////////
  void define_type(typer& t)
  {
    t.member(m_size);
    t.member(m_locs);
  }
};

} // namespace dist_view_impl

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing the set of locations
/// in the system.
///
/// @return A system_view that represents the locations of the system
///         used in the current computation.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, location_type, location_type>*
system_view(location_type nlocs = get_num_locations());

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing a set of locations
/// in the system.
///
/// @return A system_view that represents the subset of the system whose
///         location ids are specified in the vector provided.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, location_type, location_type>*
system_view(std::vector<location_type> const& locs);

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing a set of locations
/// in the system.
///
/// @return A system_view that represents the subset of the system whose
///         locations are contained in the next lower level of the
///         system hierarchy.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, location_type, location_type>*
system_view(level const lvl);

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing a set of locations
/// in the system.
///
/// @tparam ContainerGID Type of GID in the container to be constructed
///         using a view-based specification of which the system_view
///         being constructed is the base.
/// @tparam ContainerCID Type of the partition id in the container to be
///         constructed using a view-based specification of which the
///         system_view being constructed is the base.
///
/// @return A system_view that represents the subset of the system whose
///         locations are contained in the next lower level of the
///         system hierarchy.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
template <typename ContainerGID, typename ContainerCID>
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, ContainerGID, ContainerCID>*
system_view(level const lvl)
{
  typedef distribution_spec_view<dist_view_impl::system_container,
            indexed_domain<location_type>, ContainerGID, ContainerCID> type;
  return new type(new dist_view_impl::system_container(lvl),
    indexed_domain<location_type>(0),
    identity_map<location_type, location_type, ContainerGID, ContainerCID>());
}

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing a set of locations
/// in the system using an external @ref system_container.
///
/// @tparam ContainerGID Type of GID in the container to be constructed
///         using a view-based specification of which the system_view
///         being constructed is the base.
/// @tparam ContainerCID Type of the partition id in the container to be
///         constructed using a view-based specification of which the
///         system_view being constructed is the base.
/// @param  sys_ct System container initialized and managed externally.
///
/// @return A system_view that represents the subset of the system whose
///         locations are contained in the next lower level of the
///         system hierarchy.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
template <typename ContainerGID, typename ContainerCID>
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, ContainerGID, ContainerCID>*
system_view(dist_view_impl::system_container& sys_ct)
{
  typedef distribution_spec_view<dist_view_impl::system_container,
            indexed_domain<location_type>, ContainerGID, ContainerCID> type;
  return new type(
    sys_ct,
    indexed_domain<location_type>(0),
    identity_map<location_type, location_type, ContainerGID, ContainerCID>());
}

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view representing a set of @c locs
/// in the system using an external @ref system_container.
///
/// @tparam ContainerGID Type of GID in the container to be constructed
///         using a view-based specification of which the system_view
///         being constructed is the base.
/// @tparam ContainerCID Type of the partition id in the container to be
///         constructed using a view-based specification of which the
///         system_view being constructed is the base.
/// @param  locs Explicit set of locations in the system view.
///
/// @return A system_view that represents the subset of the system whose
///         locations are contained in the next lower level of the
///         system hierarchy.
/// @ingroup system_view
//////////////////////////////////////////////////////////////////////
template <typename ContainerGID, typename ContainerCID>
distribution_spec_view<dist_view_impl::system_container,
  indexed_domain<location_type>, ContainerGID, ContainerCID>*
system_view(std::vector<location_type> const& locs)
{
  typedef distribution_spec_view<dist_view_impl::system_container,
            indexed_domain<location_type>, ContainerGID, ContainerCID> type;
  return new type(
    new dist_view_impl::system_container(locs),
    indexed_domain<location_type>(0),
    identity_map<location_type, location_type, ContainerGID, ContainerCID>());
}

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view whose domain is the identifiers of
/// the partitions of an element partitioning and whose mapping function
/// maps from partition id to location id.
///
/// The view is defined over a view that represents the set of
/// locations on to which this view will map partitions.  The set of
/// operations on the view is currently empty and will be expanded as useful
/// operations are identified.
///
/// @note The system view must have a lifetime that exceeds the mapping_view.
///
/// @param system view over the system on to which parition ids will be mapped
/// @param domain set of partition ids to be mapped to locations.
/// @param mapping_func Unary functor that receives a partition id
/// and returns a location id.
/// @return A view capable of mapping partition ids to locations.
/// @ingroup mapping_view
//////////////////////////////////////////////////////////////////////
template <typename SysView, typename PartitionIds, typename MappingFunction,
          typename GID = typename PartitionIds::index_type>
distribution_spec_view<SysView, PartitionIds, GID,
  typename PartitionIds::index_type>
mapping_view(SysView& system, PartitionIds const& domain,
             MappingFunction const& mapping_func)
{
  return
    distribution_spec_view<SysView, PartitionIds, GID,
      typename PartitionIds::index_type>(&system, domain, mapping_func);
}

//////////////////////////////////////////////////////////////////////
/// @brief Creates a read-only view whose domain is the identifiers of
/// a container's elements and whose mapping function maps from element
/// GID to partition id.
///
/// When the partitioning_view is defined over a mapping_view the result
/// is a view which defines the distribution of container elements across
/// the locations of the system.
///
/// @param mapping_view view that will map the partition ids to locations
/// @param domain Set of element GIDs to be partitioned.
/// @param partitioning_func Unary functor that receives a GID and returns
/// a partition id.
/// @return A view capable of partitioning GIDs of container elements.
/// @ingroup partitioning_view
//////////////////////////////////////////////////////////////////////
template <typename MappingView, typename GIDDomain,
          typename PartitioningFunction>
distribution_spec_view<MappingView, GIDDomain, typename GIDDomain::index_type,
  typename MappingView::index_type>
partitioning_view(MappingView& mapping_view, GIDDomain const& domain,
                  PartitioningFunction const& partitioning_func)
{
  return distribution_spec_view<
    MappingView, GIDDomain,
    typename GIDDomain::index_type, typename MappingView::index_type
   >(&mapping_view, domain, partitioning_func);
}

} // namespace stapl

#endif // STAPL_VIEWS_SYSTEM_VIEW_HPP

ext: hpp


text:
/*
 * Copyright (C) 2012-2014 Open Source Robotics Foundation
 *
 * 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 <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/iostreams/filter/bzip2.hpp>
#include <boost/iostreams/filter/zlib.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/copy.hpp>
#include <boost/archive/iterators/base64_from_binary.hpp>
#include <boost/archive/iterators/binary_from_base64.hpp>
#include <boost/archive/iterators/remove_whitespace.hpp>
#include <boost/archive/iterators/istream_iterator.hpp>
#include <boost/archive/iterators/transform_width.hpp>

#include "gazebo/math/Rand.hh"

#include "gazebo/common/Exception.hh"
#include "gazebo/common/Console.hh"
#include "gazebo/common/Base64.hh"
#include "gazebo/util/LogRecord.hh"
#include "gazebo/util/LogPlay.hh"

using namespace gazebo;
using namespace util;

/////////////////////////////////////////////////
LogPlay::LogPlay()
{
  this->logStartXml = NULL;
}

/////////////////////////////////////////////////
LogPlay::~LogPlay()
{
}

/////////////////////////////////////////////////
void LogPlay::Open(const std::string &_logFile)
{
  boost::filesystem::path path(_logFile);
  if (!boost::filesystem::exists(path))
    gzthrow("Invalid logfile[" + _logFile + "]. Does not exist.");

  // Parse the log file
  if (!this->xmlDoc.LoadFile(_logFile))
    gzthrow("Unable to parse log file[" << _logFile << "]");

  // Get the gazebo_log element
  this->logStartXml = this->xmlDoc.FirstChildElement("gazebo_log");

  if (!this->logStartXml)
    gzthrow("Log file is missing the <gazebo_log> element");

  // Store the filename for future use.
  this->filename = _logFile;

  // Read in the header.
  this->ReadHeader();

  this->logCurrXml = this->logStartXml;
  this->encoding.clear();
}

/////////////////////////////////////////////////
std::string LogPlay::GetHeader() const
{
  std::ostringstream stream;
  stream << "<?xml version='1.0'?>\n"
         << "<gazebo_log>\n"
         << "<header>\n"
         << "<log_version>" << this->logVersion << "</log_version>\n"
         << "<gazebo_version>" << this->gazeboVersion << "</gazebo_version>\n"
         << "<rand_seed>" << this->randSeed << "</rand_seed>\n"
         << "</header>\n";

  return stream.str();
}

/////////////////////////////////////////////////
void LogPlay::ReadHeader()
{
  this->randSeed = math::Rand::GetSeed();
  TiXmlElement *headerXml, *childXml;

  this->logVersion.clear();
  this->gazeboVersion.clear();

  // Get the header element
  headerXml = this->logStartXml->FirstChildElement("header");
  if (!headerXml)
    gzthrow("Log file has no header");

  // Get the log format version
  childXml = headerXml->FirstChildElement("log_version");
  if (!childXml)
    gzerr << "Log file header is missing the log version.\n";
  else
    this->logVersion = childXml->GetText();

  // Get the gazebo version
  childXml = headerXml->FirstChildElement("gazebo_version");
  if (!childXml)
    gzerr << "Log file header is missing the gazebo version.\n";
  else
    this->gazeboVersion = childXml->GetText();

  // Get the random number seed.
  childXml = headerXml->FirstChildElement("rand_seed");
  if (!childXml)
    gzerr << "Log file header is missing the random number seed.\n";
  else
    this->randSeed = boost::lexical_cast<uint32_t>(childXml->GetText());

  if (this->logVersion != GZ_LOG_VERSION)
    gzwarn << "Log version[" << this->logVersion << "] in file["
           << this->filename
           << "] does not match Gazebo's log version["
           << GZ_LOG_VERSION << "]\n";

  /// Set the random number seed for simulation
  math::Rand::SetSeed(this->randSeed);
}

/////////////////////////////////////////////////
bool LogPlay::IsOpen() const
{
  return this->logStartXml != NULL;
}

/////////////////////////////////////////////////
std::string LogPlay::GetLogVersion() const
{
  return this->logVersion;
}

/////////////////////////////////////////////////
std::string LogPlay::GetGazeboVersion() const
{
  return this->gazeboVersion;
}

/////////////////////////////////////////////////
uint32_t LogPlay::GetRandSeed() const
{
  return this->randSeed;
}

/////////////////////////////////////////////////
bool LogPlay::Step(std::string &_data)
{
  std::string startMarker = "<sdf ";
  std::string endMarker = "</sdf>";
  size_t start = this->currentChunk.find(startMarker);
  size_t end = this->currentChunk.find(endMarker);

  if (start == std::string::npos || end == std::string::npos)
  {
    this->currentChunk.clear();

    if (this->logCurrXml == this->logStartXml)
      this->logCurrXml = this->logStartXml->FirstChildElement("chunk");
    else if (this->logCurrXml)
    {
      this->logCurrXml = this->logCurrXml->NextSiblingElement("chunk");
    }
    else
      return false;

    // Stop if there are no more chunks
    if (!this->logCurrXml)
      return false;

    if (!this->GetChunkData(this->logCurrXml, this->currentChunk))
    {
      gzerr << "Unable to decode log file\n";
      return false;
    }

    start = this->currentChunk.find(startMarker);
    end = this->currentChunk.find(endMarker);
  }

  _data = this->currentChunk.substr(start, end+endMarker.size()-start);

  this->currentChunk.erase(0, end + endMarker.size());

  return true;
}

/////////////////////////////////////////////////
bool LogPlay::GetChunk(unsigned int _index, std::string &_data)
{
  unsigned int count = 0;
  TiXmlElement *xml = this->logStartXml->FirstChildElement("chunk");

  while (xml && count < _index)
  {
    count++;
    xml = xml->NextSiblingElement("chunk");
  }

  if (xml && count == _index)
    return this->GetChunkData(xml, _data);
  else
    return false;
}

/////////////////////////////////////////////////
bool LogPlay::GetChunkData(TiXmlElement *_xml, std::string &_data)
{
  // Make sure we have valid xml pointer
  if (!_xml)
    return false;

  /// Get the chunk's encoding
  this->encoding = _xml->Attribute("encoding");

  // Make sure there is an encoding value.
  if (this->encoding.empty())
    gzthrow("Enconding missing for a chunk in log file[" +
        this->filename + "]");

  if (this->encoding == "txt")
    _data = _xml->GetText();
  else if (this->encoding == "bz2")
  {
    std::string data = _xml->GetText();
    std::string buffer;

    // Decode the base64 string
    buffer = Base64Decode(data);

    // Decompress the bz2 data
    {
      boost::iostreams::filtering_istream in;
      in.push(boost::iostreams::bzip2_decompressor());
      in.push(boost::make_iterator_range(buffer));

      // Get the data
      std::getline(in, _data, '\0');
      _data += '\0';
    }
  }
  else if (this->encoding == "zlib")
  {
    std::string data = _xml->GetText();
    std::string buffer;

    // Decode the base64 string
    buffer = Base64Decode(data);

    // Decompress the zlib data
    {
      boost::iostreams::filtering_istream in;
      in.push(boost::iostreams::zlib_decompressor());
      in.push(boost::make_iterator_range(buffer));

      // Get the data
      std::getline(in, _data, '\0');
      _data += '\0';
    }
  }
  else
  {
    gzerr << "Inavlid encoding[" << this->encoding << "] in log file["
      << this->filename << "]\n";
    return false;
  }

  return true;
}

/////////////////////////////////////////////////
std::string LogPlay::GetEncoding() const
{
  return this->encoding;
}

/////////////////////////////////////////////////
unsigned int LogPlay::GetChunkCount() const
{
  unsigned int count = 0;
  TiXmlElement *xml = this->logStartXml->FirstChildElement("chunk");

  while (xml)
  {
    count++;
    xml = xml->NextSiblingElement("chunk");
  }

  return count;
}

ext: cc

max_stars_repo_name: hyunoklee/Gazebo

max_issues_repo_name: hyunoklee/Gazebo

max_forks_repo_name: hyunoklee/Gazebo


text:
//
// $Id$
//
//
// Original author: Matt Chambers <matt.chambers .@. vanderbilt.edu>
//
// Copyright 2009 Vanderbilt University - Nashville, TN 37232
//
// 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.
//

#define PWIZ_SOURCE

#include "SpectrumList_Bruker.hpp"
#include "pwiz/utility/chemistry/Chemistry.hpp"

#ifdef PWIZ_READER_BRUKER
#include "pwiz/data/msdata/MSData.hpp"
#include "pwiz/utility/misc/automation_vector.h"
#include "pwiz/utility/misc/Filesystem.hpp"
#include "pwiz/utility/misc/Std.hpp"
#include "pwiz/utility/misc/IntegerSet.hpp"
#include "pwiz/utility/misc/SHA1Calculator.hpp"
#include "pwiz/utility/minimxml/XMLWriter.hpp"
#include <boost/range/algorithm/find_if.hpp>
#include <boost/spirit/include/karma.hpp>

using namespace pwiz::util;
using namespace pwiz::minimxml;
using namespace pwiz::vendor_api::Bruker;

namespace pwiz {
namespace msdata {
namespace detail {

using namespace Bruker;

PWIZ_API_DECL
SpectrumList_Bruker::SpectrumList_Bruker(MSData& msd,
                                         const string& rootpath,
                                         Reader_Bruker_Format format,
                                         CompassDataPtr compassDataPtr,
                                         const Reader::Config& config)
:   msd_(msd), rootpath_(rootpath), format_(format),
    compassDataPtr_(compassDataPtr), config_(config),
    size_(0)
{
    fillSourceList();
    createIndex();
}

PWIZ_API_DECL size_t SpectrumList_Bruker::size() const
{
    return size_;
}

PWIZ_API_DECL const SpectrumIdentity& SpectrumList_Bruker::spectrumIdentity(size_t index) const
{
    if (index >= size_)
        throw runtime_error(("[SpectrumList_Bruker::spectrumIdentity()] Bad index: "
                            + lexical_cast<string>(index)).c_str());
    return index_[index];
}

PWIZ_API_DECL size_t SpectrumList_Bruker::find(const string& id) const
{
    boost::container::flat_map<string, size_t>::const_iterator scanItr = idToIndexMap_.find(id);
    if (scanItr == idToIndexMap_.end())
    {
        if (format_ == Reader_Bruker_Format_TDF)
        {
            try
            {
                int frame = msdata::id::valueAs<int>(id, "frame");
                int scan = msdata::id::valueAs<int>(id, "scan");
                return compassDataPtr_->getSpectrumIndex(frame, scan);
            }
            catch (exception&)
            {
                // fall through and return size_ (id not found)
            }
        }

        return checkNativeIdFindResult(size_, id);
    }
    return scanItr->second;
}

MSSpectrumPtr SpectrumList_Bruker::getMSSpectrumPtr(size_t scan, vendor_api::Bruker::DetailLevel detailLevel) const
{
    if (format_ == Reader_Bruker_Format_FID)
    {
        compassDataPtr_ = CompassData::create(sourcePaths_[scan].string());
        return compassDataPtr_->getMSSpectrum(1, detailLevel);
    }

    return compassDataPtr_->getMSSpectrum(scan, detailLevel);
}

PWIZ_API_DECL SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, bool getBinaryData) const
{
    return spectrum(index, getBinaryData, MSLevelsNone);
}

PWIZ_API_DECL SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, DetailLevel detailLevel) const
{
    return spectrum(index, detailLevel, MSLevelsNone);
}

PWIZ_API_DECL SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, bool getBinaryData, const pwiz::util::IntegerSet& msLevelsToCentroid) const
{
    return spectrum(index, getBinaryData ? DetailLevel_FullData : DetailLevel_FullMetadata, msLevelsToCentroid);
}

PWIZ_API_DECL SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, DetailLevel detailLevel, const pwiz::util::IntegerSet& msLevelsToCentroid) const
{
    if (index >= size_)
        throw runtime_error(("[SpectrumList_Bruker::spectrum()] Bad index: "
                            + lexical_cast<string>(index)).c_str());

    // allocate a new Spectrum
    SpectrumPtr result(new Spectrum);
    if (!result.get())
        throw runtime_error("[SpectrumList_Bruker::spectrum()] Allocation error.");

    vendor_api::Bruker::DetailLevel brukerDetailLevel;
    switch (detailLevel)
    {
        case DetailLevel_InstantMetadata:
        case DetailLevel_FastMetadata:
            brukerDetailLevel = vendor_api::Bruker::DetailLevel_InstantMetadata;
            break;

        case DetailLevel_FullMetadata:
            brukerDetailLevel = vendor_api::Bruker::DetailLevel_FullMetadata;
            break;

        default:
        case DetailLevel_FullData:
            brukerDetailLevel = vendor_api::Bruker::DetailLevel_FullData;
            break;
    }

    const IndexEntry& si = index_[index];
    result->index = si.index;
    result->id = si.id;

    // the scan element may be empty for MALDI spectra, but it's required for a valid file
    result->scanList.scans.push_back(Scan());
    Scan& scan = result->scanList.scans[0];
    scan.instrumentConfigurationPtr = msd_.run.defaultInstrumentConfigurationPtr;

    //try
    //{
        // is this spectrum from the LC interface?
        if (si.collection > -1)
        {
            // fill the spectrum from the LC interface
            LCSpectrumSourcePtr source = compassDataPtr_->getLCSource(si.source);
            LCSpectrumPtr spectrum = compassDataPtr_->getLCSpectrum(si.source, si.scan);

            if (source->getXAxisUnit() == LCUnit_NanoMeter)
                result->set(MS_EMR_spectrum);
            else
                throw runtime_error("[SpectrumList_Bruker::spectrum()] unexpected XAxisUnit");

            double scanTime = spectrum->getTime();
            if (scanTime > 0)
                scan.set(MS_scan_start_time, scanTime, UO_minute);

            result->set(MS_profile_spectrum);

            vector<double> lcX;
            source->getXAxis(lcX);

            if (detailLevel == DetailLevel_FullData)
            {
                result->setMZIntensityArrays(vector<double>(), vector<double>(), MS_number_of_detector_counts);
                result->defaultArrayLength = lcX.size();

                BinaryDataArrayPtr mzArray = result->getMZArray();
                vector<CVParam>::iterator itr = boost::range::find_if(mzArray->cvParams, CVParamIs(MS_m_z_array));
                *itr = CVParam(MS_wavelength_array);

                swap(mzArray->data, lcX);

                vector<double> lcY;
                spectrum->getData(lcY);
                swap(result->getIntensityArray()->data, lcY);
            }
            else
                result->defaultArrayLength = lcX.size();

            return result;
        }

        // get the spectrum from MS interface; for FID formats scan is 0-based, else it's 1-based
        MSSpectrumPtr spectrum = getMSSpectrumPtr(si.scan, brukerDetailLevel);

        result->cvParams.reserve(8); // Anticipate about this many CVParams

        int msLevel = spectrum->getMSMSStage();
        result->set(MS_ms_level, msLevel);

        if (msLevel == 1)
            result->set(MS_MS1_spectrum);
        else
            result->set(MS_MSn_spectrum);

        double scanTime = spectrum->getRetentionTime();
        if (scanTime > 0)
            scan.set(MS_scan_start_time, scanTime, UO_second);

        pair<double, double> scanRange = spectrum->getScanRange();
        if (scanRange.first > 0 && scanRange.second > 0)
            scan.scanWindows.push_back(ScanWindow(scanRange.first, scanRange.second, MS_m_z));

        IonPolarity polarity = spectrum->getPolarity();
        switch (polarity)
        {
            case IonPolarity_Negative:
                result->set(MS_negative_scan);
                break;
            case IonPolarity_Positive:
                result->set(MS_positive_scan);
                break;
            default:
                break;
        }

        //sd.set(MS_base_peak_m_z, pScanStats_->BPM);
        if (spectrum->getTIC() > 0)
        {
            result->set(MS_base_peak_intensity, spectrum->getBPI());
            result->set(MS_total_ion_current, spectrum->getTIC());
        }

        double oneOverK0 = spectrum->oneOverK0();
        if (oneOverK0 > 0)
        {
            scan.set(MS_inverse_reduced_ion_mobility, oneOverK0, MS_Vs_cm_2);
            int windowGroup = spectrum->getWindowGroup();
            if (windowGroup > 0)
                scan.userParams.push_back(UserParam("windowGroup", lexical_cast<string>(windowGroup))); // diaPASEF data
        }

        if (detailLevel == DetailLevel_InstantMetadata)
            return result;

        if (msLevel > 1)
        {
            Precursor precursor;

            vector<double> fragMZs, isolMZs;
            vector<FragmentationMode> fragModes;
            vector<IsolationMode> isolModes;

            spectrum->getFragmentationData(fragMZs, fragModes);

            if (!fragMZs.empty())
            {
                spectrum->getIsolationData(isolMZs, isolModes);

                for (size_t i = 0; i < fragMZs.size(); ++i)
                {
                    if (fragMZs[i] > 0)
                    {
                        SelectedIon selectedIon(fragMZs[i]);

                        int charge = spectrum->getChargeState();
                        if (charge > 0)
                            selectedIon.set(MS_charge_state, charge);

                        switch (fragModes[i])
                        {
                            case FragmentationMode_CID:
                                precursor.activation.set(MS_CID);
                                break;
                            case FragmentationMode_ETD:
                                precursor.activation.set(MS_ETD);
                                break;
                            case FragmentationMode_CIDETD_CID:
                                precursor.activation.set(MS_CID);
                                precursor.activation.set(MS_ETD);
                                break;
                            case FragmentationMode_CIDETD_ETD:
                                precursor.activation.set(MS_CID);
                                precursor.activation.set(MS_ETD);
                                break;
                            case FragmentationMode_ISCID:
                                precursor.activation.set(MS_in_source_collision_induced_dissociation);
                                break;
                            case FragmentationMode_ECD:
                                precursor.activation.set(MS_ECD);
                                break;
                            case FragmentationMode_IRMPD:
                                precursor.activation.set(MS_IRMPD);
                                break;
                            case FragmentationMode_PTR:
                                break;
                        }
                        //precursor.activation.set(MS_collision_energy, pExScanStats_->CollisionEnergy);

                        precursor.selectedIons.push_back(selectedIon);
                    }

                    if (isolMZs[i] > 0)
                    {
                        precursor.isolationWindow.set(MS_isolation_window_target_m_z, isolMZs[i], MS_m_z);

                        double isolationWidth = spectrum->getIsolationWidth();
                        if (isolationWidth > 0)
                        {
                            precursor.isolationWindow.set(MS_isolation_window_lower_offset, isolationWidth / 2, MS_m_z);
                            precursor.isolationWindow.set(MS_isolation_window_upper_offset, isolationWidth / 2, MS_m_z);
                        }
                    }
                }
            }

            if (precursor.selectedIons.size() > 0 || !precursor.isolationWindow.empty())
                result->precursors.push_back(precursor);
        }

        if (detailLevel == DetailLevel_FastMetadata)
            return result;

        // Enumerating merged scan numbers is not instant.
        IntegerSet scanNumbers = spectrum->getMergedScanNumbers();
        if (config_.combineIonMobilitySpectra)
        {
            result->scanList.set(MS_sum_of_spectra);
            if (scanNumbers.size() < 100)
            {
                using namespace boost::spirit::karma;
                auto frameScanPair = compassDataPtr_->getFrameScanPair(si.scan);
                generate(std::back_insert_iterator<std::string>(scan.spectrumID),
                         "frame=" << int_ << " scan=" << int_,
                         frameScanPair.frame, *scanNumbers.begin());

                vector<Scan>& scans = result->scanList.scans;
                scans.reserve(scanNumbers.size());
                for (auto itr = ++scanNumbers.begin(); itr != scanNumbers.end(); ++itr)
                {
                    scans.push_back(Scan());
                    scans.back().instrumentConfigurationPtr = msd_.run.defaultInstrumentConfigurationPtr;
                    generate(std::back_insert_iterator<std::string>(scans.back().spectrumID),
                             "frame=" << int_ << " scan=" << int_,
                             frameScanPair.frame, *itr);

                    // CONSIDER: do we need this? all scan times will be the same and it's rather verbose
                    //if (scanTime > 0)
                    //    scans.back().set(MS_scan_start_time, scanTime, UO_second);

                    // CONSIDER: do we need this? all scan ranges will be the same and it's very verbose!
                    //if (scanRange.first > 0 && scanRange.second > 0)
                    //    scans.back().scanWindows.push_back(ScanWindow(scanRange.first, scanRange.second, MS_m_z));
                }
            }
        }
        else
            result->scanList.set(MS_no_combination);

        //sd.set(MS_lowest_observed_m_z, minObservedMz);
        //sd.set(MS_highest_observed_m_z, maxObservedMz);

        if (detailLevel == DetailLevel_FullData)
        {
            bool getLineData = msLevelsToCentroid.contains(msLevel);

            result->setMZIntensityArrays(vector<double>(), vector<double>(), MS_number_of_detector_counts);

            if (config_.combineIonMobilitySpectra)
            {
                auto& mz = result->getMZArray()->data;
                auto& intensity = result->getIntensityArray()->data;

                result->set(MS_centroid_spectrum); // TIMS is always centroided

                BinaryDataArrayPtr mobility(new BinaryDataArray);
                result->binaryDataArrayPtrs.push_back(mobility);
                CVParam arrayType(MS_mean_inverse_reduced_ion_mobility_array);
                arrayType.units = MS_Vs_cm_2;
                mobility->cvParams.emplace_back(arrayType);

                spectrum->getCombinedSpectrumData(mz, intensity, mobility->data, config_.sortAndJitter);
                result->defaultArrayLength = mz.size();
            }
            else
            {
                automation_vector<double> mzArray, intensityArray;
                if (!getLineData)
                {
                    spectrum->getProfileData(mzArray, intensityArray);
                    if (mzArray.size() == 0)
                        getLineData = true;  // We preferred profile, but there isn't any - try centroided
                    else
                        result->set(MS_profile_spectrum);
                }

                if (getLineData)
                {
                    result->set(MS_centroid_spectrum); // Declare this as centroided data even if scan is empty
                    spectrum->getLineData(mzArray, intensityArray);
                    if (mzArray.size() > 0 && msLevelsToCentroid.contains(msLevel))
                    {
                        result->set(MS_profile_spectrum); // let SpectrumList_PeakPicker know this was probably also a profile spectrum, but doesn't need conversion (actually checking for profile data is crazy slow)
                    }
                }
                result->getMZArray()->data.assign(mzArray.begin(), mzArray.end());
                result->getIntensityArray()->data.assign(intensityArray.begin(), intensityArray.end());
                result->defaultArrayLength = mzArray.size();
            }
        }
        else if (detailLevel == DetailLevel_FullMetadata)
        {
            if (config_.combineIonMobilitySpectra)
            {
                result->defaultArrayLength = spectrum->getCombinedSpectrumDataSize();
                result->set(MS_centroid_spectrum); // TIMS is always centroided
            }
            else
            {
                // N.B.: just getting the data size from the Bruker API is quite expensive.
                if (msLevelsToCentroid.contains(msLevel) || ((result->defaultArrayLength = spectrum->getProfileDataSize())==0))
                {
                    result->defaultArrayLength = spectrum->getLineDataSize();
                    result->set(MS_centroid_spectrum);
                }
                else
                {
                    result->set(MS_profile_spectrum);
                }
            }
        }
    /*}
    catch (_com_error& e) // not caught by either std::exception or '...'
    {
        throw runtime_error(string("[SpectrumList_Bruker::spectrum()] COM error: ") +
                            (const char*)e.Description());
    }*/

    return result;
}

namespace {

void recursivelyEnumerateFIDs(vector<bfs::path>& fidPaths, const bfs::path& rootpath)
{
    const static bfs::directory_iterator endItr = bfs::directory_iterator();

    if (rootpath.leaf() == "fid")
        fidPaths.push_back(rootpath.branch_path());
    else if (bfs::is_directory(rootpath))
    {
        for (bfs::directory_iterator itr(rootpath); itr != endItr; ++itr)
            recursivelyEnumerateFIDs(fidPaths, itr->path());
    }
}

void addSource(MSData& msd, const bfs::path& sourcePath, const bfs::path& rootPath)
{
    SourceFilePtr sourceFile(new SourceFile);
    sourceFile->id = sourcePath.string();
    sourceFile->name = BFS_STRING(sourcePath.leaf());

    // sourcePath: <source>\Analysis.yep|<source>\Analysis.baf|<source>\fid
    // rootPath: c:\path\to\<source>[\Analysis.yep|Analysis.baf|fid]
    bfs::path location = rootPath.has_branch_path() ?
                         BFS_COMPLETE(rootPath.branch_path() / sourcePath) :
                         BFS_COMPLETE(sourcePath); // uses initial path
    sourceFile->location = "file://" + location.branch_path().string();

    msd.fileDescription.sourceFilePtrs.push_back(sourceFile);
}

} // namespace

#if BOOST_FILESYSTEM_VERSION == 2
# define NATIVE_PATH_SLASH "/"
#else
# define NATIVE_PATH_SLASH "\\"
#endif

PWIZ_API_DECL void SpectrumList_Bruker::fillSourceList()
{
    switch (format_)
    {
        case Reader_Bruker_Format_FID:
            recursivelyEnumerateFIDs(sourcePaths_, rootpath_);

            // each fid's source path is a directory but the source file is the fid
            for (size_t i=0; i < sourcePaths_.size(); ++i)
            {
                // in "/foo/bar/1/1SRef/fid", replace "/foo/bar/" with "" so relativePath is "1/1SRef/fid"
                bfs::path relativePath = sourcePaths_[i] / "fid";
                if (rootpath_.has_branch_path())
                    relativePath = bal::replace_first_copy(relativePath.string(), rootpath_.branch_path().string() + NATIVE_PATH_SLASH, "");
                relativePath = bal::replace_all_copy(relativePath.string(), NATIVE_PATH_SLASH, "/");
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_FID_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_FID_format);
            }
            break;

        // a YEP's source path is the same as the source file
        case Reader_Bruker_Format_YEP:
            {
                sourcePaths_.push_back(rootpath_ / "Analysis.yep");
                // strip parent path to get "bar.d/Analysis.yep"
                bfs::path relativePath = bfs::path(rootpath_.filename()) / "Analysis.yep";
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_Agilent_YEP_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_Agilent_YEP_format);
                msd_.run.defaultSourceFilePtr = msd_.fileDescription.sourceFilePtrs.back();
            }
            break;

        // a BAF's source path is the same as the source file
        case Reader_Bruker_Format_BAF:
            {
                sourcePaths_.push_back(rootpath_ / "Analysis.baf");
                // strip parent path to get "bar.d/Analysis.baf"
                bfs::path relativePath = bfs::path(rootpath_.filename()) / "Analysis.baf";
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_BAF_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_BAF_format);
                msd_.run.defaultSourceFilePtr = msd_.fileDescription.sourceFilePtrs.back();
            }
            break;

        // a TDF's source path is the same as the source file
        case Reader_Bruker_Format_TDF:
            {
                sourcePaths_.push_back(rootpath_ / "Analysis.tdf");
                // strip parent path to get "bar.d/Analysis.tdf"
                bfs::path relativePath = bfs::path(rootpath_.filename()) / "Analysis.tdf";
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_TDF_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_TDF_format);
                msd_.run.defaultSourceFilePtr = msd_.fileDescription.sourceFilePtrs.back();
            }

            {
                sourcePaths_.push_back(rootpath_ / "Analysis.tdf_bin");
                // strip parent path to get "bar.d/Analysis.tdf_bin"
                bfs::path relativePath = bfs::path(rootpath_.filename()) / "Analysis.tdf_bin";
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_TDF_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_TDF_format);
            }
            break;

        // a BAF/U2 combo has two sources, with different nativeID formats
        case Reader_Bruker_Format_BAF_and_U2:
            {
                sourcePaths_.push_back(rootpath_ / "Analysis.baf");
                // strip parent path to get "bar.d/Analysis.baf"
                bfs::path relativePath = bfs::path(rootpath_.filename()) / "Analysis.baf";
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_BAF_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_BAF_format);
                msd_.run.defaultSourceFilePtr = msd_.fileDescription.sourceFilePtrs.back();
            }

            {
                sourcePaths_.push_back(bfs::change_extension(rootpath_, ".u2"));
                // in "/foo/bar.d/bar.u2", replace "/foo/" with "" so relativePath is "bar.d/bar.u2"
                bfs::path relativePath = sourcePaths_.back();
                if (rootpath_.has_branch_path())
                    relativePath = bal::replace_first_copy(relativePath.string(), rootpath_.branch_path().string() + NATIVE_PATH_SLASH, "");
                relativePath = bal::replace_all_copy(relativePath.string(), NATIVE_PATH_SLASH, "/");
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_U2_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_U2_format);
            }
            break;

        case Reader_Bruker_Format_U2:
            {
                sourcePaths_.push_back(bfs::change_extension(rootpath_, ".u2"));
                // in "/foo/bar.d/bar.u2", replace "/foo/" with "" so relativePath is "bar.d/bar.u2"
                bfs::path relativePath = sourcePaths_.back();
                if (rootpath_.has_branch_path())
                    relativePath = bal::replace_first_copy(relativePath.string(), rootpath_.branch_path().string() + NATIVE_PATH_SLASH, "");
                relativePath = bal::replace_all_copy(relativePath.string(), NATIVE_PATH_SLASH, "/");
                addSource(msd_, relativePath, rootpath_);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_U2_nativeID_format);
                msd_.fileDescription.sourceFilePtrs.back()->set(MS_Bruker_U2_format);
                msd_.run.defaultSourceFilePtr = msd_.fileDescription.sourceFilePtrs.back();
            }
            break;
    }
}

PWIZ_API_DECL void SpectrumList_Bruker::createIndex()
{
    using namespace boost::spirit::karma;
    map<std::string, size_t> idToIndexTempMap;

    if (format_ == Reader_Bruker_Format_U2 ||
        format_ == Reader_Bruker_Format_BAF_and_U2)
    {
        msd_.fileDescription.fileContent.set(MS_EMR_spectrum);

        for (size_t i=0, end=compassDataPtr_->getLCSourceCount(); i < end; ++i)
        {
            LCSpectrumSourcePtr source = compassDataPtr_->getLCSource(i);

            for (size_t scan=0, end=compassDataPtr_->getLCSpectrumCount(i); scan < end; ++scan)
            {
                index_.push_back(IndexEntry());
                IndexEntry& si = index_.back();
                si.source = i;
                si.collection = source->getCollectionId();
                si.scan = scan;
                si.index = index_.size()-1;
                si.id = "collection=" + lexical_cast<std::string>(si.collection) +
                        " scan=" + lexical_cast<std::string>(si.scan);
                idToIndexTempMap[si.id] = si.index;
            }
        }
    }

    if (format_ == Reader_Bruker_Format_FID)
    {
        int scan = -1;
        BOOST_FOREACH(const SourceFilePtr& sf, msd_.fileDescription.sourceFilePtrs)
        {
            index_.push_back(IndexEntry());
            IndexEntry& si = index_.back();
            si.source = si.collection = -1;
            si.index = index_.size()-1;
            si.scan = ++scan;
            si.id = "file=" + encode_xml_id_copy(sf->id);
            idToIndexTempMap[si.id] = si.index;
        }
    }
    else if (format_ == Reader_Bruker_Format_TDF)
    {
        index_.reserve(compassDataPtr_->getMSSpectrumCount());
        for (size_t scan = 1, end = compassDataPtr_->getMSSpectrumCount(); scan <= end; ++scan)
        {
            index_.emplace_back(IndexEntry());
            IndexEntry& si = index_.back();
            si.source = si.collection = -1;
            si.index = index_.size() - 1;
            si.scan = scan;
            auto frameScanPair = compassDataPtr_->getFrameScanPair(scan);
            std::back_insert_iterator<std::string> sink(si.id);
            if (config_.combineIonMobilitySpectra)
            {
                generate(sink,
                         "merged=" << int_ << " frame=" << int_ << " scanStart=" << int_ << " scanEnd=" << int_,
                         si.index, frameScanPair.frame, frameScanPair.scanStart, frameScanPair.scanEnd);
                idToIndexTempMap[si.id] = si.index;
            }
            else // not inserting into idToIndexTempMap (instead uses on-the-fly logic in find())
                generate(sink,
                         "frame=" << int_ << " scan=" << int_,
                         frameScanPair.frame, frameScanPair.scanStart);
        }
    }
    else if (format_ != Reader_Bruker_Format_U2)
    {
        for (size_t scan=1, end=compassDataPtr_->getMSSpectrumCount(); scan <= end; ++scan)
        {
            index_.push_back(IndexEntry());
            IndexEntry& si = index_.back();
            si.source = si.collection = -1;
            si.index = index_.size()-1;
            si.scan = scan;
            std::back_insert_iterator<std::string> sink(si.id);
            generate(sink,
                     "scan=" << int_,
                     si.scan);
            idToIndexTempMap[si.id] = si.index;
        }
    }

    idToIndexMap_.reserve(idToIndexTempMap.size());
    idToIndexMap_.insert(boost::container::ordered_unique_range, idToIndexTempMap.begin(), idToIndexTempMap.end());
    size_ = index_.size();
}

PWIZ_API_DECL bool SpectrumList_Bruker::hasIonMobility() const
{
    return format_ == Reader_Bruker_Format_TDF;
}

PWIZ_API_DECL bool SpectrumList_Bruker::hasPASEF() const
{
    return compassDataPtr_->hasPASEFData();
}

PWIZ_API_DECL bool SpectrumList_Bruker::canConvertIonMobilityAndCCS() const
{
    return format_ == Reader_Bruker_Format_TDF;
}

// Per email thread Aug 22 2017 bpratt, mattc, Bruker's SvenB:
// The gas is nitrogen(14.0067 AMU) and the temperature is(according to Sven) assumed to be 305K.
static const double ccs_conversion_factor = 18509.863216340458;
static const double MolWeightGas = 14.0067;
static const double Temperature = 305;

PWIZ_API_DECL double SpectrumList_Bruker::ionMobilityToCCS(double inverseK0, double mz, int charge) const
{
    double MolWeight = mz * abs(charge) + chemistry::Electron * charge;
    double ReducedMass = MolWeight * MolWeightGas / (MolWeight + MolWeightGas);
    double K0 = (inverseK0 == 0) ? 0 : (1.0 / inverseK0);
    double ccs = ccs_conversion_factor * abs(charge) / (sqrt(ReducedMass * Temperature) * K0);
    return ccs;    // in Angstrom^2
}

PWIZ_API_DECL double SpectrumList_Bruker::ccsToIonMobility(double ccs, double mz, int charge) const
{
    double MolWeight = mz * abs(charge) + chemistry::Electron * charge;
    double ReducedMass = MolWeight * MolWeightGas / (MolWeight + MolWeightGas);
    double K0 = ccs_conversion_factor * abs(charge) / (sqrt(ReducedMass * Temperature) * ccs);
    return K0 == 0 ? 0 : 1 / K0;    // in Vs/cm^2
}

} // detail
} // msdata
} // pwiz

#else // PWIZ_READER_BRUKER

//
// non-MSVC implementation
//

namespace pwiz {
namespace msdata {
namespace detail {

namespace {const SpectrumIdentity emptyIdentity;}

size_t SpectrumList_Bruker::size() const {return 0;}
const SpectrumIdentity& SpectrumList_Bruker::spectrumIdentity(size_t index) const {return emptyIdentity;}
size_t SpectrumList_Bruker::find(const std::string& id) const {return 0;}
SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, bool getBinaryData) const {return SpectrumPtr();}
SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, DetailLevel detailLevel) const {return SpectrumPtr();}
SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, bool getBinaryData, const pwiz::util::IntegerSet& msLevelsToCentroid) const {return SpectrumPtr();}
SpectrumPtr SpectrumList_Bruker::spectrum(size_t index, DetailLevel detailLevel, const pwiz::util::IntegerSet& msLevelsToCentroid) const {return SpectrumPtr();}
bool SpectrumList_Bruker::hasIonMobility() const { return false; }
bool SpectrumList_Bruker::hasPASEF() const { return false; }
bool SpectrumList_Bruker::canConvertIonMobilityAndCCS() const { return false; }
double SpectrumList_Bruker::ionMobilityToCCS(double ionMobility, double mz, int charge) const {return 0;}
double SpectrumList_Bruker::ccsToIonMobility(double ccs, double mz, int charge) const {return 0;}
} // detail
} // msdata
} // pwiz

#endif // PWIZ_READER_BRUKER

ext: cpp

max_stars_repo_name: rfellers/pwiz

max_issues_repo_name: rfellers/pwiz

max_forks_repo_name: rfellers/pwiz


text:
#include <iostream>
#include <boost/coroutine2/all.hpp>

typedef boost::coroutines2::coroutine<int>::pull_type pull_type;
typedef boost::coroutines2::coroutine<int>::push_type push_type;

void log(const std::string &s) {
    std::cout << "n = " << s << std::endl;
}

void log(int n) {
    log(std::to_string(n));
}

void f(push_type &yield) {
    log(1);
    yield(3);
    log(5);
}

int main() {
    log(0);
    pull_type source(f);
    log(2);
    log(source.get());
    log(4);
    source();
}

max_stars_repo_name: dkandalov/katas

max_issues_repo_name: dkandalov/katas

max_forks_repo_name: dkandalov/katas


text:
// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// FlexibleSUSY is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================

#ifndef MATHLINK_UTILS_H
#define MATHLINK_UTILS_H

#include "mathlink.h"

#include <complex>
#include <string>
#include <vector>
#include <Eigen/Core>

/********************* put types *********************/

inline void MLPut(MLINK link, const std::string& s)
{
   MLPutSymbol(link, s.c_str());
}

inline void MLPut(MLINK link, int c)
{
   MLPutInteger(link, c);
}

inline void MLPut(MLINK link, double c)
{
   MLPutReal(link, c);
}

inline void MLPut(MLINK link, std::complex<double> c)
{
   if (std::imag(c) == 0.) {
      MLPutReal(link, std::real(c));
   } else {
      MLPutFunction(link, "Complex", 2);
      MLPutReal(link, std::real(c));
      MLPutReal(link, std::imag(c));
   }
}

template <int M>
void MLPut(MLINK link, const Eigen::Array<double,M,1>& a)
{
   double v[M];
   for (int i = 0; i < M; i++)
      v[i] = a(i);
   MLPutRealList(link, v, M);
}

template <int M>
void MLPut(MLINK link, const Eigen::Matrix<double,M,1>& m)
{
   const Eigen::Array<double,M,1> a(m.array());
   MLPut(link, a);
}

template <int M, int N>
void MLPut(MLINK link, const Eigen::Matrix<double,M,N>& m)
{
   double mat[M][N];
   for (int i = 0; i < M; i++)
      for (int k = 0; k < N; k++)
         mat[i][k] = m(i, k);

   long dims[] = { M, N };
   MLPutDoubleArray(link, (double*)mat, dims, NULL, 2);
}

template <int M>
void MLPut(MLINK link, const Eigen::Array<std::complex<double>,M,1>& a)
{
   MLPutFunction(link, "List", M);
   for (int i = 0; i < M; i++)
      MLPut(link, a(i));
}

template <int M>
void MLPut(MLINK link, const Eigen::Matrix<std::complex<double>,M,1>& m)
{
   const Eigen::Array<std::complex<double>,M,1> a(m.array());
   MLPut(link, a);
}

template <int M, int N>
void MLPut(MLINK link, const Eigen::Matrix<std::complex<double>,M,N>& m)
{
   MLPutFunction(link, "List", M);
   for (int i = 0; i < M; i++) {
      MLPutFunction(link, "List", N);
      for (int k = 0; k < N; k++)
         MLPut(link, m(i,k));
   }
}

/********************* put single heads *********************/

inline void MLPutHeads(MLINK link, const std::vector<std::string>& heads)
{
   for (const auto& h: heads)
      MLPutFunction(link, h.c_str(), 1);
}

/********************* put rules to types *********************/

inline void MLPutRule(MLINK link, const std::string& name, const std::vector<std::string>& heads = {})
{
   MLPutFunction(link, "Rule", 2);
   MLPutHeads(link, heads);
   MLPutUTF8Symbol(link, reinterpret_cast<const unsigned char*>(name.c_str()), name.size());
}

inline void MLPutRule(MLINK link, int number, const std::vector<std::string>& heads = {})
{
   MLPutFunction(link, "Rule", 2);
   MLPutHeads(link, heads);
   MLPutInteger(link, number);
}

inline void MLPutRule(MLINK link, long number, const std::vector<std::string>& heads = {})
{
   MLPutFunction(link, "Rule", 2);
   MLPutHeads(link, heads);
   MLPutLongInteger(link, number);
}

template <class T1, class T2>
void MLPutRuleTo(MLINK link, T1 t, const T2& name, const std::vector<std::string>& heads = {})
{
   MLPutRule(link, name, heads);
   MLPut(link, t);
}

/********************* get types *********************/

inline void MLGet(MLINK link, int *c)
{
   MLGetInteger(link, c);
}

inline void MLGet(MLINK link, long *c)
{
   MLGetLongInteger(link, c);
}

inline void MLGet(MLINK link, short *c)
{
   MLGetShortInteger(link, c);
}

#endif // MATHLINK_UTILS_H


text:
// Boost.Geometry
// Unit Test

// Copyright (c) 2017, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_TEST_SRS_CHECK_GEOMETRY_HPP
#define BOOST_GEOMETRY_TEST_SRS_CHECK_GEOMETRY_HPP

#include <geometry_test_common.hpp>

#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/tag.hpp>
#include <boost/geometry/core/tags.hpp>

#include <boost/geometry/io/wkt/read.hpp>

#include <boost/geometry/views/detail/indexed_point_view.hpp>

#include <boost/range/begin.hpp>
#include <boost/range/end.hpp>
#include <boost/range/size.hpp>
#include <boost/range/value_type.hpp>

namespace test
{

struct check_point
{
    template <typename Point, typename T>
    static void apply(Point const& point1, Point const& point2, T tol)
    {
        typename bg::coordinate_type<Point>::type
            x1 = bg::get<0>(point1),
            y1 = bg::get<1>(point1),
            x2 = bg::get<0>(point2),
            y2 = bg::get<1>(point2);

        BOOST_CHECK_CLOSE(x1, x2, tol);
        BOOST_CHECK_CLOSE(y1, y2, tol);
    }
};

template <typename Policy = check_point>
struct check_range
{
    template <typename Range, typename T>
    static void apply(Range const& range1, Range const& range2, T tol)
    {
        size_t range1_count = boost::size(range1);
        size_t range2_count = boost::size(range2);
        BOOST_CHECK_EQUAL(range1_count, range2_count);
        if (range1_count == range2_count)
        {
            apply(boost::begin(range1), boost::end(range1),
                  boost::begin(range2), tol);
        }
    }
    template <typename It, typename T>
    static void apply(It first1, It last1, It first2, T tol)
    {
        for ( ; first1 != last1 ; ++first1, ++first2)
            Policy::apply(*first1, *first2, tol);
    }
};

template <typename Geometry, typename Tag = typename bg::tag<Geometry>::type>
struct check_geometry_impl
{};

template <typename Point>
struct check_geometry_impl<Point, bg::point_tag>
    : check_point
{};

template <typename Segment>
struct check_geometry_impl<Segment, bg::segment_tag>
{
    template <typename T>
    static void apply(Segment const& g1, Segment const& g2, T tol)
    {
        bg::detail::indexed_point_view<Segment const, 0> p1(g1);
        bg::detail::indexed_point_view<Segment const, 1> p2(g1);
        bg::detail::indexed_point_view<Segment const, 0> q1(g2);
        bg::detail::indexed_point_view<Segment const, 1> q2(g2);

        check_point::apply(p1, q1, tol);
        check_point::apply(p2, q2, tol);
    }
};

template <typename MultiPoint>
struct check_geometry_impl<MultiPoint, bg::multi_point_tag>
    : check_range<>
{};

template <typename Linestring>
struct check_geometry_impl<Linestring, bg::linestring_tag>
    : check_range<>
{};

template <typename MultiLinestring>
struct check_geometry_impl<MultiLinestring, bg::multi_linestring_tag>
    : check_range< check_range<> >
{};

template <typename Ring>
struct check_geometry_impl<Ring, bg::ring_tag>
    : check_range<>
{};

template <typename Polygon>
struct check_geometry_impl<Polygon, bg::polygon_tag>
{
    template <typename T>
    static void apply(Polygon const& g1, Polygon const& g2, T tol)
    {
        check_range<>::apply(bg::exterior_ring(g1), bg::exterior_ring(g2), tol);
        check_range< check_range<> >::apply(bg::interior_rings(g1), bg::interior_rings(g2), tol);
    }
};

template <typename MultiPolygon>
struct check_geometry_impl<MultiPolygon, bg::multi_polygon_tag>
    : check_range
        <
            check_geometry_impl
                <
                    typename boost::range_value<MultiPolygon>::type,
                    bg::polygon_tag
                >
        >
{};

template <typename Geometry, typename T>
inline void check_geometry(Geometry const& g1, Geometry const& g2, T tol)
{
    check_geometry_impl<Geometry>::apply(g1, g2, tol);
}

template <typename Geometry, typename T>
inline void check_geometry(Geometry const& g1, std::string const& wkt2, T tol)
{
    Geometry g2;
    bg::read_wkt(wkt2, g2);
    check_geometry_impl<Geometry>::apply(g1, g2, tol);
}

} // namespace test

#endif // BOOST_GEOMETRY_TEST_SRS_CHECK_GEOMETRY_HPP


text:
#define BOOST_TEST_DYN_LINK
#include <boost/test/auto_unit_test.hpp>
#include <boost/test/unit_test_suite.hpp>
#include <boost/test/unit_test.hpp>

#include <luacxx/binding/class.hpp>
#include <luacxx/binding/function.hpp>
#include <luacxx/binding/engine.hpp>

#include <config.h>

namespace
{
  class Cls1;
  class ClassA;
  class ClassB;
  class ClassC;
  class ClassD;
  class ClassE;

  Cls1*   last_cls1;
  ClassA* last_clsA;
  ClassB* last_clsB;
  ClassC* last_clsC;
  ClassD* last_clsD;
  ClassE* last_clsE;

  int counter_A;
  int counter_B;
  int counter_C;
  int counter_D;
  int counter_E;

  void init()
  {
    last_cls1 = nullptr;
    last_clsA = nullptr;
    last_clsB = nullptr;
    last_clsC = nullptr;
    last_clsD = nullptr;
    last_clsE = nullptr;
    counter_A = 0;
    counter_B = 0;
    counter_C = 0;
    counter_D = 0;
    counter_E = 0;
  }

  class Cls1
  {
    public:
      Cls1()
      {
        last_cls1 = this;
      }

      bool meth1(int a, const std::string& b)
      {
        return (a == 5 && b == "blabla" && this == last_cls1);
      }

      bool meth2(int a, const std::string& b) const
      {
        return (a == 7 && b == "blabla" && this == last_cls1);
      }
  };

  bool extern_meth1(Cls1* cls, int a, const std::string& b)
  {
    return (a ==  9 && b == "blabla"  && cls == last_cls1);
  }

  class Cls2
  {
    public:
      Cls2()
      {
      }

      static bool mth1(const std::string&a, const std::string& b)
      {
        return a == b;
      }
  };

  class ClassA
  {
    public:
      ClassA()
      {
        last_clsA = this;
      }

      int methdA()
      {
        return ((last_clsA == this) ?  counter_A++ : -1);
      }
  };

  class ClassB : public ClassA
  {
    public:
      ClassB()
      {
        last_clsB = this;
      }

      int methdB()
      {
        return ((last_clsB == this) ?  counter_B++ : -1);
      }
  };

  class ClassC: public ClassB
  {
    public:
      ClassC()
      {
        last_clsC = this;
      }

      int methdC()
      {
        return ((last_clsC == this) ?  counter_C++ : -1);
      }
  };

  class ClassD
  {
    public:
      ClassD()
      {
        last_clsD = this;
      }

      int methdD()
      {
        return ((last_clsD == this) ?  counter_D++ : -1);
      }
  };

  class ClassE : public ClassD, public ClassC
  {
    public:
      ClassE()
      {
        last_clsE = this;
      }

      int methdE()
      {
        return ((last_clsE == this) ?  counter_E++ : -1);
      }
  };
}

BOOST_AUTO_TEST_CASE( method_class_01 )
{
  init();
  {
    luacxx::engine engine;
    luacxx::make_class<Cls1>(engine, "Cls1")
        .ctor()
        .method("meth1", &Cls1::meth1)
        .const_method("meth2", &Cls1::meth2)
        .extern_method("meth3", extern_meth1)
        .extern_method("meth3", extern_meth1)
        .extern_method("meth4", std::function< bool (Cls1*, int, const std::string&)>(extern_meth1));
    BOOST_CHECK(engine.bind());

    std::pair<bool, std::string> ret = engine.do_file(LUACXX_TESTS_DIR "method_class_01.lua");
    if(!ret.first)
    {
      BOOST_FAIL(ret.second);
    }
  }
}

BOOST_AUTO_TEST_CASE( method_class_02 )
{
  init();
  {
    luacxx::engine engine;
    luacxx::make_class<ClassA>(engine, "ClassA")
        .ctor()
        .method("methdA", &ClassA::methdA);
    luacxx::make_class<ClassB>(engine, "ClassB")
        .ctor()
        .inheritance<ClassA>()
        .method("methdB", &ClassB::methdB);
    luacxx::make_class<ClassC>(engine, "ClassC")
        .ctor()
        .inheritance<ClassB>()
        .method("methdC", &ClassC::methdC);
    BOOST_CHECK(engine.bind());

    std::pair<bool, std::string> ret = engine.do_file(LUACXX_TESTS_DIR "method_class_02.lua");
    if(!ret.first)
    {
      BOOST_FAIL(ret.second);
    }
    else
    {
      BOOST_CHECK_EQUAL(counter_A, 3);
      BOOST_CHECK_EQUAL(counter_B, 2);
      BOOST_CHECK_EQUAL(counter_C, 1);
    }
  }
}

BOOST_AUTO_TEST_CASE( method_class_03 )
{
  init();
  {
    luacxx::engine engine;
    luacxx::make_class<ClassA>(engine, "ClassA")
        .ctor()
        .method("methdA", &ClassA::methdA);
    luacxx::make_class<ClassB>(engine, "ClassB")
        .ctor()
        .inheritance<ClassA>()
        .method("methdB", &ClassB::methdB);
    luacxx::make_class<ClassC>(engine, "ClassC")
        .ctor()
        .inheritance<ClassB>()
        .method("methdC", &ClassC::methdC);
    luacxx::make_class<ClassD>(engine, "ClassD")
        .ctor()
        .method("methdD", &ClassD::methdD);
    luacxx::make_class<ClassE>(engine, "ClassE")
        .ctor()
        .inheritance<ClassD>()
        .inheritance<ClassC>()
        .method("methdE", &ClassE::methdE);
    BOOST_CHECK(engine.bind());

    std::pair<bool, std::string> ret = engine.do_file(LUACXX_TESTS_DIR "method_class_03.lua");
    if(!ret.first)
    {
      BOOST_FAIL(ret.second);
    }
    else
    {
      BOOST_CHECK_EQUAL(counter_A, 1);
      BOOST_CHECK_EQUAL(counter_B, 1);
      BOOST_CHECK_EQUAL(counter_C, 1);
      BOOST_CHECK_EQUAL(counter_D, 1);
      BOOST_CHECK_EQUAL(counter_E, 1);
    }
  }
}

BOOST_AUTO_TEST_CASE( method_class_04 )
{
  init();
  {
    luacxx::engine engine;
    luacxx::make_class<Cls2>(engine, "Cls2")
        .ctor()
        .static_method("mth1", &Cls2::mth1);
    BOOST_CHECK(engine.bind());

    std::pair<bool, std::string> ret = engine.do_file(LUACXX_TESTS_DIR "method_class_04.lua");
    if(!ret.first)
    {
      BOOST_FAIL(ret.second);
    }
  }
}

max_stars_repo_name: EVaillant/luacxx

max_issues_repo_name: EVaillant/luacxx

max_forks_repo_name: EVaillant/luacxx