Skip to content
Snippets Groups Projects
Commit feb171c6 authored by pospelov's avatar pospelov
Browse files

Functional test for Convolve class, code cleanup.

parent 1570eb84
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
App/App.pro
+ 2
− 2
View file @ feb171c6
......@@ -17,7 +17,7 @@ SOURCES += \
src/CommandLine.cpp \
src/TestFactory.cpp \
src/TestDiffuseReflection.cpp \
src/TestInstrument.cpp
src/TestConvolution.cpp
HEADERS += \
inc/DrawHelper.h \
......@@ -32,7 +32,7 @@ HEADERS += \
inc/CommandLine.h \
inc/TestFactory.h \
inc/TestDiffuseReflection.h \
inc/TestInstrument.h
inc/TestConvolution.h
INCLUDEPATH += ./inc
DEPENDPATH += ./inc
......
App/inc/TestInstrument.h → App/inc/TestConvolution.h
+ 14
− 5
View file @ feb171c6
#ifndef TESTINSTRUMENT_H
#define TESTINSTRUMENT_H
#ifndef TESTCONVOLUTION_H
#define TESTCONVOLUTION_H
// ********************************************************************
// * The BornAgain project *
// * Simulation of neutron and x-ray scattering at grazing incidence *
......@@ -24,10 +24,10 @@
//! @class TestConvolution
//! @brief Testing Convolve class for instrumental effects studies
//- -------------------------------------------------------------------
class TestInstrument : public IFunctionalTest
class TestConvolution : public IFunctionalTest
{
public:
TestInstrument();
TestConvolution();
void execute();
......@@ -38,8 +38,17 @@ public:
void test_convolve2d();
private:
//! test function with many gaus'es on top of flat background for convolution studies
double fpeaks(double *x, double *par);
//! typedef of pair for the description of convolution mode
typedef std::pair<std::string, MathFunctions::Convolve::Mode> mode_pair_t;
//! vector of pairs defined above
std::vector<mode_pair_t> m_modes;
//! number of peaks for convolution test function
int m_npeaks;
};
#endif // TESTINSTRUMENT_H
#endif // TESTCONVOLUTION_H
App/src/TestInstrument.cpp → App/src/TestConvolution.cpp
+ 128
− 148
View file @ feb171c6
#include "TestInstrument.h"
#include "TestConvolution.h"
#include "MathFunctions.h"
#include "Convolve.h"
#include "DrawHelper.h"
......@@ -16,29 +16,14 @@
#include "TGraph.h"
#include "TBenchmark.h"
#include "TStyle.h"
#include "TText.h"
std::vector<double > Convolve_MyFFT1D_B(int mode, const std::vector<double> &data1, const std::vector<double> &data2);
int npeaks(30);
Double_t fpeaks(Double_t *x, Double_t *par) {
Double_t result = par[0] + par[1]*x[0];
for (Int_t p=0;p<npeaks;p++) {
Double_t norm = par[3*p+2];
Double_t mean = par[3*p+3];
Double_t sigma = par[3*p+4];
result += norm*TMath::Gaus(x[0],mean,sigma);
}
return result;
}
TestInstrument::TestInstrument()
TestConvolution::TestConvolution()
{
std::cout << "TestInstrument::TestInstrument() -> Info." << std::endl;
std::cout << "TestConvolution::TestConvolution() -> Info." << std::endl;
m_modes.push_back( mode_pair_t("LINEAR_FULL", MathFunctions::Convolve::FFTW_LINEAR_FULL));
m_modes.push_back( mode_pair_t("LINEAR_SAME_UNPADDED", MathFunctions::Convolve::FFTW_LINEAR_SAME_UNPADDED));
......@@ -48,107 +33,120 @@ TestInstrument::TestInstrument()
}
void TestInstrument::execute()
void TestConvolution::execute()
{
test_convolve1d();
// results are about
// LINEAR_FULL: Real Time = 0.37 seconds Cpu Time = 0.36 seconds
// LINEAR_SAME_UNPADDED: Real Time = 1.03 seconds Cpu Time = 1.03 seconds
// LINEAR_SAME: Real Time = 0.23 seconds Cpu Time = 0.23 seconds
// CIRCULAR_SAME: Real Time = 0.50 seconds Cpu Time = 0.50 seconds
// CIRCULAR_SAME_SHIFTED: Real Time = 0.51 seconds Cpu Time = 0.51 seconds
//test_convolve2d();
// LINEAR_FULL: Real Time = 1.19 seconds Cpu Time = 1.18 seconds
// LINEAR_SAME_UNPADDED: Real Time = 0.72 seconds Cpu Time = 0.72 seconds
// LINEAR_SAME: Real Time = 0.72 seconds Cpu Time = 0.73 seconds
// CIRCULAR_SAME: Real Time = 0.53 seconds Cpu Time = 0.54 seconds
// CIRCULAR_SAME_SHIFTED: Real Time = 0.61 seconds Cpu Time = 0.61 seconds
}
test_convolve2d();
/*
--- TestConvolution::test_convolve_1d() -> Peformance.
LINEAR_FULL: Real Time = 0.37 seconds Cpu Time = 0.37 seconds
LINEAR_SAME_UNPADDED: Real Time = 0.98 seconds Cpu Time = 0.98 seconds
LINEAR_SAME: Real Time = 0.22 seconds Cpu Time = 0.22 seconds
CIRCULAR_SAME: Real Time = 0.51 seconds Cpu Time = 0.50 seconds
CIRCULAR_SAME_SHIFTED: Real Time = 0.51 seconds Cpu Time = 0.51 seconds
--- TestConvolution::test_convolve_2d() -> Peformance.
LINEAR_FULL: Real Time = 1.26 seconds Cpu Time = 1.26 seconds
LINEAR_SAME_UNPADDED: Real Time = 0.76 seconds Cpu Time = 0.76 seconds
LINEAR_SAME: Real Time = 0.73 seconds Cpu Time = 0.73 seconds
CIRCULAR_SAME: Real Time = 0.57 seconds Cpu Time = 0.57 seconds
CIRCULAR_SAME_SHIFTED: Real Time = 0.62 seconds Cpu Time = 0.63 seconds
*/
}
void TestInstrument::test_convolve1d()
/* ************************************************************************* */
// test of convolution of two arrays in 1d
/* ************************************************************************* */
void TestConvolution::test_convolve1d()
{
npeaks = TMath::Abs(10);
double xmin(0.0);
double xmax(1000.);
size_t npoints = 501;
double dx = (xmax-xmin)/(npoints-1);
double sigma0=(xmax-xmin)*0.001;
double dx = (xmax-xmin)/double(npoints-1);
double sigmax=(xmax-xmin)*0.01; // precision of the measurement of x-value
TH1F *href = new TH1F("h","test", npoints-1, xmin, xmax);
href->SetMinimum(0.0);
href->SetMaximum(2.0);
//generate n peaks at random
Double_t par[3000];
// generate signal (n gaus peaks at random on top of flat background)
std::vector<double > signal;
m_npeaks = 10;
double *par = new double[m_npeaks*3+2];
par[0] = 1.0;
par[1] = -0.6/(xmax-xmin);
Int_t p;
for (p=0;p<npeaks;p++) {
par[3*p+2] = 1;
par[3*p+3] = (xmax-xmin)*0.01+gRandom->Rndm()*(xmax-xmin);
par[3*p+4] = sigma0+sigma0*5*gRandom->Rndm();
for (int i_par=0; i_par<m_npeaks;i_par++) {
par[3*i_par+2] = 1;
par[3*i_par+3] = (xmax-xmin)*0.01+gRandom->Rndm()*(xmax-xmin);
par[3*i_par+4] = sigmax/5.+sigmax*gRandom->Rndm();
}
TF1 *fspect = new TF1("f",fpeaks,0,1000,2+3*npeaks);
fspect->SetNpx(10000);
fspect->SetParameters(par);
TCanvas *c1 = new TCanvas("c1_test_convolve1d","c1_test_convolve1d",1024, 768);
c1->Divide(3,3);
c1->cd(1);
//href->FillRandom("f",20000);
href->DrawCopy();
fspect->Draw("pl same");
c1->cd(2);
double sigma_measure = sigma0*5;
TH1F *hmeasured = new TH1F("hp","hp",500, xmin, xmax);
// double sigma_measure=sigma0*3;
for(int i=0; i<1000000; i++) {
hmeasured->Fill(fspect->GetRandom(xmin, xmax) + gRandom->Gaus(0.0,sigma_measure));
}
hmeasured->Scale(1./hmeasured->GetEntries());
hmeasured->SetLineColor(kRed);
hmeasured->Draw();
std::vector<double> xspect;
xspect.resize(npoints);
TF1 *fsignal = new TF1("f",this, &TestConvolution::fpeaks, xmin, xmax,2+3*m_npeaks, "TestConvolution","fpeaks");
fsignal->SetNpx(10000);
fsignal->SetParameters(par);
signal.resize(npoints);
for(size_t i=0; i<npoints; i++) {
double x = xmin + i*dx;
xspect[i] = fspect->Eval(x);
std::cout << i << " " << x << " " << xspect[i] << std::endl;
signal[i] = fsignal->Eval(x);
}
// building responce function
std::vector<double> xresp;
xresp.resize(npoints);
// building kernel
std::vector<double> kernel;
kernel.resize(npoints);
TGraph *gr_kernel = new TGraph(npoints);
// for(size_t i=0; i<npoints/2; i++) {
// double x = i*dx;
// double y = std::exp(-x*x/sigma_measure/sigma_measure);
// xresp[i] = y;
// xresp[npoints-1-i] = xresp[i];
// kernel[i] = y;
// kernel[npoints-1-i] = kernel[i];
// }
for(size_t i=0; i<npoints; i++) {
double x = -(xmax-xmin)/2.+i*dx;
double y = std::exp(-x*x/sigma_measure/sigma_measure);
xresp[i] = y;
double y = std::exp(-x*x/sigmax/sigmax);
kernel[i] = y;
gr_kernel->SetPoint(i,x,kernel[i]);
}
TGraph *gr_resp = new TGraph(npoints);
for(size_t i=0; i<xresp.size(); i++){
double x = xmin + i*dx;
gr_resp->SetPoint(i,x,xresp[i]);
}
c1->cd(3);
// --------------
// drawing
// --------------
TCanvas *c1 = new TCanvas("c1_test_convolve1d","c1_test_convolve1d",1024, 768);
DrawHelper::instance().SetMagnifier(c1);
c1->Divide(3,3);
// drawing signal
c1->cd(1);
gPad->SetTopMargin(0.1);
TH1F *href = new TH1F("h","test", npoints-1, xmin, xmax);
href->SetMinimum(0.0);
href->SetMaximum(2.0);
href->SetStats(0);
href->DrawCopy();
gr_resp->Draw("pl same");
href->SetBit(TH1::kNoTitle, false);
fsignal->SetTitle("signal");
fsignal->Draw("pl same");
// drawing measured signal
c1->cd(2);
gPad->SetTopMargin(0.1);
TH1F *hmeasured = new TH1F("hmeasured","hmeasured",500, xmin, xmax);
for(int i=0; i<1000000; i++) {
hmeasured->Fill(fsignal->GetRandom(xmin, xmax) + gRandom->Gaus(0.0,sigmax));
}
hmeasured->Scale(1./hmeasured->GetEntries());
hmeasured->SetLineColor(kRed);
hmeasured->SetStats(0);
hmeasured->SetBit(TH1::kNoTitle, false);
hmeasured->Draw();
// drawing kernel
c1->cd(3);
gPad->SetTopMargin(0.1);
gr_kernel->SetTitle("kernel");
gr_kernel->Draw("apl");
// calculation of convolution
TBenchmark benchmark;
for(size_t i_mode=0; i_mode<m_modes.size(); i_mode++) {
std::string sname = m_modes[i_mode].first;
......@@ -160,34 +158,34 @@ void TestInstrument::test_convolve1d()
MathFunctions::Convolve cv;
cv.setMode(mode);
std::vector<double> xresult;
std::vector<double> result;
for(int i=0; i<1000; i++) {
//int mm = (int)modes[i_mode];
// xresult = Convolve_MyFFT1D_B(mm, xspect, xresp);
cv.fftconvolve(xspect, xresp, xresult);
cv.fftconvolve(signal, kernel, result);
}
benchmark.Stop(sname.c_str());
benchmark.Show(sname.c_str());
for(size_t i=0; i<xresult.size(); i++) {
//benchmark.Show(sname.c_str());
for(size_t i=0; i<result.size(); i++) {
double x = xmin + i*dx;
gr_result->SetPoint(i,x,xresult[i]);
gr_result->SetPoint(i,x,result[i]);
}
c1->cd(4+i_mode);
gPad->SetTopMargin(0.1);
href->SetMaximum(20.0);
href->SetTitle(sname.c_str());
href->DrawCopy();
gr_result->Draw("pl same");
}
float_t rp, cp;
std::cout << "--------------" << std::endl;
std::cout << "--- TestConvolution::test_convolve_1d() -> Peformance." << std::endl;
benchmark.Summary(rp, cp);
}
/* ************************************************************************* */
//
// test of convolution of two arrays in 2d
/* ************************************************************************* */
void TestInstrument::test_convolve2d()
void TestConvolution::test_convolve2d()
{
double xmin(0.0), xmax(100.);
size_t npx = 100;
......@@ -211,13 +209,12 @@ void TestInstrument::test_convolve2d()
// spikes
if(ix>1 && iy!=0) {
if( (ix%20==0 || (ix-1)%20==0) && iy%8==0) {
std::cout << ix << " " << iy << std::endl;
source[ix][iy]+=5;
}
}
} // iy
} // ix
// signal histogram
// filling signal histogram with values from signal array
TH2F *h2_signal = new TH2F("h2_signal","h2_signal", npx, xmin, xmax, npy, ymin, ymax);
h2_signal->SetStats(0);
h2_signal->SetMinimum(0);
......@@ -236,6 +233,7 @@ void TestInstrument::test_convolve2d()
// preparing resolution kernel and histogram
std::vector<std::vector<double> > kernel;
TH2F *h2_kernel = new TH2F("h2_kernel","h2_kernel", npx, xmin, xmax, npy, ymin, ymax);
h2_kernel->SetStats(0);
kernel.resize(npx);
for(size_t ix=0; ix<npx; ix++) {
kernel[ix].resize(npy);
......@@ -249,6 +247,7 @@ void TestInstrument::test_convolve2d()
// simulation of measurements with smearing measurement with resolution function
TH2F *h2_measured = new TH2F("h2_measured","h2_measured", npx, xmin, xmax, npy, ymin, ymax);
h2_measured->SetStats(0);
TRandom mr;
for(int i_meas=0; i_meas<1000000; i_meas++) {
// generating distribution according to our signal
......@@ -264,16 +263,24 @@ void TestInstrument::test_convolve2d()
gStyle->SetPalette(1);
c1->Divide(3,3);
c1->cd(1); gPad->SetRightMargin(0.15);
c1->cd(1);
gPad->SetRightMargin(0.12);
gPad->SetTopMargin(0.1);
h2_signal->Draw("colz");
c1->cd(2); gPad->SetRightMargin(0.15);
c1->cd(2);
gPad->SetRightMargin(0.12);
gPad->SetTopMargin(0.1);
h2_signal->Draw("lego");
c1->cd(3); gPad->SetRightMargin(0.15);
c1->cd(3);
gPad->SetRightMargin(0.12);
gPad->SetTopMargin(0.1);
h2_kernel->Draw("lego");
c1->cd(4); gPad->SetRightMargin(0.15);
c1->cd(4);
gPad->SetRightMargin(0.12);
gPad->SetTopMargin(0.1);
h2_measured->Draw("colz");
// calculation of convolution
......@@ -294,6 +301,8 @@ void TestInstrument::test_convolve2d()
// drawing results of last convolution
c1->cd(5+i_mode);
gPad->SetRightMargin(0.12);
gPad->SetTopMargin(0.1);
TH2F h2_result("h2_result","h2_result", npx, xmin, xmax, npy, ymin, ymax);
h2_result.SetStats(0);
for(size_t ix=0; ix<npx; ix++) {
......@@ -304,61 +313,32 @@ void TestInstrument::test_convolve2d()
}
}
h2_result.SetTitle(sname.c_str());
h2_result.SetBit(TH1::kNoTitle, false);
h2_result.DrawCopy("colz");
}
// benchmark summary
float_t rp, cp;
std::cout << "--------------" << std::endl;
std::cout << "--- TestConvolution::test_convolve_2d() -> Peformance." << std::endl;
benchmark.Summary(rp, cp);
}
std::vector<double> Convolve_MyFFT1D_B(int mode, const std::vector<double> &data1, const std::vector<double> &data2)
/* ************************************************************************* */
// test function for convolution
/* ************************************************************************* */
double TestConvolution::fpeaks(double *x, double *par)
{
size_t Ns = data1.size();
double * src = new double[Ns];
for(size_t i=0; i<Ns; i++) {
src[i] = data1[i];
Double_t result = par[0] + par[1]*x[0];
for (Int_t p=0;p<m_npeaks;p++) {
Double_t norm = par[3*p+2];
Double_t mean = par[3*p+3];
Double_t sigma = par[3*p+4];
result += norm*TMath::Gaus(x[0],mean,sigma);
}
size_t Nk = data2.size();
double * kernel = new double[Nk];
for(size_t i=0; i<Nk; i++) {
kernel[i] = data2[i];
}
FFTW_Workspace ws;
// init_workspace_fftw(ws, FFTW_LINEAR_FULL, 1,Ns,1,Nk);
// init_workspace_fftw(ws, FFTW_LINEAR_SAME_UNPADDED, 1,Ns,1,Nk);
// init_workspace_fftw(ws, FFTW_LINEAR_SAME, 1,Ns,1,Nk);
// init_workspace_fftw(ws, FFTW_LINEAR_VALID, 1,Ns,1,Nk);
// init_workspace_fftw(ws, FFTW_CIRCULAR_SAME, 1,Ns,1,Nk);
// init_workspace_fftw(ws, FFTW_CIRCULAR_SAME_SHIFTED, 1,Ns,1,Nk);
FFTW_Convolution_Mode xm = (FFTW_Convolution_Mode)mode;
init_workspace_fftw(ws, xm, 1,Ns,1,Nk);
fftw_convolve(ws, src, kernel);
// std::cout << "!!! " << ws.w_dst << std::endl;
// for(int j = 0 ; j < ws.w_dst ; ++j) std::cout << ws.dst[j] << " ";
std::vector<double > fft_result;
fft_result.resize(ws.w_dst);
for(size_t i=0; i<data1.size(); i++){
fft_result[i] = ws.dst[i];
}
clear_workspace_fftw(ws);
return fft_result;
return result;
}
App/src/TestFactory.cpp
+ 7
− 7
View file @ feb171c6
......@@ -4,7 +4,7 @@
#include "TestFormFactor.h"
#include "TestDWBAFormFactor.h"
#include "TestDiffuseReflection.h"
#include "TestInstrument.h"
#include "TestConvolution.h"
#include "TBenchmark.h"
......@@ -15,12 +15,12 @@ TestFactory::TestFactory() : m_benchmark(0)
setStoreObjects(true);
setDeleteObjects(true);
registerItem("roughness", IFactoryCreateFunction<TestRoughness, IFunctionalTest> );
registerItem("fresnel", IFactoryCreateFunction<TestFresnelCoeff, IFunctionalTest> );
registerItem("formfactor", IFactoryCreateFunction<TestFormFactor, IFunctionalTest> );
registerItem("dwba", IFactoryCreateFunction<TestDWBAFormFactor, IFunctionalTest> );
registerItem("diffuse", IFactoryCreateFunction<TestDiffuseReflection, IFunctionalTest> );
registerItem("instrument", IFactoryCreateFunction<TestInstrument, IFunctionalTest> );
registerItem("roughness", IFactoryCreateFunction<TestRoughness, IFunctionalTest> );
registerItem("fresnel", IFactoryCreateFunction<TestFresnelCoeff, IFunctionalTest> );
registerItem("formfactor", IFactoryCreateFunction<TestFormFactor, IFunctionalTest> );
registerItem("dwba", IFactoryCreateFunction<TestDWBAFormFactor, IFunctionalTest> );
registerItem("diffuse", IFactoryCreateFunction<TestDiffuseReflection, IFunctionalTest> );
registerItem("convolution", IFactoryCreateFunction<TestConvolution, IFunctionalTest> );
m_benchmark = new TBenchmark();
}
......
Core/inc/Convolve.h
+ 3
− 49
View file @ feb171c6
......@@ -52,6 +52,7 @@ public:
typedef std::vector<double1d_t> double2d_t;
//! convolution modes
//! use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance
enum Mode { FFTW_LINEAR_FULL, FFTW_LINEAR_SAME_UNPADDED, FFTW_LINEAR_SAME, FFTW_LINEAR_VALID, FFTW_CIRCULAR_SAME, FFTW_CIRCULAR_SAME_SHIFTED, FFTW_UNDEFINED };
//! convolution in 1D
......@@ -67,12 +68,7 @@ public:
void setMode(Mode mode) { m_mode = mode; }
private:
//! compute convolution of source and kernel using fast fourier transformation
void fftw_convolve(const double2d_t &source, const double2d_t &kernel);
//void fftw_convolve(double * src,double * kernel);
//! compute circual convolution of source and kernel using fast fourier transformation
//void fftw_circular_convolution(double * src, double * kernel);
void fftw_circular_convolution(const double2d_t &source, const double2d_t &kernel);
//! find closest number X>n, which can be factorised according to fftw3 favorite factorisation
......@@ -82,7 +78,7 @@ private:
bool is_optimal(int n);
//! Workspace contains input (source and kernel), intermediate and output arrays to run convolution via fft
//! 'source' it is our signal, 'kernel' it is our resolution (or delta-responce) function
//! 'source' it is our signal, 'kernel' it is our resolution (also known as delta-responce) function
//! Sizes of input arrays are adjusted; output arrays are alocated via fftw3 allocation for maximum performance
class Workspace
{
......@@ -102,6 +98,7 @@ private:
double *dst_fft; // result of production of FFT(source) and FFT(kernel)
int h_dst, w_dst; // size of resulting array
double *dst; // The array containing the result
int h_offset, w_offset; // offsets to copy result into output arrays
fftw_plan p_forw_src;
fftw_plan p_forw_kernel;
fftw_plan p_back;
......@@ -116,47 +113,4 @@ private:
} // namespace MathFunctions
typedef enum
{
FFTW_LINEAR_FULL,
FFTW_LINEAR_SAME_UNPADDED,
FFTW_LINEAR_SAME,
FFTW_LINEAR_VALID,
FFTW_CIRCULAR_SAME,
FFTW_CIRCULAR_SAME_SHIFTED,
FFTW_CIRCULAR_CUSTOM
} FFTW_Convolution_Mode;
typedef struct FFTW_Workspace
{
//fftw_complex * in_src, *out_src;
double * in_src, *out_src, *in_kernel, *out_kernel;
int h_src, w_src, h_kernel, w_kernel;
int w_fftw, h_fftw;
FFTW_Convolution_Mode mode;
double * dst_fft;
double * dst; // The array containing the result
int h_dst, w_dst; // its size ; This is automatically set by init_workspace, despite you use FFTW_CIRCULAR_CUSTOM
fftw_plan p_forw_src;
fftw_plan p_forw_kernel;
fftw_plan p_back;
} FFTW_Workspace;
extern void factorize (const int n,
int *n_factors,
int factors[],
int * implemented_factors);
extern bool is_optimal(int n, int * implemented_factors);
extern int find_closest_factor(int n, int * implemented_factor);
extern void init_workspace_fftw(FFTW_Workspace & ws, FFTW_Convolution_Mode mode, int h_src, int w_src, int h_kernel, int w_kernel);
extern void clear_workspace_fftw(FFTW_Workspace & ws);
extern void fftw_convolve(FFTW_Workspace &ws, double * src,double * kernel);
extern void fftw_circular_convolution(FFTW_Workspace &ws, double * src, double * kernel);
#endif // CONVOLVE_H
Core/src/Convolve.cpp
+ 108
− 448
View file @ feb171c6
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment