data plot with filters

CMakeLists.txt

@@ -13,7 +13,7 @@ include(FindPackageMessage)
 include(GNUInstallDirs)
 
 # C++ compile flags
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_CXX_EXTENSIONS OFF)
 check_cxx_compiler_flag("-march=native" COMPILER_SUPPORTS_MARCH_NATIVE)
@@ -29,7 +29,7 @@ set(CMAKE_CXX_FLAGS_RELWITHDEBINFO
     "-O3 -g -DNDEBUG ${MARCH_FLAG} ${MTUNE_FLAG}")
 set(gcc_like_cxx "$<COMPILE_LANG_AND_ID:CXX,ARMClang,AppleClang,Clang,GNU,LCC>")
 add_library(compiler_flags INTERFACE)
-target_compile_features(compiler_flags INTERFACE cxx_std_14)
+target_compile_features(compiler_flags INTERFACE cxx_std_17)
 target_compile_options(
   compiler_flags
   INTERFACE "$<${gcc_like_cxx}:-Wall;-Wextra;-Wshadow;-Wformat=2;-Wunused>")

lib/LatAnalyze/Core/Plot.cpp

@@ -206,6 +206,42 @@ PlotData::PlotData(const XYStatData &data, const Index i, const Index j, const b
     setCommand("'" + tmpFileName + "' " + usingCmd);
 }
 
