From 471c553f959e2f687251c45d093ac108fd8dab65 Mon Sep 17 00:00:00 2001
From: pospelov <pospelov@fz-juelich.de>
Date: Tue, 20 Nov 2012 18:29:15 +0100
Subject: [PATCH] possibility to make 1D slices from 2D OutputData in
IsGISAXSTools and to make 1D plots
---
App/inc/IsGISAXSTools.h | 15 +-
App/inc/TestIsGISAXS12.h | 24 ++-
App/src/FitSuiteHelper.cpp | 77 ++++-----
App/src/IsGISAXSTools.cpp | 212 +++++++++++++++++-------
App/src/TestFittingModule2.cpp | 6 +-
App/src/TestIsGISAXS12.cpp | 161 ++++++++++++++++--
Core/Algorithms/src/GISASExperiment.cpp | 1 -
Core/Tools/inc/FitSuite.h | 3 +
Core/Tools/inc/OutputData.h | 2 +-
Core/Tools/src/FitSuite.cpp | 17 ++
Core/Tools/src/FitSuiteKit.cpp | 3 +
Core/Tools/src/FitSuiteStrategy.cpp | 2 +-
12 files changed, 394 insertions(+), 129 deletions(-)
diff --git a/App/inc/IsGISAXSTools.h b/App/inc/IsGISAXSTools.h
index 2740bc04496..5f14ecba501 100644
--- a/App/inc/IsGISAXSTools.h
+++ b/App/inc/IsGISAXSTools.h
@@ -17,6 +17,7 @@
#include <string>
+class TH1;
class TH2D;
//- -------------------------------------------------------------------
@@ -25,6 +26,12 @@ class TH2D;
//- -------------------------------------------------------------------
class IsGISAXSTools {
public:
+ struct AxisStructure {
+ int nbins;
+ std::vector<double> xbins;
+ std::string name;
+ };
+
//! draw 2D histogram representing logarithm of OutputData (in new canvas)
static void drawLogOutputData(const OutputData<double> &output, const std::string &canvas_name,
const std::string &canvas_title, const std::string &draw_options,
@@ -83,9 +90,15 @@ public:
, std::vector<std::vector<double > > &v_axis0
, std::vector<std::vector<double > > &v_axis1);
- //! create TH2D from OutputData
+ //! create two-dimensional TH2D from OutputData
static TH2D *getOutputDataTH2D(const OutputData<double>& output, const std::string &histo_name);
+ //! create one, thwo, three-dimensional histograms from OutputData
+ static TH1 *getOutputDataTH123D(const OutputData<double>& output, const std::string &histo_name);
+
+ //! slice 2D output data into 1D having fixed axes value on one of the axes
+ static OutputData<double> *sliceOutputData(const OutputData<double > &data, const std::string &fixed_axes, double fixed_axis_value);
+
private:
static double m_hist_min; // minimum value of y-axis (for 1D histograms), or z-axis (for 2D histograms)
static double m_hist_max; // maximum value of y-axis (for 1D histograms), or z-axis (for 2D histograms)
diff --git a/App/inc/TestIsGISAXS12.h b/App/inc/TestIsGISAXS12.h
index 3403c49990f..062d7b293d0 100644
--- a/App/inc/TestIsGISAXS12.h
+++ b/App/inc/TestIsGISAXS12.h
@@ -18,10 +18,11 @@
#include "OutputData.h"
#include "ISample.h"
#include "ISampleBuilder.h"
-#include "GISASExperiment.h"
-
#include <string>
+class GISASExperiment;
+class FitSuite;
+
//- -------------------------------------------------------------------
//! @class TestIsGISAXS12
//! @brief Comparison with IsGISAXS ex-12: constrained fit example
@@ -82,13 +83,32 @@ private:
double getFixedAngle() { return (fixed_phif ? isgiDataVector.at(0).phif : isgiDataVector.at(0).alphaf); }
};
+ //! initialize experiment
+ void initialiseExperiment();
+
+ //! run standard isgisaxs comparison for the sample of that kind
+ void run_isgisaxs_comparison();
+
+ //! run test fit
+ void run_test_fit();
+
+ //! run isgisaxs ex-12 style fit
+ void run_isgisaxs_fit();
+
+ //! generate test real data
+ OutputData<double > *createNoisyData(const OutputData<double> &exact_data, double noise_factor = 0.1);
+
//! read special isgisaxs *.dat file with data to fit
void read_isgisaxs_datfile(const std::string &filename);
+
std::string m_data_path;
//! represent content of isgisaxs *.dat file (in given ex-12 case, two scans - one with fixed alpha_f and another with fixed phi_f, which called in isgisaxs "2theta_f")
std::vector<IsgiScan > m_isgiCrossSections;
GISASExperiment *m_experiment;
+ ISampleBuilder *m_sample_builder;
+ FitSuite *m_fitSuite;
+
};
diff --git a/App/src/FitSuiteHelper.cpp b/App/src/FitSuiteHelper.cpp
index 80f9b998e65..18af3aa7c1a 100644
--- a/App/src/FitSuiteHelper.cpp
+++ b/App/src/FitSuiteHelper.cpp
@@ -27,30 +27,31 @@ void FitSuiteObserverPrint::update(IObservable *subject)
FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject);
if( !fitSuite ) throw NullPointerException("FitSuiteObserverPrint::update() -> Error! Can't cast FitSuite");
- // printing parameter values
- std::cout << "FitSuiteObserverPrint::update() -> Info."
- << " NumberOfVariables:" << fitSuite->getMinimizer()->getNumberOfVariables()
- << " NCall:" << fitSuite->getNCall()
- << " NStrategy:" << fitSuite->getNStrategy()
- << " Chi2:" << std::scientific << std::setprecision(8) << fitSuite->getSuiteKit()->getChiSquaredModule()->getValue() << std::endl;
- timeval call_time;
- gettimeofday(&call_time, 0);
- clock_t call_clock = clock();
- double timediff = (call_time.tv_sec - m_last_call_time.tv_sec) + double(call_time.tv_usec - m_last_call_time.tv_usec)/1e+06;
- std::cout << "Time spend since last call, cpu:" << std::fixed << std::setprecision(2) << double(call_clock-m_last_call_clock)/CLOCKS_PER_SEC << " " << "sec, wall time " << timediff << "sec" << std::endl;
- m_last_call_clock = call_clock;
- m_last_call_time = call_time;
-
- int npar(0);
- for(FitSuiteParameters::iterator it = fitSuite->getFitParameters()->begin(); it!=fitSuite->getFitParameters()->end(); ++it, ++npar) {
- std::cout << " # "<< npar << " " << (*(*it)) << std::endl;
- }
-
if(fitSuite->isLastIteration()) {
std::cout << std::endl;
std::cout << "FitSuiteObserverPrint::update() -> Info. Printing results" << std::endl;
fitSuite->getMinimizer()->printResults();
+ } else {
+ // printing parameter values
+ std::cout << "FitSuiteObserverPrint::update() -> Info."
+ << " NumberOfVariables:" << fitSuite->getMinimizer()->getNumberOfVariables()
+ << " NCall:" << fitSuite->getNCall()
+ << " NStrategy:" << fitSuite->getNStrategy()
+ << " Chi2:" << std::scientific << std::setprecision(8) << fitSuite->getSuiteKit()->getChiSquaredModule()->getValue() << std::endl;
+ timeval call_time;
+ gettimeofday(&call_time, 0);
+ clock_t call_clock = clock();
+ double timediff = (call_time.tv_sec - m_last_call_time.tv_sec) + double(call_time.tv_usec - m_last_call_time.tv_usec)/1e+06;
+ std::cout << "Time spend since last call, cpu:" << std::fixed << std::setprecision(2) << double(call_clock-m_last_call_clock)/CLOCKS_PER_SEC << " " << "sec, wall time " << timediff << "sec" << std::endl;
+ m_last_call_clock = call_clock;
+ m_last_call_time = call_time;
+
+ int npar(0);
+ for(FitSuiteParameters::iterator it = fitSuite->getFitParameters()->begin(); it!=fitSuite->getFitParameters()->end(); ++it, ++npar) {
+ std::cout << " # "<< npar << " " << (*(*it)) << std::endl;
+ }
}
+
}
@@ -62,7 +63,8 @@ void FitSuiteObserverDraw::update(IObservable *subject)
FitSuite *fitSuite = dynamic_cast<FitSuite *>(subject);
if( !fitSuite ) throw NullPointerException("FitSuiteObserverDraw::update() -> Error! Can't cast FitSuite");
- if( !fitSuite->isLastIteration() && (fitSuite->getNCall() % m_draw_every_nth != 0) ) return;
+ if( fitSuite->isLastIteration()) return;
+ if( (fitSuite->getNCall() % m_draw_every_nth != 0) && fitSuite->getNCall()!=0) return;
TCanvas *c1 = dynamic_cast<TCanvas *>( gROOT->FindObject(m_canvas_name.c_str()) );
if(!c1) {
@@ -94,6 +96,12 @@ void FitSuiteObserverDraw::update(IObservable *subject)
IsGISAXSTools::drawOutputDataInPad(*diff_map_relative, "COLZ", "relative difference map");
delete diff_map_relative;
+ std::cout << "FitSuiteObserverDraw::update() -> debug " << std::endl;
+ const OutputData<double > *real = fitSuite->getSuiteKit()->getRealData();
+ const OutputData<double > *simul = fitSuite->getSuiteKit()->getSimulatedData();
+ std::cout << " real: " << real->getAllocatedSize() << " " << real->getNdimensions() << " " << std::endl;
+ std::cout << " real: " << simul->getAllocatedSize() << " " << simul->getNdimensions() << " " << std::endl;
+
// chi2 difference
c1->cd(4);
gPad->SetLogz();
@@ -122,35 +130,6 @@ void FitSuiteObserverDraw::update(IObservable *subject)
}
m_ptext->Draw();
- if(fitSuite->isLastIteration()) {
-// TCanvas *c2 = new TCanvas("FitSuiteObserverDraw_c2", "FitSuiteObserverDraw_c2", 1024, 768);
-// c2->Divide(2,2);
-// ROOT::Math::Minimizer *minim = (dynamic_cast<ROOTMinimizer *>(fitSuite->getMinimizer()))->getROOTMinimizer();
-// if(!minim) {
-// throw NullPointerException("FitSuiteObserverDraw::update() -> Can't get right minimizer from FitSuite");
-// }
-// int npad=1;
-// for(unsigned int i=0; i<minim->NDim(); ++i) {
-// for(unsigned int j=0; j<i; ++j) {
-// if(i != j) {
-// c2->cd(npad++);
-// unsigned int np=50;
-// std::vector<double> x;
-// std::vector<double> y;
-// x.resize(np, 0);
-// y.resize(np, 0);
-// minim->Contour(i,j, np,&x[0], &y[0]);
-// x.push_back(x.front());
-// y.push_back(y.front());
-// TGraph *gr = new TGraph(np+1, &x[0], &y[0]);
-// gr->Draw("apl");
-// }
-// }
-// }
-// std::cout << "Last iteration " << std::endl;
- }
-
-
}
diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index f512bd8df67..47cefdd1da9 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -3,7 +3,9 @@
#include "Exceptions.h"
#include "TCanvas.h"
+#include "TH1.h"
#include "TH2.h"
+#include "TH3.h"
#include "TStyle.h"
#include <iostream>
@@ -48,7 +50,7 @@ void IsGISAXSTools::drawOutputData(const OutputData<double>& output,
}
/* ************************************************************************* */
-// draw 2D histogram representing OutputData (in current gPad)
+// draw 1D or 2D histograms representing OutputData (in current gPad)
/* ************************************************************************* */
void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output,
const std::string& draw_options, const std::string &histogram_title)
@@ -57,65 +59,20 @@ void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output,
throw NullPointerException("IsGISAXSTools::drawOutputDataInPad() -> Error! No canvas exists.");
}
-// if (output.getDimension() != 2) return;
-
-// // creation of 2D histogram from calculated intensities
-// output.resetIndex();
-// const NamedVector<double> *p_y_axis = reinterpret_cast<const NamedVector<double>*>(output.getAxes()[0]);
-// const NamedVector<double> *p_z_axis = reinterpret_cast<const NamedVector<double>*>(output.getAxes()[1]);
-// std::string y_axis_name = p_y_axis->getName();
-// std::string z_axis_name = p_z_axis->getName();
-// size_t y_size = p_y_axis->getSize();
-// size_t z_size = p_z_axis->getSize();
-// double y_start = (*p_y_axis)[0]/Units::degree;
-// double y_end = (*p_y_axis)[y_size-1]/Units::degree;
-// double z_start = (*p_z_axis)[0]/Units::degree;
-// double z_end = (*p_z_axis)[z_size-1]/Units::degree;
-// std::string histo_name = histogram_title;
-// if (histo_name.empty()) {
-// histo_name = gPad->GetTitle();
-// }
-// TH2D p_hist2D("p_hist2D", histo_name.c_str(), y_size, y_start, y_end, z_size, z_start, z_end);
-// //p_hist2D->UseCurrentStyle();
-// p_hist2D.GetXaxis()->SetTitle(y_axis_name.c_str());
-// p_hist2D.GetYaxis()->SetTitle(z_axis_name.c_str());
-
-// while (output.hasNext())
-// {
-// size_t index_y = output.getCurrentIndexOfAxis(y_axis_name.c_str());
-// size_t index_z = output.getCurrentIndexOfAxis(z_axis_name.c_str());
-// //std::cout << "!!! " << index_y << " " << index_z << std::endl;
-// double x_value = (*p_y_axis)[index_y]/Units::degree;
-// double y_value = (*p_z_axis)[index_z]/Units::degree;
-// double z_value = output.next();
-// p_hist2D.Fill(x_value, y_value, z_value);
-// }
-// p_hist2D.SetContour(50);
-// p_hist2D.SetStats(0);
-// p_hist2D.GetYaxis()->SetTitleOffset(1.3);
-
-// gStyle->SetPalette(1);
-// gStyle->SetOptStat(0);
-//// gPad->SetLogz();
-// gPad->SetRightMargin(0.115);
-// gPad->SetLeftMargin(0.115);
-// if( hasMinimum() ) p_hist2D.SetMinimum(m_hist_min);
-// if( hasMaximum() ) p_hist2D.SetMaximum(m_hist_max);
-
-// p_hist2D.DrawCopy(draw_options.c_str());
-// //delete p_output;
-
- TH2D *h2 = IsGISAXSTools::getOutputDataTH2D(output, "p_hist2D");
- h2->SetContour(99);
+ TH1 *hist(0);
+ if(output.getNdimensions() == 2 ) {
+ hist = IsGISAXSTools::getOutputDataTH2D(output, "p_hist1D");
+ } else {
+ hist = IsGISAXSTools::getOutputDataTH123D(output, "p_hist1D");
+ }
+ hist->SetContour(99);
gStyle->SetPalette(1);
gStyle->SetOptStat(0);
-// gPad->SetRightMargin(0.115);
-// gPad->SetLeftMargin(0.115);
- if( hasMinimum() ) h2->SetMinimum(m_hist_min);
- if( hasMaximum() ) h2->SetMaximum(m_hist_max);
- h2->SetTitle(histogram_title.c_str());
- h2->DrawCopy(draw_options.c_str());
- delete h2;
+ if( hasMinimum() ) hist->SetMinimum(m_hist_min);
+ if( hasMaximum() ) hist->SetMaximum(m_hist_max);
+ hist->SetTitle(histogram_title.c_str());
+ hist->DrawCopy(draw_options.c_str());
+ delete hist;
}
@@ -178,6 +135,86 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
}
+
+/* ************************************************************************* */
+// create TH1D, TH2D or TH3D from OutputData
+/* ************************************************************************* */
+TH1 *IsGISAXSTools::getOutputDataTH123D(const OutputData<double>& output, const std::string &histo_name)
+{
+ if (output.getNdimensions() >3) {
+ std::cout << "IsGISAXSTools::getOutputDataTH123D() -> Warning! Expected number of dimensions should be not more than 3" << std::endl;
+ return 0;
+ }
+
+ std::vector<AxisStructure > haxises;
+ haxises.resize(output.getNdimensions());
+
+ // we assume variable bin size and prepare [nbins+1] array of low edges of each bin
+ for(size_t i_axis=0; i_axis<output.getNdimensions(); ++i_axis) {
+ const NamedVector<double> *axis = reinterpret_cast<const NamedVector<double>*>(output.getAxes()[i_axis]);
+ double dx(0);
+ haxises[i_axis].nbins = axis->getSize();
+ haxises[i_axis].name = axis->getName();
+ for(size_t i_bin=0; i_bin<axis->getSize(); ++i_bin) {
+ if(i_bin == 0) {
+ dx = (*axis)[1]-(*axis)[0];
+ } else {
+ dx = (*axis)[i_bin] - (*axis)[i_bin-1];
+ }
+ haxises[i_axis].xbins.push_back( (*axis)[i_bin] - dx/2.);
+ }
+ haxises[i_axis].xbins.push_back((*axis)[axis->getSize()-1] + dx/2.); // right bin edge of last bin, so for 100 bins size of vector will be 101
+ }
+
+// for(size_t i_axis=0; i_axis<2; ++i_axis) {
+// std::cout << "axis " << i_axis << " size:" << histo_axises[i_axis].size() << std::endl;
+// for(size_t i_bin=0; i_bin<histo_axises[i_axis].size(); ++i_bin) {
+// std::cout << histo_axises[i_axis][i_bin] << " ";
+// }
+// std::cout << std::endl;
+// }
+
+ // creation of 2D histogram with variable bin size
+ TH1 *hist;
+ if(output.getNdimensions() == 1) {
+ hist = new TH1D(histo_name.c_str(), histo_name.c_str(), haxises[0].nbins, &haxises[0].xbins[0]);
+ hist->GetXaxis()->SetTitle( haxises[0].name.c_str() );
+ } else if(output.getNdimensions() == 2) {
+ hist = new TH2D(histo_name.c_str(), histo_name.c_str(), haxises[0].nbins, &haxises[0].xbins[0], haxises[1].nbins, &haxises[1].xbins[0]);
+ hist->GetXaxis()->SetTitle( haxises[0].name.c_str() );
+ hist->GetYaxis()->SetTitle( haxises[1].name.c_str() );
+ } else if(output.getNdimensions() == 3) {
+ hist = new TH3D(histo_name.c_str(), histo_name.c_str(), haxises[0].nbins, &haxises[0].xbins[0], haxises[1].nbins, &haxises[1].xbins[0], haxises[1].nbins, &haxises[1].xbins[0]);
+ hist->GetXaxis()->SetTitle( haxises[0].name.c_str() );
+ hist->GetYaxis()->SetTitle( haxises[1].name.c_str() );
+ hist->GetZaxis()->SetTitle( haxises[2].name.c_str() );
+ } else {
+ throw LogicErrorException("IsGISAXSTools::getOutputDataTH123D() -> Warning! Wrong number of dimensions.");
+ }
+
+ OutputData<double>::const_iterator it = output.begin();
+ while (it != output.end())
+ {
+ std::vector<double > xyz;
+ for(size_t i_axis=0; i_axis<haxises.size(); ++i_axis) {
+ xyz.push_back(output.getValueOfAxis<double>( haxises[i_axis].name, it.getIndex() ) );
+ }
+ double value = *it++;
+ if(output.getNdimensions() == 1) dynamic_cast<TH1D *>(hist)->Fill(xyz[0], value);
+ if(output.getNdimensions() == 2) dynamic_cast<TH2D *>(hist)->Fill(xyz[0], xyz[1], value);
+ if(output.getNdimensions() == 3) dynamic_cast<TH3D *>(hist)->Fill(xyz[0], xyz[1], xyz[2], value);
+ }
+ hist->SetContour(50);
+ hist->SetStats(0);
+ hist->GetYaxis()->SetTitleOffset(1.1);
+
+ gStyle->SetPalette(1);
+ gStyle->SetOptStat(0);
+ return hist;
+}
+
+
+
/* ************************************************************************* */
// draw 1D distribution over values stored in OutputData
// binning of histogram is calculated on the fly
@@ -447,3 +484,62 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output
}
+// TODO re-implement sliceOutputData 1) using iterators 2) as a projection along of one axis having x1,x2 range defined for another axes 3) have sum (or sum averaged) within this [x1,x2] range
+// TODO implement creation of OutputData as a subset of bigger OutputData using [alphamin, alphamax][phimin,phimax] window
+/* ************************************************************************* */
+// slice 2D output data into 1D having fixed axes value on one of the axes
+/* ************************************************************************* */
+OutputData<double> *IsGISAXSTools::sliceOutputData(const OutputData<double > &data, const std::string &fixed_axis_name, double fixed_axis_value)
+{
+ if (data.getNdimensions() != 2) {
+ throw LogicErrorException("IsGISAXSTools::sliceOutputData() -> Error. Number of dimensions should be 2");
+ }
+ if( !data.getAxis(fixed_axis_name) ) {
+ throw LogicErrorException("IsGISAXSTools::sliceOutputData() -> Error. No axis with name "+fixed_axis_name);
+ }
+
+ OutputData<double > *sliced_data = new OutputData<double >;
+
+ const NamedVector<double> *fixed_axis(0);
+ int fixed_axis_index(-1);
+ for(size_t i_axis=0; i_axis<data.getNdimensions(); i_axis++) {
+ const NamedVector<double> *axis = dynamic_cast<const NamedVector<double>*>(data.getAxes()[i_axis]);
+ if( axis->getName() != fixed_axis_name ) {
+ sliced_data->addAxis(axis->clone());
+ } else {
+ fixed_axis = axis;
+ fixed_axis_index = i_axis;
+ }
+ }
+
+ // finding bin number on fixed_axis which is closest to fixed_axis_value
+ std::vector<double > buff;
+ buff.resize(fixed_axis->getSize(), 0);
+ for(size_t i=0; i<fixed_axis->getSize(); ++i) {
+ buff[i]= (*fixed_axis)[i];
+ }
+ std::vector<double >::iterator before = std::lower_bound(buff.begin(), buff.end(), fixed_axis_value);
+ int nbin(0);
+ if(before == buff.end() ) --before;
+ if(before == buff.begin() ) ++before;
+ std::vector<double >::iterator after = before;
+ --before;
+ ( *after-fixed_axis_value) < (fixed_axis_value - *before) ? nbin = std::distance(buff.begin(), after) : nbin = std::distance(buff.begin(), before);
+
+ // filling sliced data structure
+ OutputData<double>::const_iterator it_data = data.begin();
+ OutputData<double>::iterator it_sliced = sliced_data->begin();
+ while (it_data != data.end())
+ {
+ size_t fixed_axis_nbin = data.toCoordinates(it_data.getIndex())[fixed_axis_index];
+ if((int)fixed_axis_nbin == nbin) {
+ *it_sliced = *it_data;
+ ++it_sliced;
+ }
+ ++it_data;
+ }
+
+ return sliced_data;
+}
+
+
diff --git a/App/src/TestFittingModule2.cpp b/App/src/TestFittingModule2.cpp
index 5c649ee9be5..b764728dcb6 100644
--- a/App/src/TestFittingModule2.cpp
+++ b/App/src/TestFittingModule2.cpp
@@ -115,9 +115,9 @@ void TestFittingModule2::execute()
// Applying fit strategy: resizing real data
m_fitSuite->addFitStrategy(new FitSuiteStrategyAdjustData(3));
-// m_fitSuite->addFitStrategy(new FitSuiteStrategyAdjustData(2));
-// m_fitSuite->addFitStrategy(new FitSuiteStrategyAdjustData(1));
-// m_fitSuite->addFitStrategy(new FitSuiteStrategyDefault());
+ m_fitSuite->addFitStrategy(new FitSuiteStrategyAdjustData(2));
+ m_fitSuite->addFitStrategy(new FitSuiteStrategyAdjustData(1));
+ m_fitSuite->addFitStrategy(new FitSuiteStrategyDefault());
// Applying fit strategy: disturning data to get out of local minima
//m_fitSuite->addFitStrategy(new FitSuiteStrategyBootstrap());
diff --git a/App/src/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp
index eb47d0998fc..e15ac1e109a 100644
--- a/App/src/TestIsGISAXS12.cpp
+++ b/App/src/TestIsGISAXS12.cpp
@@ -17,6 +17,11 @@
#include "InterferenceFunction1DParaCrystal.h"
#include "GISASExperiment.h"
#include "DrawHelper.h"
+#include "FitSuite.h"
+#include "FitSuiteHelper.h"
+#include "ResolutionFunction2DSimple.h"
+#include "MathFunctions.h"
+#include "ROOTMinimizer.h"
#include <iostream>
#include <fstream>
@@ -30,7 +35,10 @@
-TestIsGISAXS12::TestIsGISAXS12() : IFunctionalTest("TestIsGISAXS12"), m_experiment(0)
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+TestIsGISAXS12::TestIsGISAXS12() : IFunctionalTest("TestIsGISAXS12"), m_experiment(0), m_sample_builder(0), m_fitSuite(0)
{
std::cout << "TestIsGISAXS12::TestIsGISAXS12() -> Info" << std::endl;
m_data_path = std::string(Utils::FileSystem::GetHomePath()+"./Examples/IsGISAXS_examples/ex-12/");
@@ -40,35 +48,48 @@ TestIsGISAXS12::TestIsGISAXS12() : IFunctionalTest("TestIsGISAXS12"), m_experime
TestIsGISAXS12::~TestIsGISAXS12()
{
delete m_experiment;
+ delete m_sample_builder;
+ delete m_fitSuite;
}
+/* ************************************************************************* */
+// execute
+/* ************************************************************************* */
void TestIsGISAXS12::execute()
{
- read_isgisaxs_datfile(m_data_path+"isgi_fitconstraints.dat");
+ // initializing experiment and sample builder
+ initialiseExperiment();
- m_experiment = new GISASExperiment(mp_options);
- m_experiment->setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
- m_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+ // run standard isgisaxs comparison for the sample we have
+ //run_isgisaxs_comparison();
- TestSampleBuilder builder;
+ // run test of our style fit
+ //run_test_fit();
- // run simulation for default sample parameters
- m_experiment->setSampleBuilder(&builder);
- m_experiment->runSimulation();
+ // run isgisaxs ex-12 style fit
+ run_isgisaxs_fit();
-
- IsGISAXSTools::writeOutputDataToFile(*(m_experiment->getOutputData()), m_data_path+"this_fitconstraints.ima");
}
-void TestIsGISAXS12::finalise()
+
+/* ************************************************************************* */
+// standard ixgisaxs comparison
+/* ************************************************************************* */
+void TestIsGISAXS12::run_isgisaxs_comparison()
{
- // ---------------------------------------
+ // run simulation for default sample parameters
+ m_experiment->runSimulation();
+ IsGISAXSTools::writeOutputDataToFile(*(m_experiment->getOutputData()), m_data_path+"this_fitconstraints.ima");
+
// plotting results of comparison we/isgisaxs for the sample with default parameters
std::string isgi_file(m_data_path+"isgi_fitconstraints.ima");
std::string this_file(m_data_path+"this_fitconstraints.ima");
+ // -------------
+ // plot results
+ // -------------
OutputData<double> *isgi_data = IsGISAXSTools::readOutputDataFromFile(isgi_file);
OutputData<double> *our_data = IsGISAXSTools::readOutputDataFromFile(this_file);
@@ -100,6 +121,83 @@ void TestIsGISAXS12::finalise()
delete isgi_data;
delete our_data;
+}
+
+
+/* ************************************************************************* */
+// run test fitting
+/* ************************************************************************* */
+void TestIsGISAXS12::run_test_fit()
+{
+ m_experiment->setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
+ m_experiment->setBeamIntensity(1e10);
+ m_experiment->runSimulation();
+ m_experiment->normalize();
+ OutputData<double > *test_real_data = createNoisyData( *m_experiment->getOutputData() );
+
+ // creating fit suite
+ m_fitSuite = new FitSuite();
+ m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
+
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius1", 0.2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius2", 0.2*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio1", 0.8*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio2", 0.8*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 6*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_probability", 12*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius1", 4*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+ m_fitSuite->addFitParameter("*SampleBuilder/particle_radius2", 4*Units::nanometer, 1*Units::nanometer, AttLimits::lowerLimited(0.01) );
+
+ m_fitSuite->attachObserver( new FitSuiteObserverPrint() );
+ m_fitSuite->attachObserver( new FitSuiteObserverDraw() );
+
+ m_fitSuite->addExperimentAndRealData(m_experiment, test_real_data);
+
+ m_fitSuite->runFit();
+
+ delete test_real_data;
+}
+
+
+/* ************************************************************************* */
+// run isgisaxs ex-12 style fit
+/* ************************************************************************* */
+void TestIsGISAXS12::run_isgisaxs_fit()
+{
+ // read_isgisaxs_datfile(m_data_path+"isgi_fitconstraints.dat");
+
+ // run simulation for default sample parameters
+ m_experiment->runSimulation();
+
+ TCanvas *c2 = DrawHelper::instance().createAndRegisterCanvas("TestIsGISAXS12_c2", "ex-12: Mixture of cylindrical particles with different size distribution");
+
+ c2->Divide(2,2);
+
+ IsGISAXSTools::setMinimum(1.);
+ // our calculations
+ c2->cd(1); gPad->SetLogz();
+ TH2D *hist1 = dynamic_cast<TH2D* >(IsGISAXSTools::getOutputDataTH123D(*m_experiment->getOutputData(), "hist1"));
+ hist1->SetMinimum(1.0);
+ hist1->Draw("CONT4 Z");
+
+ OutputData<double > *sliced_data = IsGISAXSTools::sliceOutputData(*m_experiment->getOutputData(), "phi_f", 0.01);
+ std::cout << "XXX " << sliced_data->getAllocatedSize() << std::endl;
+ c2->cd(2);
+ TH1D *hist2 = dynamic_cast<TH1D* >(IsGISAXSTools::getOutputDataTH123D(*sliced_data, "hist2"));
+ hist2->SetMinimum(1.0);
+ hist2->Draw();
+
+}
+
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+void TestIsGISAXS12::finalise()
+{
+ return;
// ----------------------------------------------------------
// plotting isgisaxs "experimental data" together with projections of our output_data for default sample parameters
TCanvas *c2 = DrawHelper::instance().createAndRegisterCanvas("TestIsGISAXS12_c2", "ex-12: experimental data to fit");
@@ -176,6 +274,43 @@ void TestIsGISAXS12::finalise()
}
+/* ************************************************************************* */
+// initialize experiment
+/* ************************************************************************* */
+void TestIsGISAXS12::initialiseExperiment()
+{
+ delete m_sample_builder;
+ m_sample_builder = new TestSampleBuilder();
+ delete m_experiment;
+ m_experiment = new GISASExperiment(mp_options);
+ m_experiment->setSampleBuilder(m_sample_builder);
+ m_experiment->setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
+ m_experiment->setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+ // no resolution function defined yet, since we want to run comparison with isgisaxs fitst
+}
+
+
+/* ************************************************************************* */
+// add noise to data
+/* ************************************************************************* */
+OutputData<double > *TestIsGISAXS12::createNoisyData(const OutputData<double> &exact_data, double noise_factor)
+{
+ OutputData<double > *real_data = exact_data.clone();
+
+ OutputData<double>::iterator it = real_data->begin();
+ while (it != real_data->end()) {
+ double current = *it;
+ double sigma = noise_factor*std::sqrt(current);
+ double random = MathFunctions::GenerateNormalRandom(current, sigma);
+ if (random<0.0) random = 0.0;
+ *it = random;
+ ++it;
+ }
+
+ return real_data;
+}
+
+
/* ************************************************************************* */
// read special isgisaxs *.dat file with data to fit
/* ************************************************************************* */
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index 736ecc7d7ac..dae5c6c339f 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -83,7 +83,6 @@ void GISASExperiment::runSimulation()
}
}
- std::cout << "OOO " << getOutputData()->getNdimensions() << " " << getOutputData()->getAllocatedSize() << std::endl;
m_detector.applyDetectorResolution(&m_intensity_map);
}
diff --git a/Core/Tools/inc/FitSuite.h b/Core/Tools/inc/FitSuite.h
index 4f60b019bb1..789809ac485 100644
--- a/Core/Tools/inc/FitSuite.h
+++ b/Core/Tools/inc/FitSuite.h
@@ -92,6 +92,9 @@ private:
FitSuite &operator=(const FitSuite &);
FitSuite(const FitSuite &);
+ //! check if all prerequisites to run fit fit are filled
+ bool check_prerequisites();
+
FitSuiteKit m_suite_kit; //! kit which contains pairs of <experiment,real_data> to fit
FitSuiteParameters m_fit_parameters; //! collection of fit parameters
fitstrategies_t m_fit_strategies; //! collection of strategies which are executed before every minimization round
diff --git a/Core/Tools/inc/OutputData.h b/Core/Tools/inc/OutputData.h
index b88fcdc8a7c..07980a51f48 100644
--- a/Core/Tools/inc/OutputData.h
+++ b/Core/Tools/inc/OutputData.h
@@ -474,7 +474,7 @@ bool OutputData<T>::hasSameShape(const OutputData<T> &right) const
{
if(!hasSameDimensions(right)) return false;
- // TODO: move check of consistency between dimensions of LLData and NamedVector into UNIT test
+ // TODO: move check of consistency between dimensions of LLData and NamedVector into UNIT test, if it's not there
if( (mp_ll_data->getRank() != m_value_axes.size()) || (right.mp_ll_data->getRank() != right.m_value_axes.size()) ) {
throw LogicErrorException("OutputData<T>::hasSameShape() -> Panic! Inconsistent dimensions in LLData and axes");
}
diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp
index f8c2022397c..b2aeb43de66 100644
--- a/Core/Tools/src/FitSuite.cpp
+++ b/Core/Tools/src/FitSuite.cpp
@@ -70,6 +70,7 @@ void FitSuite::link_fit_parameters()
{
// loop over all experiments defined
for(size_t i_exp = 0; i_exp<m_suite_kit.size(); ++i_exp) {
+ std::cout << "XXX 1.1" << i_exp << std::endl;
m_fit_parameters.link_to_experiment(m_suite_kit.getExperiment(i_exp));
}
}
@@ -87,17 +88,33 @@ void FitSuite::minimize()
for(size_t i_par = 0; i_par<m_fit_parameters.size(); i_par++) {
m_minimizer->setVariable(i_par, m_fit_parameters[i_par] );
}
+ if( m_fit_parameters.size() != m_minimizer->getNumberOfVariables()) std::cout << "FitSuite::minimize() -> Warning. Something unexpected" << std::endl;
// minimizing
m_minimizer->minimize();
}
+/* ************************************************************************* */
+// run fit
+/* ************************************************************************* */
+bool FitSuite::check_prerequisites()
+{
+ if( !m_minimizer ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No minimizer found.");
+ if( !m_suite_kit.size() ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No experiment defined");
+ if( !m_fit_parameters.size() ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No fit parameters defined");
+ return true;
+}
+
+
/* ************************************************************************* */
// run fit
/* ************************************************************************* */
void FitSuite::runFit()
{
+ // check if all prerequisites are fullfilled before starting minimization
+ check_prerequisites();
+
m_is_last_iteration = false;
// linking fit parameters with parameters defined in the experiment
diff --git a/Core/Tools/src/FitSuiteKit.cpp b/Core/Tools/src/FitSuiteKit.cpp
index 6e8dd295fd2..f8c11d71e7c 100644
--- a/Core/Tools/src/FitSuiteKit.cpp
+++ b/Core/Tools/src/FitSuiteKit.cpp
@@ -12,6 +12,9 @@ FitSuiteKit::KitItem::KitItem(Experiment *experiment, const OutputData<double >
, m_real_data(0)
, m_chi2_module(0)
{
+ if(!m_experiment) {
+ throw LogicErrorException("FitSuiteKit::KitItem::KitItem -> Error! Experiment can't be 0");
+ }
if(real_data) m_real_data = real_data->clone();
if(chi2_module) {
m_chi2_module = chi2_module->clone();
diff --git a/Core/Tools/src/FitSuiteStrategy.cpp b/Core/Tools/src/FitSuiteStrategy.cpp
index 4be0bc5c119..4d86b33e8d0 100644
--- a/Core/Tools/src/FitSuiteStrategy.cpp
+++ b/Core/Tools/src/FitSuiteStrategy.cpp
@@ -27,9 +27,9 @@ void FitSuiteStrategyDefault::execute()
/* ************************************************************************* */
// adjust (rebin) data before running fit suite minimization round
/* ************************************************************************* */
+// TODO: refactor this all
void FitSuiteStrategyAdjustData::execute()
{
- // TODO: refactor this all
if( !m_fit_suite ) throw NullPointerException("FitSuiteStrategyAdjustData::execute() -> FitSuite doesn't exists");
// if no data rediction was requested, just call FitSuite's minimization
--
GitLab