Merge pull request #12 from GooFit/experimental_integrator2

Experimental integrator2

.ci/build_root.sh

@@ -0,0 +1,20 @@
+pushd $DEPS_DIR 
+set -evx
+
+os=$1
+if [[ $os == "osx" ]] ; then 
+  wget -nv http://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda_${os}.sh
+elif [[ $os == "linux" ]] ; then 
+  wget -nv http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda_${os}.sh
+fi
+
+bash miniconda_${os}.sh -b -p $DEPS_DIR/miniconda
+export PATH="$DEPS_DIR/miniconda/bin:$PATH"
+hash -r
+conda config --add channels conda-forge
+# conda config --set channel_priority strict
+
+conda create --yes -n env_${os} root doxygen -c conda-forge
+
+set +evx 
+popd

.ci/build_root_linux.sh

@@ -1,12 +1,20 @@
 pushd $DEPS_DIR 
 
-wget -nv http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
+if [[ $1 == "osx" ]] ; then 
+  wget -nv http://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda.sh
+elif [[ $1 == "linux" ]] ; then 
+  wget -nv http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
+fi
+
 bash miniconda.sh -b -p $DEPS_DIR/miniconda
 export PATH="$DEPS_DIR/miniconda/bin:$PATH"
 hash -r
 conda config --add channels conda-forge
-conda install --quiet --yes -c conda-forge/label/gcc8 root 
+conda config --set channel_priority strict
 
-source "$DEPS_DIR/miniconda/bin/thisroot.sh"
-export CXX="$DEPS_DIR/miniconda/bin/g++"
+# conda install --quiet --yes -c conda-forge/label/gcc8 root_base doxygen
+conda create --quiet --yes -n my_root_env root doxygen zstd=1.3.7 -c conda-forge
+# conda init bash
+# source "$DEPS_DIR/miniconda/bin/thisroot.sh"
+# export CXX="$DEPS_DIR/miniconda/bin/g++"
 popd

.ci/travis_linux.sh

@@ -1,24 +1,23 @@
+#!/bin/bash
+
+set -evx
+# from https://stackoverflow.com/questions/55342122/conda-activate-on-travis-ci
+export PATH="$DEPS_DIR/miniconda/bin:$PATH"
+. $(conda info --root)/etc/profile.d/conda.sh
+conda activate env_${TRAVIS_OS_NAME}
+
 echo -en 'travis_fold:start:script.build\\r'
 echo "Building..."
 echo "Building under OS: $TRAVIS_OS_NAME, CXX =$CXX"
-set -evx
 
-mkdir -p build
-cd build
-cmake .. -DCMAKE_CXX_COMPILER=$CXX
+mkdir -p build.conda
+cd build.conda
+cmake .. 
 cmake --build . -- -j2
 
 set +evx
 
 cd ..
 
-./build/bin/Generator options/example_b2kstarll.opt --CompilerWrapper::Verbose --nEvents 1000
-
-
-# echo -e 'travis_fold:end:script.build\\r'
-# echo -en 'travis_fold:start:script.test\\r'
-# echo "Testing..."
-# set -evx
-
-# ctest --output-on-failure
+./build.conda/bin/Generator options/example_b2kstarll.opt --CompilerWrapper::Verbose --nEvents 10000
 

.ci/travis_osx.sh

@@ -1,22 +1,18 @@
+#!/bin/bash
+
+. $(conda info --root)/etc/profile.d/conda.sh
+conda activate env_${TRAVIS_OS_NAME}
+
 echo -en 'travis_fold:start:script.build\\r'
 echo "Building..."
 echo "Building under OS: $TRAVIS_OS_NAME"
 
 mkdir -p build
 cd build
 echo "CMake-ing, CXX = $CXX"
-cmake .. -DCMAKE_CXX_COMPILER=clang
+cmake .. -DCMAKE_CXX_COMPILER=clang -DUSE_SIMD="" -DUSE_OPENMP=0
 echo "Building ..."
 cmake --build . -- -j2 
 cd ..
 echo "Running test job ..."
-./build/bin/Generator options/example_b2kstarll.opt --CompilerWrapper::Verbose --nEvents 1000
-
-
-# echo -e 'travis_fold:end:script.build\\r'
-# echo -en 'travis_fold:start:script.test\\r'
-# echo "Testing..."
-# set -evx
-
-# ctest --output-on-failure
-
+./build/bin/Generator options/example_b2kstarll.opt --CompilerWrapper::Verbose --nEvents 10000

.travis.yml

@@ -10,17 +10,17 @@ matrix:
           packages:
             - libomp
     - os: linux
-      addons:
-       apt:
-         sources:
-           - ubuntu-toolchain-r-test
-         packages:
-           - g++-8
-           - doxygen
-           - doxygen-doc
-           - doxygen-gui
-           - graphviz
-           - libtbb-dev
+      #      addons:
+      #       apt:
+      #         sources:
+      #           - ubuntu-toolchain-r-test
+      #         packages:
+      #           - g++-8
+      #           - doxygen
+      #           - doxygen-doc
+      #           - doxygen-gui
+      #           - graphviz
+      #           - libtbb-dev
       env:
         - MATRIX_EVAL="CC=gcc-8 && CXX=g++-8"
 
@@ -48,8 +48,8 @@ env:
 
 before_install:
   - eval "${MATRIX_EVAL}"
-  - chmod +x .ci/build_root_${TRAVIS_OS_NAME}.sh
-  - source .ci/build_root_${TRAVIS_OS_NAME}.sh 
+  - chmod +x .ci/build_root.sh 
+  - source .ci/build_root.sh ${TRAVIS_OS_NAME} 
 
 # Build your code e.g. by calling make
 script:

AmpGen/ASTResolver.h

@@ -36,14 +36,16 @@ namespace AmpGen {
       {
         auto it = m_cacheFunctions.find(name);
         if( it != m_cacheFunctions.end() ) return it->second->address();
-        auto cacheFunction = std::make_shared<TYPE>(m_nParameters, args... );
+        auto cacheFunction = std::make_shared<TYPE>(m_nParameters, name, args... );
         m_cacheFunctions[name] = cacheFunction;
         m_nParameters += cacheFunction->size();
         return m_nParameters - cacheFunction->size();
       }
       size_t nParams() const { return m_nParameters ; }       
       bool enableCuda() const { return m_enable_cuda ; }
+      bool enableAVX()  const { return m_enable_avx; }
       bool enableCompileConstants() const { return m_enable_compileTimeConstants ;} 
+      void setEnableAVX(){ m_enable_avx = true ; }
       std::map<std::string, std::shared_ptr<CacheTransfer>> cacheFunctions() const;
       void addResolvedParameter(const IExpression* param, const std::string& thing);
       void addResolvedParameter(const IExpression* param, const size_t& address, const size_t& arg=0);
@@ -58,10 +60,12 @@ namespace AmpGen {
       std::map<std::string, unsigned>                       m_evtMap;                      /// Event specification 
       std::map<std::string, std::string>                    m_parameterMapping;            /// Mapping of parameters to compile parameters
       const MinuitParameterSet*                             m_mps;                         /// Set of MinuitParameters 
-      std::map<const SubTree*, uint64_t>                    m_tempTrees;                   /// temporary store of sub-trees for performing cse reduction 
+      std::map<const IExpression*, const SubTree*>          m_tempTrees;                   /// temporary store of sub-trees for performing cse reduction 
       unsigned int                                          m_nParameters;                 /// Number of parameters
-      bool                                                  m_enable_cuda;                 /// flag to generate CUDA code <<experimental>>
-      bool                                                  m_enable_compileTimeConstants; /// flag to enable compile time constants <<experimental>> 
+      bool                                                  m_enable_cuda                 {false}; /// flag to generate CUDA code <<experimental>>
+      bool                                                  m_enable_compileTimeConstants {false}; /// flag to enable compile time constants <<experimental>> 
+      bool                                                  m_enable_avx                  {false}; /// flag to generate code using AVX instructions <<experimental>>
+      bool                                                  m_check_hashes                {false}; /// flag to check that hashes are unique 
   };
 
   template <> void ASTResolver::resolve<Parameter>( const Parameter& obj );

AmpGen/AmplitudeRules.h

@@ -14,6 +14,9 @@
 #include "AmpGen/Event.h"
 #include "AmpGen/Particle.h"
 #include "AmpGen/ExpressionParser.h"
+#if ENABLE_AVX 
+  #include "AmpGen/EventListSIMD.h"
+#endif
 
 namespace AmpGen
 {
@@ -73,7 +76,7 @@ namespace AmpGen
       AmplitudeRules( const MinuitParameterSet& mps );
       std::vector<Coupling> rulesForDecay(const std::string& head, const std::string& prefix="");
       bool hasDecay( const std::string& head );
-      std::map<std::string, std::vector<Coupling>> rules();
+      const std::map<std::string, std::vector<Coupling>>& rules() const;
       std::vector<std::pair<Particle, TotalCoupling>> getMatchingRules( 
           const EventType& type, const std::string& prefix="" );
       std::vector<Coupling> processesThatProduce(const Particle& particle) const; 
@@ -82,39 +85,48 @@ namespace AmpGen
       std::map<std::string, std::vector<Coupling>> m_rules;
   };
 
-  template <class RT> struct TransitionMatrix 
+  template <class RT> struct TransitionMatrix : public CompiledExpression<RT(const real_t*, const float_v*)>  
   {
+    using amp_type = CompiledExpression<RT(const real_t*, const float_v*)>; 
     TransitionMatrix() = default;
     TransitionMatrix(const Particle& dt, 
                      const TotalCoupling& coupling, 
-                     const CompiledExpression<RT, const real_t*, const real_t*> & amp) : 
+                     const amp_type& amp) :
+          amp_type(amp),
           decayTree(dt), 
-          coupling(coupling), 
-          amp(amp) {}
+          coupling(coupling) {}
 
-    TransitionMatrix(Particle& dt, 
+    TransitionMatrix(const Particle& dt, 
                      const TotalCoupling& coupling, 
                      const MinuitParameterSet& mps,
                      const std::map<std::string, unsigned>& evtFormat, 
                      const bool& debugThis=false) :
+      amp_type(Particle(dt).getExpression(debugThis ? &db : nullptr ), dt.decayDescriptor(), evtFormat, db, &mps ),
       decayTree(dt),
-      coupling(coupling),
-      amp(decayTree.getExpression(debugThis ? &db : nullptr ), decayTree.decayDescriptor(), evtFormat, db, &mps ) {}
+      coupling(coupling) {}
 
-    const RT operator()(const Event& event) const { return amp(event.address() ); }
-    const RT operator()(const Event& event, const size_t& cacheOffset) const { return amp(event.address() + cacheOffset); }
+    #if ENABLE_AVX
+    const RT operator()(const Event& event) const { return amp_type::operator()(EventListSIMD::makeEvent(event).data()); }
+    void debug( const Event& event ) const {               amp_type::debug(EventListSIMD::makeEvent(event).data() ) ; } 
+
+    #else
+    const RT operator()(const Event& event) const { return amp_type::operator()(event.address()) ; }
+    void debug( const Event& event )        const {        amp_type::debug(event.address()) ; }
+    #endif
+    template <class... arg_types> auto operator()(arg_types... args ) const { return amp_type::operator()(args...) ; }
+
+    const RT operator()(const float_v* t) const     { return amp_type::operator()(t) ; }
+    void debug( const float_v* t )        const     {        amp_type::debug(t) ; } 
     const std::string decayDescriptor() const { return decayTree.decayDescriptor() ; }  
 
     Particle                                            decayTree;
     TotalCoupling                                       coupling;
     complex_t                                           coefficient;
     DebugSymbols                                        db; 
-    CompiledExpression<RT,const real_t*,const real_t*>  amp; 
-    size_t                                              addressData = {999};
+    bool                                                workToDo    = {false};
   };
 
-  template <class RT>  
-  std::vector<size_t> processIndex(const std::vector<TransitionMatrix<RT>>& tm, const std::string& label)
+  template <class RT> std::vector<size_t> processIndex(const std::vector<TransitionMatrix<RT>>& tm, const std::string& label)
   {
     std::vector<size_t> indices;
     for ( size_t i = 0; i < tm.size(); ++i ) {
@@ -123,8 +135,7 @@ namespace AmpGen
     return indices;
   }
 
-  template <class RT>
-  size_t findIndex(const std::vector<TransitionMatrix<RT>>& tm, const std::string& decayDescriptor)
+  template <class RT> size_t findIndex(const std::vector<TransitionMatrix<RT>>& tm, const std::string& decayDescriptor)
   {
     for ( size_t i = 0; i < tm.size(); ++i ) {
       if ( tm[i].decayDescriptor() == decayDescriptor ) return i;
@@ -143,43 +154,51 @@ namespace AmpGen
     return rt;
   }
 
-  template <> struct TransitionMatrix<void>
+  template <> struct TransitionMatrix<void> : public CompiledExpression<void(complex_v*, const size_t&, const real_t*, const float_v*)>
   {
+    using amp_type    = CompiledExpression<void(complex_v*, const size_t&, const real_t*, const float_v*)>;
     TransitionMatrix() = default;
     TransitionMatrix(const Particle& dt, 
                      const TotalCoupling& coupling, 
-                     const CompiledExpression<void, complex_t*, const real_t*, const real_t*> & amp) : 
+                     const amp_type& amp) : 
+          amp_type(amp),
           decayTree(dt), 
-          coupling(coupling), 
-          amp(amp) {}
+          coupling(coupling) {}
 
     TransitionMatrix(const Particle& dt, 
                      const TotalCoupling& coupling, 
                      const MinuitParameterSet& mps,
                      const std::map<std::string, unsigned>& evtFormat, 
                      const bool& debugThis=false) :
+      amp_type(Particle(dt).getExpression(debugThis ? &db : nullptr ), dt.decayDescriptor(), evtFormat, db, &mps ),
       decayTree(dt),
-      coupling(coupling),
-      amp(decayTree.getExpression(debugThis ? &db : nullptr ), decayTree.decayDescriptor(), evtFormat, db, &mps ) { amp.use_rto();}
+      coupling(coupling)
+      { use_rto();}
 
-    const std::vector<complex_t> operator()(const Event& event) const { 
-      std::vector<complex_t> rt; 
-      amp(rt.data(), amp.externBuffer().data(), event.address() ); 
-      return rt;
-    }
-    const std::vector<complex_t> operator()(const Event& event, const size_t& cacheOffset) const { 
-      std::vector<complex_t> rt; 
-      amp(rt.data(), amp.externBuffer().data(), event.address() + cacheOffset); 
+    const std::vector<complex_v> operator()(const Event& event) const 
+    { 
+      std::vector<complex_v> rt(size); 
+      #if ENABLE_AVX 
+      amp_type::operator()(rt.data(), 1, externBuffer().data(), EventListSIMD::makeEvent(event).data());
+      #else
+      amp_type::operator()(rt.data(), 1, externBuffer().data(), event.address()); 
+      #endif
       return rt;
     }
+    template <class... arg_types> auto operator()(arg_types... args ) const { return amp_type::operator()(args...) ; }
+    #if ENABLE_AVX
+    void debug( const Event& event ) const { amp_type::debug(EventListSIMD::makeEvent(event).data() ) ; } 
+    #else 
+    void debug( const Event& event ) const { amp_type::debug(event.address()) ; }
+    #endif 
     const std::string decayDescriptor() const { return decayTree.decayDescriptor() ; }  
 
     Particle                                            decayTree;
     TotalCoupling                                       coupling;
     complex_t                                           coefficient;
     DebugSymbols                                        db; 
-    CompiledExpression<void, complex_t*, const real_t*, const real_t*>  amp; 
-    size_t                                              addressData = {999};
+    bool                                                workToDo    = {false};
+    unsigned                                            size        = {0};
   };
 
 } // namespace AmpGen

AmpGen/ArgumentPack.h

@@ -12,7 +12,7 @@ namespace AmpGen
   #define DECLARE_ARGUMENT(X, Y)                          \
   struct X : public AmpGen::Argument<Y> {                 \
     template<class Z>                                     \
-    explicit X(Z val) : AmpGen::Argument<Y>(val){}        \
+    explicit X(Z val = Z()) : AmpGen::Argument<Y>(val){}  \
     X() : AmpGen::Argument<Y>(){}                         \
   }
   /** @class IArgument 
@@ -76,14 +76,20 @@ namespace AmpGen
           std::tuple<ARGS...> argTuple( args... );
           for_each(argTuple, [this](const auto& f){ this->addArgument(f) ; } );
         }
-      template <typename ARG, typename DEFAULT_TYPE=ARG> 
-        ARG getArg( const DEFAULT_TYPE& default_argument = DEFAULT_TYPE() ) const
+      template <typename arg_type> arg_type* get() const 
       {
-        for ( auto param : m_parameters ) {
-          auto ptr = dynamic_cast<ARG*>( param.get() );
-          if ( ptr != nullptr ) return *ptr;
+        for( const auto& param : m_parameters )
+        {
+          auto ptr = dynamic_cast<arg_type*>(param.get());
+          if( ptr != nullptr ) return ptr;
         }
-        return ARG(default_argument);
+        return nullptr;
+      }
+      template <typename arg_type, typename default_arg_type=arg_type> 
+        arg_type getArg( const default_arg_type& default_argument = default_arg_type() ) const
+      {
+        auto p = get<arg_type>();
+        return p == nullptr ? arg_type(default_argument) : *p;
       }
     private:
       std::vector<std::shared_ptr<IArgument>> m_parameters;

AmpGen/Array.h

@@ -31,11 +31,11 @@ namespace AmpGen
       complex_t operator()() const override;
       Expression operator[]( const Expression& address ) const;
       Expression top() const { return m_top ; } 
-
+      unsigned size() const { return m_size; }
     private:
       Expression m_top; 
       Expression m_address;
-      size_t     m_size; 
+      unsigned   m_size; 
   };
 } // namespace AmpGen