CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Model1D.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2014 by Federico Marulli and Alfonso Veropalumbo *
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 "Model1D.h"
35 
36 using namespace std;
37 
38 using namespace cbl;
39 
40 
41 // ======================================================================================
42 
43 
45 {
46  m_function = [this, function](vector<vector<double>> xx, shared_ptr<void> inputs, vector<double> &parameters)
47  {
48  vector<vector<double>> res(1, function(xx[0], inputs, parameters));
49  return res;
50  };
51 }
52 
53 
54 // ======================================================================================
55 
56 
57 void cbl::statistics::Model1D::set_function (const std::vector<double> (*function)(const std::vector<double> xx, std::vector<double> &val))
58 {
59  m_function = [this, function](vector<vector<double>> xx, shared_ptr<void> inputs, vector<double> &parameters)
60  {
61  (void)inputs;
62  vector<vector<double>> res(1, function(xx[0], parameters));
63  return res;
64  };
65 }
66 
67 
68 // ======================================================================================
69 
70 
71 void cbl::statistics::Model1D::stats_from_chains (const vector<double> xx, vector<double> &median_model, vector<double> &low_model, vector<double> &up_model, const int start, const int thin)
72 {
73  vector<vector<double>> _median, _low, _up;
74 
75  Model::stats_from_chains({1, xx}, _median, _low, _up, start, thin);
76 
77  median_model = _median[0];
78  low_model = _low[0];
79  up_model = _up[0];
80 }
81 
82 
83 // ======================================================================================
84 
85 
86 void cbl::statistics::Model1D::write (const string output_dir, const string output_file, const vector<double> xx, const vector<double> parameters)
87 {
88  vector<double> pp = parameters;
89  vector<double> xx_unique = different_elements(xx);
90 
91  if (xx.size() % xx_unique.size()!=0)
92  ErrorCBL("model.size() is not a multiple of xx.size()!", "write", "Model1D.cpp");
93 
94  const int nmodels = xx.size()/xx_unique.size();
95  vector<double> model = this->operator()(xx, pp);
96 
97  const string mkdir = "mkdir -p "+output_dir;
98  if (system(mkdir.c_str())) {}
99 
100  const string file = output_dir+output_file;
101 
102  ofstream fout(file.c_str()); checkIO(fout, file);
103 
104  fout << "### [1] x";
105  for (int nn=0; nn<nmodels; nn++) fout << " # [" << nn+2 << "] y_" << nn << "(x)";
106  fout << " ###" << endl;
107 
108  for (size_t i=0; i<xx_unique.size(); ++i) {
109  fout << setprecision(5) << setw(10) << right << xx_unique[i];
110  for (int nn=0; nn<nmodels; nn++)
111  fout << " " << setprecision(5) << setw(10) << right << model[i+nn*xx_unique.size()];
112  fout << endl;
113  }
114 
115  fout.clear(); fout.close(); coutCBL << "I wrote the file: " << file << endl << endl;
116 }
117 
118 
119 // ======================================================================================
120 
121 
122 void cbl::statistics::Model1D::write_at_bestfit (const string output_dir, const string output_file, const vector<double> xx)
123 {
124  vector<double> bf = m_parameters->bestfit_value();
125  write(output_dir, output_file, xx, bf);
126 }
127 
128 
129 // ======================================================================================
130 
131 
132 void cbl::statistics::Model1D::write_from_chains (const string output_dir, const string output_file, const vector<double> xx, const int start, const int thin)
133 {
134  vector<double> xx_unique = different_elements(xx);
135 
136  if (xx.size() % xx_unique.size()!=0)
137  ErrorCBL("model.size() is not a multiple of xx.size()!", "write_from_chains", "Model1D.cpp");
138 
139  const int nmodels = xx.size()/xx_unique.size();
140 
141  vector<double> median_model, low_model, up_model;
142  stats_from_chains(xx, median_model, low_model, up_model, start, thin);
143 
144  const string mkdir = "mkdir -p "+output_dir;
145  if (system(mkdir.c_str())) {}
146 
147  const string file = output_dir+output_file;
148 
149  ofstream fout(file.c_str()); checkIO(fout, file);
150 
151  fout << "### [1] x";
152  for (int nn=0; nn<nmodels; nn++) fout << " # [" << nn+2 << "] median y(x) # [" << nn+3 << "] 16% percentile y(x)# [" << nn+4 << "] 84% percentile y(x)";
153  fout << " ###" << endl;
154 
155  for (size_t i=0; i<xx_unique.size(); i++) {
156  fout << setprecision(5) << setw(10) << right << xx_unique[i] << " ";
157  for (int nn=0; nn<nmodels; nn++)
158  fout << setprecision(5) << setw(10) << right << median_model[i+nn*xx_unique.size()] << " "
159  << setprecision(5) << setw(10) << right << low_model[i+nn*xx_unique.size()] << " "
160  << setprecision(5) << setw(10) << right << up_model[i+nn*xx_unique.size()] << " ";
161  fout << endl;
162  }
163 
164  fout.clear(); fout.close(); coutCBL << "I wrote the file: " << file << endl;
165 }
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class Model1D.
void write(const std::string output_dir, const std::string output_file, const std::vector< double > xx, const std::vector< double > parameters) override
write the model at xx for given parameters
Definition: Model1D.cpp:86
void write_from_chains(const std::string output_dir, const std::string output_file, const std::vector< double > xx, const int start=0, const int thin=1) override
write the model at xx computing 16th, 50th and 84th percentiles from the chains.
Definition: Model1D.cpp:132
void write_at_bestfit(const std::string output_dir, const std::string output_file, const std::vector< double > xx) override
write the model at xx with best-fit parameters obtained from likelihood maximization
Definition: Model1D.cpp:122
void set_function(const model_function_1D function)
set the model function
Definition: Model1D.cpp:44
void stats_from_chains(const std::vector< double > xx, std::vector< double > &median_model, std::vector< double > &low_model, std::vector< double > &up_model, const int start=0, const int thin=1) override
compute the median and percentiles of the model from MCMC chains
Definition: Model1D.cpp:71
std::function< std::vector< double >std::vector< double >, std::shared_ptr< void >, std::vector< double > &)> model_function_1D
1D function: the inputs are a vector of values at which the function is computed, a pointer to a set ...
Definition: Model.h:58
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
std::vector< T > different_elements(const std::vector< T > vect_input)
get the unique elements of a std::vector
Definition: Kernel.h:1349
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