CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Data2D_extra.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2016 by Federico Marulli *
3  * federico.marulli3@unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
34 #include "Data2D_extra.h"
35 
36 using namespace std;
37 
38 using namespace cbl;
39 using namespace data;
40 using namespace glob;
41 
42 
43 // ======================================================================================
44 
45 
46 void cbl::data::Data2D_extra::read (const string input_file, const int skip_nlines, const vector<int> column, const vector<int> column_data, const vector<int> column_errors, const std::vector<int> column_edges)
47 {
48  if (column_data.size()!=column_errors.size()) ErrorCBL("the sizes of column_data and columt_errors must be equal!", "read", "Data2D_extra.cpp");
49 
50  // default column of input data values
51  int column_data_default = column[1]+1;
52 
53  // default column of input error values
54  int column_error_default = column[1]+2;
55 
56  // default first column of input edge values
57  int column_edges_default = column[1]+3;
58 
59  ifstream fin(input_file.c_str()); checkIO(fin, input_file);
60  string line;
61 
62  // loop on the data (and, possibly, error) vectors to be read; if
63  // more than one data vectors are read, the data (and, possibly,
64  // error) vectors will be added one after the other, in order to
65  // construct a single data object
66  const int cl_max = (column_data.size()<2) ? 1 : column_data.size();
67  std::vector<double> dummy_edges_xx, dummy_edges_yy;
68  for (int cl=0; cl<cl_max; ++cl) {
69 
70  // skip the first nlines (in case of header lines)
71  if (skip_nlines>0)
72  for (int i=0; i<skip_nlines; ++i)
73  getline(fin, line);
74 
75 
76  // read the file lines
77 
78  double last_bin_edge_xx = par::defaultDouble;
79  double last_bin_edge_yy = par::defaultDouble;
80 
81  while (getline(fin, line)) {
82 
83  stringstream ss(line);
84  vector<double> num; double NUM;
85  while (ss>>NUM) num.emplace_back(NUM);
86 
87  // if column[0] is larger than 0, the column of x coordinates is
88  // the one specified in input
89  const int ind_x = max(0, column[0]-1);
90  checkDim(num, ind_x+1, "num", false);
91  m_x.emplace_back(num[ind_x]);
92 
93  // if column[1] is larger than 0, the column of y coordinates is
94  // the one specified in input
95  const int ind_y = max(0, column[1]-1);
96  checkDim(num, ind_y+1, "num", false);
97  m_y.emplace_back(num[ind_y]);
98 
99  // if the size of the column_data vector is 0, the columns of
100  // data values are the ones specified in input
101  const int ind_data = (column_data.size()==0) ? column_data_default-1 : column_data[cl]-1;
102  checkDim(num, ind_data+1, "num", false);
103  m_data.emplace_back(num[ind_data]);
104 
105  // if the size of the column_errors vector is 0, the columns of
106  // error values are the ones specified in input; if the error
107  // column is not present, the errors will be set to 1, by
108  // default
109  const size_t ind_error = (column_errors.size()==0) ? column_error_default-1 : column_errors[cl]-1;
110  if (num.size()<ind_error+1 && column_errors.size()>1)
111  WarningMsgCBL("the errors cannot be retrieved from the provided input file, and will be set to 1", "read", "Data2D_extra.cpp");
112  m_error.emplace_back((num.size()<ind_error+1) ? 1 : num[ind_error]);
113 
114  // if the size of the column_edges vector is 0, the columns of
115  // bin edges values are the ones specified in input
116  const size_t ind_edges = (column_edges.size()==0) ? column_edges_default-1 : column_edges[cl]-1;
117  if (num.size()>ind_edges+1) {
118  dummy_edges_xx.emplace_back(num[ind_edges]);
119  last_bin_edge_xx = num[ind_edges+1];
120  }
121 
122  // if the size of the column_edges vector is 0, the columns of
123  // bin edges values are the ones specified in input
124  const size_t ind_edges_yy = (column_edges.size()==0) ? column_edges_default-1+2 : column_edges[cl]-1+2;
125  if (num.size()>ind_edges_yy+1) {
126  dummy_edges_yy.emplace_back(num[ind_edges_yy]);
127  last_bin_edge_yy = num[ind_edges_yy+1];
128  }
129 
130  // get the extra info
131  int ind = 0;
132  for (size_t ex=num.size()-m_extra_info.size(); ex<num.size(); ++ex)
133  m_extra_info[ind++].emplace_back(num[ex]);
134  }
135 
136  std::sort( dummy_edges_xx.begin(), dummy_edges_xx.end() );
137  dummy_edges_xx.erase( std::unique( dummy_edges_xx.begin(), dummy_edges_xx.end() ), dummy_edges_xx.end() );
138  std::sort( dummy_edges_yy.begin(), dummy_edges_yy.end() );
139  dummy_edges_yy.erase( std::unique( dummy_edges_yy.begin(), dummy_edges_yy.end() ), dummy_edges_yy.end() );
140 
141  for (size_t i=0; i<dummy_edges_xx.size(); i++)
142  m_edges_xx.emplace_back(dummy_edges_xx[i]);
143  m_edges_xx.emplace_back(last_bin_edge_xx);
144  for (size_t i=0; i<dummy_edges_yy.size(); i++)
145  m_edges_yy.emplace_back(dummy_edges_yy[i]);
146  m_edges_yy.emplace_back(last_bin_edge_yy);
147 
148  // go back at the beginning of the file
149  fin.clear(); fin.seekg(ios::beg);
150 
151  column_data_default += 2;
152  column_error_default += 2;
153  }
154  fin.clear(); fin.close();
155 
156 
157  // set the data size and diagonal covariance
158 
159  unique_unsorted(m_x);
160  unique_unsorted(m_y);
161 
162  m_xsize = m_x.size();
163  m_ysize = m_y.size();
164  m_ndata = m_data.size();
165 
166  m_covariance.resize(m_ndata, vector<double>(m_ndata, 0));
167  for (int i=0; i<m_ndata; i++)
168  m_covariance[i][i] = pow(m_error[i], 2);
169 }
170 
171 
172 // ======================================================================================
173 
174 
175 void cbl::data::Data2D_extra::Print (const int precision) const
176 {
177  for (int i=0; i<m_xsize; ++i)
178  for (int j=0; j<m_ysize; ++j) {
179  int index = j+m_ysize*i;
180  coutCBL << setprecision(precision) << setw(8) << right << m_x[i]
181  << " " << setprecision(precision) << setw(8) << right << m_y[j]
182  << " " << setprecision(precision) << setw(8) << right << m_data[index]
183  << " " << setprecision(precision) << setw(8) << right << m_error[index];
184 
185  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
186  coutCBL << " " << setprecision(precision) << setw(8) << right << m_extra_info[ex][index];
187  cout << endl;
188  }
189 }
190 
191 
192 // ======================================================================================
193 
194 
195 void cbl::data::Data2D_extra::write (const string dir, const string file, const string header, const bool full, const int prec, const int ww, const int rank) const
196 {
197  (void)rank;
198 
199  string file_out = dir+file;
200  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
201 
202  if (header!=par::defaultString)
203  fout << "### " << header << " ###" << endl;
204 
205  const int bp = std::cout.precision();
206 
207  for (int i=0; i<m_xsize; ++i)
208  for (int j=0; j<m_ysize; ++j) {
209  int index = j+m_ysize*i;
210  cbl::Print(m_x[i], prec, ww, "", " ", false, fout);
211  cbl::Print(m_y[j], prec, ww, "", " ", false, fout);
212  cbl::Print(m_data[index], prec, ww, "", " ", false, fout);
213  cbl::Print(m_error[index], prec, ww, "", " ", false, fout);
214  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
215  cbl::Print(m_extra_info[ex][index], prec, ww, "", " ", false, fout);
216  fout << endl;
217  }
218 
219 
220  if (full) { // duplicate the information in the other 3 quadrants
221 
222  for (int i=0; i<m_xsize; ++i)
223  for (int j=0; j<m_ysize; ++j) {
224  int index = j+m_ysize*i;
225  cbl::Print(m_x[i], prec, ww, "", " ", false, fout);
226  cbl::Print(-m_y[j], prec, ww, "", " ", false, fout);
227  cbl::Print(m_data[index], prec, ww, "", " ", false, fout);
228  cbl::Print(m_error[index], prec, ww, "", " ", false, fout);
229  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
230  cbl::Print(m_extra_info[ex][index], prec, ww, "", " ", false, fout);
231  fout << endl;
232  }
233 
234  for (int i=0; i<m_xsize; ++i)
235  for (int j=0; j<m_ysize; ++j) {
236  int index = j+m_ysize*i;
237  cbl::Print(-m_x[i], prec, ww, "", " ", false, fout);
238  cbl::Print(-m_y[j], prec, ww, "", " ", false, fout);
239  cbl::Print(m_data[index], prec, ww, "", " ", false, fout);
240  cbl::Print(m_error[index], prec, ww, "", " ", false, fout);
241  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
242  cbl::Print(m_extra_info[ex][index], prec, ww, "", " ", false, fout);
243  fout << endl;
244  }
245 
246  for (int i=0; i<m_xsize; ++i)
247  for (int j=0; j<m_ysize; ++j) {
248  int index = j+m_ysize*i;
249  cbl::Print(-m_x[i], prec, ww, "", " ", false, fout);
250  cbl::Print(m_y[j], prec, ww, "", " ", false, fout);
251  cbl::Print(m_data[index], prec, ww, "", " ", false, fout);
252  cbl::Print(m_error[index], prec, ww, "", " ", false, fout);
253  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
254  cbl::Print(m_extra_info[ex][index], prec, ww, "", " ", false, fout);
255  fout << endl;
256  }
257  }
258 
259  std::cout.precision(bp);
260  fout.close(); cout << endl; coutCBL << "I wrote the file: " << file_out << endl << endl;
261 }
262 
263 
264 // ======================================================================================
265 
266 
267 shared_ptr<Data> cbl::data::Data2D_extra::cut (const double xmin, const double xmax, const double ymin, const double ymax) const
268 {
269  vector<bool> mask(m_ndata, true);
270  vector<double> xx, yy, edges_xx, edges_yy;
271 
272  for (int i=0; i<m_xsize; i++)
273  for (int j=0; j<m_ysize; j++)
274  if (((m_x[i]<xmin) || (m_x[i]>xmax)) || ((m_y[j]<ymin) || (m_y[j]>ymax)))
275  mask[j+m_ysize*i] = false;
276 
277  int last_index = 0;
278  for (int i=0; i<m_xsize; i++){
279  if ( !((m_x[i] < xmin) || (m_x[i]>xmax))){
280  xx.push_back(m_x[i]);
281  if (m_edges_xx.size()>0)
282  edges_xx.push_back(m_edges_xx[i]);
283  last_index = i;
284  }
285  }
286  if (m_edges_xx.size()>0)
287  edges_xx.push_back(m_edges_xx[last_index+1]);
288 
289  last_index = 0;
290  for (int i=0; i<m_ysize; i++){
291  if (!((m_y[i] < ymin) || (m_y[i]>ymax))){
292  yy.push_back(m_y[i]);
293  if (m_edges_yy.size()>0)
294  edges_yy.push_back(m_edges_yy[i]);
295  last_index = i;
296  }
297  }
298  if (m_edges_yy.size()>0)
299  edges_yy.push_back(m_edges_yy[last_index+1]);
300 
301  vector<double> data, error;
302  vector<vector<double>> covariance;
303  vector<vector<double>> extra_info;
304 
305  int ndata_eff = 0;
306 
307  for (int i=0; i<m_ndata; i++)
308  if (mask[i])
309  ndata_eff ++;
310 
311  if (ndata_eff<1)
312  ErrorCBL("no elements left!", "cut", "Data2D_extra.cpp");
313 
314  data.resize(ndata_eff, 0);
315  error.resize(ndata_eff, 0);
316  covariance.resize(ndata_eff, vector<double>(ndata_eff, 0));
317  extra_info.resize(m_extra_info.size(), vector<double>(ndata_eff, 0));
318 
319  int index1 = 0;
320  for (int i=0; i<m_ndata; i++) {
321  if (mask[i]) {
322  data[index1] = m_data[i];
323  error[index1] = m_error[i];
324 
325  for (size_t j=0; j<m_extra_info.size(); j++)
326  extra_info[j][index1] = m_extra_info[j][i];
327 
328  int index2 = 0;
329  for (int j=0; j<m_ndata; j++) {
330  if (mask[j]) {
331  covariance[index1][index2] = m_covariance[i][j];
332  index2 ++;
333  }
334  }
335  index1 ++;
336  }
337  }
338 
339  shared_ptr<Data> dd = make_shared<Data2D_extra>(Data2D_extra(xx, yy, data, covariance, extra_info, edges_xx, edges_yy));
340 
341  return dd;
342 }
343 
The class Data2D_extra.
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class Data2D_extra.
Definition: Data2D_extra.h:52
virtual void read(const std::string input_file, const int skip_nlines=0, const std::vector< int > column={1, 2}, const std::vector< int > column_data={}, const std::vector< int > column_errors={}, const std::vector< int > column_edges={}) override
read the data
std::shared_ptr< Data > cut(const double xmin, const double xmax, const double ymin, const double ymax) const
cut the data
virtual void write(const std::string dir, const std::string file, const std::string header, const bool full, const int prec=4, const int ww=8, const int rank=0) const override
write the data
virtual void Print(const int precision=4) const override
print the data on screen
static const std::string defaultString
default std::string value
Definition: Constants.h:336
static const double defaultDouble
default double value
Definition: Constants.h:348
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
void Print(const T value, const int prec, const int ww, const std::string header="", const std::string end="\n", const bool use_coutCBL=true, std::ostream &stream=std::cout, const std::string colour=cbl::par::col_default)
function to print values with a proper homegenised format
Definition: Kernel.h:1142
void checkDim(const std::vector< T > vect, const int val, const std::string vector, bool equal=true)
check if the dimension of a std::vector is equal/lower than an input value
Definition: Kernel.h:1532
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
void unique_unsorted(std::vector< int > &vv)
erase all the equal elements of the input std::vector
Definition: Kernel.cpp:284
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747