text:
// Copyright (c) 2016 Matt Overby
//
// MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause)
// Redistribution and use in source and binary forms, with or without modification, are
// permitted provided that the following conditions are met:
// 1. Redistributions of source code must retain the above copyright notice, this list of
// conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list
// of conditions and the following disclaimer in the documentation and/or other materials
// provided with the distribution.
// THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR  A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
// OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
// IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef INTERP_HPP
#define INTERP_HPP 1

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

namespace mpm
{

	// Cubic B-spline
	static inline double cspline(double x)
	{
		if (x==0.0) { x=1e-12; }
		x = fabs(x);
		if (x < 1.0) { return x*x*(x*0.5 - 1.0) + 2.0/3.0; }
		else if (x < 2.0) { return x*(x*(-x/6.0 + 1.0) - 2.0) + 4.0/3.0; }
		return 0.0;
	}

	// Slope of cubic spline
	static inline double d_cspline(double x)
	{
		if (x==0.0) { x=1e-12; }
		double abs_x = fabs(x);
		if (abs_x < 1.0) { return 1.5*x*abs_x - 2.0*x; }
		else if (x < 2.0) { return -x*abs_x*0.5 + 2.0*x - 2.0*x/abs_x; }
		return 0.0;
	}

	// Quadratic spline
	static inline double qspline( double x )
	{
		double fx = fabs(x);
		if (fx < 0.5) { return ( 0.75 - x*x ); }
		else if (fx < 1.5) { return ( 0.5*x*x - 1.5*fx + 9.0/8.0 ); }
		return 0.0;
	}

	// Slope of quadratic spline
	static inline double d_qspline( double x )
	{
		double fx = fabs(x);
		if (fx < 0.5) { return ( -2.0 * x ); }
		else if (fx < 1.5)
		{
			if (x < 0.0) { return fabs( fx - 1.5 ); }
			else { return (fx - 1.5); }
		}
		return 0.0;
	}

	static inline bool isreal(double n)
	{
		return !std::isnan(n) && !std::isinf(n);
	}

} // end namespace mpm

#endif

ext: hpp

max_stars_repo_name: NTForked/mpm-optimization

max_issues_repo_name: NTForked/mpm-optimization

max_forks_repo_name: NTForked/mpm-optimization


text:
///////////////////////////////////////////////////////////////////
//  Copyright Eduardo Quintana 2021
//  Copyright Janek Kozicki 2021
//  Copyright Christopher Kormanyos 2021
//  Distributed under the Boost Software License,
//  Version 1.0. (See accompanying file LICENSE_1_0.txt
//  or copy at http://www.boost.org/LICENSE_1_0.txt)

// This file provides two traits:
//  boost::multiprecision::complex<Real>                 - creates a coplex type supported by Boost
//  boost::multiprecision::is_boost_complex<T>::value    - a bool telling whether the type is a complex number supported by Boost
//                                                         it has nothing to do with boost::is_complex<T> which recognizes only std::complex<T>
//
// TODO: this file should be split apart and each specializetion should go to the end of the relevant file in Boost.Multiprecision
//   --> specialization for float128      --> complex128    goes to <boost/multiprecision/complex128.hpp>
//   --> specialization for cpp_bin_float --> cpp_complex   goes to <boost/multiprecision/cpp_complex.hpp>
//   --> specialization for mpfr_float    --> mpc_complex   goes to <boost/multiprecision/mpc.hpp>

#ifndef BOOST_MATH_FFT_MULTIPRECISION_COMPLEX
  #define BOOST_MATH_FFT_MULTIPRECISION_COMPLEX
  #include <boost/config.hpp>
  #include <boost/cstdfloat.hpp>
  #if defined(BOOST_MATH_USE_FLOAT128)
  #include <boost/multiprecision/float128.hpp>
  #endif
  #include <boost/multiprecision/cpp_bin_float.hpp>
  #include <boost/multiprecision/cpp_dec_float.hpp>
  #if defined(__GNUC__)
  //#include <boost/multiprecision/mpfr.hpp>
  #endif
  #include <boost/core/demangle.hpp>
  #include <boost/core/enable_if.hpp>
  // now include complex types
  #if defined(BOOST_MATH_USE_FLOAT128)
  #include <boost/multiprecision/complex128.hpp>
  #endif
  #include <boost/multiprecision/cpp_complex.hpp>
  #if defined(__GNUC__)
  //#include <boost/multiprecision/mpc.hpp>
  #endif

  namespace boost { namespace multiprecision { namespace detail {

    template<typename T>
    struct make_boost_complex {
      using type = std::complex<T>;
    };

    // float128 --> <boost/multiprecision/complex128.hpp>
    #if defined(BOOST_MATH_USE_FLOAT128)
    template<>
    struct make_boost_complex<boost::multiprecision::float128> {
      using type = boost::multiprecision::complex128;
    };
    #endif

    // cpp_bin_float --> <boost/multiprecision/cpp_complex.hpp>
    template <unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent, expression_template_option ExpressionTemplates>
    struct make_boost_complex< number< cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>, ExpressionTemplates> > {
      using type = number<complex_adaptor<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>>, ExpressionTemplates>;
    };

    #if defined(__GNUC__)
    // mpc_complex --> <boost/multiprecision/mpc.hpp>
    //template <unsigned Digits, mpfr_allocation_type AllocationType, expression_template_option ExpressionTemplates>
    //struct make_boost_complex< number< backends::mpfr_float_backend<Digits, AllocationType>, ExpressionTemplates> > {
    //  using type = number<mpc_complex_backend<Digits>, ExpressionTemplates>;
    //};
    #endif

    #if __cplusplus < 201700L
    template<typename... Ts> struct make_void { typedef void type;};
    template<typename... Ts> using void_t = typename make_void<Ts...>::type;
    #else
    using std::void_t;
    #endif

    template< class, class = void > struct has_value_type : std::false_type { };
    template< class T > struct has_value_type<T, void_t<typename T::value_type>> : std::true_type { };

    template <typename T, bool = has_value_type<T>::value > struct is_boost_complex {
      static constexpr bool value = false;
    };

    template <typename T> struct is_boost_complex<T, true> {
      static constexpr bool value = std::is_same<
                                      typename make_boost_complex<typename T::value_type>::type
                                    , typename std::decay<T>::type
                                    >::value;
    };

  }  // boost::multiprecision::detail
  template<typename T> using complex          = typename detail::make_boost_complex<T>::type;
  template<typename T> using is_boost_complex = detail::is_boost_complex<T>; // may be a duplicate of https://github.com/BoostGSoC21/math/pull/8#discussion_r660549430

  }} // boost::multiprecision

#endif


