log scale basis in plots

.github/workflows/build-macos.yml

@@ -1,6 +1,6 @@
 name: Build macOS
 
-on: [push]
+on: [push, workflow_dispatch]
 
 jobs:
   build:

examples/exRand.cpp

@@ -7,7 +7,7 @@
 using namespace std;
 using namespace Latan;
 
-constexpr Index  size          = 8;
+constexpr Index  n             = 8;
 constexpr Index  nDraw         = 20000;
 constexpr Index  nSample       = 2000;
 const     string stateFileName = "exRand.seed";
@@ -40,14 +40,14 @@ int main(void)
     p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);
     p.display();
 
-    DMat       var(size, size);
-    DVec       mean(size);
-    DMatSample sample(nSample, size, 1);
+    DMat       var(n, n);
+    DVec       mean(n);
+    DMatSample sample(nSample, n, 1);
 
     cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
-    var   = DMat::Random(size, size);
+    var   = DMat::Random(n, n);
     var  *= var.adjoint();
-    mean  = DVec::Random(size); 
+    mean  = DVec::Random(n); 
     RandomNormal mgauss(mean, var, rd());
     sample[central] = mgauss();
     FOR_STAT_ARRAY(sample, s)

lib/Core/Math.cpp

@@ -18,6 +18,7 @@
  */
 
 #include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Numerical/GslFFT.hpp>
 #include <LatAnalyze/includes.hpp>
 #include <gsl/gsl_cdf.h>
 
@@ -48,16 +49,42 @@ DMat MATH_NAMESPACE::corrToVar(const DMat &corr, const DVec &varDiag)
     return res;
 }
 
-double MATH_NAMESPACE::svdDynamicRange(const DMat &mat)
+double MATH_NAMESPACE::conditionNumber(const DMat &mat)
 {
     DVec s = mat.singularValues();
 
     return s.maxCoeff()/s.minCoeff();
 }
 
