Merge pull request #210 from aglowacki/master

Bug fixes

Author: aglowacki

.github/workflows/cmake.yml

@@ -32,7 +32,7 @@ jobs:
         key: ${{ runner.os }}-vcpkg
     - name: Run vcpkg
       if: ${{ steps.cache-vcpkg-restore.outputs.cache-hit != 'true' }}
-      run: ./vcpkg/vcpkg install hdf5 netcdf-c yaml-cpp zeromq jsoncpp
+      run: ./vcpkg/vcpkg install hdf5 netcdf-c yaml-cpp zeromq cppzmq jsoncpp nlopt
     - name: Save vcpkg cache
       id: cache-vcpkg-save
       uses: actions/cache/save@v3
@@ -87,7 +87,7 @@ jobs:
         rsync -RavzP bin install
         cd install
         echo '{"branchName":"${{ github.ref }}","commitHash":"${{ github.sha }}"}' >> version.json
-    - uses: actions/upload-artifact@v3
+    - uses: actions/upload-artifact@v4
       with:
         name: XRF-Maps_${{ matrix.os }}_x86_64
         path: ${{github.workspace}}/install/ 

azure-pipelines.yml

@@ -24,7 +24,7 @@ jobs:
         path: $(VCPKG_CONFIG_CACHE)
         cacheHitVar: CACHE_RESTORED
         displayName: Cache vcpkg packages
-    - script: .\vcpkg\vcpkg.exe install hdf5 netcdf-c yaml-cpp zeromq jsoncpp --triplet x64-windows
+    - script: .\vcpkg\vcpkg.exe install hdf5 netcdf-c yaml-cpp zeromq cppzmq jsoncpp nlopt --triplet x64-windows
       condition: ne(variables.CACHE_RESTORED, 'true')
     - script: mkdir build
     - task: CMake@1

src/core/process_whole.cpp

@@ -473,7 +473,7 @@ bool perform_quantification(data_struct::Analysis_Job<double>* analysis_job, boo
                     logI << Fitting_Routine_To_Str.at(fit_itr.first) << " " << quant_itr.first << "\n";
                     Fit_Parameters<double> fit_params;
                     // min, and max values doen't matter because we are free fitting amplitude only
-                    fit_params.add_parameter(Fit_Param<double>("quantifier", 0.0, std::numeric_limits<double>::max()/2.0, 1.0, 0.0001, E_Bound_Type::FIT));
+                    fit_params.add_parameter(Fit_Param<double>("quantifier", 0.0000000001, std::numeric_limits<double>::max()*.00001, 1.0, 0.0001, E_Bound_Type::LIMITED_LO_HI));
                     //initial guess: parinfo_value[0] = 100000.0 / factor
                     fit_params["quantifier"].value = (double)100000.0 / quant_itr.second;
                     optimizer->minimize_quantification(&fit_params, &detector->all_element_quants[fit_itr.first][quant_itr.first], &quantification_model);

src/fitting/optimizers/nlopt_optimizer.cpp

@@ -86,10 +86,14 @@ double residuals_nlopt(const std::vector<double> &x, std::vector<double> &grad,
 
     double sum = 0.0;
     double dy = 0.;
+    double dt = 0.;
     //Calculate residuals
     for (int i=0; i<ud->spectra.size(); i++)
     {
-        sum += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0) * ud->weights[i];
+        // orig
+        //sum += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0) * ud->weights[i];
+        dy += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0);
+        dt += pow((ud->spectra[i] - ud->spec_avg), 2.0);
         /*
     	if (std::isfinite(dy) == false)
 		{
@@ -111,6 +115,7 @@ double residuals_nlopt(const std::vector<double> &x, std::vector<double> &grad,
 		}
         */
     }
+    sum = dy / dt;
     //logI << "f = " << sum << "\n";
     ud->cur_itr++;
     if (ud->status_callback != nullptr)
@@ -145,11 +150,16 @@ double gen_residuals_nlopt(const std::vector<double> &x, std::vector<double> &gr
     ud->spectra_model += ud->spectra_background;
 
     double sum = 0.0;
+    double dy = 0.;
+    double dt = 0.;
     // Calculate residuals
     for (int i=0; i<ud->spectra.size(); i++)
     {
-        sum += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0) * ud->weights[i];
+        //sum += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0) * ud->weights[i];
+        dy += pow((ud->spectra[i] - ud->spectra_model[i]), 2.0);
+        dt += pow((ud->spectra[i] - ud->spec_avg), 2.0);
     }
+    sum = dy / dt;
 
     return sum;
 }
