CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Data1D_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 "Data1D_extra.h"
35 
36 using namespace std;
37 
38 using namespace cbl;
39 using namespace data;
40 
41 
42 // ======================================================================================
43 
44 
45 void cbl::data::Data1D_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 vector<int> column_edges)
46 {
47  if (column_data.size()!=column_errors.size()) ErrorCBL("the sizes of column_data and column_errors must be equal!", "read", "Data1D_extra.cpp");
48 
49  // default column of input data values
50  int column_data_default = column[0]+1;
51 
52  // default column of input error values
53  int column_error_default = column[0]+2;
54 
55  // default first column of input edge values
56  int column_edges_default = column[0]+3;
57 
58  ifstream fin(input_file.c_str()); checkIO(fin, input_file);
59  string line;
60 
61  // loop on the data (and, possibly, error) vectors to be read; if
62  // more than one data vectors are read, the data (and, possibly,
63  // error) vectors will be added one after the other, in order to
64  // construct a single data object
65  const int cl_max = (column_data.size()<2) ? 1 : column_data.size();
66  for (int cl=0; cl<cl_max; ++cl) {
67 
68  // skip the first nlines (in case of header lines)
69  if (skip_nlines>0)
70  for (int i=0; i<skip_nlines; ++i)
71  getline(fin, line);
72 
73  // read the file lines
74  double last_bin_edge = par::defaultDouble;
75  while (getline(fin, line)) {
76 
77  stringstream ss(line);
78  vector<double> num; double NUM;
79  while (ss>>NUM) num.emplace_back(NUM);
80 
81  // if column[0] is larger than 0, the column of x coordinates is
82  // the one specified in input
83  const int ind_x = max(0, column[0]-1);
84  checkDim(num, ind_x+1, "num", false);
85  m_x.emplace_back(num[ind_x]);
86 
87  // if the size of the column_data vector is 0, the columns of
88  // data values are the ones specified in input
89  const int ind_data = (column_data.size()==0) ? column_data_default : column_data[cl]-1;
90  checkDim(num, ind_data+1, "num", false);
91  m_data.emplace_back(num[ind_data]);
92 
93  // if the size of the column_errors vector is 0, the columns of
94  // error values are the ones specified in input; if the error
95  // column is not present, the errors will be set to 1, by
96  // default
97  const size_t ind_error = (column_errors.size()==0) ? column_error_default : column_errors[cl]-1;
98  if (num.size()<ind_error+1 && column_errors.size()>1)
99  WarningMsgCBL("the errors cannot be retrieved from the provided input file, and will be set to 1", "read", "Data1D_extra.cpp");
100  m_error.emplace_back((num.size()<ind_error+1) ? 1 : num[ind_error]);
101 
102  // if the size of the column_edges vector is 0, the columns of
103  // bin edges values are the ones specified in input
104  const size_t ind_edges = (column_edges.size()==0) ? column_edges_default-1 : column_edges[cl]-1;
105  if (num.size()>ind_edges+1) {
106  m_edges_xx.emplace_back(num[ind_edges]);
107  last_bin_edge = num[ind_edges+1];
108  }
109 
110  // get the extra info
111  int ind = 0;
112  for (size_t ex=num.size()-m_extra_info.size(); ex<num.size(); ++ex)
113  m_extra_info[ind++].emplace_back(num[ex]);
114  }
115 
116  // add the last bin edge to m_edges_xx
117  if (last_bin_edge>par::defaultDouble)
118  m_edges_xx.emplace_back(last_bin_edge);
119 
120  // go back at the beginning of the file
121  fin.clear(); fin.seekg(ios::beg);
122 
123  column_data_default += 2;
124  column_error_default += 2;
125  }
126 
127  fin.clear(); fin.close();
128 
129 
130  // set the data size and diagonal covariance
131 
132  m_ndata = m_data.size();
133 
134  m_covariance.resize(m_ndata, vector<double>(m_ndata, 0));
135  for (int i=0; i<m_ndata; i++)
136  m_covariance[i][i] = pow(m_error[i], 2);
137 
138 }
139 
140 
141 // ======================================================================================
142 
143 
144 void cbl::data::Data1D_extra::Print (const int precision) const
145 {
146  for (size_t i=0; i<m_x.size(); ++i) {
147  coutCBL << setprecision(precision) << setw(8) << right << m_x[i]
148  << " " << setprecision(precision) << setw(8) << right << m_data[i]
149  << " " << setprecision(precision) << setw(8) << right << m_error[i];
150 
151  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
152  coutCBL << " " << setprecision(precision) << setw(8) << m_extra_info[ex][i];
153  cout << endl;
154  }
155 }
156 
157 
158 // ======================================================================================
159 
160 
161 void cbl::data::Data1D_extra::write (const string dir, const string file, const string header, const int prec, const int ww, const int rank) const
162 {
163  (void)rank;
164 
165  string file_out = dir+file;
166  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
167 
168  fout << "### "<< header <<" ###" << endl;
169 
170  for (size_t i=0; i<m_x.size(); ++i) {
171  cbl::Print(m_x[i], prec, ww, "", " ", false, fout);
172  cbl::Print(m_data[i], prec, ww, "", " ", false, fout);
173  cbl::Print(m_error[i], prec, ww, "", " ", false, fout);
174  for (size_t ex=0; ex<m_extra_info.size(); ++ex)
175  cbl::Print(m_extra_info[ex][i], prec, ww, "", " ", false, fout);
176  fout << endl;
177  }
178 
179  fout.close(); cout << endl ; coutCBL << "I wrote the file: " << file_out << endl;
180 }
181 
182 
183 // ======================================================================================
184 
185 
186 std::shared_ptr<Data> cbl::data::Data1D_extra::cut (const std::vector<bool> mask) const
187 {
188  checkDim(mask, m_ndata, "mask");
189 
190  vector<double> xx, edges_xx;
191  vector<double> data, error;
192  vector<vector<double>> covariance;
193  vector<vector<double>> extra_info;
194 
195  int ndata_eff = 0;
196 
197  for (int i=0; i<m_ndata; i++)
198  if (mask[i])
199  ndata_eff +=1;
200 
201  if (ndata_eff <1)
202  ErrorCBL("no elements left!", "cut", "Data1D_extra.cpp");
203 
204  xx.resize(ndata_eff, 0);
205  if (m_edges_xx.size() > 0)
206  edges_xx.resize(ndata_eff, 0);
207  data.resize(ndata_eff, 0);
208  error.resize(ndata_eff, 0);
209  covariance.resize(ndata_eff, vector<double>(ndata_eff, 0));
210  extra_info.resize(m_extra_info.size(), vector<double>(ndata_eff, 0));
211 
212  int index1 = 0;
213  int last_index1 = 0;
214  for (int i=0; i<m_ndata; i++) {
215  if (mask[i]) {
216  xx[index1] = m_x[i];
217  if (m_edges_xx.size() > 0)
218  edges_xx[index1] = m_edges_xx[i];
219  data[index1] = m_data[i];
220  error[index1] = m_error[i];
221 
222  for (size_t j=0; j<m_extra_info.size(); j++)
223  extra_info[j][index1] = m_extra_info[j][i];
224 
225  int index2 = 0;
226  for (int j=0; j<m_ndata; j++) {
227  if (mask[j]) {
228  covariance[index1][index2] = m_covariance[i][j];
229  index2++;
230  }
231  }
232  index1++;
233  last_index1 = i;
234  }
235  }
236  if (m_edges_xx.size() > 0)
237  edges_xx.push_back(m_edges_xx[last_index1+1]);
238 
239  shared_ptr<Data> dd = make_shared<Data1D_extra>(Data1D_extra(xx, data, covariance, extra_info, edges_xx));
240 
241  return dd;
242 }
243 
244 
245 // ======================================================================================
246 
247 
248 shared_ptr<Data> cbl::data::Data1D_extra::cut (const double xmin, const double xmax) const
249 {
250  vector<bool> mask(m_ndata, true);
251  for (int i=0; i<m_ndata; i++)
252  if ((m_x[i] < xmin) || (m_x[i]>xmax))
253  mask[i] = false;
254 
255  return cut(mask);
256 }
The class Data1D_extra.
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class Data1D_extra.
Definition: Data1D_extra.h:52
std::shared_ptr< Data > cut(const double xmin, const double xmax) const override
cut the data
virtual void Print(const int precision=4) const override
print the data on screen
virtual void read(const std::string input_file, const int skip_nlines=0, const std::vector< int > column={1}, const std::vector< int > column_data={}, const std::vector< int > column_errors={}, const std::vector< int > column_edges={}) override
read the data
virtual void write(const std::string dir, const std::string file, const std::string header, const int prec=4, const int ww=8, const int rank=0) const override
write the data
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 WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747