-double MATH_NAMESPACE::svdDynamicRangeDb(const DMat &mat)
+double MATH_NAMESPACE::cdr(const DMat &mat)
 {
-    return 10.*log10(svdDynamicRange(mat));
+    return 10.*log10(conditionNumber(mat));
+}
+
+template <typename FFT>
+double nsdr(const DMat &m)
+{
+    Index  n = m.rows();
+    FFT    fft(n);
+    CMat   buf(n, 1);
+
+    FOR_VEC(buf, i)
+    {
+        buf(i) = 0.;
+        for (Index j = 0; j < n; ++j)
+        {
+            buf(i) += m(j, (i+j) % n);
+        }
+        buf(i) /= n;
+    }
+    fft(buf, FFT::Forward);
+
+    return 10.*log10(buf.real().maxCoeff()/buf.real().minCoeff());
+}
+
+double MATH_NAMESPACE::nsdr(const DMat &mat)
+{
+    return ::nsdr<GslFFT>(mat);
 }
 
 /******************************************************************************

lib/Core/Math.hpp

@@ -73,8 +73,9 @@ namespace MATH_NAMESPACE
     DMat corrToVar(const DMat &corr, const DVec &varDiag);
 
     // matrix SVD dynamic range
-    double svdDynamicRange(const DMat &mat);
-    double svdDynamicRangeDb(const DMat &mat);
+    double conditionNumber(const DMat &mat);
+    double cdr(const DMat &mat);
+    double nsdr(const DMat &mat);
     
     // Constants
     constexpr double pi  = 3.1415926535897932384626433832795028841970;

lib/Core/Plot.cpp

@@ -515,14 +515,16 @@ void Dash::operator()(PlotOptions &option) const
 }
 
 // LogScale constructor ////////////////////////////////////////////////////////
-LogScale::LogScale(const Axis axis)
+LogScale::LogScale(const Axis axis, const double basis)
 : axis_(axis)
+, basis_(basis)
 {}
 
 // Logscale modifier ///////////////////////////////////////////////////////////
 void LogScale::operator()(PlotOptions &option) const
 {
-    option.scaleMode[static_cast<int>(axis_)] |= Plot::Scale::log;
+    option.scaleMode[static_cast<int>(axis_)]     |= Plot::Scale::log;
+    option.logScaleBasis[static_cast<int>(axis_)]  = basis_;
 }
 
 // PlotRange constructors //////////////////////////////////////////////////////
@@ -915,11 +917,11 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
     out << "unset log" << endl;
     if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
-        out << "set log x" << endl;
+        out << "set log x " << plot.options_.logScaleBasis[x] << endl;;
     }
     if (plot.options_.scaleMode[y] & Plot::Scale::log)
     {
-        out << "set log y" << endl;
+        out << "set log y " << plot.options_.logScaleBasis[y] << endl;
     }
     if (!plot.options_.label[x].empty())
     {

lib/Core/Plot.hpp

@@ -227,6 +227,7 @@ struct PlotOptions
     std::string              caption;
     std::string              title;
     unsigned int             scaleMode[2];
+    double                   logScaleBasis[2];
     Range                    scale[2];
     std::string              label[2];
     std::string              lineColor;
@@ -314,13 +315,14 @@ class LogScale: public PlotModifier
 {
 public:
     // constructor
-    explicit LogScale(const Axis axis);
+    explicit LogScale(const Axis axis, const double basis = 10);
     // destructor
     virtual ~LogScale(void) = default;
     // modifier
     virtual void operator()(PlotOptions &option) const;
 private:
     const Axis axis_;
+    const double basis_;
 };
 
 class PlotRange: public PlotModifier

lib/Numerical/DWT.cpp

@@ -135,3 +135,26 @@ DVec DWT::backward(const std::vector<DWTLevel>& dwt) const
 
     return res;
 }
+
+// concatenate levels //////////////////////////////////////////////////////////
+DVec DWT::concat(const std::vector<DWTLevel> &dwt, const int maxLevel, const bool dropLow)
+{
+    unsigned int level = ((maxLevel >= 0) ? (maxLevel + 1) : dwt.size());
+    Index nlast = dwt[level - 1].first.size();
+    Index n = 2*dwt.front().first.size() - ((dropLow) ? nlast : 0);
+    Index pt = n, nl;
+    DVec  res(n);
+
+    for (unsigned int l = 0; l < level; ++l)
+    {
+        nl  = dwt[l].second.size();
+        pt -= nl;
+        res.segment(pt, nl) = dwt[l].second;
+    }
+    if (!dropLow)
+    {
+        res.segment(0, nl) = dwt[level-1].first;
+    }
+
+    return res;
+}

lib/Numerical/DWT.hpp

@@ -46,6 +46,8 @@ public:
     // DWT
     std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
     DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+    // concatenate levels
+    static DVec concat(const std::vector<DWTLevel>& dwt, const int maxLevel = -1, const bool dropLow = false);
 private:
     DWTFilter filter_;
 };

lib/Statistics/XYSampleData.cpp

@@ -343,7 +343,7 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
     result.nPar_       = sampleResult.getNPar();
     result.nDof_       = sampleResult.nDof_;
     result.parName_    = sampleResult.parName_;
-    result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
+    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
     
     return result;
 }

lib/Statistics/XYStatData.cpp

@@ -358,7 +358,7 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
         result    = (*m)(chi2);
         totalInit = result;
     }
-    result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
+    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
     result.chi2_       = chi2(result);
     result.nPar_       = nPar;
     result.nDof_       = layout.totalYSize - nPar;

physics/2pt-fit.cpp

@@ -24,7 +24,7 @@ int main(int argc, char *argv[])
 {
     // parse arguments /////////////////////////////////////////////////////////
     OptParser            opt;
-    bool                 parsed, doLaplace, doPlot, doHeatmap, doCorr, fold, doScan;
+    bool                 parsed, doLaplace, doPlot, doHeatmap, doCorr, fold, doScan, noGuess;
     string               corrFileName, model, outFileName, outFmt, savePlot;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
@@ -59,6 +59,8 @@ int main(int argc, char *argv[])
                   "show the fit plot");
     opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
                   "show the fit correlation heatmap");
+    opt.addOption("", "no-guess"  , OptParser::OptType::trigger, true,
+                  "do not try to guess fit parameters");
     opt.addOption("", "save-plot", OptParser::OptType::value, true,
                     "saves the source and .pdf", "");
     opt.addOption("", "scan", OptParser::OptType::trigger, true,
@@ -87,6 +89,7 @@ int main(int argc, char *argv[])
     fold         = opt.gotOption("fold");
     doPlot       = opt.gotOption("p");
     doHeatmap    = opt.gotOption("h");
+    noGuess      = opt.gotOption("no-guess");
     savePlot     = opt.optionValue("save-plot");
     doScan       = opt.gotOption("scan");
     switch (opt.optionValue<unsigned int>("v"))
@@ -167,13 +170,14 @@ int main(int argc, char *argv[])
     fitter.setThinning(thinning);
 
     // set initial values ******************************************************
-    if (modelPar.type != CorrelatorType::undefined)
+    if ((modelPar.type != CorrelatorType::undefined) and !noGuess)
     {
         init = CorrelatorModels::parameterGuess(corr, modelPar);
     }
     else
     {
-        init.fill(0.1);
+        init.fill(1.);
+        init(0) = 0.2;
     }
 
     // set limits for minimisers ***********************************************

utils/sample-dwt.cpp

@@ -17,6 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/Core/OptParser.hpp>
 #include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/Io/Io.hpp>
@@ -53,6 +54,12 @@ int main(int argc, char *argv[])
         cerr << "usage: " << argv[0];
         cerr << " <options> <input file>" << endl;
         cerr << endl << "Possible options:" << endl << opt << endl;
+        cerr << "Available DWT filters:" << endl;
+        for (auto &fv: DWTFilters::fromName)
+        {
+            cerr << fv.first << " ";
+        }
+        cerr << endl << endl;
 
         return EXIT_FAILURE;
     }
@@ -68,22 +75,45 @@ int main(int argc, char *argv[])
     DMatSample            in = Io::load<DMatSample>(inFilename), res;
     Index                 nSample = in.size(), n = in[central].rows();
     vector<DMatSample>    out(ss ? 1 : level, DMatSample(nSample)), 
-                          outh(ss ? 0 : level, DMatSample(nSample));
+                          outh(ss ? 0 : level, DMatSample(nSample)),
+                          concath(ss ? 0 : level, DMatSample(nSample));
+    DMatSample            concat(nSample, n, 1);
     DWT                   dwt(*DWTFilters::fromName.at(filterName));
     vector<DWT::DWTLevel> dataDWT(level);
 
+    FOR_STAT_ARRAY(in, s)
+    {
+        in[s].conservativeResize(n, 1);
+    }
     if (!ss)
     {
+        DMatSample buf(nSample);
+
         cout << "-- compute discrete wavelet transform" << endl;
         cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
         FOR_STAT_ARRAY(in, s)
         {
-            dataDWT = dwt.forward(in[s].col(0), level);
+            dataDWT = dwt.forward(in[s], level);
             for (unsigned int l = 0; l < level; ++l)
             {
-                out[l][s]  = dataDWT[l].first;
-                outh[l][s] = dataDWT[l].second;
+                out[l][s]     = dataDWT[l].first;
+                outh[l][s]    = dataDWT[l].second;
+                concath[l][s] = DWT::concat(dataDWT, l, true);
             }
+            concat[s] = DWT::concat(dataDWT);
+        }
+        cout << "Data CDR " << Math::cdr(in.correlationMatrix()) << " dB" << endl;
+        cout << "DWT CDR " << Math::cdr(concat.correlationMatrix()) << " dB" << endl;
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT level " << l << " CDR: L= ";
+            cout << Math::cdr(out[l].correlationMatrix()) << " dB / H= ";
+            cout << Math::cdr(outh[l].correlationMatrix()) << " dB" << endl;
+        }
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT detail level " << l << " CDR: ";
+            cout << Math::cdr(concath[l].correlationMatrix()) << " dB" << endl;
         }
     }
     else
@@ -102,7 +132,7 @@ int main(int argc, char *argv[])
         }
         FOR_STAT_ARRAY(in, s)
         {
-            dataDWT.back().first = in[s].col(0);
+            dataDWT.back().first = in[s];
             out[0][s] = dwt.backward(dataDWT);
         }
     }
@@ -115,7 +145,9 @@ int main(int argc, char *argv[])
             {
                 Io::save<DMatSample>(out[l], outFilename + "/L" + strFrom(l) + ".h5");
                 Io::save<DMatSample>(outh[l], outFilename + "/H" + strFrom(l) + ".h5");
+                Io::save<DMatSample>(concath[l], outFilename + "/concatH" + strFrom(l) + ".h5");
             }
+            Io::save<DMatSample>(concat, outFilename + "/concat.h5");
         }
         else
         {

utils/sample-plot-corr.cpp

@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
     var    = sample.varianceMatrix();
     corr   = sample.correlationMatrix();
 
-    cout << "dynamic range " << Math::svdDynamicRangeDb(corr) << " dB" << endl;
+    cout << "dynamic range " << Math::cdr(corr) << " dB" << endl;
     p << PlotCorrMatrix(corr);
     p.display();
     if (!outVarName.empty())