+PlotData::PlotData(const XYStatData & data, XYStatData::CoordFilter f, Index i, const Index j, const bool abs)
+{
+    string usingCmd, tmpFileName;
+    
+    if (!abs)
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
+    }
+    else
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:(abs($3)):4 w yerr" : "u 1:(abs($3)):2:4 w xyerr";
+    }
+   
+    tmpFileName = dumpToTmpFile(data.getTable(i, j, f));
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' " + usingCmd);
+}
+
+PlotData::PlotData(const XYStatData & data, XYStatData::PointFilter f, Index i, const Index j, const bool abs)
+{
+    string usingCmd, tmpFileName;
+    
+    if (!abs)
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
+    }
+    else
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:(abs($3)):4 w yerr" : "u 1:(abs($3)):2:4 w xyerr";
+    }
+   
+    tmpFileName = dumpToTmpFile(data.getTable(i, j, f));
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' " + usingCmd);
+}
+
 // PlotPoint constructor ///////////////////////////////////////////////////////
 PlotPoint::PlotPoint(const double x, const double y)
 {
@@ -775,7 +811,7 @@ void Plot::display(void)
         ostringstream scriptBuf;
 
         getProgramPath();
-        command     = gnuplotPath_ + "/" + gnuplotBin_ + " 2>/dev/null";
+        command     = gnuplotPath_ + "/" + gnuplotBin_;
         gnuplotPipe = popen(command.c_str(), "w");
         if (!gnuplotPipe)
         {

lib/LatAnalyze/Core/Plot.hpp

@@ -94,6 +94,10 @@ public:
     PlotData(const DMatSample &x, const DVec       &y, const bool abs = false);
     PlotData(const XYStatData &data, const Index i = 0, const Index j = 0, 
              const bool abs = false);
+    PlotData(const XYStatData &data, XYStatData::CoordFilter f, Index i = 0, 
+             const Index j = 0, const bool abs = false);
+    PlotData(const XYStatData &data, XYStatData::PointFilter f, Index i = 0, 
+             const Index j = 0, const bool abs = false);
     // destructor
     virtual ~PlotData(void) = default;
 };

lib/LatAnalyze/Numerical/CompositeMinimizer.cpp

@@ -0,0 +1,71 @@
+/*
+ * CompositeMinimizer.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2024 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/CompositeMinimizer.hpp>
+#include <LatAnalyze/includes.hpp>
+#include "CompositeMinimizer.hpp"
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                   CompositeMinimizer implementation                        *
+ ******************************************************************************/
+// constructor ////////////////////////////////////////////////////////////////
+CompositeMinimizer::CompositeMinimizer(const std::vector<Minimizer *> &minimizer)
+: minVec_(minimizer)
+{
+    if (minVec_.size() == 0)
+    {
+        LATAN_ERROR(Size, "minimizer vector has size 0");
+    }
+}
+
+// access /////////////////////////////////////////////////////////////////////
+bool CompositeMinimizer::supportLimits(void) const
+{
+    return false;
+}
+
+// minimization ///////////////////////////////////////////////////////////////
+const DVec &CompositeMinimizer::operator()(const DoubleFunction &f)
+{
+    DVec &x = getState();
+    DVec tmp;
+    unsigned int i = 1, nMin = minVec_.size();
+
+    // resize minimizer state to match function number of arguments
+    if (f.getNArg() != x.size())
+    {
+        resize(f.getNArg());
+    }
+
+    for (auto &m: minVec_)
+    {
+        if (getVerbosity() >= Verbosity::Normal)
+        {
+            cout << "********** Composite minimizer step " << i << "/" << nMin << endl;
+        }
+        m->setInit(x);
+        x = (*m)(f);
+        i++;
+    }
+
+    return x;
+}

lib/LatAnalyze/Numerical/CompositeMinimizer.hpp

@@ -0,0 +1,68 @@
+/*
+ * CompositeMinimizer.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2024 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_CompositeMinimizer_hpp_
+#define Latan_CompositeMinimizer_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
+#include <LatAnalyze/Numerical/Minimizer.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                       Class for chaining minimizers                        *
+ ******************************************************************************/
+class CompositeMinimizer: public Minimizer
+{
+public:
+    // constructors
+    explicit CompositeMinimizer(const std::vector<Minimizer *> &minimizer);
+    template <typename... Ts>
+    CompositeMinimizer(Minimizer &min, Ts & ... minimisers);
+    // destructor
+    virtual ~CompositeMinimizer(void) = default;
+    // access
+    virtual bool supportLimits(void) const;
+    // minimization
+    virtual const DVec & operator()(const DoubleFunction &f);
+private:
+    std::vector<Minimizer *> minVec_;
+};
+
+/******************************************************************************
+ *              CompositeMinimizer template implementation                    *
+ ******************************************************************************/
+template <typename... Ts>
+CompositeMinimizer::CompositeMinimizer(Minimizer &min, Ts & ... minimisers)
+{
+    static_assert(static_or<std::is_assignable<Minimizer &, Ts>::value...>::value,
+                  "model arguments are not compatible with Minimizer");
+
+
+    minVec_ = {&min, &minimisers...};
+    if (minVec_.size() == 0)
+    {
+        LATAN_ERROR(Size, "minimizer vector has size 0");
+    }
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_CompositeMinimizer_hpp_

lib/LatAnalyze/Numerical/MinuitMinimizer.cpp

@@ -164,8 +164,7 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
     {
         if (getVerbosity() >= Verbosity::Normal)
         {
-            cout << "========== Minuit minimization, pass #" << n + 1;
-            cout << " =========" << endl;
+            cout << "========== Minuit minimization, pass #" << n + 1 << endl;
         }
         min.Minimize();
         status = min.Status();

lib/LatAnalyze/Numerical/NloptMinimizer.cpp

@@ -94,8 +94,7 @@ const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
     {
         if (getVerbosity() >= Verbosity::Normal)
         {
-            cout << "========== NLopt minimization, pass #" << n + 1;
-            cout << " ==========" << endl;
+            cout << "========== NLopt minimization, pass #" << n + 1 << endl;
             cout << "Algorithm: " << min.get_algorithm_name() << endl;
             cout << "Max eval.= " << min.get_maxeval();
             cout << " -- Precision= " << min.get_xtol_rel() << endl;

lib/LatAnalyze/Physics/DataFilter.hpp

@@ -54,6 +54,9 @@ private:
  ******************************************************************************/
 class LaplaceDataFilter: public DataFilter
 {
+public:
+    template <typename MatType, Index o>
+    using ObjectiveFunction = std::function<double(const StatArray<MatType, o> &)>;
 public:
     // constructor
     LaplaceDataFilter(const bool downsample = false);
@@ -63,10 +66,11 @@ public:
     template <typename MatType, Index o>
     void operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in, 
                     const double lambda = 0.);
-    // correlation optimisation
+    // metric optimisation
     template <typename MatType, Index o>
-    double optimiseCdr(const StatArray<MatType, o> &data, Minimizer &min,
-                       const unsigned int nPass = 3);
+    double optimiseFunction(const StatArray<MatType, o> &data, 
+                            ObjectiveFunction<MatType, o> &fn,
+                            Minimizer &min, const unsigned int nPass = 3);
 };
 
 /******************************************************************************
@@ -98,38 +102,25 @@ void LaplaceDataFilter::operator()(StatArray<MatType, o> &out,
 
 // correlation optimisation ///////////////////////////////////////////////////
 template <typename MatType, Index o>
-double LaplaceDataFilter::optimiseCdr(const StatArray<MatType, o> &data, 
-                                      Minimizer &min, const unsigned int nPass)
+double LaplaceDataFilter::optimiseFunction(const StatArray<MatType, o> &data, 
+                                           ObjectiveFunction<MatType, o> &fn, 
+                                           Minimizer &min,
+                                           const unsigned int nPass)
 {
     StatArray<MatType, o> fdata(data.size());
     DVec init(1);
     double reg, prec;
-    DoubleFunction cdr([&data, &fdata, this](const double *x)
+    DoubleFunction fnReg([&data, &fdata, &fn, this](const double *x)
     {
         double res;
         (*this)(fdata, data, x[0]);
-        res = Math::cdr(fdata.correlationMatrix());
+        res = fn(fdata);
         
         return res;
     }, 1);
-
-    min.setLowLimit(0., -0.1);
-    min.setHighLimit(0., 100000.);
-    init(0) = 0.1;
-    min.setInit(init);
-    prec = 0.1;
+    
     min.setPrecision(prec);
-    reg = min(cdr)(0);
-    for (unsigned int pass = 0; pass < nPass; pass++)
-    {
-        min.setLowLimit(0., (1.-10.*prec)*reg);
-        min.setHighLimit(0., (1.+10.*prec)*reg);
-        init(0) = reg;
-        min.setInit(init);
-        prec *= 0.1;
-        min.setPrecision(prec);
-        reg = min(cdr)(0);
-    }
+    reg = min(fnReg)(0);
 
     return reg;
 }

lib/LatAnalyze/Statistics/XYStatData.cpp

@@ -184,7 +184,7 @@ void XYStatData::setXError(const Index i, const DVec &err)
 
 void XYStatData::setYError(const Index j, const DVec &err)
 {
-    checkXDim(j);
+    checkYDim(j);
     checkErrVec(err, yyVar_(j, j));
     yyVar_(j, j).diagonal() = err.cwiseProduct(err);
     scheduleFitVarMatInit();
@@ -251,6 +251,62 @@ DMat XYStatData::getTable(const Index i, const Index j) const
     return table;
 }
 
+DMat XYStatData::getTable(const Index i, const Index j, CoordFilter &coordFilter) const
+{
+    checkXDim(i);
+    checkYDim(j);
+ 
+    DMat  table(getYSize(j), 4);
+    Index row = 0;
+    
+    for (auto &p: yData_[j])
+    {
+        Index k = p.first;
+        auto  c = dataCoord(k);
+        
+        if (coordFilter(c))
+        {
+            Index r = c[i];
+            table(row, 0) = x(k)(i);
+            table(row, 2) = p.second;
+            table(row, 1) = xxVar_(i, i).diagonal().cwiseSqrt()(r);
+            table(row, 3) = yyVar_(j, j).diagonal().cwiseSqrt()(row);
+            row++;
+        }
+    }
+    table.conservativeResize(row, 4);
+    
+    return table;
+}
+
+DMat XYStatData::getTable(const Index i, const Index j, PointFilter &ptFilter) const
+{
+    checkXDim(i);
+    checkYDim(j);
+ 
+    DMat  table(getYSize(j), 4);
+    Index row = 0;
+    
+    for (auto &p: yData_[j])
+    {
+        Index k = p.first;
+        auto  c = dataCoord(k);
+        
+        if (ptFilter(x(k)))
+        {
+            Index r = c[i];
+            table(row, 0) = x(k)(i);
+            table(row, 2) = p.second;
+            table(row, 1) = xxVar_(i, i).diagonal().cwiseSqrt()(r);
+            table(row, 3) = yyVar_(j, j).diagonal().cwiseSqrt()(row);
+            row++;
+        }
+    }
+    table.conservativeResize(row, 4);
+    
+    return table;
+}
+
 // get total fit variance matrix ///////////////////////////////////////////////
 const DMat & XYStatData::getFitVarMat(void)
 {

lib/LatAnalyze/Statistics/XYStatData.hpp

@@ -65,6 +65,9 @@ private:
  ******************************************************************************/
 class XYStatData: public FitInterface
 {
+public:
+    typedef std::function<bool(const std::vector<Index> &)> CoordFilter;
+    typedef std::function<bool(const DVec &)> PointFilter;
 public:
     // constructor
     XYStatData(void) = default;
@@ -89,6 +92,8 @@ public:
     DVec           getXError(const Index i) const;
     DVec           getYError(const Index j) const;
     DMat           getTable(const Index i, const Index j) const;
+    DMat           getTable(const Index i, const Index j, CoordFilter &coordFilter) const;
+    DMat           getTable(const Index i, const Index j, PointFilter &ptFilter) const;
     // get total fit variance & correlation matrices and their pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);

lib/source-list.cmake

@@ -25,6 +25,7 @@ set(LATAN_HEADERS
     LatAnalyze/Io/IoObject.hpp
     LatAnalyze/Io/Xml/tinyxml2.hpp
     LatAnalyze/Io/XmlReader.hpp
+    LatAnalyze/Numerical/CompositeMinimizer.hpp
     LatAnalyze/Numerical/DWT.hpp
     LatAnalyze/Numerical/DWTFilters.hpp
     LatAnalyze/Numerical/Derivative.hpp
@@ -77,6 +78,7 @@ set(LATAN_SOURCES
     LatAnalyze/Io/Hdf5File.cpp
     LatAnalyze/Io/Io.cpp
     LatAnalyze/Io/XmlReader.cpp
+    LatAnalyze/Numerical/CompositeMinimizer.cpp
     LatAnalyze/Numerical/DWT.cpp
     LatAnalyze/Numerical/DWTFilters.cpp
     LatAnalyze/Numerical/Derivative.cpp

utils/sample-read.cpp

@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
             Io::save(s, copy, File::Mode::write, name);
         }
     }
-    catch (Exceptions::Definition)
+    catch (Exceptions::Definition &)
     {
         DSample s    = Io::load<DSample>(fileName);
         string  name = Io::getFirstName(fileName);