text:
/*==============================================================================
  Copyright (c) 2017, 2018 Matt Calabrese

  Distributed under the Boost Software License, Version 1.0. (See accompanying
  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/

#include "future_like_testing.hpp"

#include <argot/concepts/future.hpp>
#include <argot/concepts/same_type.hpp>
#include <argot/executor/immediate.hpp>
#include <argot/gen/concept_ensure.hpp>
#include <argot/fut_traits/destructive_then.hpp>
#include <argot/fut_traits/persistent_then.hpp>
#include <argot/fut_traits/value_type.hpp>
#include <argot/packager/boost_future.hpp>
#include <argot/void_.hpp>

#include <argot/gen/concept_ensure.hpp>

#include <boost/thread/barrier.hpp>
#include <boost/thread/future.hpp>

#include <thread>
#include <type_traits>
#include <utility>

#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>

namespace {

namespace fut_traits = argot::fut_traits;
namespace executor = argot::executor;
namespace packager = argot::packager;

using executor::immediate;
using argot::Future;
using argot::SameType;
using argot::void_;
using argot::void_to_regular_void_t;
using fut_traits::destructive_then;
using fut_traits::value_type_t;

using argot_future_like_testing::wrapped;
using argot_future_like_testing::wrap;

using packager_t = packager::boost_future;

template< class T >
using hesitate = argot_future_like_testing::hesitate< boost::promise, T >;

template< class T >
using promise = boost::promise< T >;

template< class T >
using future = boost::future< T >;

template< class T >
using then_result = boost::future< T >;

ARGOT_CONCEPT_ENSURE( Future< future< int > > );
ARGOT_CONCEPT_ENSURE( Future< future< void > > );
ARGOT_CONCEPT_ENSURE( Future< future< int& > > );

ARGOT_CONCEPT_ENSURE( SameType< value_type_t< future< int > >, int > );
ARGOT_CONCEPT_ENSURE( SameType< value_type_t< future< void > >, void_ > );
ARGOT_CONCEPT_ENSURE( SameType< value_type_t< future< int& > >, int& > );

BOOST_AUTO_TEST_CASE( test_destructive_then_object_single_threaded )
{
  promise< int > prom;
  future< int > fut = prom.get_future();

  prom.set_value( 5 );

  decltype( auto ) result
    = destructive_then< packager_t >
      ( std::move( fut ), immediate, wrap< int >() );

  ARGOT_CONCEPT_ENSURE
  ( SameType< decltype( result ), then_result< wrapped< int > > > );

  BOOST_TEST( result.get().value() == 5 );
}

BOOST_AUTO_TEST_CASE( test_destructive_then_void_single_threaded )
{
  promise< void > prom;
  future< void > fut = prom.get_future();

  prom.set_value();

  decltype( auto ) result
    = destructive_then< packager_t >
    ( std::move( fut ), immediate, wrap< void_ >() );

  ARGOT_CONCEPT_ENSURE
  ( SameType< decltype( result ), then_result< wrapped< void_ > > > );

  result.get();
}

BOOST_AUTO_TEST_CASE( test_destructive_then_reference_single_threaded )
{
  promise< int& > prom;
  future< int& > fut = prom.get_future();
  int val = 5;

  prom.set_value( val );

  decltype( auto ) result
    = destructive_then< packager_t >
    ( std::move( fut ), immediate, wrap< int& >() );

  ARGOT_CONCEPT_ENSURE
  ( SameType< decltype( result ), then_result< wrapped< int& > > > );

  BOOST_TEST( &result.get().value() == &val );
}

BOOST_AUTO_TEST_CASE( test_destructive_then_object_multi_threaded )
{
  promise< int > prom;
  future< int > fut = prom.get_future();
  boost::barrier wall( 2 );

  std::thread other_thread( hesitate< int >{ std::move( prom ), &wall, 5 } );

  decltype( auto ) result
    = destructive_then< packager_t >
    ( std::move( fut ), immediate, wrap< int >() );

  wall.wait();

  BOOST_TEST( result.get().value() == 5 );

  other_thread.join();
}

BOOST_AUTO_TEST_CASE( test_destructive_then_void_multi_threaded )
{
  promise< void > prom;
  future< void > fut = prom.get_future();
  boost::barrier wall( 2 );

  std::thread other_thread( hesitate< void >{ std::move( prom ), &wall, {} } );

  decltype( auto ) result
    = destructive_then< packager_t >
      ( std::move( fut ), immediate, wrap< void_ >() );

  wall.wait();

  result.get();

  other_thread.join();
}

BOOST_AUTO_TEST_CASE( test_destructive_then_reference_multi_threaded )
{
  promise< int& > prom;
  future< int& > fut = prom.get_future();
  int val = 5;
  boost::barrier wall( 2 );

  std::thread other_thread( hesitate< int& >{ std::move( prom ), &wall, val } );

  decltype( auto ) result
    = destructive_then< packager_t >
      ( std::move( fut ), immediate, wrap< int& >() );

  wall.wait();

  BOOST_TEST( &result.get().value() == &val );

  other_thread.join();
}

}  // namespace

ext: cpp

max_stars_repo_name: mattcalabrese/argot

max_issues_repo_name: mattcalabrese/argot

max_forks_repo_name: mattcalabrese/argot


text:
// All content Copyright (C) 2018 Genomics plc
#include "caller/job.hpp"
#include "caller/params.hpp"
#include "caller/annotation.hpp"
#include "caller/diploid/diploid.hpp"
#include "caller/regionUtils.hpp"
#include "caller/diploid/referenceCalling.hpp"
#include "caller/alignPhasing.hpp"

#include "variant/variantGenerator.hpp"
#include "variant/variantFilter.hpp"
#include "variant/variantContainer.hpp"
#include "variant/type/variant.hpp"
#include "variant/clustering.hpp"
#include "variant/haplotypeGenerator.hpp"
#include "variant/type/variant.hpp"
#include "variant/breakpointVariantGenerator.hpp"

#include "io/bamFile.hpp"
#include "io/readRange.hpp"
#include "vcf/reader.hpp"

#include "readrecalibration/readRecalibration.hpp"

#include "utils/exceptions.hpp"
#include "utils/logging.hpp"
#include "utils/flatten.hpp"

#include <algorithm>
#include <limits>
#include <boost/filesystem/operations.hpp>

using namespace wecall::io;
using wecall::variant::varPtr_t;

namespace wecall
{
namespace caller
{
    //-----------------------------------------------------------------------------------------

    Job::Job( caller::params::Application applicationParams,
              caller::params::Data dataParams,
              caller::params::System systemParams,
              caller::params::PrivateSystem privateSystemParams,
              caller::params::Filters filterParams,
              caller::params::PrivateCalling privateCallingParams,
              caller::params::Calling callingParams,
              caller::params::PrivateData privateDataParams )
        : m_applicationParams( applicationParams ),
          m_dataParams( dataParams ),
          m_systemParams( systemParams ),
          m_privateSystemParams( privateSystemParams ),
          m_filterParams( filterParams ),
          m_privateCallingParams( privateCallingParams ),
          m_callingParams( callingParams ),
          m_privateDataParams( privateDataParams ),
          m_candidateVCFFile(
              privateDataParams.m_candidateVariantsFile.empty()
                  ? nullptr
                  : std::make_shared< io::TabixVCFFile >( privateDataParams.m_candidateVariantsFile,
                                                          privateDataParams.m_candidateVariantsFile + ".tbi" ) ),
          m_genotypeAllelesVCFFile(
              privateDataParams.genotypingMode()
                  ? std::make_shared< io::TabixVCFFile >( privateDataParams.genotypeAllelesFile(),
                                                          privateDataParams.genotypeAllelesFile() + ".tbi" )
                  : nullptr ),
          m_intermediateOutputWriter( dataParams.inputDataSources(), privateDataParams.m_intermediateRecalibFileStem ),
          m_variantSoftFilterBank( privateCallingParams.m_varFilterIDs,
                                   callingParams.m_minAlleleBiasP,
                                   callingParams.m_minStrandBiasP,
                                   callingParams.m_minAllelePlusStrandBiasP,
                                   callingParams.m_minRMSMappingQ,
                                   callingParams.m_minSNPQOverDepth,
                                   callingParams.m_minIndelQOverDepth,
                                   callingParams.m_minBadReadsScore,
                                   callingParams.m_minCallQual ),
          m_readDataReader( systemParams,
                            privateSystemParams,
                            filterParams,
                            privateSystemParams.m_biteSize,
                            dataParams.inputDataSources() ),
          m_vcOut( dataParams.outputDataSink(), dataParams.outputRefCalls(), callingParams.m_outputPhasedGenotypes ),
          m_ref( dataParams.refFile() ),
          // TODO(ES): Tie together contig, calling and output regions together into nice container.
          m_outputRegions( utils::functional::flatten( dataParams.dataRegions() ) ),
          m_callingRegions( utils::functional::flatten(
              caller::padPartitionedRegions( dataParams.dataRegions(), m_privateCallingParams.m_regionPadding ) ) ),
          m_model( callingParams.m_badReadsWindowSize,
                   privateCallingParams.m_maxHaplotypesPerCluster,
                   m_readDataReader.getSampleNames() ),
          m_variantCallBuilder( dataParams.outputRefCalls(),
                                callingParams.m_outputPhasedGenotypes,
                                privateCallingParams.m_allVariants,
                                privateDataParams.genotypingMode(),
                                m_readDataReader.getSampleNames().size() )
    {
    }

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

    void Job::process()
    {
        WECALL_LOG( DEBUG, "Writing VCF header" );

        std::vector< vcf::FilterDesc > filterDescs = m_variantSoftFilterBank.getFilterDescs();
        filterDescs.emplace_back( vcf::filter::NC_key,
                                  "Not called: Indicates a variant that was not positively genotyped in any sample." );

        m_vcOut.writeHeader( m_dataParams.outputFormat(), m_applicationParams, m_dataParams.refFile(),
                             m_readDataReader.getSampleNames(), filterDescs,
                             caller::getContigsFromRegions( m_outputRegions, m_ref.indexFile().contigs() ) );

        // Currently defaulting to same ploidy all samples.
        const std::vector< std::size_t > ploidyPerSample( m_readDataReader.getSampleNames().size(),
                                                          m_privateCallingParams.m_ploidy );

        for ( const Region & region : m_callingRegions )
        {
            const auto contig = region.contig();
            const auto refContigs = m_ref.indexFile().contigs();
            const auto contigLookup = refContigs.find( contig );

            if ( contigLookup == refContigs.end() )
            {
                continue;
            }

            const auto contigRegion = Region( contig, contigLookup->second );

            if ( not contigRegion.overlaps( region ) )
            {
                continue;
            }

            const auto callingRegion = contigRegion.getIntersect( region );

            if ( callingRegion.size() == 0 )
            {
                continue;
            }

            m_vcOut.contig( contig );

            // Get required reference.
            const auto maxReadLength =
                1000;  // TODO made sure Read construction doesn't fail it it goes byond this limit
            const auto assemblePadding = 3 * params::defaults::maxBreakpointKmerSize;
            const caller::Region paddedRefRegion = callingRegion.getPadded( maxReadLength + assemblePadding );
            m_ref.cacheSequence( paddedRefRegion );
            const auto referenceSequence =
                std::make_shared< utils::ReferenceSequence >( m_ref.getSequence( paddedRefRegion ) );

            auto blockIterator = m_readDataReader.readRegion( callingRegion, referenceSequence );

            while ( auto readDataset = blockIterator.getReadDatasetForNextBlock() )
            {
                const auto block = readDataset->region();

                if ( not readDataset->isEmpty() or m_dataParams.outputRefCalls() or
                     m_privateDataParams.genotypingMode() )
                {
                    const auto actualEnd = processBlock( readDataset, blockIterator.isLastBlock(), ploidyPerSample );

                    if ( actualEnd < block.end() )
                    {
                        blockIterator.chopCurrentBlock( actualEnd );
                    }
                }
                else
                {
                    WECALL_LOG( INFO, "Skipping " << block << " due to no read data." );
                }
            }
        }

        WECALL_LOG( INFO, "Job completed successfully." );
    }

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

    std::vector< variant::VariantCluster > Job::generateVariantClustersInBlock(
        const caller::Region & blockRegion,
        io::perSampleRegionsReads_t allReads,
        utils::referenceSequencePtr_t referenceSequence )
    {
        const variant::VariantGenerator varGen( referenceSequence, m_privateCallingParams.m_minBaseQual,
                                                m_filterParams.m_readMappingFilterQ );

        auto varContainer = varGen.generateVariantsFromReads( allReads );

        if ( m_candidateVCFFile != nullptr )
        {
            for ( const auto & record : m_candidateVCFFile->fetch( blockRegion ) )
            {
                const auto variants = record.getVariants( referenceSequence );
                const auto priors = caller::getPriorsFromInfo( record.m_info );
                for ( std::size_t index = 0; index < variants.size(); ++index )
                {
                    varContainer.addCandidateVariant( variants[index], priors[index] );
                }
            }
        }

        if ( m_genotypeAllelesVCFFile != nullptr )
        {
            for ( const auto & record : m_genotypeAllelesVCFFile->fetch( blockRegion ) )
            {
                const auto variants = record.getVariants( referenceSequence );
                for ( std::size_t index = 0; index < variants.size(); ++index )
                {
                    varContainer.addGenotypingVariant( variants[index] );
                }
            }
        }

        varContainer.computeCoverage( blockRegion, allReads );

        const auto varFilter = variant::VariantFilter( m_privateCallingParams.m_minReadsPerVar,
                                                       m_privateCallingParams.m_perSamPercentReadsPerVar );

        auto variants = varFilter.getSortedFilteredVariants( blockRegion, varContainer );

        if ( m_privateCallingParams.m_turnOnLargeVariantCalls )
        {
            const auto breakpoints = varContainer.getBreakpoints();

            const variant::MultiLocusBreakpointVariantGenerator multiLocusBreakpointVariantGenerator(
                blockRegion, referenceSequence, params::defaults::defaultBreakpointKmerSize,
                params::defaults::maxBreakpointKmerSize, m_privateCallingParams.m_largeVariantSizeDefinition );
            const auto bpVariants = multiLocusBreakpointVariantGenerator.getVariantCandidates( breakpoints, allReads );

            for ( const auto & bpVar : bpVariants )
            {
                // This is a problem that small CNVs know no boundaries...
                if ( blockRegion.contains( bpVar->region() ) )
                {
                    variants.insert( bpVar );
                }
            }
        }

        return variant::generateMergedClusters( variants, m_privateCallingParams.m_minClusterDist,
                                                m_privateCallingParams.m_maxClusterDist,
                                                m_privateCallingParams.m_maxClusterSize, blockRegion, referenceSequence,
                                                m_privateCallingParams.m_maxClusterVariantCombinations,
                                                m_privateCallingParams.m_minReadsToMakeCombinationClaim );
    }

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

    int64_t Job::processBlock( io::readDataset_t readDataset,
                               bool lastBlock,
                               const std::vector< std::size_t > & ploidyPerSample )
    {
        auto blockRegion = readDataset->region();
        WECALL_LOG( INFO, "Processing:\t" << blockRegion );
        // TODO: Calibrate min value below.
        const auto allReads = readDataset->getAllReads( m_filterParams.m_readMappingFilterQ );

        // Get required reference.
        const auto maxReadLength = io::perSampleMaxAlignedReadLength( allReads );
        const auto assemblePadding = 3 * params::defaults::maxBreakpointKmerSize;
        const caller::Region paddedRefRegion = blockRegion.getPadded( maxReadLength + assemblePadding );
        m_ref.cacheSequence( paddedRefRegion );

        const auto referenceSequence =
            std::make_shared< utils::ReferenceSequence >( m_ref.getSequence( paddedRefRegion ) );

        auto clusters = this->generateVariantClustersInBlock( blockRegion, allReads, referenceSequence );
        if ( not lastBlock )
        {
            variant::removeClustersNearBlockEnd( blockRegion, clusters, m_privateCallingParams.m_maxClusterDist );
        }

        WECALL_LOG( INFO, "Calling variants:\t" << blockRegion );

        if ( clusters.empty() )
        {
            this->callReference( readDataset->region().contig(), blockRegion.start(), blockRegion.end(), readDataset,
                                 ploidyPerSample );
        }
        else
        {
            // init with first cluster
            std::vector< variant::VariantCluster >::iterator clusterIterator = clusters.begin();
            auto firstCluster = *clusterIterator;
            auto callsPrevCluster =
                this->processBigCluster( firstCluster, readDataset, referenceSequence, ploidyPerSample );
            this->writeCallsForCluster( callsPrevCluster, firstCluster, blockRegion.start(), readDataset,
                                        ploidyPerSample );
            Region regionPrevCluster = firstCluster.region();

            // iterate over remaining clusters
            clusterIterator++;
            for ( ; clusterIterator != clusters.cend(); ++clusterIterator )
            {
                // call variants in cluster
                const auto cluster = *clusterIterator;
                auto calls = this->processBigCluster( cluster, readDataset, referenceSequence, ploidyPerSample );

                callVector_t aligned_calls;

                // align phasing
                if ( callsPrevCluster.empty() )
                {
                    aligned_calls = calls;
                }
                else
                {
                    // get padded reference sequence for phasing
                    const Region combinedRegion =
                        Region( cluster.region().contig(), regionPrevCluster.start(), cluster.region().end() );
                    const io::perSampleRegionsReads_t combinedRegionReads =
                        readDataset->getRegionsReads( combinedRegion, m_filterParams.m_readMappingFilterQ );
                    const variant::VariantCluster combinedCluster = variant::VariantCluster( {}, combinedRegion );
                    const auto paddedReferenceSequence =
                        this->getReferenceForCluster( combinedCluster, combinedRegionReads, referenceSequence );

                    // attempt to phase clusters
                    alignPhasingBetweenClusters( calls, callsPrevCluster, cluster.region(), regionPrevCluster,
                                                 combinedRegion, paddedReferenceSequence, ploidyPerSample,
                                                 combinedRegionReads, readDataset->getSampleNames() );
                    aligned_calls = calls;
                }

                this->writeCallsForCluster( aligned_calls, cluster, regionPrevCluster.end(), readDataset,
                                            ploidyPerSample );

                // move one cluster forward
                callsPrevCluster = calls;
                regionPrevCluster = cluster.region();
            }

            // call block at the end
            this->callReference( readDataset->region().contig(), regionPrevCluster.end(), blockRegion.end(),
                                 readDataset, ploidyPerSample );
        }
        return blockRegion.end();
    }

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

    void Job::writeCallsForCluster( const callVector_t calls,
                                    const variant::VariantCluster cluster,
                                    const int64_t refStart,
                                    io::readDataset_t readDataset,
                                    const std::vector< std::size_t > & ploidyPerSample )
    {
        callVector_t outputCalls = this->filterOutputCalls( readDataset->region().contig(), calls );
        // output reference calls between previous and current cluster
        this->callReference( readDataset->region().contig(), refStart, cluster.zeroIndexedVcfStart( outputCalls ),
                             readDataset, ploidyPerSample );

        // write calls (variant + ref calls) for current cluster
        m_variantSoftFilterBank.applyFilterAnnotation( outputCalls );
        m_vcOut.writeCallSet( m_ref, outputCalls );
    }

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

    void Job::callReference( const std::string & contig,
                             int64_t start,
                             int64_t end,
                             io::readDataset_t readDataset,
                             const std::vector< std::size_t > & ploidyPerSample )
    {
        if ( m_dataParams.outputRefCalls() and start < end )
        {
            const caller::Region refInterval( contig, start, end );
            const auto reads = readDataset->getRegionsReads( refInterval, m_filterParams.m_readMappingFilterQ );

            callVector_t calls;
            const auto maxUncalledVaraintQuality = 0.0;
            const int64_t maxRefCallSize = m_dataParams.maxRefCallSize();

            for ( auto begBlock = refInterval.start(); begBlock < refInterval.end(); begBlock += maxRefCallSize )
            {
                const int64_t endBlock = std::min( begBlock + maxRefCallSize, refInterval.end() );
                const caller::Region region( contig, utils::Interval( begBlock, endBlock ) );
                const auto chunkedRefCalls =
                    caller::model::buildRefCall( region, reads, maxUncalledVaraintQuality, ploidyPerSample,
                                                 m_privateCallingParams.m_referenceCallQualityDeltaThreshold );
                calls.insert( calls.end(), chunkedRefCalls.begin(), chunkedRefCalls.end() );
            }

            const auto outputCalls = this->filterOutputCalls( contig, calls );

            m_vcOut.writeCallSet( m_ref, outputCalls );
        }
    }

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

    void Job::recalibrateBaseQualities( const io::perSampleRegionsReads_t & readRangesPerSample,
                                        const caller::Region & region ) const
    {
        const auto errorCorrectionParameters = corrector::ErrorCorrectionParameters();
        corrector::recalibrateReads( readRangesPerSample, m_ref, region, errorCorrectionParameters );
        m_intermediateOutputWriter.writeReads( readRangesPerSample, region.contig() );
    }

    std::vector< varPtr_t > getCandidateVariants( const std::vector< varPtr_t > & variants,
                                                  const variant::variantSet_t & newVariants )
    {
        if ( newVariants.size() == 0 )
        {
            return variants;
        }
        else
        {
            variant::variantSet_t variantsWithMNPs( variants.cbegin(), variants.cend() );
            variantsWithMNPs.insert( newVariants.cbegin(), newVariants.cend() );
            return std::vector< varPtr_t >( variantsWithMNPs.cbegin(), variantsWithMNPs.cend() );
        }
    }

    std::vector< std::size_t > getReducedPloidies( const callVector_t & calls,
                                                   const std::vector< std::size_t > & defaultPloidies,
                                                   const caller::Region & thisRegion )
    {
        auto reducedPloidies = defaultPloidies;
        for ( std::size_t sampleIndex = 0; sampleIndex < defaultPloidies.size(); ++sampleIndex )
        {
            auto & currentPloidy = reducedPloidies[sampleIndex];
            for ( const auto & call : calls )
            {
                if ( call.var->region().contains( thisRegion ) )
                {
                    const auto & genotypes = call.samples[sampleIndex].genotypeCalls;
                    for ( const auto & genotype : genotypes )
                    {
                        if ( genotype == caller::Call::VAR )
                        {
                            --currentPloidy;
                        }
                    }
                }
            }
        }

        return reducedPloidies;
    }

    callVector_t Job::processBigCluster( const variant::VariantCluster & cluster,
                                         io::readDataset_t readDataset,
                                         const utils::referenceSequencePtr_t & blockReferenceSequence,
                                         const std::vector< std::size_t > & ploidyPerSample )
    {
        variant::setDefaultPriors( cluster.variants() );

        bool hasLargeVariant = false;
        if ( m_privateCallingParams.m_turnOnLargeVariantCalls )
        {
            const auto largeVariantThreshold = m_privateCallingParams.m_maxClusterDist * 2;
            const auto isLargeVar = [largeVariantThreshold]( const varPtr_t & var )
            {
                return var->sequenceLengthInRef() >= largeVariantThreshold;
            };
            hasLargeVariant = std::any_of( cluster.variants().cbegin(), cluster.variants().cend(), isLargeVar );
        }

        const auto bigClusterReads =
            readDataset->getRegionsReads( cluster.region(), m_filterParams.m_readMappingFilterQ );
        const auto allReads = readDataset->getAllReads( m_filterParams.m_readMappingFilterQ );

        if ( not hasLargeVariant )
        {
            return processCluster( cluster, bigClusterReads, allReads, blockReferenceSequence, ploidyPerSample );
        }
        else
        {
            // TODO fix ref calling for large variant caller: see Regeneron ticket 11
            const auto data = cluster.buildSubClusters( m_privateCallingParams.m_maxClusterDist );
            auto lvcCluster = std::get< 0 >( data );
            auto smallVariantClustersNotTouchingLargeVariants = std::get< 1 >( data );

            lvcCluster.computeVariantCombinations(
                m_privateCallingParams.m_minReadsToMakeCombinationClaim, m_privateCallingParams.m_maxClusterDist,
                m_privateCallingParams.m_maxClusterVariantCombinations, blockReferenceSequence );
            for ( auto & miniCluster : smallVariantClustersNotTouchingLargeVariants )
            {
                miniCluster.computeVariantCombinations(
                    m_privateCallingParams.m_minReadsToMakeCombinationClaim, m_privateCallingParams.m_maxClusterDist,
                    m_privateCallingParams.m_maxClusterVariantCombinations, blockReferenceSequence );
            }

            const auto lvcReads = io::reduceRegionSet( bigClusterReads, lvcCluster.readRegions() );
            const auto lvcCalls =
                processCluster( lvcCluster, lvcReads, allReads, blockReferenceSequence, ploidyPerSample );

            callVector_t allCalls = lvcCalls;
            for ( const auto & leftOverCluster : smallVariantClustersNotTouchingLargeVariants )
            {
                const auto clusterReads = io::reduceRegionSet( bigClusterReads, leftOverCluster.readRegions() );
                const auto thisAreasPloidy = getReducedPloidies( lvcCalls, ploidyPerSample, leftOverCluster.region() );
                const auto leftOverCalls =
                    processCluster( leftOverCluster, clusterReads, allReads, blockReferenceSequence, thisAreasPloidy );
                allCalls.insert( allCalls.end(), leftOverCalls.begin(), leftOverCalls.end() );
            }

            std::sort( allCalls.begin(), allCalls.end(), caller::CallComp() );

            return allCalls;
        }
    }

    utils::referenceSequencePtr_t Job::getReferenceForCluster(
        const variant::VariantCluster & cluster,
        const io::perSampleRegionsReads_t & regionReads,
        const utils::referenceSequencePtr_t & blockReferenceSequence ) const
    {
        const int64_t minRefPadding = 40L;  // TODO: Handle this more safely. Needed for kmer building later in code.
        const auto maxReadLength = io::perSampleMaxReadCigarLength( regionReads );
        const auto refPaddingSize = std::max( minRefPadding, maxReadLength + constants::needlemanWunschPadding );

        return std::make_shared< utils::ReferenceSequence >(
            blockReferenceSequence->getPaddedSequence( cluster.region(), cluster.paddedRegion(), refPaddingSize ) );
    }

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

    callVector_t Job::processCluster( const variant::VariantCluster & cluster,
                                      const io::perSampleRegionsReads_t & regionReads,
                                      const io::perSampleRegionsReads_t & allReads,
                                      const utils::referenceSequencePtr_t & blockReferenceSequence,
                                      const std::vector< std::size_t > & ploidyPerSample )
    {
        WECALL_LOG( DEBUG, "Processing: " << cluster.toString() );

        if ( m_callingParams.m_recalibrateBaseQs )
        {
            this->recalibrateBaseQualities( regionReads, cluster.region() );
        }

        const auto paddedRefSequence = this->getReferenceForCluster( cluster, regionReads, blockReferenceSequence );

        variant::HaplotypeVector haplotypes( cluster.readRegions(), paddedRefSequence );

        if ( cluster.allCombinationsComputed() )
        {
            for ( const auto & varCombo : cluster.variantCombinations() )
            {
                haplotypes.push_back( variant::variantSet_t( varCombo.cbegin(), varCombo.cend() ) );
            }
            haplotypes.sort();
            haplotypes.merge();
        }
        else
        {
            WECALL_LOG( DEBUG, "Re-clustering " << cluster.region()
                                                 << " with nVariants: " << cluster.variants().size() );
            const variant::AlignmentHaplotypeGenerator hapGen(
                cluster.variants(), cluster.readRegions(), regionReads, paddedRefSequence,
                m_privateCallingParams.m_maxHaplotypesPerCluster,
                m_privateCallingParams.m_minReadsToMakeCombinationClaim );

            haplotypes = hapGen.generateHaplotypes();
            if ( haplotypes.size() <= 1 )
            {
                WECALL_LOG( WARNING, "Skipping region: " << cluster.region() << " with " << cluster.variants().size()
                                                          << " variants." );
                return {};
            }
        }
        WECALL_LOG( DEBUG, "Processing " << haplotypes.size() << " haplotypes" );

        std::vector< varPtr_t > candidateVariants;
        if ( m_privateCallingParams.m_normalizeVariantCalls )
        {
            const auto variants = haplotypes.withNormalizedVariants();
            candidateVariants = std::vector< varPtr_t >( variants.cbegin(), variants.cend() );
            variant::setDefaultPriors( candidateVariants );
        }
        else if ( m_callingParams.m_allowMNPCalls )
        {
            const auto newVariants = haplotypes.withMNPs();
            candidateVariants = getCandidateVariants( cluster.variants(), newVariants );

            variant::setDefaultPriors( std::vector< varPtr_t >( newVariants.cbegin(), newVariants.cend() ) );
        }
        else
        {
            candidateVariants = cluster.variants();
        }

        const auto results = m_model.getResults( regionReads, haplotypes, candidateVariants, ploidyPerSample );

        const auto calls = m_variantCallBuilder.getAnnotatedVariantCalls(
            cluster.region().start(), haplotypes.region(), regionReads, allReads, candidateVariants,
            results.getVariantQualities(), results.getVariantMetadata(), results.getCalledGenotypes(),
            results.getGenotypeMetadata(), ploidyPerSample, haplotypes,
            m_privateCallingParams.m_referenceCallQualityDeltaThreshold );

        return calls;
    }

    callVector_t Job::filterOutputCalls( const std::string & contig, const callVector_t & calls ) const
    {
        callVector_t outputCalls;
        for ( const auto & call : calls )
        {
            const Region callRegion( contig, call.interval );  // One bp before variant region for indels.
            const auto overlapCall = [callRegion]( const Region & outputRegion )
            {
                return outputRegion.overlaps( callRegion );
            };

            // TODO: Replace with better handling than linear scan.
            if ( std::any_of( m_outputRegions.begin(), m_outputRegions.end(), overlapCall ) )
            {
                outputCalls.push_back( call );
            }
        }

        return outputCalls;
    }
}
}

max_stars_repo_name: dylex/wecall

max_issues_repo_name: dylex/wecall

max_forks_repo_name: dylex/wecall


text:
// Copyright 2005 The Trustees of Indiana University.

// 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)

//  Authors: Douglas Gregor
//           Andrew Lumsdaine
#include <boost/graph/breadth_first_search.hpp>
#include <boost/graph/python/queue.hpp>
#include <boost/python.hpp>

namespace boost { namespace graph { namespace python {

using boost::python::object;

#define BGL_PYTHON_VISITOR maybe_bfs_visitor
#define BGL_PYTHON_EVENTS_HEADER <boost/graph/python/bfs_events.hpp>
#include <boost/graph/python/visitor.hpp>
#undef BGL_PYTHON_EVENTS_HEADER
#undef BGL_PYTHON_VISITOR

template<typename Graph>
void
breadth_first_search
  (const Graph& g,
   typename graph_traits<Graph>::vertex_descriptor s,
   object in_Q,
   object in_visitor,
   vector_property_map<
     default_color_type,
     typename property_map<Graph, vertex_index_t>::const_type>* in_color)
{
  typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

  typedef vector_property_map<
            default_color_type,
            typename property_map<Graph, vertex_index_t>::const_type>
    ColorMap;

  ColorMap color(in_color? *in_color
                 /*default*/: ColorMap(num_vertices(g), get(vertex_index, g)));

  bool has_default_queue = in_Q == object();
  if (!has_default_queue || in_visitor != object()) {
    typedef queue<vertex_descriptor> Queue;
    maybe_python_queue<Queue> Q(in_Q, Queue());
    maybe_bfs_visitor vis(in_visitor);
    boost::breadth_first_search(g, s, color_map(color).visitor(vis).buffer(Q));
  } else {
    boost::breadth_first_search(g, s, color_map(color));
  }
}

} } } // end namespace boost::graph::python