@@ -513,13 +523,17 @@ OPTIMIZER_OUTCOME NLOPT_Optimizer<T_real>::minimize_quantification(Fit_Parameter
         return OPTIMIZER_OUTCOME::STOPPED;
     }
 
-    nlopt::opt opt(_algo, fitp_arr.size());
+    //nlopt::opt opt(nlopt::algorithm::GN_ESCH, fitp_arr.size());
+    //nlopt::opt opt(nlopt::algorithm::GN_ISRES, fitp_arr.size());
+    //nlopt::opt opt(nlopt::algorithm::LN_SBPLX, fitp_arr.size());
+    nlopt::opt opt(nlopt::algorithm::GN_CRS2_LM, fitp_arr.size());
     opt.set_lower_bounds(lb_arr);
     opt.set_upper_bounds(ub_arr);
     opt.set_default_initial_step(step_arr);
     opt.set_min_objective(quantification_residuals_nlopt<T_real>, (void*)&ud);
     opt.set_xtol_rel(_options.at(STR_OPT_XTOL));
-    opt.set_maxeval(_options.at(STR_OPT_MAXITER));
+    //opt.set_maxeval(_options.at(STR_OPT_MAXITER));
+    opt.set_maxeval(1000000);
 
     double minf;
     nlopt::result result;

src/fitting/optimizers/optimizer.h

@@ -114,8 +114,8 @@ struct User_Data
     Range energy_range;
     Spectra<T_real>  spectra_model;
     const Spectra<T_real> *orig_spectra;
-    ArrayTr<T_real> norm_arr;
-    T_real normalizer;
+    T_real spec_avg;
+    //T_real normalizer;
     Callback_Func_Status_Def* status_callback;
     size_t cur_itr;
     size_t total_itr;
@@ -137,7 +137,8 @@ struct Gen_User_Data
     Range energy_range;
 	Gen_Func_Def<T_real> func;
 	Spectra<T_real>  spectra_model;
-    T_real normalizer;
+    T_real spec_avg;
+    //T_real normalizer;
 };
 
 template<typename T_real>
@@ -171,8 +172,7 @@ void fill_user_data(User_Data<T_real> &ud,
     ud.fit_model = (Base_Model<T_real>*)model;
     // set spectra to fit
     ud.spectra = spectra->sub_spectra(energy_range.min, energy_range.count());
-    ud.norm_arr = ud.spectra.pow(2.0); // square the spectra and sum it
-    ud.normalizer = ud.norm_arr.sum();
+    ud.spec_avg = ud.spectra.sum() / energy_range.count();
     //not allocating memory. see https://eigen.tuxfamily.org/dox/group__TutorialMapClass.html
     //new (&(ud.spectra)) Eigen::Map<const ArrayTr<T_real>>(spectra->data() + energy_range.min, energy_range.count());
     ud.orig_spectra = spectra;
@@ -243,8 +243,9 @@ void fill_gen_user_data(Gen_User_Data<T_real>& ud,
     ud.fit_model = (Base_Model<T_real>*)model;
     // set spectra to fit
     ud.spectra = spectra->sub_spectra(energy_range.min, energy_range.count());
-    ArrayTr<T_real> norm_arr = ud.spectra.pow(2.0); // square the spectra and sum it
-    ud.normalizer = norm_arr.sum();
+    ud.spec_avg = ud.spectra.sum() / energy_range.count();
+    //ArrayTr<T_real> norm_arr = ud.spectra.pow(2.0); // square the spectra and sum it
+    //ud.normalizer = norm_arr.sum();
     //not allocating memory. see https://eigen.tuxfamily.org/dox/group__TutorialMapClass.html
     //new (&ud.spectra) Eigen::Map<const ArrayTr<T_real>>(spectra->data() + energy_range.min, energy_range.count());
     ud.fit_parameters = fit_params;

src/io/file/hdf5_io.cpp

@@ -57,6 +57,7 @@ POSSIBILITY OF SUCH DAMAGE.
 #include <mutex>
 #include <algorithm>
 #include <cctype>
+#include <filesystem>
 
 
 const std::vector<std::string> xrf_analysis_save_names = { STR_FIT_ROI,
@@ -381,7 +382,23 @@ bool HDF5_IO::start_save_seq(const std::string filename, bool force_new_file, bo
     }
 
     if(false == force_new_file)
+    {
         _cur_file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+    }
+    else
+    {
+        try 
+        {
+            if (std::filesystem::remove(filename))
+            {
+               logI << "file " << filename << " deleted.\n";
+            }
+        }
+        catch(const std::filesystem::filesystem_error& err) 
+        {
+            logE << "filesystem error: " << err.what() << '\n';
+        }
+    }
     if (_cur_file_id < 1)
     {
         if (open_file_only == false)