Dear co-rooters,
I am writing a function/macro that will be either used standalone or will be called from another function/macro.
This function makes some matrix computations using ublas and should either return a matrix or nothing (for instance a zero matrix).
The way I thought it would work, was to include a bool argument and if it’s true a matrix will be returned otherwise another zero matrix would be the output using the following check at the end of the function
if (return_TMatrix) return COV_USER;
else {zero_matrix<float> m_ZERO(1, 1); return m_ZERO;}
The thing is that when I run the code it breaks… I guess I am trying to do something strange?
Any idea or advice?
The code can be found below
//#include "/afs/cern.ch/user/a/astamato/public/lib_thanos/c_header.h"
//#include "/afs/cern.ch/user/a/astamato/public/lib_thanos/root_header.h"
//ROOT6
// source /afs/cern.ch/sw/lcg/external/gcc/4.8/x86_64-slc6/setup.csh
// source /afs/cern.ch/sw/lcg/app/releases/ROOT/6.02.02/x86_64-slc6-gcc48-opt/root/bin/thisroot.csh
//ROOT6 on lxplus7
// source /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.12.06/x86_64-centos7-gcc48-opt/root/bin/thisroot.csh
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>// std::max()
#include "TROOT.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include "TMatrixF.h"
#include "TMatrixFSym.h"
#include "TMatrixFSparse.h"
#include "TMatrixTSym.h"
#include "TString.h"
#include "TRandom.h"
#include </usr/include/boost/numeric/ublas/matrix.hpp>
#include </usr/include/boost/numeric/ublas/matrix_sparse.hpp>
#include </usr/include/boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include </usr/include/boost/numeric/ublas/io.hpp>
#include "covariance_lib.cc"
using namespace std;
using namespace boost::numeric::ublas;
// covariance_linear_interpolation("ENDF_B-VII_1_reduced.dat", -2, 0, 100, 0)
boost::numeric::ublas::symmetric_matrix<float> covariance_linear_interpolation(TString filename, int E_MIN, int E_MAX, int N_BPDEC, bool return_TMatrix){
// (1) Define the isolethargic binning - MeV
int ndec = E_MAX - E_MIN;
int nbins = (int) ndec*N_BPDEC;
float step = (float) ndec / nbins;
float* En_USER = new float[nbins];
for(int i=0; i < nbins; i++) {
En_USER[i] = (float) pow(10., step * (float) i + E_MIN);
}
// (2) Read the evaluation and interpolate
int nbins_EVAL = count_lines(filename);
float* xs_USER = new float[nbins];
float* En_EVAL = new float[nbins_EVAL];
float* xs_EVAL = new float[nbins_EVAL];
int* idx = new int[nbins];
input_ascii(filename, En_EVAL, xs_EVAL);
linear_interpolation(nbins, nbins_EVAL, En_USER, xs_USER, En_EVAL, xs_EVAL, idx);
// (3) Define the size for the covariance matrices
const int size_EVAL = nbins_EVAL;
const int size_USER = nbins;
// (4) Define the types of T-objects that will be used
float *covariance_EVAL = new float[size_EVAL*(size_EVAL+1)/2];
for (int i=0; i<size_EVAL*(size_EVAL+1)/2; ++i)
covariance_EVAL[i] = 1.e-3*i;
symmetric_matrix<float, upper> COV_EVAL (size_EVAL, size_EVAL);
COV_EVAL = fill_symmetric( COV_EVAL, covariance_EVAL, size_EVAL);
delete[] covariance_EVAL;
delete[] xs_USER;
delete[] xs_EVAL;
// (5) Define the non-zero elements of the SParse TMatrix
int* icol = new int[2*(nbins)];
int* irow = new int[2*(nbins)];
float* coeff = new float[2*(nbins)];
for (int i=0; i<(nbins); ++i){
icol[2*i] = idx[i];
icol[2*i+1] = idx[i] + 1;
irow[2*i] = i;
irow[2*i+1] = i;
coeff[2*i] = 1. + ( En_USER[i] - En_EVAL[ idx[i] ] )/( En_EVAL[idx[i]]-En_EVAL[idx[i+1]] );
coeff[2*i+1] = coeff[2*i] - 1;
if ( !isfinite(coeff[2*i]) ){// If En_USER == En_EVAL then sensitivity coefficient is 1
coeff[2*i] = 1;
coeff[2*i+1] = 0;
}
}
delete[] En_USER;
delete[] En_EVAL;
delete[] idx;
mapped_matrix<float> G(size_USER, size_EVAL,2*nbins);
G = fill_sparse( G, 2*nbins, irow, icol, coeff);
delete[] icol;
delete[] irow;
delete[] coeff;
// (6) Calculate the covariance matrix - will be symmetric
matrix<float> COV_USER;
COV_USER = prod( G,COV_EVAL );
COV_USER = prod( COV_USER, trans(G) );
if (return_TMatrix) return COV_USER;
else {zero_matrix<float> m_ZERO(1, 1); return m_ZERO;}
}
covariance_lib.cc
// Inputs : -name = name of the ascii file to be read - two column file with x and y values, rescpectively assumed
// -X_USER = array that contains the user defined energy binning
// -output = if output = TRUE exports the interpolation array
// if output = FALSE exports the index of the input values used in the interpolation
void linear_interpolation(int size_USER, int size_EVAL, float* X_USER, float* Y_USER, float* X_EVAL, float* Y_EVAL, int* index){
float x1, x2,y1, y2, a, b;
bool boo;
x1 = 0.; x2 = 0.; y1 = 0; y2 = 0.; a = 0.; b = 0.;
for(int i=0; i<=size_USER; ++i){
boo = true;
for (int j=0; j<=size_EVAL; j++){
if (boo==true){
if(X_USER[i] == X_EVAL[j]){
Y_USER[i] = Y_EVAL[j];
index[i] = j;
boo = false;
}
else if(X_USER[i]>X_EVAL[j]){
x1 = X_EVAL[j];
y1 = Y_EVAL[j];
index[i] = j;
}
else if( X_USER[i]<X_EVAL[j] && boo == true ){
x2 = X_EVAL[j];
y2 = Y_EVAL[j];
a = (y2-y1)/(x2-x1);
b = y2-a*x2;
Y_USER[i] = a*X_USER[i]+b;
//index[i] = j;
boo = false;
}
}
}//end of loop over EVAL
}//end of loop over USER
}//____linear_interpolation()
void input_ascii(TString filename, float* X_EVAL, float* Y_EVAL){
std::ifstream myfile;
myfile.open(filename);
float x = 0, y = 0;
int i = 0;
while(1){
myfile >> x >> y;
if (!myfile.good()) break;
X_EVAL[i] = x;
Y_EVAL[i] = y;
i++;
}
myfile.close();
}//___input_ascii()
int count_lines(TString filename){
float x, y;
int nlines;
// (1) Open the file to count its lines
std::ifstream myfile;
myfile.open(filename);
if(!myfile){
cout << "File not found" << endl;
}
nlines = 0;
while(1){
myfile >> x >> y;
if (!myfile.good()) break;
nlines++;
}
myfile.close();
return nlines;
}//___count_lines()
std::vector<float> isolethargic_binning(int min, int max, int bpd){
int ndec = max - min;
int nbins = (int) ndec*bpd;
float step = (float) ndec / nbins;
std::vector<float> xbins(nbins+1);
for(int i=0; i <= nbins; i++) {
xbins[i] = (float) pow(10., step * (float) i + min);
}
return xbins;
}//___isolethargic_binning()
boost::numeric::ublas::symmetric_matrix<float> fill_symmetric (boost::numeric::ublas::symmetric_matrix<float> m_sym, float* filler, int size_FILLER){
float* in = filler;
for (size_t i=0; i<m_sym.size1(); ++ i)
for (size_t j = 0; j <= i && in != &filler[size_FILLER+1]; ++ j)
m_sym (i, j) = *in++;
return m_sym;
}//___fill_symmetric()
boost::numeric::ublas::mapped_matrix<float> fill_sparse (boost::numeric::ublas::mapped_matrix<float> m_spar, int size_DATA, int* row, int* col, float* data){
for (int i=0; i<size_DATA; ++ i){
m_spar( row[i], col[i] ) = data[i];}
return m_spar;
}//___fill_sparse()
boost::numeric::ublas::diagonal_matrix<float> fill_diagonal (boost::numeric::ublas::diagonal_matrix<float> m_diag, float* data) {
for (int i = 0; i < (int) m_diag.size1(); ++ i)
m_diag (i, i) = data[i];
return m_diag;
}//___fill_diagonal()
A sample file can be found here. I am running ROOT6.12.04 on lxplus7 and compiling the code in ACLiC using
.L covariance_linear-interpolation.C++
The code can be run using the following
covariance_linear_interpolation("ENDF_B-VII_1_reduced.dat", -2, 0, 100, 0)