From 2daace08d3cc69f769f038070c9c5ab4e2064db4 Mon Sep 17 00:00:00 2001
From: Bertram Kopf <bertram@ep1.rub.de>
Date: Mon, 21 Aug 2023 17:36:08 +0200
Subject: [PATCH] added autodiff tool Adept in test app ReverseCollections
---
.../Tutorial/LineShapes/BarrierFactorApp.cc | 3 +-
Examples/Tutorial/LineShapes/BwShape.cc | 3 +-
JamrootAL9AutoDiff | 5 +
PwaDynamics/BarrierFactor.cc | 4 +
PwaDynamics/BreitWignerFunction.cc | 4 +
PwaDynamics/Jamfile | 2 +
SetEnv_rub_AutoDiff | 2 +-
autodiffExamples/ReverseCollections.cc | 137 +++++++++++++-----
autodiffExamples/ReverseCollections.hh | 11 +-
qft++Extension/PawianUtils.cc | 5 +
10 files changed, 132 insertions(+), 44 deletions(-)
diff --git a/Examples/Tutorial/LineShapes/BarrierFactorApp.cc b/Examples/Tutorial/LineShapes/BarrierFactorApp.cc
index b2189fe3..09a741c5 100644
--- a/Examples/Tutorial/LineShapes/BarrierFactorApp.cc
+++ b/Examples/Tutorial/LineShapes/BarrierFactorApp.cc
@@ -359,7 +359,8 @@ int main(int __argc, char *__argv[]) {
double qProd = PawianQFT::breakupMomQDefault(motherMass, recoilMass, dHist->GetBinCenter(i)).real();
double qNorm=PawianQFT::breakupMomQDefault(motherMass, recoilMass, decayParticle1Mass+decayParticle2Mass).real();
std::complex<double> valueProd = BarrierFactor::BlattWeisskopf(lp, qProd, qRprod)/BarrierFactor::BlattWeisskopf(lp, qNorm, qRprod);
- std::complex<double> breitWigner = BreitWignerFunction::BlattWRel(ld, dHist->GetBinCenter(i), resonanceMass,
+ double currentBinCenter=dHist->GetBinCenter(i);
+ std::complex<double> breitWigner = BreitWignerFunction::BlattWRel(ld, currentBinCenter, resonanceMass,
resonanceWidth,decayParticle1Mass, decayParticle2Mass, qRdec );
dHist->SetBinContent(i, std::norm(breitWigner * valueProd));
diff --git a/Examples/Tutorial/LineShapes/BwShape.cc b/Examples/Tutorial/LineShapes/BwShape.cc
index 5365a064..e350d7f3 100644
--- a/Examples/Tutorial/LineShapes/BwShape.cc
+++ b/Examples/Tutorial/LineShapes/BwShape.cc
@@ -95,7 +95,8 @@ BwShape::BwShape(double MassRes, double MassWidth, double MassDec1, double MassD
Vector4<double> res4V(massIt, 0., 0., 0.);
// complex<double> currentBW=BreitWignerBlattW(res4V, MassDec1, MassDec2, MassRes, MassWidth, lIt);
- complex<double> currentBWnew=BreitWignerFunction::BlattWRel(lIt, res4V.Mass(), MassRes, MassWidth, MassDec1, MassDec2);
+ double currentres4VMass=res4V.Mass();
+ complex<double> currentBWnew=BreitWignerFunction::BlattWRel(lIt, currentres4VMass, MassRes, MassWidth, MassDec1, MassDec2);
// double weight=norm(currentBW);
// currentHist->Fill(massIt,weight);
currentHistNew->Fill(massIt,norm(currentBWnew));
diff --git a/JamrootAL9AutoDiff b/JamrootAL9AutoDiff
index e48f8978..5e14b79e 100644
--- a/JamrootAL9AutoDiff
+++ b/JamrootAL9AutoDiff
@@ -32,6 +32,7 @@ project :
<include>/data/jollyj/bertram/autoDiff/autodiff
<include>/data/jollyj/bertram/eigen/eigen3-d0923f82df9296ac8d880a4d7f661c3b6f12077b/
<include>/data/jollyj/bertram/CoDiPack/CoDiPack/include
+ <include>/data/jollyj/bertram/Adept/AdeptInst/include
<cxxflags>--std=c++17
<cxxflags>-ftemplate-depth=256
<cxxflags>-DBOOST_BIND_GLOBAL_PLACEHOLDERS
@@ -40,6 +41,10 @@ project :
<linkflags>$(ROOTLIBS)
<linkflags>-L$(BOOSTLIBPATH)
<linkflags>$(BOOSTLIBS)
+ <linkflags>-L/data/jollyj/bertram/Adept/AdeptInst/lib
+ <linkflags>-ladept
+# <linkflags>-L/data/jollyj/bertram/Adept/OpenBlasInst/lib/
+# <linkflags>-lopenblas
<linkflags>-lgsl
<linkflags>-lgslcblas
<cxxflags>-fPIC
diff --git a/PwaDynamics/BarrierFactor.cc b/PwaDynamics/BarrierFactor.cc
index 91a62795..0cafc1f6 100644
--- a/PwaDynamics/BarrierFactor.cc
+++ b/PwaDynamics/BarrierFactor.cc
@@ -25,6 +25,7 @@
//#include <autodiff/forward/dual.hpp>
#include <autodiff/reverse/var.hpp>
#include <codi.hpp>
+#include "adept.h"
#include "PwaDynamics/BarrierFactor.hh"
#include "ErrLogger/ErrLogger.hh"
@@ -54,6 +55,7 @@ std::complex<FloatP> BarrierFactor::BlattWeisskopf(std::complex<FloatP> q){
template std::complex<double> BarrierFactor::BlattWeisskopf(std::complex<double>);
template std::complex<autodiff::var> BarrierFactor::BlattWeisskopf(std::complex<autodiff::var>);
template std::complex<codi::RealReverse> BarrierFactor::BlattWeisskopf(std::complex<codi::RealReverse>);
+template std::complex<adept::adouble> BarrierFactor::BlattWeisskopf(std::complex<adept::adouble>);
template<typename FloatP>
std::complex<FloatP> BarrierFactor::BlattWeisskopf(int l, std::complex<FloatP> q, double qR){
@@ -64,6 +66,7 @@ std::complex<FloatP> BarrierFactor::BlattWeisskopf(int l, std::complex<FloatP> q
}
template std::complex<double> BarrierFactor::BlattWeisskopf(int, std::complex<double>, double);
template std::complex<autodiff::var> BarrierFactor::BlattWeisskopf(int, std::complex<autodiff::var>, double);
+template std::complex<adept::adouble> BarrierFactor::BlattWeisskopf(int, std::complex<adept::adouble>, double);
std::complex<double> BarrierFactor::BlattWeisskopfTensor(int l, std::complex<double> q, double qR){
double z(std::norm(q)/std::norm(qR));
@@ -135,6 +138,7 @@ std::complex<FloatP> BarrierFactor::BlattWeisskopfRatio(int l, std::complex<Floa
template std::complex<double> BarrierFactor::BlattWeisskopfRatio(int, std::complex<double>, std::complex<double>, double);
template std::complex<autodiff::var> BarrierFactor::BlattWeisskopfRatio(int, std::complex<autodiff::var>, std::complex<autodiff::var>, double);
template std::complex<codi::RealReverse> BarrierFactor::BlattWeisskopfRatio(int, std::complex<codi::RealReverse>, std::complex<codi::RealReverse>, double);
+template std::complex<adept::adouble> BarrierFactor::BlattWeisskopfRatio(int, std::complex<adept::adouble>, std::complex<adept::adouble>, double);
std::complex<double> BarrierFactor::BlattWeisskopfTensorRatio(int l, std::complex<double> q, std::complex<double> q0, double qR){
std::complex<double> denom=BlattWeisskopfTensor(l, q0, qR);
diff --git a/PwaDynamics/BreitWignerFunction.cc b/PwaDynamics/BreitWignerFunction.cc
index 4f17b64f..58880163 100644
--- a/PwaDynamics/BreitWignerFunction.cc
+++ b/PwaDynamics/BreitWignerFunction.cc
@@ -29,6 +29,7 @@
#include <autodiff/reverse/var.hpp>
//#include <autodiff/reverse/var/eigen.hpp>
#include <codi.hpp>
+#include "adept.h"
template<typename FloatP>
std::complex<FloatP> BreitWignerFunction::NonRel(FloatP currentMass,FloatP mass0, FloatP width){
@@ -42,6 +43,7 @@ std::complex<FloatP> BreitWignerFunction::NonRel(FloatP currentMass,FloatP mass0
template std::complex<double> BreitWignerFunction::NonRel(double currentMass,double mass0, double width);
template std::complex<autodiff::var> BreitWignerFunction::NonRel(autodiff::var currentMass,autodiff::var mass0, autodiff::var width);
template std::complex<codi::RealReverse> BreitWignerFunction::NonRel(codi::RealReverse currentMass, codi::RealReverse mass0, codi::RealReverse width);
+template std::complex<adept::adouble> BreitWignerFunction::NonRel(adept::adouble currentMass, adept::adouble mass0, adept::adouble width);
template<typename FloatP>
@@ -57,6 +59,7 @@ complex<FloatP> BreitWignerFunction::Rel(FloatP currentMass, FloatP mass0, Float
template std::complex<double> BreitWignerFunction::Rel(double, double, double, double, double);
template std::complex<autodiff::var> BreitWignerFunction::Rel(autodiff::var, autodiff::var, autodiff::var, autodiff::var, autodiff::var);
template std::complex<codi::RealReverse> BreitWignerFunction::Rel(codi::RealReverse, codi::RealReverse, codi::RealReverse, codi::RealReverse, codi::RealReverse);
+template std::complex<adept::adouble> BreitWignerFunction::Rel(adept::adouble, adept::adouble, adept::adouble, adept::adouble, adept::adouble);
template<typename FloatP>
complex<FloatP> BreitWignerFunction::BlattWRel(int orbMom, FloatP currentMass, FloatP mass0, FloatP width, FloatP massA, FloatP massB, double qR){
@@ -76,6 +79,7 @@ complex<FloatP> BreitWignerFunction::BlattWRel(int orbMom, FloatP currentMass,
template std::complex<double> BreitWignerFunction::BlattWRel(int, double, double, double, double, double, double);
template std::complex<autodiff::var> BreitWignerFunction::BlattWRel(int, autodiff::var, autodiff::var, autodiff::var, autodiff::var, autodiff::var, double);
template std::complex<codi::RealReverse> BreitWignerFunction::BlattWRel(int, codi::RealReverse, codi::RealReverse, codi::RealReverse, codi::RealReverse, codi::RealReverse, double);
+template std::complex<adept::adouble> BreitWignerFunction::BlattWRel(int, adept::adouble, adept::adouble, adept::adouble, adept::adouble, adept::adouble, double);
complex<double> BreitWignerFunction::BlattWTensorRel(int orbMom, double currentMass,double mass0, double width, double massA, double massB, double qR){
diff --git a/PwaDynamics/Jamfile b/PwaDynamics/Jamfile
index 1c6953ad..fb62590b 100644
--- a/PwaDynamics/Jamfile
+++ b/PwaDynamics/Jamfile
@@ -47,6 +47,8 @@ exe LUTFileApp : LUTFileApp.cc PwaDynamics
unit-test test_PwaDynamics
: [ glob test/*.cc ]
$(TOP)//boost_test
+ $(TOP)/qft++//qft++
+ $(TOP)/qft++Extension//qft++Extension
$(TOP)/ErrLogger//ErrLogger
PwaDynamics
;
\ No newline at end of file
diff --git a/SetEnv_rub_AutoDiff b/SetEnv_rub_AutoDiff
index 8cc7b58f..f587b93d 100755
--- a/SetEnv_rub_AutoDiff
+++ b/SetEnv_rub_AutoDiff
@@ -1,7 +1,7 @@
source /opt/root/6-28.04-AL9.2-gcc12.2.0/bin/thisroot.csh
setenv TOP_DIR `pwd | sed -e 's/\/nfs//'`
setenv extern /data/sleipnir/PANDA/PWA
-setenv LD_LIBRARY_PATH ${TOP_DIR}/lib:${extern}/libgcc920-opt
+setenv LD_LIBRARY_PATH ${TOP_DIR}/lib:${extern}/libgcc920-opt:/data/jollyj/bertram/Adept/AdeptInst/lib
setenv BOOST_BUILD_PATH /usr/share/boost-build
setenv PATH ${TOP_DIR}/bin:${PATH}
setenv KMAT_DIR /data/duldul/bertram/KMatStore/
diff --git a/autodiffExamples/ReverseCollections.cc b/autodiffExamples/ReverseCollections.cc
index ae3eee52..f2189a52 100644
--- a/autodiffExamples/ReverseCollections.cc
+++ b/autodiffExamples/ReverseCollections.cc
@@ -32,6 +32,8 @@
#include <boost/random.hpp>
#include <chrono>
//#include <boost/timer/timer.hpp>
+//#include "adept.h"
+
boost::timer::cpu_timer theTimer;
using Tape = typename codi::RealReverse::Tape;
@@ -68,6 +70,7 @@ FloatP ReverseCollections::dynFunc(std::vector<FloatP>& currentMasses, ReverseCo
template autodiff::var ReverseCollections::dynFunc(std::vector<autodiff::var>& currentMasses, ReverseCollection::BWParams<autodiff::var>& bwParams, std::string dynType);
template codi::RealReverse ReverseCollections::dynFunc(std::vector<codi::RealReverse>& currentMasses, ReverseCollection::BWParams<codi::RealReverse>& bwParams, std::string dynType);
+//template adept::adouble ReverseCollections::dynFunc(std::vector<adept::adouble>& currentMasses, ReverseCollection::BWParams<adept::adouble>& bwParams, std::string dynType);
template<typename FloatP>
FloatP ReverseCollections::dynFuncCombine(std::vector<std::vector<FloatP>>& currentMassesCollect, std::vector<ReverseCollection::BWParams<FloatP>>& bwParamsCollect, std::vector<std::string>& dynTypeCollect){
@@ -215,31 +218,48 @@ void ReverseCollections::evalDerivBWDynamics(int noOfEvts){
ReverseCollection::bwParamsDouble bwParamsBWDouble(1.9, 0.3);
ReverseCollection::bwParamsAutoD bwParamsBWAutoD(1.9, 0.3);
ReverseCollection::bwParamsValDual bwParamsBWValDual(1.9, 0.3);
- ReverseCollection:: bwParamsRealReverse bwParamsBWValRealReverse(1.9, 0.3);
+ ReverseCollection::bwParamsRealReverse bwParamsBWValRealReverse(1.9, 0.3);
+ // ReverseCollection::bwParamsAdept bwParamsBWAdep(1.9, 0.3);
+ std::vector<double> currentMassesBWDouble;
+ generateMasses(currentMassesBWDouble, noOfEvts, bwParamsBWDouble, "BW");
std::vector<ReverseCollection::ValDual> currentMassesBWValDual;
std::vector<autodiff::var> currentMassesBWAutoVar;
std::vector<codi::RealReverse> currentMassesBWRealReverse;
- generateMasses(currentMassesBWValDual, currentMassesBWAutoVar, currentMassesBWRealReverse, noOfEvts, bwParamsBWDouble, "BW");
+ // std::vector<adept::adouble> currentMassesBWAdep;
+ fillOtherMassType(currentMassesBWDouble, currentMassesBWValDual);
+ fillOtherMassType(currentMassesBWDouble, currentMassesBWAutoVar);
+ fillOtherMassType(currentMassesBWDouble, currentMassesBWRealReverse);
+ //fillOtherMassType(currentMassesBWDouble, currentMassesBWAdep);
- std::vector<ReverseCollection::ValDual> currentMassesBWrelValDual;
- std::vector<autodiff::var> currentMassesBWrelAutoVar;
- std::vector<codi::RealReverse> currentMassesBWrelRealReverse;
ReverseCollection::bwParamsDouble bwParamsBWrelDouble(1.9, 0.3, 0.8, 0.9);
ReverseCollection::bwParamsAutoD bwParamsBWrelAutoD(1.9, 0.3, 0.8, 0.9);
ReverseCollection::bwParamsRealReverse bwParamsBWrelRealReverse(1.9, 0.3, 0.8, 0.9);
ReverseCollection::bwParamsValDual bwParamsBWrelValDual(1.9, 0.3, 0.8, 0.9);
- generateMasses(currentMassesBWrelValDual, currentMassesBWrelAutoVar, currentMassesBWrelRealReverse, noOfEvts, bwParamsBWrelDouble, "BWrel");
+ std::vector<double> currentMassesBWrelDouble;
+ generateMasses(currentMassesBWrelDouble, noOfEvts, bwParamsBWrelDouble, "BWrel");
+ std::vector<ReverseCollection::ValDual> currentMassesBWrelValDual;
+ std::vector<autodiff::var> currentMassesBWrelAutoVar;
+ std::vector<codi::RealReverse> currentMassesBWrelRealReverse;
+ fillOtherMassType(currentMassesBWrelDouble, currentMassesBWrelValDual);
+ fillOtherMassType(currentMassesBWrelDouble, currentMassesBWrelAutoVar);
+ fillOtherMassType(currentMassesBWrelDouble, currentMassesBWrelRealReverse);
- std::vector<ReverseCollection::ValDual> currentMassesBWBlattWValDual;
- std::vector<autodiff::var> currentMassesBWBlattWAutoVar;
- std::vector<codi::RealReverse> currentMassesBWBlattWRealReverse;
ReverseCollection::bwParamsDouble bwParamsBWBlattWrelDouble(1.9, 0.3, 0.8, 0.9, 2, 0.166);
ReverseCollection::bwParamsAutoD bwParamsBWBlattWrelAutoD(1.9, 0.3, 0.8, 0.9, 2, 0.166);
ReverseCollection::bwParamsRealReverse bwParamsBWBlattWrelRealReverse(1.9, 0.3, 0.8, 0.9, 2, 0.166);
ReverseCollection::bwParamsValDual bwParamsBWBlattWrelValDual(1.9, 0.3, 0.8, 0.9, 2, 0.166);
- generateMasses(currentMassesBWBlattWValDual, currentMassesBWBlattWAutoVar, currentMassesBWBlattWRealReverse, noOfEvts, bwParamsBWBlattWrelDouble, "BWBlattWrel");
-
- theTimer.start();
+ std::vector<double> currentMassesBWBlattWDouble;
+ generateMasses(currentMassesBWBlattWDouble, noOfEvts, bwParamsBWBlattWrelDouble, "BWBlattWrel");
+ std::vector<ReverseCollection::ValDual> currentMassesBWBlattWValDual;
+ std::vector<autodiff::var> currentMassesBWBlattWAutoVar;
+ std::vector<codi::RealReverse> currentMassesBWBlattWRealReverse;
+ // std::vector<adept::adouble> currentMassesBWBlattWAdep;
+ fillOtherMassType(currentMassesBWBlattWDouble, currentMassesBWBlattWValDual);
+ fillOtherMassType(currentMassesBWBlattWDouble, currentMassesBWBlattWAutoVar);
+ fillOtherMassType(currentMassesBWBlattWDouble, currentMassesBWBlattWRealReverse);
+ //fillOtherMassType(currentMassesBWBlattWDouble, currentMassesBWBlattWAdep);
+
+ theTimer.start();
autodiff::var u = dynFunc(currentMassesBWAutoVar, bwParamsBWAutoD, "BW");
auto [dudmass0, dudwidth] = derivatives(u, wrt(bwParamsBWAutoD.mass0, bwParamsBWAutoD.width));
theTimer.stop();
@@ -272,7 +292,15 @@ void ReverseCollections::evalDerivBWDynamics(int noOfEvts){
InfoMsg << "d uRealReverse / d width= " << bwParamsBWValRealReverse.width.getGradient() << endmsg;
tape.reset();
dumpTimerInfo(theTimer, "CoDi derivative BW");
-
+
+ // adept::Stack stack;
+ // stack.new_recording(); // Start recording
+ // adouble uAdep = dynFunc(currentMassesBWAdep, bwParamsBWAdep, "BW");
+ //uAdep.set_gradient(1.0); // Defines y as the objective function
+ //stack.compute_adjoint();
+ //bwParamsBWAdep.mass0.get_gradient();
+
+
theTimer.start();
double BWNumResult=dynFuncNumDerivative(currentMassesBWValDual, bwParamsBWValDual, "BW");
theTimer.stop();
@@ -308,26 +336,50 @@ void ReverseCollections::evalDerivBWDynamics(int noOfEvts){
theTimer.stop();
InfoMsg << "uRealReverseBWRel = " << uRealReverseBWRel // print the evaluated output u
- << "\ndudmass0 = " << bwParamsBWrelRealReverse.mass0.getGradient()
- << "\ndudwidth = " << bwParamsBWrelRealReverse.width.getGradient()
- << "\ndudmassA = " << bwParamsBWrelRealReverse.massA.getGradient()
- << "\ndudmassB = " << bwParamsBWrelRealReverse.massB.getGradient() << endmsg;
+ << "\ndudmass0 = " << bwParamsBWrelRealReverse.mass0.getGradient()
+ << "\ndudwidth = " << bwParamsBWrelRealReverse.width.getGradient()
+ << "\ndudmassA = " << bwParamsBWrelRealReverse.massA.getGradient()
+ << "\ndudmassB = " << bwParamsBWrelRealReverse.massB.getGradient() << endmsg;
tape.reset();
dumpTimerInfo(theTimer, "CoDi derivative BWRel");
+
+ //Adept BWRel
+ adept::Stack stack;
+ ReverseCollection::bwParamsAdept bwParamsBWrelAdep(1.9, 0.3, 0.8, 0.9);
+ std::vector<adept::adouble> currentMassesBWrelAdep;
+ fillOtherMassType(currentMassesBWrelDouble, currentMassesBWrelAdep);
+ theTimer.start();
+ stack.new_recording(); // Start recording
+ adept::adouble uAdeptBWRel = dynFunc(currentMassesBWrelAdep, bwParamsBWrelAdep, "BWrel");
+ uAdeptBWRel.set_gradient(1.0);
+ //stack.compute_adjoint();
+ stack.reverse();
+ theTimer.stop();
+
+ InfoMsg << "BWrelAdept: " << uAdeptBWRel
+ << "\ndudmass0 (num) = " << bwParamsBWrelAdep.mass0.get_gradient()
+ << "\ndudwidth (num) = " << bwParamsBWrelAdep.width.get_gradient()
+ << "\ndudmassA (num) = " << bwParamsBWrelAdep.massA.get_gradient()
+ << "\ndudmassB (num) = " << bwParamsBWrelAdep.massB.get_gradient()
+ << endmsg;
+ dumpTimerInfo(theTimer, "Adep derivative BWRel");
+
+
theTimer.start();
double BWrelNumResult=dynFuncNumDerivative(currentMassesBWrelValDual, bwParamsBWrelValDual, "BWrel");
theTimer.stop();
InfoMsg << "BWrelNumResult: " << BWrelNumResult
- << "\ndudmass0 (num) = " << bwParamsBWrelValDual.mass0.grad
- << "\ndudwidth (num) = " << bwParamsBWrelValDual.width.grad
- << "\ndudmassA (num) = " << bwParamsBWrelValDual.massA.grad
- << "\ndudmassB (num) = " << bwParamsBWrelValDual.massB.grad << endmsg;
+ << "\ndudmass0 (num) = " << bwParamsBWrelValDual.mass0.grad
+ << "\ndudwidth (num) = " << bwParamsBWrelValDual.width.grad
+ << "\ndudmassA (num) = " << bwParamsBWrelValDual.massA.grad
+ << "\ndudmassB (num) = " << bwParamsBWrelValDual.massB.grad << endmsg;
dumpTimerInfo(theTimer, "numerical derivative BWRel");
+
theTimer.start();
autodiff::var uBWBlattWRel = dynFunc( currentMassesBWBlattWAutoVar, bwParamsBWBlattWrelAutoD, "BWBlattWrel");
std::array<double, 4> gBWBlattWRel = derivatives(uBWBlattWRel, wrt(bwParamsBWBlattWrelAutoD.mass0, bwParamsBWBlattWrelAutoD.width, bwParamsBWBlattWrelAutoD.massA, bwParamsBWBlattWrelAutoD.massB));
@@ -370,10 +422,10 @@ void ReverseCollections::evalDerivBWDynamics(int noOfEvts){
theTimer.stop();
InfoMsg << "BWBlattWrelNumResult: " << BWBlattWrelNumResult
- << "\ndudmass0 (num) = " << bwParamsBWBlattWrelValDual.mass0.grad
- << "\ndudwidth (num) = " << bwParamsBWBlattWrelValDual.width.grad
- << "\ndudmassA (num) = " << bwParamsBWBlattWrelValDual.massA.grad
- << "\ndudmassB (num) = " << bwParamsBWBlattWrelValDual.massB.grad << endmsg;
+ << "\ndudmass0 (num) = " << bwParamsBWBlattWrelValDual.mass0.grad
+ << "\ndudwidth (num) = " << bwParamsBWBlattWrelValDual.width.grad
+ << "\ndudmassA (num) = " << bwParamsBWBlattWrelValDual.massA.grad
+ << "\ndudmassB (num) = " << bwParamsBWBlattWrelValDual.massB.grad << endmsg;
dumpTimerInfo(theTimer, "numerical derivative BWBlattWRel");
@@ -425,17 +477,17 @@ void ReverseCollections::evalDerivBWDynamics(int noOfEvts){
theTimer.stop();
InfoMsg << "uBWRealReverseCombine = " << uBWRealReverseCombine // print the evaluated output u
- << "\ndudmass0 BW0 = " << bwParamsBWRealReverseCollection.at(0).mass0.getGradient()
- << "\ndudwidth BW0 = " << bwParamsBWRealReverseCollection.at(0).width.getGradient()
- << "\ndudmass0 BW1 = " << bwParamsBWRealReverseCollection.at(1).mass0.getGradient()
- << "\ndudwidth BW1 = " << bwParamsBWRealReverseCollection.at(1).width.getGradient()
+ << "\ndudmass0 BW0 = " << bwParamsBWRealReverseCollection.at(0).mass0.getGradient()
+ << "\ndudwidth BW0 = " << bwParamsBWRealReverseCollection.at(0).width.getGradient()
+ << "\ndudmass0 BW1 = " << bwParamsBWRealReverseCollection.at(1).mass0.getGradient()
+ << "\ndudwidth BW1 = " << bwParamsBWRealReverseCollection.at(1).width.getGradient()
<< "\ndudmassA BW1 = " << bwParamsBWRealReverseCollection.at(1).massA.getGradient()
<< "\ndudmassB BW1 = " << bwParamsBWRealReverseCollection.at(1).massB.getGradient()
<< "\ndudmass0 BW2 = " << bwParamsBWRealReverseCollection.at(2).mass0.getGradient()
- << "\ndudwidth BW2 = " << bwParamsBWRealReverseCollection.at(2).width.getGradient()
+ << "\ndudwidth BW2 = " << bwParamsBWRealReverseCollection.at(2).width.getGradient()
<< "\ndudmassA BW2 = " << bwParamsBWRealReverseCollection.at(2).massA.getGradient()
<< "\ndudmassB BW2 = " << bwParamsBWRealReverseCollection.at(2).massB.getGradient()
- << endmsg;
+ << endmsg;
tape.reset();
@@ -474,7 +526,7 @@ void ReverseCollections::dumpTimerInfo(boost::timer::cpu_timer& theCPUTimer, std
}
}
-void ReverseCollections::generateMasses(std::vector<ReverseCollection::ValDual>& massVecDual, std::vector<autodiff::var>& massVecAutoD, std::vector<codi::RealReverse>& massRealReverse, int NoEvts, ReverseCollection::bwParamsDouble& bwParams, std::string dynType){
+void ReverseCollections::generateMasses(std::vector<double>& massVecDouble, int NoEvts, ReverseCollection::bwParamsDouble& bwParams, std::string dynType){
double maxVal=0.;
if(dynType=="BW"){
@@ -509,13 +561,20 @@ void ReverseCollections::generateMasses(std::vector<ReverseCollection::ValDual>&
currentBWVal=std::norm(BreitWignerFunction::BlattWRel(bwParams.orbMom, currentMass, bwParams.mass0, bwParams.width, bwParams.massA, bwParams.massB, bwParams.qR));
}
if (currentBWVal>currentRndNo || currentBWVal>0.01){ //accept evt
- ReverseCollection::ValDual currentMassDual(currentBWVal);
- massVecDual.push_back(currentMassDual);
- autodiff::var currentMassAutoD(currentBWVal);
- massVecAutoD.push_back(currentMassAutoD);
- codi::RealReverse currentMassRealReverse(currentBWVal);
- massRealReverse.push_back(currentMassRealReverse);
- noAcceptedEvts++;
+ massVecDouble.push_back(currentMass);
+ ++noAcceptedEvts;
}
}
}
+
+
+template<typename FloatP>
+void ReverseCollections::fillOtherMassType(std::vector<double>& inputMasses, std::vector<FloatP>& outputMasses){
+ for(unsigned int i=0; i<inputMasses.size(); ++i){
+ outputMasses.push_back((FloatP) inputMasses.at(i));
+ }
+}
+template void ReverseCollections::fillOtherMassType(std::vector<double>&, std::vector<autodiff::var>&);
+template void ReverseCollections::fillOtherMassType(std::vector<double>&, std::vector<codi::RealReverse>&);
+template void ReverseCollections::fillOtherMassType(std::vector<double>&, std::vector<ReverseCollection::ValDual>&);
+template void ReverseCollections::fillOtherMassType(std::vector<double>&, std::vector<adept::adouble>&);
diff --git a/autodiffExamples/ReverseCollections.hh b/autodiffExamples/ReverseCollections.hh
index 2be0a276..c0f506a3 100644
--- a/autodiffExamples/ReverseCollections.hh
+++ b/autodiffExamples/ReverseCollections.hh
@@ -32,6 +32,9 @@
#include <boost/timer/timer.hpp>
#include <codi.hpp>
+//#include <adept_source.h>
+#include "adept.h"
+
#include "ErrLogger/ErrLogger.hh"
namespace ReverseCollection{
@@ -101,6 +104,7 @@ namespace ReverseCollection{
typedef BWParams<autodiff::var> bwParamsAutoD;
typedef BWParams<ValDual> bwParamsValDual;
typedef BWParams<codi::RealReverse> bwParamsRealReverse;
+ typedef BWParams<adept::adouble> bwParamsAdept;
}
@@ -125,7 +129,10 @@ class ReverseCollections {
void dumpTimerInfo(boost::timer::cpu_timer& theCPUTimer, std::string logInfo);
- void generateMasses(std::vector<ReverseCollection::ValDual>& massVecDual, std::vector<autodiff::var>& massVecAutoD, std::vector<codi::RealReverse>& massRealReverse, int NoEvts, ReverseCollection::bwParamsDouble& bwParams, std::string dynType);
+ void generateMasses(std::vector<double>& massVecDouble, int NoEvts, ReverseCollection::bwParamsDouble& bwParams, std::string dynType);
+
+ template<typename FloatP>
+ void fillOtherMassType(std::vector<double>& inputMasses, std::vector<FloatP>& outputMasses);
- const double _epsilon;
+ const double _epsilon;
};
diff --git a/qft++Extension/PawianUtils.cc b/qft++Extension/PawianUtils.cc
index 48b3a684..e6e7c255 100644
--- a/qft++Extension/PawianUtils.cc
+++ b/qft++Extension/PawianUtils.cc
@@ -29,6 +29,7 @@
//#include <autodiff/forward/dual.hpp>
#include <autodiff/reverse/var.hpp>
#include <codi.hpp>
+#include "adept.h"
vector<LS> PawianQFT::GetValidLSWeak(const Spin &__j, const Spin &__s1, const Spin &__s2){
vector<LS> valid_ls;
@@ -341,6 +342,7 @@ template std::complex<double> PawianQFT::phaseSpaceFacDefault(double, double, do
//template std::complex<autodiff::dual> PawianQFT::phaseSpaceFacDefault(autodiff::dual, autodiff::dual, autodiff::dual);
template std::complex<autodiff::var> PawianQFT::phaseSpaceFacDefault(autodiff::var, autodiff::var, autodiff::var);
template std::complex<codi::RealReverse> PawianQFT::phaseSpaceFacDefault(codi::RealReverse, codi::RealReverse, codi::RealReverse);
+template std::complex<adept::adouble> PawianQFT::phaseSpaceFacDefault(adept::adouble, adept::adouble, adept::adouble);
template<typename FloatP>
std::complex<FloatP> PawianQFT::phaseSpaceFacDefault(std::complex<FloatP> mass, FloatP massDec1, FloatP massDec2){
@@ -366,6 +368,7 @@ template std::complex<double> PawianQFT::phaseSpaceFacDefault(std::complex<doubl
//template std::complex<autodiff::dual> PawianQFT::phaseSpaceFacDefault(std::complex<autodiff::dual>, autodiff::dual, autodiff::dual);
template std::complex<autodiff::var> PawianQFT::phaseSpaceFacDefault(std::complex<autodiff::var>, autodiff::var, autodiff::var);
template std::complex<codi::RealReverse> PawianQFT::phaseSpaceFacDefault(std::complex<codi::RealReverse>, codi::RealReverse, codi::RealReverse);
+template std::complex<adept::adouble> PawianQFT::phaseSpaceFacDefault(std::complex<adept::adouble>, adept::adouble, adept::adouble);
complex<double> PawianQFT::phaseSpaceFacAS(double mass, double massDec1, double massDec2){
@@ -423,6 +426,7 @@ complex<FloatP> PawianQFT::breakupMomQDefault(FloatP mass, FloatP massDec1, Floa
template complex<double> PawianQFT::breakupMomQDefault(double, double, double);
template complex<autodiff::var> PawianQFT::breakupMomQDefault(autodiff::var, autodiff::var, autodiff::var);
template complex<codi::RealReverse> PawianQFT::breakupMomQDefault(codi::RealReverse, codi::RealReverse, codi::RealReverse);
+template complex<adept::adouble> PawianQFT::breakupMomQDefault(adept::adouble, adept::adouble, adept::adouble);
template<typename FloatP>
complex<FloatP> PawianQFT::breakupMomQDefault(complex<FloatP> mass, FloatP massDec1, FloatP massDec2){
@@ -432,6 +436,7 @@ complex<FloatP> PawianQFT::breakupMomQDefault(complex<FloatP> mass, FloatP massD
template complex<double> PawianQFT::breakupMomQDefault(complex<double>, double, double);
template complex<autodiff::var> PawianQFT::breakupMomQDefault(complex<autodiff::var>, autodiff::var, autodiff::var);
template complex<codi::RealReverse> PawianQFT::breakupMomQDefault(complex<codi::RealReverse>, codi::RealReverse, codi::RealReverse);
+template complex<adept::adouble> PawianQFT::breakupMomQDefault(complex<adept::adouble>, adept::adouble, adept::adouble);
template<typename MassType>
--
GitLab