max_stars_repo_name: erwinvaneijk/bgl-python

max_issues_repo_name: erwinvaneijk/bgl-python

max_forks_repo_name: erwinvaneijk/bgl-python


text:
/**
 * Copyright 2021 Scott Brauer
 *
 * 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  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.
 */

/**
 * @file   GenerateSineCandlesticksProgramOptions.hpp
 * @author Scott Brauer
 * @date   11-29-2021
 */
#ifndef OFFCENTER_TRADING_GENERATESINECANDLESTICKS_GENERATESINECANDLESTICKSPROGRAMOPTIONS__H_
#define OFFCENTER_TRADING_GENERATESINECANDLESTICKS_GENERATESINECANDLESTICKSPROGRAMOPTIONS__H_

#include <functional>
#include <iostream>
#include <sstream>

#include <boost/bind/bind.hpp>

#include "offcenter/common/program_options/ProgramOptionsGroup.hpp"
#include "offcenter/trading/generatesinecandlesticks/GenerateSineCandlesticksOptions.hpp"

namespace offcenter {
namespace trading {
namespace generatesinecandlesticks {

class GenerateSineCandlesticksProgramOptions:
		public offcenter::common::program_options::ProgramOptionsGroup<offcenter::trading::generatesinecandlesticks::GenerateSineCandlesticksOptions>
{
public:
	GenerateSineCandlesticksProgramOptions():
		ProgramOptionsGroup("Generate Sine Candlestick Options")
	{
		optionsDescription().add_options()
				("start,s", po::value<std::string>()->required()->notifier(boost::bind(&GenerateSineCandlesticksOptions::setStart, &optionConfig(), boost::placeholders::_1)), "Start time (2004-06-09T20:59:55.000000000Z)")
				("end,e", po::value<std::string>()->required()->notifier(boost::bind(&GenerateSineCandlesticksOptions::setEnd, &optionConfig(), boost::placeholders::_1)), "End time (2004-08-09T20:59:55.000000000Z)")
				("duration,d", po::value<int>(&optionConfig().m_duration)->required()->default_value(5), "Duration in seconds")

				("sine0.amplitude", po::value<double>(&optionConfig().m_sine0.m_amplitude)->required(), "Sine 0 amplitude")
				("sine0.wavelength", po::value<int>(&optionConfig().m_sine0.m_wavelength)->required(), "Sine 0 wavelength")
				("sine0.amplitude_shift", po::value<double>(&optionConfig().m_sine0.m_amplitudeShift)->required(), "Sine 0 amplitude shift")
				("sine0.wavelength_shift", po::value<int>(&optionConfig().m_sine0.m_wavelengthShift)->required(), "Sine 0 wavelength shift")

				("sine1.amplitude", po::value<double>(&optionConfig().m_sine1.m_amplitude), "Sine 1 amplitude")
				("sine1.wavelength", po::value<int>(&optionConfig().m_sine1.m_wavelength), "Sine 1 wavelength")
				("sine1.amplitude_shift", po::value<double>(&optionConfig().m_sine1.m_amplitudeShift), "Sine 1 amplitude shift")
				("sine1.wavelength_shift", po::value<int>(&optionConfig().m_sine1.m_wavelengthShift), "Sine 1 wavelength shift")

				("sine2.amplitude", po::value<double>(&optionConfig().m_sine2.m_amplitude), "Sine 2 amplitude")
				("sine2.wavelength", po::value<int>(&optionConfig().m_sine2.m_wavelength), "Sine 2 wavelength")
				("sine2.amplitude_shift", po::value<double>(&optionConfig().m_sine2.m_amplitudeShift), "Sine 2 amplitude shift")
				("sine2.wavelength_shift", po::value<int>(&optionConfig().m_sine2.m_wavelengthShift), "Sine 2 wavelength shift")

				("sine3.amplitude", po::value<double>(&optionConfig().m_sine3.m_amplitude), "Sine 3 amplitude")
				("sine3.wavelength", po::value<int>(&optionConfig().m_sine3.m_wavelength), "Sine 3 wavelength")
				("sine3.amplitude_shift", po::value<double>(&optionConfig().m_sine3.m_amplitudeShift), "Sine 3 amplitude shift")
				("sine3.wavelength_shift", po::value<int>(&optionConfig().m_sine3.m_wavelengthShift), "Sine 3 wavelength shift")

				;
	}

	virtual ~GenerateSineCandlesticksProgramOptions() = default;

	virtual void validateOptions(po::variables_map& vm) {

	}
};

} /* namespace generatesinecandlesticks */
} /* namespace trading */
} /* namespace offcenter */

