diff --git a/inc/LauParameter.hh b/inc/LauParameter.hh index aab9e11..debd293 100644 --- a/inc/LauParameter.hh +++ b/inc/LauParameter.hh @@ -1,603 +1,588 @@ /* Copyright 2006 University of Warwick 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. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauParameter.hh \brief File containing declaration of LauParameter class. */ #ifndef LAU_PARAMETER #define LAU_PARAMETER #include "LauAbsRValue.hh" #include "LauBlind.hh" #include "TObject.h" #include "TString.h" #include #include #include #include /*! \class LauParameter \brief Class for defining the fit parameter objects. Holds all relevant information for the parameters for both generation and fitting step: current, initial and generated value, maximum and minimum range, error, asymmetric error, fix and float and etc. */ class LauParameter final : public TObject, public LauAbsRValue { public: //! Default constructor LauParameter() = default; //! Constructor for named parameter /*! \param [in] parName the parameter name */ explicit LauParameter(const TString& parName); - //! Constructor for parameter value - /*! - \param [in] parValue the parameter value - */ - explicit LauParameter(const Double_t parValue); - - //! Constructor double limit parameter - /*! - \param [in] parValue the parameter value - \param [in] min the minimum value of the parameter - \param [in] max the maximum value of the parameter - */ - LauParameter(const Double_t parValue, const Double_t min, const Double_t max); - - //! Constructor double limit fixed parameter - /*! - \param [in] parValue the parameter value - \param [in] min the minimum value of the parameter - \param [in] max the maximum value of the parameter - \param [in] parFixed boolean flag to fix or float parameter - */ - LauParameter(const Double_t parValue, const Double_t min, const Double_t max, const Bool_t parFixed); - - //! Constructor for double error and limit parameter - /*! - \param [in] parValue the parameter value - \param [in] parError the parameter error - \param [in] min the minimum value of the parameter - \param [in] max the maximum value of the parameter - */ - LauParameter(const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max); - - //! Constructor for parameter value and name + //! Constructor for named parameter with (fixed or floated) value /*! \param [in] parName the parameter name \param [in] parValue the parameter value + \param [in] parFixed boolean flag to fix (kTRUE) or float (kFALSE) the parameter */ - LauParameter(const TString& parName, const Double_t parValue); + LauParameter(const TString& parName, const Double_t parValue, const Bool_t parFixed = kTRUE); - //! Constructor double limit parameter and name + //! Constructor for named parameter with (fixed or floated) value and error /*! \param [in] parName the parameter name \param [in] parValue the parameter value - \param [in] min the minimum value of the parameter - \param [in] max the maximum value of the parameter + \param [in] parError the parameter error + \param [in] parFixed boolean flag to fix (kTRUE) or float (kFALSE) the parameter */ - LauParameter(const TString& parName, const Double_t parValue, const Double_t min, const Double_t max); + LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Bool_t parFixed = kTRUE); - //! Constructor double limit fixed parameter and name + //! Constructor for named parameter with (fixed or floated) value and limits /*! \param [in] parName the parameter name \param [in] parValue the parameter value \param [in] min the minimum value of the parameter \param [in] max the maximum value of the parameter \param [in] parFixed boolean flag to fix (kTRUE) or float (kFALSE) the parameter */ - LauParameter(const TString& parName, const Double_t parValue, const Double_t min, const Double_t max, const Bool_t parFixed); + LauParameter(const TString& parName, const Double_t parValue, const Double_t min, const Double_t max, const Bool_t parFixed = kTRUE); - //! Constructor double error and limit parameter and name + //! Constructor for named parameter with (fixed or floated) value, error and limits /*! \param [in] parName the parameter name \param [in] parValue the parameter value \param [in] parError the parameter error \param [in] min the minimum value of the parameter \param [in] max the maximum value of the parameter + \param [in] parFixed boolean flag to fix (kTRUE) or float (kFALSE) the parameter */ - LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max); + LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max, const Bool_t parFixed = kTRUE); // Destructor virtual ~LauParameter() noexcept; //! Copy constructor /*! \param [in] rhs the parameter to be copied */ LauParameter(const LauParameter& rhs) = default; //! Copy assignment operator /*! \param [in] rhs the parameter to be copied */ LauParameter& operator=(const LauParameter& rhs) = default; //! Move constructor /*! \param [in] rhs the parameter to be moved from */ LauParameter(LauParameter&& rhs) = default; //! Move assignment operator /*! \param [in] rhs the parameter to be moved from */ LauParameter& operator=(LauParameter&& rhs) = default; // the simple accessor functions //! The parameter name /*! \return the name of the parameter */ const TString& name() const override {return name_;} //! The blinding state /*! \return the blinding state: kTRUE means that it is blinded, kFALSE that it is not blinded */ Bool_t blind() const override {return blinder_.active();} //! Access the blinder object /*! \return the blinder */ const LauBlind* blinder() const {return blinder_.active() ? &blinder_ : nullptr;} //! The value of the parameter /*! \return the value of the parameter */ Double_t value() const override {return value_;} //! The unblinded value of the parameter /*! \return the unblinded value of the parameter */ Double_t unblindValue() const override {return blinder_.active() ? blinder_.unblind(value_) : value_;} //! The error on the parameter /*! \return the error on the parameter */ Double_t error() const {return error_;} //! The lower error on the parameter /*! \return the lower error on the parameter */ Double_t negError() const {return negError_;} //! The upper error on the parameter /*! \return the upper error on the parameter */ Double_t posError() const {return posError_;} //! The value generated for the parameter /*! \return the value generated for the parameter */ Double_t genValue() const override {return genValue_;} //! The initial value of the parameter /*! \return the initial value of the parameter given to the fitter */ Double_t initValue() const override {return initValue_;} //! The minimum value allowed for the parameter /*! \return the minimum value allowed for the parameter */ Double_t minValue() const {return minValue_;} //! The maximum value allowed for the parameter /*! \return the maximum value allowed for the parameter */ Double_t maxValue() const {return maxValue_;} + //! The limits for the parameter + /*! + \return the minimum and maximum value allowed for the parameter + */ + std::pair limits() const {return std::make_pair(minValue_,maxValue_);} + //! The range allowed for the parameter /*! \return the range allowed for the parameters, defined as the difference between the max and min value */ Double_t range() const {return this->maxValue() - this->minValue();} + //! Check whether the parameter has limits + /*! + \return the boolean flag true/false whether the parameter is limited + */ + Bool_t limited() const {return minValue_ != 0.0 || maxValue_ != 0.0;} + //! Check whether the parameter is fixed or floated /*! \return the boolean flag true/false whether the parameter is fixed */ Bool_t fixed() const override {return fixed_;} //! Check whether the parameter should be floated only in the second stage of a two stage fit /*! \return the boolean flag true/false whether it floats only in the second stage */ Bool_t secondStage() const {return secondStage_;} //! Check whether a Gaussian constraints is applied /*! \return the boolean flag true/false whether a Gaussian constraint is applied */ Bool_t gaussConstraint() const override {return gaussConstraint_;} //! The penalty term from the Gaussian constraint /*! \return the penalty term from the Gaussian constraint */ Double_t constraintPenalty() const override; //! The parameter global correlation coefficient /*! \return the global correlation coefficient */ Double_t globalCorrelationCoeff() const {return gcc_;} //! The bias in the parameter /*! \return the bias in the parameter, defined as the difference between the value and the generated value */ Double_t bias() const {return bias_;} //! The pull value for the parameter /*! \return the pull value for the parameter, defined as the bias divided by the error */ Double_t pull() const {return pull_;} //! Boolean to say it is an L value /*! \return kTRUE, LauParameters are L values */ Bool_t isLValue() const override {return kTRUE;} //! Get the LauParameter itself /*! \return a vector of the LauParameter */ std::vector getPars() override; // the simple "setter" functions //! Set the parameter name /*! \param [in] newName the name of the parameter */ void name(const TString& newName) override; //! Set the value of the parameter /*! \param [in] newValue the value of the parameter */ void value(const Double_t newValue); //! Set the error on the parameter /*! \param [in] newError the error on the parameter */ void error(const Double_t newError); //! Set the lower error on the parameter /*! \param [in] newNegError the lower error on the parameter */ void negError(const Double_t newNegError); //! Set the upper error on the parameter /*! \param [in] newPosError the upper error on the parameter */ void posError(const Double_t newPosError); //! Set the error values on the parameter /*! \param [in] newError the error on the parameter \param [in] newNegError the lower error on the parameter \param [in] newPosError the upper error on the parameter */ void errors(const Double_t newError, const Double_t newNegError, const Double_t newPosError); //! Set the value and errors on the parameter /*! \param [in] newValue the value of the parameter \param [in] newError the error on the parameter \param [in] newNegError the lower error on the parameter (default set to zero) \param [in] newPosError the upper error on the parameter (default set to zero) */ void valueAndErrors(const Double_t newValue, const Double_t newError, const Double_t newNegError = 0.0, const Double_t newPosError = 0.0); //! Set the global correlation coefficient /*! \param [in] newGCCValue the value of the coefficient */ void globalCorrelationCoeff(const Double_t newGCCValue); //! Set the generated value for the parameter /*! \param [in] newGenValue the generated value for the parameter */ void genValue(const Double_t newGenValue); //! Set the inital value for the parameter /*! \param [in] newInitValue the initial value for the parameter */ void initValue(const Double_t newInitValue); //! Set the minimum value for the parameter /*! \param [in] newMinValue the minimum value for the parameter */ void minValue(const Double_t newMinValue); //! Set the maximum value for the parameter /*! \param [in] newMaxValue the maximum value for the parameter */ void maxValue(const Double_t newMaxValue); //! Set the range for the parameter /*! \param [in] newMinValue the minimum value for the parameter \param [in] newMaxValue the maximum value for the parameter */ void range(const Double_t newMinValue, const Double_t newMaxValue); //! Set the value and range for the parameter /*! \param [in] newValue the value of the parameter \param [in] newMinValue the minimum value for the parameter \param [in] newMaxValue the maximum value for the parameter */ void valueAndRange(const Double_t newValue, const Double_t newMinValue, const Double_t newMaxValue); + //! Remove limits from the parameter + void removeLimits() { this->range(0.0,0.0); } + //! Fix or float the given parameter /*! \param [in] parFixed boolean flag to fix or float the parameter */ void fixed(const Bool_t parFixed); //! Set parameter as second-stage or not of the fit /*! \param [in] secondStagePar boolean flag to check whether is a second-stage parameter */ void secondStage(const Bool_t secondStagePar); //! Add a Gaussian constraint (or modify an existing one) /*! \param [in] newGaussMean the new value of the Gaussian constraint mean \param [in] newGaussWidth the new value of the Gaussian constraint width */ void addGaussianConstraint(const Double_t newGaussMean, const Double_t newGaussWidth); //! Remove the Gaussian constraint void removeGaussianConstraint(); //! Generate per-experiment constraint mean void generateConstraintMean() override; //! Blind the parameter /*! See LauBlind documentation for details of blinding procedure \param [in] blindingString the unique blinding string used to seed the random number generator \param [in] width the width of the Gaussian from which the offset should be sampled \param [in] flipSign activate possible random sign flip (off by default) */ void blindParameter(const TString& blindingString, const Double_t width, const Bool_t flipSign = kFALSE); // functions for the cloning mechanism //! Check whether is a clone or not /*! \return true/false whether is a clone */ Bool_t clone() const {return (parent_ != nullptr);} //! Method to create a clone from the parent parameter using the copy constructor /*! \param [in] constFactor the optional constant factor by which the clone shold be multiplied \return the cloned parameter */ LauParameter* createClone(const Double_t constFactor = 1.0); //! Method to create a clone from the parent parameter using the copy constructor and setting a new name /*! \param [in] newName the new name of the cloned parameter \param [in] constFactor the optional constant factor by which the clone shold be multiplied \return the cloned parameter */ LauParameter* createClone(const TString& newName, const Double_t constFactor = 1.0); //! The parent parameter /*! \return the parent parameter */ LauParameter* parent() const {return parent_;} //! Call to update the bias and pull values void updatePull(); //! Randomise the value of the parameter (if it is floating). /*! The pre-defined parameter range is used as the randomisation range. */ void randomiseValue(); //! Randomise the value of the parameter (if it is floating). /*! Use the given range unless either of the given values are outside the range of the parameter, in which case that value will be altered to the current max or min. \param [in] minVal the minimum value for the parameter \param [in] maxVal the maximum value for the parameter */ void randomiseValue(const Double_t minVal, const Double_t maxVal); //! Write state to a JSON record /*! \param [in,out] j the JSON record to write to */ void serialiseToJson( nlohmann::json& j ) const override; protected: //! Method to check whether value provided is within the range and that the minimum and maximum limits make sense /*! \param [in] val the value of the parameter \param [in] minVal the minimum value allowed \param [in] maxVal the maximum value allowed */ void checkRange(const Double_t val, const Double_t minVal, const Double_t maxVal); //! Method to check whether value provided is whithin the range and that the minimum and maximum limits make sense void checkRange() { this->checkRange(this->value(),this->minValue(),this->maxValue()); } //! Mark this as a clone of the given parent /*! \param theparent the parent parameter */ void clone(LauParameter* theparent) { parent_ = theparent; } //! Method to remove a clone from the list of clones /*! This is used in the destructor to allow a clone to inform its parent it is no longer around \param [in] clone the clone to be removed from the list */ void removeFromCloneList(LauParameter* clone) { auto iter = clones_.find( clone ); if ( iter != clones_.end() ) { clones_.erase( iter ); } } //! Method to clear the clone parameters void wipeClones() {clones_.clear();} //! Method to update clone values /*! \param [in] justValue boolean flag to determine whether it is necessary to update all the parameter settings or only its value. */ void updateClones(const Bool_t justValue); private: //! LauFitNtuple is a friend class friend class LauFitNtuple; //! The output streaming operator function is also a friend friend std::ostream& operator << (std::ostream& stream, const LauParameter& par); //! The parameter name TString name_; //! The parameter value Double_t value_{0.0}; //! The error on the parameter Double_t error_{0.0}; //! The lower error on the parameter Double_t negError_{0.0}; //! The upper error on the parameter Double_t posError_{0.0}; //! Toy generation value Double_t genValue_{0.0}; //! Initial fit value Double_t initValue_{0.0}; //! Minimum value for the parameter Double_t minValue_{0.0}; //! Maximum value for the parameter Double_t maxValue_{0.0}; //! Fix/float option for parameter Bool_t fixed_{kTRUE}; //! Flag whether it is floated only in the second stage of the fit Bool_t secondStage_{kFALSE}; //! Choice to use Gaussian constraint Bool_t gaussConstraint_{kFALSE}; //! True mean of the Gaussian constraint Double_t constraintTrueMean_{0.0}; //! Mean value of the Gaussian constraint Double_t constraintMean_{0.0}; //! Width of the Gaussian constraint Double_t constraintWidth_{0.0}; //! Global correlation coefficient Double_t gcc_{0.0}; //! Parameter bias Double_t bias_{0.0}; //! Parameter pull Double_t pull_{0.0}; //! The parent parameter LauParameter* parent_{nullptr}; //! The clones of this parameter std::map clones_; //! The blinding engine LauBlind blinder_; ClassDefOverride(LauParameter, 5) }; //! Output stream operator std::ostream& operator << (std::ostream& stream, const LauParameter& par); //! Type to define an array of parameters typedef std::vector< std::vector > LauParArray; //! \cond DOXYGEN_IGNORE namespace nlohmann { template <> struct adl_serializer { static LauParameter from_json(const json& j); }; } //! \endcond #endif diff --git a/src/LauMinuit.cc b/src/LauMinuit.cc index b682b34..f5f67e2 100644 --- a/src/LauMinuit.cc +++ b/src/LauMinuit.cc @@ -1,301 +1,305 @@ /* Copyright 2013 University of Warwick 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. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauMinuit.cc \brief File containing implementation of LauMinuit methods. */ #include "LauMinuit.hh" #include "LauFitObject.hh" #include "LauFitter.hh" #include "LauParameter.hh" #include "LauParamFixed.hh" #include "TMatrixD.h" #include "TVirtualFitter.h" #include #include #include // It's necessary to define an external function that specifies the address of the function // that Minuit needs to minimise. Minuit doesn't know about any classes - therefore // use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this). // Here, we use TVirtualFitter* fitter instead of gMinuit, defined below. // Then, within the external function, invoke an object from this class (LauAllModel), // and use the member functions to access the parameters/variables. extern void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); ClassImp(LauMinuit) LauMinuit::LauMinuit( const UInt_t maxPar, const LauOutputLevel verbosity ) : LauAbsFitter(), maxPar_{maxPar}, outputLevel_{verbosity} { TVirtualFitter::SetDefaultFitter( "Minuit" ); minuit_ = TVirtualFitter::Fitter( nullptr, maxPar_ ); // Set the printout level std::array argL { static_cast(outputLevel_) }; minuit_->ExecuteCommand("SET PRINT", argL.data(), argL.size()); if ( outputLevel_ == LauOutputLevel::None ) { minuit_->ExecuteCommand("SET NOWARNINGS", argL.data(), 0); } } void LauMinuit::initialise( LauFitObject* fitObj, const std::vector& parameters ) { // Check whether we're going to use asymmetric errors if ( outputLevel_ > LauOutputLevel::Quiet ) { if (useAsymmFitErrors_ == kTRUE) { std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl; std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl; } else { std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl; } } // Store the parameters params_ = parameters; // Hook the external likelihood function to this LauFitter::fitter() and this class. minuit_->SetFCN( logLikeFun ); minuit_->SetObjectFit( fitObj ); // Clear any stored parameters etc... before using minuit_->Clear(); nParams_ = params_.size(); if ( outputLevel_ > LauOutputLevel::Quiet ) { std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl; std::cout << " : Total number of parameters = " << nParams_ << std::endl; } // Define the default relative error const Double_t defaultError(0.01); // Set-up the parameters for (UInt_t i = 0; i < nParams_; ++i) { TString name = params_[i]->name(); Double_t initVal = params_[i]->initValue(); Double_t initErr = params_[i]->error(); // If we do not have a supplied estimate of the error, we should make a reasonable guess if ( initErr == 0.0 ) { if ( initVal == 0.0 ) { initErr = defaultError; } else if ( TMath::Abs(initErr/initVal) < 1e-6 ) { initErr = TMath::Abs(defaultError * initVal); } } Double_t minVal = params_[i]->minValue(); Double_t maxVal = params_[i]->maxValue(); Bool_t secondStage = params_[i]->secondStage(); if (this->twoStageFit() && secondStage == kTRUE) { params_[i]->fixed(kTRUE); } Bool_t fixVar = params_[i]->fixed(); if ( outputLevel_ > LauOutputLevel::Quiet ) { - std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; + if ( params_[i]->limited() ) { + std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl; + } else { + std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << " and error " << initErr << std::endl; + } } minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal); // Fix parameter if required if (fixVar == kTRUE) { if ( outputLevel_ > LauOutputLevel::Quiet ) { std::cout << " : Fixing parameter " << i << std::endl; } minuit_->FixParameter(i); } } LauParamFixed pred; nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors // for maximum likelihood fit. Very important command, otherwise all // extracted errors will be too big, and pull distributions will be too narrow! // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L) std::array argL { 0.5 }; fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL.data(), argL.size()); //argL[0] = 0; //fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL.data(), argL.size()); } LauFitObject* LauMinuit::getFitObject() { return (minuit_!=nullptr) ? dynamic_cast( minuit_->GetObjectFit() ) : nullptr; } const LauAbsFitter::FitStatus& LauMinuit::minimise() { std::array arglist { 1000.0*nParams_, // maximum iterations 0.05 // tolerance -> min EDM = 0.001*tolerance (0.05) }; fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist.data(), arglist.size()); // Dummy variables - need to feed them to the function // used for getting NLL, EDM and error matrix status Double_t errdef; Int_t nvpar, nparx; if (fitStatus_.status != 0) { if ( outputLevel_ > LauOutputLevel::None ) { std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl; } } else { // Check that the error matrix is ok fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); if ( outputLevel_ > LauOutputLevel::Quiet ) { std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl; } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite // 3= full accurate covariance matrix // Fit result was OK. Now get the more precise errors. fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist.data(), 1); if (fitStatus_.status != 0) { if ( outputLevel_ > LauOutputLevel::None ) { std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl; } } else { // Check that the error matrix is ok fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); if ( outputLevel_ > LauOutputLevel::Quiet ) { std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl; } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite // 3= full accurate covariance matrix // Symmetric errors and eror matrix were OK. // Get asymmetric errors if asked for. if (useAsymmFitErrors_ == kTRUE) { LauFitObject* fitObj = this->getFitObject(); fitObj->withinAsymErrorCalc( kTRUE ); fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist.data(), 1); fitObj->withinAsymErrorCalc( kFALSE ); if (fitStatus_.status != 0) { std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl; } } } } // Print results fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx); if ( outputLevel_ > LauOutputLevel::None ) { std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl; } // 0= not calculated at all // 1= approximation only, not accurate // 2= full matrix, but forced positive-definite // 3= full accurate covariance matrix if ( outputLevel_ > LauOutputLevel::None ) { minuit_->PrintResults(3, fitStatus_.NLL); } // Retrieve the covariance matrix from the fitter // For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_ // but only the elements within the sub-matrix nFreeParams_ x nFreeParams_ have values, // the "trailing" elements are zero, so we trim them off. Double_t* covMatrix = minuit_->GetCovarianceMatrix(); covMatrix_.Clear(); covMatrix_.ResizeTo( nParams_, nParams_ ); covMatrix_.SetMatrixArray( covMatrix ); covMatrix_.ResizeTo( nFreeParams_, nFreeParams_ ); return fitStatus_; } void LauMinuit::fixSecondStageParameters() { for (UInt_t i{0}; i < nParams_; ++i) { if ( params_[i]->secondStage() ) { params_[i]->fixed(kTRUE); minuit_->FixParameter(i); } } LauParamFixed pred; nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); } void LauMinuit::releaseSecondStageParameters() { for (UInt_t i{0}; i < nParams_; ++i) { if ( params_[i]->secondStage() ) { params_[i]->fixed(kFALSE); minuit_->ReleaseParameter(i); } } LauParamFixed pred; nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred); } void LauMinuit::updateParameters() { for (UInt_t i{0}; i < nParams_; ++i) { // Get the value and errors from MINUIT Double_t value { minuit_->GetParameter(i) }; Double_t error{0.0}; Double_t negError{0.0}; Double_t posError{0.0}; Double_t globalcc{0.0}; minuit_->GetErrors(i, posError, negError, error, globalcc); params_[i]->valueAndErrors(value, error, negError, posError); params_[i]->globalCorrelationCoeff(globalcc); } } // Definition of the fitting function for Minuit void logLikeFun(Int_t& npar, Double_t* /*first_derivatives*/, Double_t& f, Double_t* par, Int_t /*iflag*/) { // Routine that specifies the negative log-likelihood function for the fit. // Used by the MINUIT minimising code. LauFitObject* theModel = LauFitter::fitter().getFitObject(); // Set the internal parameters for this model using parameters from Minuit (pars): theModel->setParsFromMinuit( par, npar ); // Set the value of f to be the total negative log-likelihood for the data sample. f = theModel->getTotNegLogLikelihood(); } diff --git a/src/LauParameter.cc b/src/LauParameter.cc index d360547..ee66e18 100644 --- a/src/LauParameter.cc +++ b/src/LauParameter.cc @@ -1,709 +1,665 @@ /* Copyright 2006 University of Warwick 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. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauParameter.cc \brief File containing implementation of LauParameter class. */ #include "LauParameter.hh" #include "LauJsonTools.hh" #include "LauRandom.hh" #include "TRandom.h" #include #include ClassImp(LauParameter) LauParameter::LauParameter(const TString& parName) : name_{parName} { } -LauParameter::LauParameter(const Double_t parValue) : - value_{parValue}, - genValue_{parValue}, - initValue_{parValue}, - minValue_{parValue-1e-6}, - maxValue_{parValue+1e-6} -{ -} - -LauParameter::LauParameter(const TString& parName, const Double_t parValue) : +LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Bool_t parFixed) : name_{parName}, value_{parValue}, genValue_{parValue}, initValue_{parValue}, - minValue_{parValue-1e-6}, - maxValue_{parValue+1e-6} -{ -} - -LauParameter::LauParameter(const Double_t parValue, const Double_t min, const Double_t max) : - value_{parValue}, - genValue_{parValue}, - initValue_{parValue}, - minValue_{min}, - maxValue_{max} -{ - this->checkRange(); -} - -LauParameter::LauParameter(const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max) : - value_{parValue}, - error_{parError}, - genValue_{parValue}, - initValue_{parValue}, - minValue_{min}, - maxValue_{max} -{ - this->checkRange(); -} - -LauParameter::LauParameter(const Double_t parValue, const Double_t min, const Double_t max, const Bool_t parFixed) : - value_{parValue}, - genValue_{parValue}, - initValue_{parValue}, - minValue_{min}, - maxValue_{max}, fixed_{parFixed} { - this->checkRange(); } -LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Double_t min, const Double_t max) : +LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Bool_t parFixed) : name_{parName}, value_{parValue}, + error_{parError}, genValue_{parValue}, initValue_{parValue}, - minValue_{min}, - maxValue_{max} + fixed_{parFixed} { - this->checkRange(); } LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Double_t min, const Double_t max, const Bool_t parFixed) : name_{parName}, value_{parValue}, genValue_{parValue}, initValue_{parValue}, minValue_{min}, maxValue_{max}, fixed_{parFixed} { this->checkRange(); } -LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max) : +LauParameter::LauParameter(const TString& parName, const Double_t parValue, const Double_t parError, const Double_t min, const Double_t max, const Bool_t parFixed) : name_{parName}, value_{parValue}, error_{parError}, genValue_{parValue}, initValue_{parValue}, minValue_{min}, - maxValue_{max} + maxValue_{max}, + fixed_{parFixed} { this->checkRange(); } LauParameter::~LauParameter() noexcept { // if we're a clone, we just need to inform our parent of our demise if ( this->clone() ) { parent_->removeFromCloneList(this); return; } // if we have no clones there's nothing to do if ( clones_.empty() ) { return; } // otherwise if we have clones we need to make one of them the new parent and inform the rest of the change // let's (arbitrarily) make the first parameter in the map the new parent auto iter = clones_.begin(); auto [ newParent, constFactor ] { *iter }; // remove that entry in the map clones_.erase(iter); // for all the other entries, we need to tell them they are clones of the new parent // and also rescale the constants so that they are relative to the new parent for ( auto& [ theClone, theFactor ] : clones_ ) { theClone->clone(newParent); theFactor /= constFactor; } // transfer the list of clones to the new parent newParent->clones_ = std::move(clones_); // finally, the new parent has to be told that it isn't a clone anymore newParent->clone(nullptr); } std::vector LauParameter::getPars() { std::vector list; list.push_back(this); return list; } void LauParameter::value(const Double_t newValue) { if (this->clone()) { parent_->value(newValue); } else { this->checkRange(newValue,this->minValue(),this->maxValue()); this->updateClones(kTRUE); } } void LauParameter::error(const Double_t newError) { if (this->clone()) { parent_->error(newError); } else { error_ = TMath::Abs(newError); this->updateClones(kFALSE); } } void LauParameter::negError(const Double_t newNegError) { if (this->clone()) { parent_->negError(newNegError); } else { negError_ = TMath::Abs(newNegError); this->updateClones(kFALSE); } } void LauParameter::posError(const Double_t newPosError) { if (this->clone()) { parent_->posError(newPosError); } else { posError_ = TMath::Abs(newPosError); this->updateClones(kFALSE); } } void LauParameter::errors(const Double_t newError, const Double_t newNegError, const Double_t newPosError) { if (this->clone()) { parent_->errors(newError,newNegError,newPosError); } else { error_ = TMath::Abs(newError); negError_ = TMath::Abs(newNegError); posError_ = TMath::Abs(newPosError); this->updateClones(kFALSE); } } void LauParameter::valueAndErrors(const Double_t newValue, const Double_t newError, const Double_t newNegError, const Double_t newPosError) { if (this->clone()) { parent_->valueAndErrors(newValue,newError,newNegError,newPosError); } else { this->checkRange(newValue,this->minValue(),this->maxValue()); error_ = TMath::Abs(newError); negError_ = TMath::Abs(newNegError); posError_ = TMath::Abs(newPosError); this->updateClones(kFALSE); } } void LauParameter::globalCorrelationCoeff(const Double_t newGCCValue) { if (this->clone()) { parent_->globalCorrelationCoeff(newGCCValue); } else { gcc_ = newGCCValue; this->updateClones(kFALSE); } } void LauParameter::genValue(const Double_t newGenValue) { if (this->clone()) { parent_->genValue(newGenValue); } else { genValue_ = newGenValue; this->updateClones(kFALSE); } } void LauParameter::initValue(const Double_t newInitValue) { if (this->clone()) { parent_->initValue(newInitValue); } else { initValue_ = newInitValue; this->updateClones(kFALSE); } } void LauParameter::minValue(const Double_t newMinValue) { if (this->clone()) { parent_->minValue(newMinValue); } else { this->checkRange(this->value(),newMinValue,this->maxValue()); this->updateClones(kFALSE); } } void LauParameter::maxValue(const Double_t newMaxValue) { if (this->clone()) { parent_->maxValue(newMaxValue); } else { this->checkRange(this->value(),this->minValue(),newMaxValue); this->updateClones(kFALSE); } } void LauParameter::range(const Double_t newMinValue, const Double_t newMaxValue) { if (this->clone()) { parent_->range(newMinValue,newMaxValue); } else { this->checkRange(this->value(),newMinValue,newMaxValue); this->updateClones(kFALSE); } } void LauParameter::valueAndRange(const Double_t newValue, const Double_t newMinValue, const Double_t newMaxValue) { if (this->clone()) { parent_->valueAndRange(newValue,newMinValue,newMaxValue); } else { this->checkRange(newValue,newMinValue,newMaxValue); this->updateClones(kFALSE); } } void LauParameter::name(const TString& newName) { // no need to update clones here // clones are allowed to have different names name_ = newName; } void LauParameter::fixed(const Bool_t parFixed) { if (this->clone()) { parent_->fixed(parFixed); } else { fixed_ = parFixed; this->updateClones(kFALSE); } } void LauParameter::secondStage(const Bool_t secondStagePar) { if (this->clone()) { parent_->secondStage(secondStagePar); } else { secondStage_ = secondStagePar; this->updateClones(kFALSE); } } void LauParameter::addGaussianConstraint(const Double_t newGaussMean, const Double_t newGaussWidth) { if (this->clone()) { parent_->addGaussianConstraint(newGaussMean,newGaussWidth); } else { gaussConstraint_ = kTRUE; constraintTrueMean_ = newGaussMean; constraintMean_ = newGaussMean; constraintWidth_ = newGaussWidth; this->updateClones(kFALSE); } } void LauParameter::removeGaussianConstraint() { if (this->clone()) { parent_->removeGaussianConstraint(); } else { gaussConstraint_ = kFALSE; this->updateClones(kFALSE); } } void LauParameter::generateConstraintMean() { if (this->clone()) { parent_->generateConstraintMean(); } else { constraintMean_ = LauRandom::randomFun()->Gaus( constraintTrueMean_, constraintWidth_ ); this->updateClones(kFALSE); } } Double_t LauParameter::constraintPenalty() const { const Double_t val { this->unblindValue() }; const Double_t diff { val - constraintMean_ }; const Double_t term { diff * diff }; return term / ( 2.0 * constraintWidth_ * constraintWidth_ ); } void LauParameter::blindParameter(const TString& blindingString, const Double_t width, const Bool_t flipSign) { if (this->clone()) { parent_->blindParameter(blindingString,width,flipSign); return; } if ( blinder_.active() ) { std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for this parameter" << std::endl; return; } blinder_ = LauBlind{ blindingString, width, flipSign }; for ( auto& [ clonePar, _ ] : clones_ ) { if ( clonePar->blinder_.active() ) { std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for a clone of this parameter - it will be replaced!" << std::endl; } clonePar->blinder_ = blinder_; } } void LauParameter::updatePull() { if (this->clone()) { parent_->updatePull(); return; } // calculate the bias bias_ = value_ - genValue_; // if we have errors calculated then calculate // the pull using the best error available if ((bias_ > 0.0) && (negError_ > 1e-10)) { pull_ = bias_ / negError_; } else if ((bias_ < 0.0) && (posError_ > 1e-10)) { pull_ = bias_ / posError_; } else if (error_ > 1e-10) { pull_ = bias_ / error_; } else { pull_ = 0.0; } this->updateClones(kFALSE); } void LauParameter::checkRange(const Double_t val, const Double_t minVal, const Double_t maxVal) { // first check that min is less than max (or they are the same - this is allowed) if (minVal > maxVal) { std::cerr<<"ERROR in LauParameter::checkRange : minValue: "< maxValue_) { minValue_ = maxValue_; std::cerr<<" : Setting both to "< maxVal)) { - if (name_ != "") { - std::cerr<<"ERROR in LauParameter::checkRange : value: "<limited() && ( val < minValue_ || val > maxValue_ ) ) { + std::cerr<<"ERROR in LauParameter::checkRange : value: "<clone()) { LauParameter* clonePar = parent_->createClone(constFactor); clonePar->name(this->name()); return clonePar; } // clone ourselves using the copy-constructor LauParameter* clonePar = new LauParameter(*this); Double_t newValue = clonePar->value() * constFactor; clonePar->value( newValue ); clonePar->wipeClones(); clonePar->clone(this); clones_.insert( std::make_pair( clonePar, constFactor ) ); return clonePar; } LauParameter* LauParameter::createClone(const TString& newName, const Double_t constFactor) { // self message to create the clone LauParameter* clonePar = this->createClone(constFactor); // set the new name clonePar->name(newName); // and return return clonePar; } void LauParameter::updateClones(const Bool_t justValue) { // if we don't have any clones then there's nothing to do if ( clones_.empty() ) { return; } // we have to set the values directly rather than using member functions because otherwise we'd get into an infinite loop if (justValue) { for ( auto& [ clonePar, constFactor ] : clones_ ) { clonePar->value_ = constFactor*value_; } } else { for ( auto& [ clonePar, constFactor ] : clones_ ) { clonePar->value_ = constFactor*value_; clonePar->error_ = constFactor*error_; clonePar->negError_ = constFactor*negError_; clonePar->posError_ = constFactor*posError_; clonePar->genValue_ = constFactor*genValue_; clonePar->initValue_ = constFactor*initValue_; clonePar->minValue_ = constFactor*minValue_; clonePar->maxValue_ = constFactor*maxValue_; clonePar->fixed_ = fixed_; clonePar->secondStage_ = secondStage_; clonePar->gaussConstraint_ = gaussConstraint_; clonePar->constraintTrueMean_ = constraintTrueMean_; clonePar->constraintMean_ = constraintMean_; clonePar->constraintWidth_ = constraintWidth_; clonePar->gcc_ = gcc_; clonePar->bias_ = bias_; clonePar->pull_ = pull_; } } } void LauParameter::randomiseValue() { this->randomiseValue(this->minValue(), this->maxValue()); } void LauParameter::randomiseValue(const Double_t minVal, const Double_t maxVal) { // if we're fixed then do nothing if (this->fixed()) { return; } // check supplied values are sensible if (maxVal < minVal) { std::cerr<<"ERROR in LauParameter::randomiseValue : Supplied maximum value smaller than minimum value."<Uniform( TMath::Max( minVal, this->minValue() ), TMath::Min( maxVal, this->maxValue() ) ) }; this->initValue(val); } void LauParameter::serialiseToJson( nlohmann::json& j ) const { using nlohmann::json; LauAbsRValue::serialiseToJson( j ); // First check whether this is a clone because there's much less to do if it is const Bool_t cloned { this->clone() }; j["clone"] = cloned; if ( cloned ) { j["parent"] = parent_->name(); auto me {const_cast(this)}; j["constFactor"] = parent_->clones_.at(me); return; } // TODO - do we need to write this? // - I guess it's maybe useful for a human reader but we don't use it at all in deserialisation if ( ! clones_.empty() ) { j["clones"] = json::array(); for ( auto& [ par, factor ] : clones_ ) { json obj = json::object( { {"name", par->name()}, {"factor", factor} } ); j["clones"].push_back( obj ); } } j["value"] = value_; j["genValue"] = genValue_; j["initValue"] = initValue_; j["error"] = error_; j["negError"] = negError_; j["posError"] = posError_; j["minValue"] = minValue_; j["maxValue"] = maxValue_; j["secondStage"] = secondStage_; j["gcc"] = gcc_; j["bias"] = bias_; j["pull"] = pull_; j["fixed"] = fixed_; const Bool_t blind { this->blind() }; j["blind"] = blind; if ( blind ) { j["blindingString"] = blinder_.blindingString(); j["blindingWidth"] = blinder_.blindingWidth(); j["blindingFlip"] = blinder_.flipSign(); } j["gaussConstraint"] = gaussConstraint_; if ( gaussConstraint_ ) { j["constraintMean"] = constraintTrueMean_; j["constraintWidth"] = constraintWidth_; } } // ostream operator std::ostream& operator<<(std::ostream& stream, const LauParameter& par) { stream << par.name() << " : "; stream << par.value(); if ( par.negError() > 1e-10 && par.posError() > 1e-10 ) { stream << " +/- (" << par.negError() << ", " << par.posError() << ")"; } else if ( par.error() > 1e-10 ) { stream << " +/- " << par.error(); } if ( par.fixed() ) { stream << " C"; } - stream << " L(" << par.minValue() << ", " << par.maxValue() << ")"; + if ( par.limited() ) { + stream << " L(" << par.minValue() << ", " << par.maxValue() << ")"; + } if ( par.gaussConstraint_ ) { stream << " G(" << par.constraintTrueMean_ << " +/- " << par.constraintWidth_ << ")"; } if ( par.blind() ) { stream << " BLIND(" << par.blinder()->blindingString() << ", " << par.blinder()->blindingWidth() << ", " << std::boolalpha << par.blinder()->flipSign() << ")"; } return stream; } //! \cond DOXYGEN_IGNORE LauParameter nlohmann::adl_serializer::from_json(const json& j) { using LauJsonTools::getValue; const auto isLValue { getValue( j, "isLValue" ) }; if ( ! isLValue ) { throw LauJsonTools::InvalidJson{"Cannot create LauParameter from non L-value"}; } const auto clone { getValue( j, "clone" ) }; if ( clone ) { throw LauJsonTools::InvalidJson{"Cannot create cloned LauParameter standalone"}; } // Get everything we need to construct the parameter const auto name { getValue( j, "name" ) }; const auto value { getValue( j, "value" ) }; const auto minValue { getValue( j, "minValue" ) }; const auto maxValue { getValue( j, "maxValue" ) }; const auto fixed { getValue( j, "fixed" ) }; LauParameter par { name, value, minValue, maxValue, fixed }; // Then get any additional information and use the appropriate functions to set the values const auto genValue { getValue( j, "genValue" ) }; par.genValue( genValue ); const auto initValue { getValue( j, "initValue" ) }; par.initValue( initValue ); const auto error { getValue( j, "error" ) }; const auto negError { getValue( j, "negError" ) }; const auto posError { getValue( j, "posError" ) }; par.errors( error, negError, posError ); par.updatePull(); const auto secondStage { getValue( j, "secondStage" ) }; par.secondStage( secondStage ); const auto gcc { getValue( j, "gcc" ) }; par.globalCorrelationCoeff( gcc ); const auto gaussCons { getValue( j, "gaussConstraint" ) }; if ( gaussCons ) { const auto mean { getValue( j, "constraintMean" ) }; const auto width { getValue( j, "constraintWidth" ) }; par.addGaussianConstraint( mean, width ); } const auto blind { getValue( j, "blind" ) }; if ( blind ) { const auto blindingStr { getValue( j, "blindingString" ) }; const auto blindingWidth { getValue( j, "blindingWidth" ) }; const auto blindingFlip { getValue( j, "blindingFlip" ) }; par.blindParameter( blindingStr, blindingWidth, blindingFlip ); } return par; } //! \endcond