#endif /* OFFCENTER_TRADING_GENERATESINECANDLESTICKS_GENERATESINECANDLESTICKSPROGRAMOPTIONS__H_ */


text:
#include <boost/test/unit_test.hpp>
#include "mod/vector.hpp"

#include "objects/numerics/float.hpp"
#include "objects/numerics/integer.hpp"
#include "objects/collections/ast.hpp"

#include "objects/invokable/effect.hpp"
#include "objects/invokable/continuation.hpp"
#include "objects/invokable/inbuilt_function.hpp"
#include "objects/invokable/inbuilt_macro.hpp"
#include "objects/invokable/generic_function.hpp"

#include "globals/core.hpp"
#include "globals/parents.hpp"
#include <iostream>

using obj::Object;

BOOST_AUTO_TEST_SUITE( invokable_suite )

BOOST_AUTO_TEST_CASE( mod_functions_inbuilt )
{
  bool func_called = false;
  Object* add =
    new obj::InbuiltFunction<>
    ([&](context& ctx){
      func_called = true;
      return core::mtrue;
    } ,"","");

  Object* out = std::get<Object*>((*add)({}, *context::empty));

  BOOST_TEST(func_called);
  BOOST_TEST(out == core::mtrue);
}

BOOST_AUTO_TEST_CASE( mod_inbuilt_function_single_arg )
{
  Object* id =
    new obj::InbuiltFunction<Object*>
    ([&](context& ctx, Object* obj){
      return obj;
    } ,"","");

  Object* out = std::get<Object*>((*id)({core::mtrue}, *context::empty));

  BOOST_TEST(out == core::mtrue);
}

BOOST_AUTO_TEST_CASE( mod_inherit_leq_native_objects_true )
{
  using obj::Float;
  Object* o1 = new Float(2.0);
  Object* o2 = parents::floatp;

  BOOST_TEST(inherit_leq(o1, o2));

  using obj::Integer;
  o1 = new Integer(2);
  o2 = parents::integerp;

  BOOST_TEST(inherit_leq(o1, o2));
}

BOOST_AUTO_TEST_CASE( mod_inherit_leq_native_objects_false )
{
  using obj::Float;
  Object* o1 = new Float(2.0);
  Object* o2 = parents::floatp;

  BOOST_TEST(!inherit_leq(o2, o1));

  using obj::Integer;
  o1 = new Integer(2);
  o2 = parents::integerp;

  BOOST_TEST(!inherit_leq(o2, o1));
}

BOOST_AUTO_TEST_CASE( mod_inherit_leq_nullptr_test )
{
  using obj::Float;
  Object* o1 = new Float(2.0);

  BOOST_TEST(!inherit_leq(nullptr, o1));

  BOOST_TEST(inherit_leq(o1, nullptr));
}

BOOST_AUTO_TEST_CASE( mod_inherit_args_native_objects_true )
{
  using obj::Float;
  using obj::Integer;
  Object* o1 = new Float(2.0);
  Object* o2 = parents::floatp;

  Object* o3 = new Integer(2);
  Object* o4 = parents::integerp;

  BOOST_TEST(args_leq({o1, o3}, {o2, o4}));
}

BOOST_AUTO_TEST_CASE( mod_inherit_args_native_objects_false )
{
  using obj::Float;
  using obj::Integer;
  Object* o1 = new Float(2.0);
  Object* o2 = parents::floatp;

  Object* o3 = new Integer(2);
  Object* o4 = parents::floatp;

  BOOST_TEST(!args_leq({o2, o4}, {o1, o3}));
}

BOOST_AUTO_TEST_CASE( effect_invocation_returns_effect_with_args )
{
  using obj::Integer;
  using obj::Effect;
  using obj::RaiseEffect;
  Object* myeffect = new Effect(2);
  Object* i1 = new Integer(1);
  Object* i2 = new Integer(2);

  obj::MaybeEffect out = (*myeffect)({i1, i2}, *context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(out));

  RaiseEffect eff = std::get<RaiseEffect>(out);

  BOOST_TEST( eff.action == myeffect );
  BOOST_TEST( eff.args.size() == 2 );

  BOOST_TEST( eff.args[0] == i1 );
  BOOST_TEST( eff.args[1] == i2 );
}

BOOST_AUTO_TEST_CASE( effect_continuation_is_identity )
{
  using obj::Integer;
  using obj::Effect;
  using obj::RaiseEffect;
  Object* myeffect = new Effect(0);
  Object* i1 = new Integer(-10);

  obj::MaybeEffect out = (*myeffect)({}, *context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(out));

  RaiseEffect eff = std::get<RaiseEffect>(out);

  BOOST_TEST( eff.action == myeffect );
  BOOST_TEST( eff.args.size() == 0 );

  obj::MaybeEffect cntcall = (*eff.continuation)({i1}, *context::empty);

  BOOST_TEST(std::holds_alternative<Object*>(cntcall));

  Object* cntres = std::get<Object*>(cntcall);

  BOOST_TEST(cntres == i1);
}

BOOST_AUTO_TEST_CASE( function_late_continuation_is_correct )
{
  // DETAILS OF TEST: we have some code which adds two numbers:
  // (+ <effect> 5). When the effect is raised, the resulting continuation
  // should be identical to the function λ x. x + 5

  using obj::Integer;
  using obj::Effect;
  using obj::RaiseEffect;
  using obj::mk_ast;

  Object* myeffect = new Effect(0);
  Object* i1 = new Integer(-10);
  Object* i2 = new Integer(5);

  Object* add =
    new obj::InbuiltFunction<Integer*, Integer*>
    ([&](context& ctx, Integer* i1, Integer* i2){
      return new Integer(i1->val + i2->val);
    } ,"","");

  Object* ast = obj::mk_ast({add, obj::mk_ast({myeffect}), i2});

  obj::MaybeEffect out = ast->eval(*context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(out));

  RaiseEffect eff = std::get<RaiseEffect>(out);

  BOOST_TEST( eff.action == myeffect );
  BOOST_TEST( eff.args.size() == 0 );

  obj::MaybeEffect cntcall = (*eff.continuation)({i1}, *context::empty);

  BOOST_TEST(std::holds_alternative<Object*>(cntcall));

  Integer* cntres = dynamic_cast<Integer*>(std::get<Object*>(cntcall));

  BOOST_TEST(cntres != nullptr);
  BOOST_TEST(cntres->val == -5);
}

BOOST_AUTO_TEST_CASE( function_early_continuation_is_correct )
{
  // DETAILS OF TEST: we have some code which adds two numbers:
  // (<effect> -10 5). When the effect is raised, the resulting continuation
  // should be identical to the function λ x. x + 5

  using obj::Integer;
  using obj::Effect;
  using obj::RaiseEffect;
  using obj::mk_ast;

  Object* myeffect = new Effect(0);
  Object* i1 = new Integer(-10);
  Object* i2 = new Integer(5);

  Object* add =
    new obj::InbuiltFunction<Integer*, Integer*>
    ([&](context& ctx, Integer* i1, Integer* i2) {
      return new Integer(i1->val + i2->val);
    } ,"test_add","test_add_doc");

  Object* ast = obj::mk_ast({obj::mk_ast({myeffect}), i1, i2});

  obj::MaybeEffect out = ast->eval(*context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(out));

  RaiseEffect eff = std::get<RaiseEffect>(out);

  BOOST_TEST( eff.action == myeffect );
  BOOST_TEST( eff.args.size() == 0 );

  obj::MaybeEffect cntcall = (*eff.continuation)({add}, *context::empty);

  BOOST_TEST(std::holds_alternative<Object*>(cntcall));

  Integer* cntres = dynamic_cast<Integer*>(std::get<Object*>(cntcall));

  BOOST_TEST(cntres != nullptr);
  BOOST_TEST(cntres->val == -5);
}

BOOST_AUTO_TEST_CASE( function_early_multiple_continuations_pos )
{
  using obj::Integer;
  using obj::Effect;
  using obj::RaiseEffect;
  using obj::mk_ast;

  Object* myeffect = new Effect(0);
  Object* inte = new Integer(5);

  Object* add =
    new obj::InbuiltFunction<Integer*, Integer*>
    ([&](context& ctx, Integer* i1, Integer* i2) {
      return new Integer(i1->val + i2->val);
    } ,"test_add","test_add_doc");

  Object* ast = obj::mk_ast({add, obj::mk_ast({myeffect}), obj::mk_ast({myeffect})});

  obj::MaybeEffect out = ast->eval(*context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(out));

  RaiseEffect eff = std::get<RaiseEffect>(out);
  obj::MaybeEffect cntcall1 = (*eff.continuation)({inte}, *context::empty);

  BOOST_TEST(std::holds_alternative<RaiseEffect>(cntcall1));

  eff = std::get<RaiseEffect>(cntcall1);
  obj::MaybeEffect cntcall2 = (*eff.continuation)({inte}, *context::empty);

  BOOST_TEST(std::holds_alternative<Object*>(cntcall2));

  Integer* result = dynamic_cast<Integer*>(std::get<Object*>(cntcall2));

  BOOST_TEST(result != nullptr);
  BOOST_TEST(result->val == 10);
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: rationalis-petra/modulus

max_issues_repo_name: rationalis-petra/modulus

max_forks_repo_name: rationalis-petra/modulus


text:
/*
 * Copyright 2021 MusicScience37 (Kenta Kabashima)
 *
 * 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.
 */
/*!
 * \file
 * \brief Test of kernel_parameter_optimizer class.
 */
#include "num_collect/interp/kernel/impl/kernel_parameter_optimizer.h"

#include <Eigen/Core>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating.hpp>

#include "eigen_approx.h"
#include "is_finite.h"
#include "num_collect/interp/kernel/calc_kernel_mat.h"
#include "num_collect/interp/kernel/euclidean_distance.h"
#include "num_collect/interp/kernel/gaussian_rbf.h"
#include "num_collect/interp/kernel/impl/auto_regularizer_wrapper.h"
#include "num_collect/interp/kernel/rbf_kernel.h"

TEST_CASE("num_collect::interp::kernel::impl::kernel_parameter_optimizer") {
    using num_collect::interp::kernel::euclidean_distance;
    using num_collect::interp::kernel::gaussian_rbf;
    using num_collect::interp::kernel::rbf_kernel;
    using num_collect::interp::kernel::impl::auto_regularizer_wrapper;
    using num_collect::interp::kernel::impl::kernel_parameter_optimizer;

    const auto vars = std::vector<double>{0.0, 0.1, 0.2, 0.4, 0.6, 1.0};
    const auto data = Eigen::VectorXd{{0.0, 0.2, 0.4, 0.7, 1.0, 2.0}};

    using kernel_type =
        rbf_kernel<euclidean_distance<double>, gaussian_rbf<double>>;
    auto kernel = kernel_type();

    auto interpolator = auto_regularizer_wrapper<double>();

    SECTION("compute") {
        auto optimizer =
            kernel_parameter_optimizer<kernel_type>(interpolator, kernel);

        optimizer.compute(vars, data);

        kernel.kernel_param(optimizer.opt_param());
        const Eigen::MatrixXd kernel_mat = calc_kernel_mat(kernel, vars);

        interpolator.compute(kernel, vars, data);
        Eigen::VectorXd coeff;
        interpolator.solve(coeff);

        const Eigen::VectorXd retrieved_data = kernel_mat * coeff;
        constexpr double tol_error = 1e-4;
        REQUIRE_THAT(retrieved_data, eigen_approx(data, tol_error));
    }
}


text:
//
//  FileSystem.cpp
//  Bump
//
//  Created by Christian Noon on 12/4/12.
//  Copyright (c) 2012 Christian Noon. All rights reserved.
//

// Boost headers
#include <boost/filesystem.hpp>
#include <boost/foreach.hpp>

// Bump headers
#include <bump/FileSystem.h>
#include <bump/FileSystemError.h>

// C++ headers
#include <fstream>

namespace bump {

namespace FileSystem {

//====================================================================================
//                               Path Coversion Methods
//====================================================================================

String convertToWindowsPath(String path)
{
	// Replace a forward slash next to a backward slash with a forward slash
	String converted = path.replace("/\\", "/");

	// Replace a backward slash next to a forward slash with a forward slash
	converted = path.replace("\\/", "/");

	// Replace all duplicate forward slashes with a single forward slash
	converted = path.replace("///", "/");
	converted = path.replace("//", "/");

	// Replace all duplicate backward slashes with a single backward slash
	converted = path.replace("\\\\\\", "\\");
	converted = path.replace("\\\\", "\\");

	// Convert all forward slashes to backward slashes
	converted = path.replace("/", "\\");

	return converted;
}

String convertToUnixPath(String path)
{
	// Replace a forward slash next to a backward slash with a forward slash
	String converted = path.replace("/\\", "/");

	// Replace a backward slash next to a forward slash with a forward slash
	converted = path.replace("\\/", "/");

	// Replace all duplicate forward slashes with a single forward slash
//	converted = path.replace("///", "/");
//	converted = path.replace("//", "/");

	// Replace all duplicate backward slashes with a single backward slash
//	converted = path.replace("\\\\\\", "\\");
//	converted = path.replace("\\\\", "\\");

	// Convert all backward slashes to forward slashes
	converted = path.replace("\\", "/");

	return converted;
}

String join(const String& path1, const String& path2)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4, const String& path5)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	joined_path /= path5.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4, const String& path5,
			const String& path6)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	joined_path /= path5.c_str();
	joined_path /= path6.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4, const String& path5,
			const String& path6, const String& path7)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	joined_path /= path5.c_str();
	joined_path /= path6.c_str();
	joined_path /= path7.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4, const String& path5,
			const String& path6, const String& path7, const String& path8)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	joined_path /= path5.c_str();
	joined_path /= path6.c_str();
	joined_path /= path7.c_str();
	joined_path /= path8.c_str();
	return convertToUnixPath(joined_path.string());
}

String join(const String& path1, const String& path2, const String& path3, const String& path4, const String& path5,
			const String& path6, const String& path7, const String& path8, const String& path9)
{
	boost::filesystem::path joined_path(path1.c_str());
	joined_path /= path2.c_str();
	joined_path /= path3.c_str();
	joined_path /= path4.c_str();
	joined_path /= path5.c_str();
	joined_path /= path6.c_str();
	joined_path /= path7.c_str();
	joined_path /= path8.c_str();
	joined_path /= path9.c_str();
	return convertToUnixPath(joined_path.string());
}

//====================================================================================
//                                System Path Methods
//====================================================================================

bool setCurrentPath(const String& path)
{
	// Set the current path
	try
	{
		boost::filesystem::current_path(boost::filesystem::path(path.c_str()));
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

String currentPath()
{
	String path = boost::filesystem::current_path().string();
	return convertToUnixPath(path);
}

String temporaryPath()
{
	String path = boost::filesystem::temp_directory_path().string();
	return convertToUnixPath(path);
}

//====================================================================================
//                               Path Query Methods
//====================================================================================

bool exists(const String& path)
{
	FileInfo path_info(path);
	return path_info.exists();
}

bool isDirectory(const String& path)
{
	FileInfo path_info(path);
	return path_info.isDirectory();
}

bool isFile(const String& path)
{
	FileInfo path_info(path);
	return path_info.isFile();
}

bool isSymbolicLink(const String& path)
{
	FileInfo path_info(path);
	return path_info.isSymbolicLink();
}

//====================================================================================
//                                 Directory Methods
//====================================================================================

bool createDirectory(const String& path)
{
	try
	{
		return boost::filesystem::create_directory(boost::filesystem::path(path.c_str()));
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool createFullDirectoryPath(const String& path)
{
	return boost::filesystem::create_directories(boost::filesystem::path(path.c_str()));
}

bool removeDirectory(const String& path)
{
	// Fail if the path is not a directory
	if (!FileInfo(path).isDirectory())
	{
		return false;
	}

	try
	{
		return boost::filesystem::remove(boost::filesystem::path(path.c_str()));
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool removeDirectoryAndContents(const String& path)
{
	// Fail if the path is not a directory
	if (!FileInfo(path).isDirectory())
	{
		return false;
	}

	try
	{
		boost::filesystem::remove_all(boost::filesystem::path(path.c_str()));
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool copyDirectory(const String& source, const String& destination)
{
	// Fail if the source path is not a directory
	if (!FileInfo(source).isDirectory())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::copy(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool copyDirectoryAndContents(const String& source, const String& destination)
{
	// Fail if the source path is not a directory
	if (!FileInfo(source).isDirectory())
	{
		return false;
	}

	// Create a FileInfo objects for the source and destination
	FileInfo source_info(source);
	FileInfo destination_info(destination);

	// Make sure the destination does not exist
	if (destination_info.exists())
	{
		return false;
	}

	// Create the destination directory
	bool created_destination_directory = createFullDirectoryPath(destination);
	if (!created_destination_directory)
	{
		return false;
	}

	// Iterate through the local source directory
	FileInfoList local_source_files = directoryInfoList(source);
	BOOST_FOREACH (const FileInfo& local_source_file, local_source_files)
	{
		bool copied_successfully;

		// Attempt to copy all types of file system objects
		if (local_source_file.isDirectory())
		{
			// Create the local destination directory path
			String filename = local_source_file.filename();
			String local_destination_directory = join(destination, filename);

			// Copy the local source directory to the local destination directory
			copied_successfully = copyDirectoryAndContents(local_source_file.absolutePath(), local_destination_directory);
		}
		else if (local_source_file.isSymbolicLink())
		{
			// Create the local destination symbolic link
			String filename = local_source_file.filename();
			String local_destination_file = join(destination, filename);

			// Copy the local source symbolic link to the local destination
			copied_successfully = copySymbolicLink(local_source_file.absolutePath(), local_destination_file);
		}
		else // local_source_info.isFile()
		{
			// Create the local destination filepath
			String filename = local_source_file.filename();
			String local_destination_file = join(destination, filename);

			// Copy the local source file to the local destination
			copied_successfully = copyFile(local_source_file.absolutePath(), local_destination_file);
		}

		// Fail if we did not copy the file system object successfully
		if (!copied_successfully)
		{
			return false;
		}
	}

	return true;
}

bool renameDirectory(const String& source, const String& destination)
{
	// Fail if the source path is not a directory
	if (!FileInfo(source).isDirectory())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::rename(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

StringList directoryList(const String& path)
{
	// Throw an exception if the path does not exist
	if (!FileInfo(path).exists())
	{
		String msg = String("The following path is not valid: %1").arg(path);
		throw FileSystemError(msg, BUMP_LOCATION);
	}

	// Throw an exception if the path is not a directory
	if (!FileInfo(path).isDirectory())
	{
		String msg = String("The following path is not a directory: %1").arg(path);
		throw FileSystemError(msg, BUMP_LOCATION);
	}

	// Create a path and two directory iterators
	boost::filesystem::path directory_path(path.c_str());
	boost::filesystem::directory_iterator iter(directory_path);
	boost::filesystem::directory_iterator end_iter;

	// Iterate through the directory collecting all the item paths as strings
	StringSet directory_set;
	BOOST_FOREACH (const boost::filesystem::path& item, std::make_pair(iter, end_iter))
	{
		String unix_item_str = convertToUnixPath(item.string());
		directory_set.insert(unix_item_str);
	}

	// Build a string list from the sorted string set
	StringList directory_list;
	BOOST_FOREACH (const String& item, directory_set)
	{
		directory_list.push_back(item);
	}

	return directory_list;
}

FileInfoList directoryInfoList(const String& path)
{
	// Throw an exception if the path does not exist
	if (!FileInfo(path).exists())
	{
		String msg = String("The following path is not valid: %1").arg(path);
		throw FileSystemError(msg, BUMP_LOCATION);
	}

	// Throw an exception if the path is not a directory
	if (!FileInfo(path).isDirectory())
	{
		String msg = String("The following path is not a directory: %1").arg(path);
		throw FileSystemError(msg, BUMP_LOCATION);
	}

	// Create a path and two directory iterators
	boost::filesystem::path directory_path(path.c_str());
	boost::filesystem::directory_iterator iter(directory_path);
	boost::filesystem::directory_iterator end_iter;

	// Iterate through the directory collecting all the item paths as strings
	StringSet directory_set;
	BOOST_FOREACH (const boost::filesystem::path& item, std::make_pair(iter, end_iter))
	{
		String unix_item_str = convertToUnixPath(item.string());
		directory_set.insert(unix_item_str);
	}

	// Build a file info list from the sorted string set
	FileInfoList directory_list;
	BOOST_FOREACH (const String& item, directory_set)
	{
		FileInfo item_info(item);
		directory_list.push_back(item_info);
	}

	return directory_list;
}

//====================================================================================
//                                   File Methods
//====================================================================================

bool createFile(const String& path)
{
	// Make sure the file doesn't already exist
	FileInfo path_info(path);
	if (path_info.exists())
	{
		return false;
	}

	// Create the file and close it
	std::ofstream stream(path.c_str());
	bool success = !stream.fail();
	stream.close();

	return success;
}

bool removeFile(const String& path)
{
	// Fail if path is not a file
	if (!FileInfo(path).isFile())
	{
		return false;
	}

	return boost::filesystem::remove(boost::filesystem::path(path.c_str()));
}

bool copyFile(const String& source, const String& destination)
{
	// Fail if source is not a file
	if (!FileInfo(source).isFile())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::copy(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool renameFile(const String& source, const String& destination)
{
	// Fail if source is not a file
	if (!FileInfo(source).isFile())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::rename(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

//====================================================================================
//                               Symbolic Link Methods
//====================================================================================

bool createDirectorySymbolicLink(const String& source, const String& destination)
{
	try
	{
		// Convert the source path to the preferred OS
		boost::filesystem::path source_path(boost::filesystem::path(source.c_str()).make_preferred().string());
		boost::filesystem::path destination_path(destination.c_str());

		// Create the directory symlink
		boost::filesystem::create_directory_symlink(source_path, destination_path);

		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool createFileSymbolicLink(const String& source, const String& destination)
{
	try
	{
		// Convert the source path to the preferred OS
		boost::filesystem::path source_path(boost::filesystem::path(source.c_str()).make_preferred().string());
		boost::filesystem::path destination_path(destination.c_str());

		// Create the file symlink
		boost::filesystem::create_symlink(source_path, destination_path);

		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool removeSymbolicLink(const String& path)
{
	// Fail if path is not a symbolic link
	if (!FileInfo(path).isSymbolicLink())
	{
		return false;
	}

	try
	{
		return boost::filesystem::remove(boost::filesystem::path(path.c_str()));
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool copySymbolicLink(const String& source, const String& destination)
{
	// Fail if source is not a symbolic link
	if (!FileInfo(source).isSymbolicLink())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::copy_symlink(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

bool renameSymbolicLink(const String& source, const String& destination)
{
	// Fail if source is not a symbolic link
	if (!FileInfo(source).isSymbolicLink())
	{
		return false;
	}

	try
	{
		boost::filesystem::path source_path(source.c_str());
		boost::filesystem::path destination_path(destination.c_str());
		boost::filesystem::rename(source_path, destination_path);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

//====================================================================================
//                                  Date Methods
//====================================================================================

bool setModifiedDate(const String& path, const std::time_t& date)
{
	// Fail if path does NOT exist
	if (!exists(path))
	{
		return false;
	}

	try
	{
		boost::filesystem::last_write_time(boost::filesystem::path(path.c_str()), date);
		return true;
	}
	catch (const boost::filesystem::filesystem_error& /*e*/)
	{
		return false;
	}
}

std::time_t modifiedDate(const String& path)
{
	FileInfo info(path);
	return info.modifiedDate();
}

}	// End of FileSystem namespace

}	// End of bump namespace

max_stars_repo_name: sourada/bump

max_issues_repo_name: sourada/bump

max_forks_repo_name: sourada/bump


text:
/*
Copyright (c) Members of the EGEE Collaboration. 2004.
See http://www.eu-egee.org/partners for details on the
copyright holders.

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.
*/

// File: common.cpp
// Author: Salvatore Monforte <Salvatore.Monforte@ct.infn.it>
// Copyright (c) 2002 EU DataGrid.

// $Id: common.cpp,v 1.8.2.4.6.2.4.3 2012/06/22 11:51:31 mcecchi Exp $

#include <boost/tokenizer.hpp>
#include <boost/bind.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/erase.hpp>
#include <boost/algorithm/string/regex.hpp>

#include "glite/wms/ism/ism.h"
#include "glite/wms/ism/purchaser/common.h"
#include "glite/wms/common/logger/logger_utils.h"
#include "glite/wmsutils/classads/classad_utils.h"
#include "glite/wms/common/configuration/Configuration.h"
#include "common/src/configuration/WMConfiguration.h"
#include "common/src/configuration/NSConfiguration.h"

#include "schema_utils.h"

using namespace std;
namespace utils = glite::wmsutils::classads;
namespace ba = boost::algorithm;

namespace glite {
namespace wms {
namespace ism {
namespace purchaser {

bdii_schema_info_type& bdii_schema_info()
{
  static bdii_schema_info_type glue_schema("GLUE");
  return glue_schema;
}

bdii_schema_info_type& bdii_schema_info_g2()
{
  static bdii_schema_info_type glue2_schema("GLUE2");
  return glue2_schema;
}

inline bool iequals(std::string const& a, std::string const& b)
{
 return ba::iequals(a,b);
}

inline bool istarts_with(std::string const& a, std::string const& b)
{
 return ba::istarts_with(a,b);
}

inline std::string strip_prefix(std::string const& prefix, std::string const& s)
{
  return boost::algorithm::istarts_with(s,prefix) ?
    boost::algorithm::erase_head_copy(s, prefix.length()) : s;
}

void insert_values(
  bool is_schema_version_20,
  std::string const& name,
  boost::shared_array<struct berval*> values,
  std::list<std::string> const& prefix,
  classad::ClassAd& ad
) {
  std::vector<std::string>::const_iterator list_b, number_b;
  std::vector<std::string>::const_iterator list_e, number_e;

  if (is_schema_version_20) {
    boost::tie(list_b, list_e) = bdii_schema_info_g2().multi_valued();
    boost::tie(number_b, number_e) = bdii_schema_info_g2().number_valued();
  } else {
    boost::tie(list_b, list_e) = bdii_schema_info().multi_valued();
    boost::tie(number_b, number_e) = bdii_schema_info().number_valued();
  }

  bool is_list = ba::iequals(name,"objectclass") || std::find_if(
    list_b, list_e, boost::bind(iequals,_1, name)
  ) != list_e;

  bool is_number = std::find_if(
    number_b, number_e, boost::bind(iequals, _1, name)
  ) != number_e;

  std::string result;
  for (size_t i = 0; values[i] != 0; ++i) {

    if (i) {
      result.append(",");
    }
    if (is_number || ba::iequals("undefined", values[i]->bv_val) ||
      ba::iequals("false", values[i]->bv_val) || ba::iequals("true", values[i]->bv_val)) {

      result.append(values[i]->bv_val);
    }
    else {
      result.append("\"").append(values[i]->bv_val).append("\"");
    }
  }
  if (is_list) {
    std::string s("{");
    s.append(result).append("}");
    result.swap(s);
  }
  classad::ExprTree* e = 0;
  classad::ClassAdParser parser;
  parser.ParseExpression(result, e);

  if (e) {
    std::list<std::string>::const_iterator it =  std::find_if(
      prefix.begin(), prefix.end(),
      boost::bind(istarts_with, name, _1)
    );
    ad.Insert(
      ( it != prefix.end() ) ? strip_prefix(*it, name) : name,
      e
    );
  }
}

classad::ClassAd*
create_classad_from_ldap_entry(
  LDAP* ld,
  LDAPMessage* lde,
  std::list<std::string> prefix,
  bool is_schema_version_20
) {
  classad::ClassAd* result(new classad::ClassAd);
  BerElement* ber = 0;
  for (char* attr = ldap_first_attribute(ld, lde, &ber);
    attr;
    attr = ldap_next_attribute(ld, lde, ber)
  ) {
    boost::shared_ptr<void> attr_guard(
      attr, ber_memfree
    );
    boost::shared_array<struct berval*> values(
      ldap_get_values_len(ld, lde, attr),
      ldap_value_free_len
    );

    if (!values) continue;
    insert_values(
      is_schema_version_20, attr, values, prefix, *result
    );
  }
  boost::shared_ptr<void> ber_guard(
    static_cast<void*>(0),boost::bind(ber_free, ber, 0)
  );
  return result;
}

namespace {

boost::shared_ptr<classad::ClassAd> glue13_mapping_ad(
   glite::wmsutils::classads::parse_classad(
    "["
   "GlueCECapability = GLUE2.Computing.Endpoint.Capability;"
   "GlueCEImplementationName = GLUE2.Computing.Endpoint.ImplementationName;"
   "GlueCEImplementationVersion = GLUE2.Computing.Endpoint.ImplementationVersion;"
   "GlueCEInfoLRMSVersion = GLUE2.Manager.ProductVersion;"
   "GlueCEPolicyMaxWallClockTime = GLUE2.Computing.Share.MaxWallTime;"
   "GlueCEPolicyMaxCPUTime = GLUE2.Computing.Share.MaxCPUTime;"
   "GlueCEPolicyMaxRunningJobs  = GLUE2.Computing.Share.MaxRunningJobs;"
   "GlueCEStateStatus = GLUE2.Computing.Share.ServingState;"
   "GlueCEStateWaitingJobs = GLUE2.Computing.Share.WaitingJobs;"
   "GlueCEStateTotalJobs = GLUE2.Computing.Share.TotalJobs;"
   "GlueCEStateFreeJobSlots  = GLUE2.Computing.Share.FreeSlots;"
   "GlueCEStateRunningJobs = GLUE2.Computing.Share.RunningJobs;"
   "GlueCEStateEstimatedResponseTime = GLUE2.Computing.Share.EstimatedAverageWaitingTime;"
   "GlueCEStateWorstResponseTime = GLUE2.Computing.Share.EstimatedWorstWaitingTime;"
   "GlueHostArchitecturePlatformType = GLUE2.ExecutionEnvironment.Platform;"
   "GlueHostArchitectureSMPSize = GLUE2.ExecutionEnvironment.OtherInfo.SmpSize;"
   "GlueHostProcessorModel = GLUE2.ExecutionEnvironment.CPUModel;"
   "GlueHostProcessorVendor = GLUE2.ExecutionEnvironment.CPUVendor;"
   "GlueHostProcessorClockSpeed = GLUE2.ExecutionEnvironment.CPUClockSpeed;"
   "GlueHostOperatingSystemName = GLUE2.ExecutionEnvironment.OSName;"
   "GlueHostMainMemoryRAMSize = GLUE2.ExecutionEnvironment.MainMemorySize;"
   "GlueHostMainMemoryVirtualSize = GLUE2.ExecutionEnvironment.VirtualMemorySize;"
   "GlueHostNetworkAdapterInboundIP = GLUE2.ExecutionEnvironment.ConnectivityIn;"
   "GlueHostNetworkAdapterOutboundIP = GLUE2.ExecutionEnvironment.ConnectivityOut;"
   "GlueSubClusterLogicalCPUs = GLUE2.ExecutionEnvironment.LogicalCPUs;"
   "GlueSubClusterPhysicalCPUs = GLUE2.ExecutionEnvironment.PhysicalCPUs;"
    "]")
);

} // {anonymous}

void populate_ism(
  vector<gluece_info_iterator>& gluece_info_container_updated_entries,
  size_t the_ism_index, // typically ce or se
  update_function_type const& uf
  )
{
  static glite::wms::common::configuration::Configuration const& config(
    *glite::wms::common::configuration::Configuration::instance()
  );
  static const time_t expiry_time(
    config.wm()->ism_ii_purchasing_rate() + config.ns()->ii_timeout()
  );

  vector<gluece_info_iterator>::const_iterator it(
    gluece_info_container_updated_entries.begin()
  );
  vector<gluece_info_iterator>::const_iterator const e(
    gluece_info_container_updated_entries.end()
  );

  // no locking is needed here (before the switch)
  int dark_side = ism::dark_side();

  bool insert = false;
  ism_type::iterator ism_entry;

  for ( ; it != e; ++it ) {
    boost::tie(ism_entry, insert) =
      get_ism(the_ism_index, dark_side).insert(
        make_ism_entry(
          (*it)->first,
          std::time(0),
          (*it)->second,
          uf,
          expiry_time
        )
      );
    if (!insert) { // existing entry (typically glue13), need to merge info
      ism_entry->second.get<2>()->Update(*(*it)->second);
      ism_entry->second.get<2>()->Update(*(glue13_mapping_ad));

      Debug((*it)->first << " updated with GLUE2.0");
    }
    else {
       Debug((*it)->first << " added to ISM ");
    }
  }
}

void apply_skip_predicate(
  gluece_info_container_type& gluece_info_container,
  vector<gluece_info_iterator>& gluece_info_container_updated_entries,
  skip_predicate_type skip,
  std::string const& purchasedby
)
{
  gluece_info_iterator it = gluece_info_container.begin();
  gluece_info_iterator const gluece_info_container_end(
    gluece_info_container.end()
  );

  for ( ; it != gluece_info_container_end; ++it) {
    if (!skip(it->first)) {
      it->second->InsertAttr("PurchasedBy", purchasedby);
      gluece_info_container_updated_entries.push_back(it);
    } else {
      Debug("Skipping " << it->first << " due to skip predicate settings");
    }
  }
}

void tokenize_ldap_dn(std::string const& s, std::vector<std::string> &v)
{
  boost::escaped_list_separator<char> ldap_dn_sep("",",","");
  boost::tokenizer<boost::escaped_list_separator<char> >
    ldap_dn_tok(s,ldap_dn_sep);

  boost::tokenizer< boost::escaped_list_separator<char> >::iterator
    ldap_dn_tok_it(
      ldap_dn_tok.begin()
    );
  boost::tokenizer< boost::escaped_list_separator<char> >::iterator const
    ldap_dn_tok_end(
       ldap_dn_tok.end()
    );

  for( ; ldap_dn_tok_it != ldap_dn_tok_end; ++ldap_dn_tok_it)
    v.push_back(
      boost::algorithm::trim_copy(*ldap_dn_tok_it)
    );

}

void cleanup_glue_info(ad_ptr ad, std::list<std::string> a)
{
  a.push_back("objectClass");
  std::for_each(
   a.begin(), a.end(),
   boost::bind(&classad::ClassAd::Delete, ad, _1)
  );
}

namespace {

   std::string const gangmatch_storage_ad_str(
    "["
    "  storage =  ["
    "     VO = parent.other.VirtualOrganisation;"
    "     CloseSEs = retrieveCloseSEsInfo( VO );"
    "  ];"
    "]"
  );

  boost::scoped_ptr<classad::ClassAd> gangmatch_storage_ad;
}

bool expand_information_service_info(gluece_info_type& gluece_info)
{
  string isURL;
  bool result = false;
  try {

    isURL = static_cast<std::string>(
      utils::evaluate_attribute(*gluece_info, "GlueInformationServiceURL")
    );
    static boost::regex  expression_gisu("\\S.*://(.*):([0-9]+)/(.*)");
    boost::smatch        pieces_gisu;

    if (boost::regex_match(isURL, pieces_gisu, expression_gisu)) {

      string ishost(pieces_gisu[1].first, pieces_gisu[1].second);
      string isport(pieces_gisu[2].first, pieces_gisu[2].second);
      string isbasedn(pieces_gisu[3].first, pieces_gisu[3].second);

      gluece_info->InsertAttr("InformationServiceDN", isbasedn);
      gluece_info->InsertAttr("InformationServiceHost", ishost);
      gluece_info->InsertAttr("InformationServicePort",
                              boost::lexical_cast<int>(isport));
      result = true;
    }
  } catch (utils::InvalidValue& e) {
    Error("Cannot evaluate GlueInformationServiceURL...");
    result = false;
  }
  return result;
}

bool insert_gangmatch_storage_ad(gluece_info_type& gluece_info)
{
 try {
    if(!gangmatch_storage_ad) {
      gangmatch_storage_ad.reset(
        utils::parse_classad(gangmatch_storage_ad_str)
      );
    }
    gluece_info->Update( *gangmatch_storage_ad);
  }
  catch(...) {
    assert(false);
  }
  return true;
}

bool expand_glueceid_info(gluece_info_type& gluece_info)
{
  string ce_str;
  ce_str.assign(utils::evaluate_attribute(*gluece_info, "GlueCEUniqueID"));
  static boost::regex  expression_ceid("(.+/[^\\-]+-([^\\-]+))-(.+)");
  boost::smatch  pieces_ceid;
  string gcrs, type, name;

  if (boost::regex_match(ce_str, pieces_ceid, expression_ceid)) {

    gcrs.assign(pieces_ceid[1].first, pieces_ceid[1].second);
    try {
      type.assign(utils::evaluate_attribute(*gluece_info, "GlueCEInfoLRMSType"));
    }
    catch(utils::InvalidValue& e) {
      // Try to fall softly in case the attribute is missing...
      type.assign(pieces_ceid[2].first, pieces_ceid[2].second);
      Warning("Cannot evaluate GlueCEInfoLRMSType using value from contact string: " << type);
    }
    // ... or in case the attribute is empty.
    if (type.length() == 0) type.assign(pieces_ceid[2].first, pieces_ceid[2].second);
    name.assign(pieces_ceid[3].first, pieces_ceid[3].second);
  }
  else {
    Warning("Cannot parse CEid=" << ce_str);
    return false;
  }
  gluece_info -> InsertAttr("GlobusResourceContactString", gcrs);
  gluece_info -> InsertAttr("LRMSType", type);
  gluece_info -> InsertAttr("QueueName", name);
  gluece_info -> InsertAttr("CEid", ce_str);

  return true;
}

bool split_information_service_url(classad::ClassAd const& ad, boost::tuple<std::string, int, std::string>& i)
{
 try {
  std::string ldap_dn;
  std::string ldap_host;
  std::string ldap_url;

  ldap_url.assign( utils::evaluate_attribute(ad, "GlueInformationServiceURL") );
  static boost::regex expression_gisu( "\\S.*://(.*):([0-9]+)/(.*)" );
  boost::smatch pieces_gisu;
  std::string port;

  if (boost::regex_match(ldap_url, pieces_gisu, expression_gisu)) {

    ldap_host.assign (pieces_gisu[1].first, pieces_gisu[1].second);
    port.assign      (pieces_gisu[2].first, pieces_gisu[2].second);
    ldap_dn.assign   (pieces_gisu[3].first, pieces_gisu[3].second);

    i = boost::make_tuple(ldap_host, std::atoi(port.c_str()), ldap_dn);
  }
  else {
  return false;
  }
 }
 catch (utils::InvalidValue& e) {
   return false;
 }
 return true;
}

}}}}

max_stars_repo_name: italiangrid/wms

max_issues_repo_name: italiangrid/wms

max_forks_repo_name: italiangrid/wms


text:
#include "RobotService.h"
#include <boost/thread.hpp>

#include "deps/easylogging++.h"
#include <iface_msg.hpp>

namespace {
  /**
   * Callback invoked when the async connect operation completes.
   */
  void connect_handler(boost::system::error_code const& error) {
    if (!error) {
      LINFO << "AsyncRobotService: Connected to remote host.";
    } else {
      LERROR << "AsyncRobotService:  Failed to connect to the remote host.";
    }
  }

  /**
   * A simple function that is used to spin up the io service event loop in a
   * dedicated thread.
   * Makes sure that the loop does not end when the service does not have any
   * work to do (temporarily).
   */
  void service_thread(boost::asio::io_service* io_service) {
    // Prevent the IO service from running out of work (and exitting when no
    // async operations are queued).
    boost::asio::io_service::work work(*io_service);
    io_service->run();
    LINFO << "AsyncRobotService: Exiting service thread...";
  }

  /**
   * A callback for write operations.
   */
  void write_handler(boost::system::error_code const& error,
               std::size_t sent) {
    if (!error) {
      LINFO << "AsyncRobotService: Send complete. "
            << "Sent " << sent << " bytes.";
    } else {
      LERROR << "AsyncRobotService: Error sending message.";
    }
  }
}

void AsyncRobotService::start() {
  // Start it up...
  boost::asio::ip::tcp::endpoint endpoint(
    boost::asio::ip::address::from_string(remote_), port_);
  LINFO << "AsyncRobotService (" << remoteName_ << "): Initiating a connection asynchronously...";
  socket_.async_connect(endpoint, connect_handler);
  // Start the service thread in the background...
  boost::thread(boost::bind(service_thread, &io_service_));
}

void AsyncRobotService::inner_send(VisionMessage const& next_message) {
  am2b_iface::MsgHeader msg_header = { am2b_iface::VISION_MESSAGE, (uint32_t)sizeof(VisionMessageHeader) + next_message.header.len };
//  LINFO << "AsyncRobotService (" << remoteName_ << "): Sending a queued message: "
//        << "msg == " << next_message;
  // Synchronously send the message, i.e. block until the send is complete.
  try {
    sendBytes(reinterpret_cast<char const*>(&msg_header), sizeof(msg_header));
    sendBytes(reinterpret_cast<char const*>(&next_message), sizeof(VisionMessageHeader));
    sendBytes(next_message.content, next_message.header.len);
  } catch (...) {
    LERROR << "AsyncRobotService (" << remoteName_ << "): Error sending message.";
  }
  // After each sent message, we want to wait a pre-defined amount of time
  // before sending the next one.
  // This is because we do not want to overwhelm the robot with a large
  // number of messages all sent in the same time.
  boost::this_thread::sleep(message_timeout_);
}

void AsyncRobotService::sendBytes(char const* data, size_t length) {
  size_t sent_bytes = 0;

  do {
    sent_bytes += socket_.send(boost::asio::buffer(&data[sent_bytes], length - sent_bytes));
  } while (sent_bytes < length);
}

void AsyncRobotService::sendMessage(VisionMessage const& msg) {
  // Just queue another message to be sent by the io_service thread.
  // Since `inner_send` makes sure to wait after it's finished sending
  // each message (i.e. the io_service thread is blocked), the queued
  // messages will all be sent with a preset time delay between subsequent
  // messages.
  io_service_.post(boost::bind(&AsyncRobotService::inner_send, this, msg));
}