Commit 3b6aa615 authored by Laurent Belcour's avatar Laurent Belcour

Adding a convertion tool between vertical segment objects.

Updating the parametrizations to be more complete.
Changing the data object to handle cosine correction.
parent b77123ed
...@@ -22,22 +22,24 @@ class data : public parametrized ...@@ -22,22 +22,24 @@ class data : public parametrized
virtual void load(const std::string& filename) = 0 ; virtual void load(const std::string& filename) = 0 ;
virtual void load(const std::string& filename, const arguments& args) = 0 ; virtual void load(const std::string& filename, const arguments& args) = 0 ;
// Save the data to a file // Save the data to a file
virtual void save(const std::string& filename) const virtual void save(const std::string& filename) const
{ {
std::ofstream file(filename.c_str(), std::ios_base::trunc); std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl; file << "#DIM " << _nX << " " << _nY << std::endl;
for(int i=0; i<size(); ++i) file << "PARAM_IN " << params::get_name(input_parametrization()) << std::endl;
{ file << "PARAM_OUT " << params::get_name(output_parametrization()) << std::endl;
vec x = this->get(i); for(int i=0; i<size(); ++i)
for(int j=0; j<_nX+_nY; ++j) {
{ vec x = this->get(i);
file << x[j] << "\t"; for(int j=0; j<_nX+_nY; ++j)
} {
file << std::endl; file << x[j] << "\t";
} }
file.close(); file << std::endl;
} }
file.close();
}
// Acces to data // Acces to data
virtual vec get(int i) const = 0 ; virtual vec get(int i) const = 0 ;
......
...@@ -430,6 +430,7 @@ bool compound_function::load(std::istream& in) ...@@ -430,6 +430,7 @@ bool compound_function::load(std::istream& in)
nb_good++; nb_good++;
} }
} }
std::cout << "<<DEBUG>> number of correct function loaded in compound: " << nb_good << " / " << fs.size() << std::endl;
return nb_good > 0; return nb_good > 0;
} }
......
...@@ -32,6 +32,7 @@ std::map<params::input, const param_info> create_map() ...@@ -32,6 +32,7 @@ std::map<params::input, const param_info> create_map()
_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_VH, param_info("RUSIN_VH", 3, "Vector representation of the Half angle only"))); _map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_VH, param_info("RUSIN_VH", 3, "Vector representation of the Half angle only")));
_map.insert(std::make_pair<params::input, const param_info>(params::SCHLICK_VK, param_info("SCHLICK_VH", 3, "Vector representation of the Back angle only"))); _map.insert(std::make_pair<params::input, const param_info>(params::SCHLICK_VK, param_info("SCHLICK_VH", 3, "Vector representation of the Back angle only")));
_map.insert(std::make_pair<params::input, const param_info>(params::ISOTROPIC_TL_TV_PROJ_DPHI, param_info("ISOTROPIC_TL_TV_PROJ_DPHI", 3, "Isoptropic projected phi parametrization."))); _map.insert(std::make_pair<params::input, const param_info>(params::ISOTROPIC_TL_TV_PROJ_DPHI, param_info("ISOTROPIC_TL_TV_PROJ_DPHI", 3, "Isoptropic projected phi parametrization.")));
_map.insert(std::make_pair<params::input, const param_info>(params::SCHLICK_TL_TK_PROJ_DPHI, param_info("SCHLICK_TL_TK_PROJ_DPHI", 3, "Isoptropic projected phi parametrization centered around the back vector.")));
/* 4D Params */ /* 4D Params */
_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_TH_PH_TD_PD, param_info("RUSIN_TH_PH_TD_PD", 4, "Complete Half angle parametrization"))); _map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_TH_PH_TD_PD, param_info("RUSIN_TH_PH_TD_PD", 4, "Complete Half angle parametrization")));
...@@ -173,6 +174,9 @@ void params::to_cartesian(const double* invec, params::input intype, ...@@ -173,6 +174,9 @@ void params::to_cartesian(const double* invec, params::input intype,
outvec[5] = cos(invec[0]); outvec[5] = cos(invec[0]);
} }
break; break;
case SCHLICK_TL_TK_PROJ_DPHI:
NOT_IMPLEMENTED();
break;
case RETRO_TL_TVL_PROJ_DPHI: case RETRO_TL_TVL_PROJ_DPHI:
{ {
const double theta = std::fabs(sqrt(invec[1]*invec[1] + invec[2]*invec[2]) - invec[0]); const double theta = std::fabs(sqrt(invec[1]*invec[1] + invec[2]*invec[2]) - invec[0]);
...@@ -354,6 +358,19 @@ void params::from_cartesian(const double* invec, params::input outtype, ...@@ -354,6 +358,19 @@ void params::from_cartesian(const double* invec, params::input outtype,
outvec[0] = theta_l; outvec[0] = theta_l;
outvec[1] = theta_v * cos(dphi); outvec[1] = theta_v * cos(dphi);
outvec[2] = theta_v * sin(dphi); outvec[2] = theta_v * sin(dphi);
}
break;
case SCHLICK_TL_TK_PROJ_DPHI:
{
const double vkx = invec[0]-invec[3];
const double vky = invec[1]-invec[4];
const double vkz = invec[2]+invec[5];
const double norm = sqrt(vkx*vkx + vky*vky + vkz*vkz);
const double theta_k = acos(vkz / norm);
const double dphi = atan2(invec[4], invec[3]) - atan2(vky/norm, vkx/norm);
outvec[0] = acos(invec[5]);
outvec[1] = theta_k * cos(dphi);
outvec[2] = theta_k * sin(dphi);
} }
break; break;
......
...@@ -52,9 +52,9 @@ class params ...@@ -52,9 +52,9 @@ class params
SCHLICK_TK_PK, /*!< Schlick's back vector parametrization */ SCHLICK_TK_PK, /*!< Schlick's back vector parametrization */
SCHLICK_VK, /*!< Schlick's back vector */ SCHLICK_VK, /*!< Schlick's back vector */
SCHLICK_TK_PROJ_DPHI, /*!< 2D Parametrization where the phi component is projected and SCHLICK_TL_TK_PROJ_DPHI,/*!< 3D Parametrization where the phi component is projected and
the parametrization is centered around the back direction. the parametrization is centered around the back direction.
\f$[x, y] = [\theta_K \cos(\phi_K), \theta_K \sin(\phi_K)]\f$*/ \f$[\theta_L, x, y] = [\theta_L, \theta_K \cos(\phi_K), \theta_K \sin(\phi_K)]\f$*/
COS_TK, /*!< Schlick's back vector dot product with the normal */ COS_TK, /*!< Schlick's back vector dot product with the normal */
...@@ -118,6 +118,15 @@ class params ...@@ -118,6 +118,15 @@ class params
//! \brief look for the string associated with a parametrization //! \brief look for the string associated with a parametrization
//! type. //! type.
static std::string get_name(const params::input param); static std::string get_name(const params::input param);
//! \brief look for the string associated with a parametrization
//! type.
//! \todo Finish this implementation. It requires another static
//! object.
static std::string get_name(const params::output param)
{
return std::string("UNKNOWN_OUTPUT");
}
//! \brief static function for input type convertion. This //! \brief static function for input type convertion. This
//! function allocate the resulting vector. //! function allocate the resulting vector.
......
...@@ -374,7 +374,7 @@ function* plugins_manager::get_function(const arguments& args) ...@@ -374,7 +374,7 @@ function* plugins_manager::get_function(const arguments& args)
} }
} }
//* /*
// Correction of the data by 1/cosine(theta_L) // Correction of the data by 1/cosine(theta_L)
if(args.is_defined("data-correct-cosine")) if(args.is_defined("data-correct-cosine"))
{ {
...@@ -382,7 +382,7 @@ function* plugins_manager::get_function(const arguments& args) ...@@ -382,7 +382,7 @@ function* plugins_manager::get_function(const arguments& args)
func = new product_function(cosine, dynamic_cast<nonlinear_function*>(func)); func = new product_function(cosine, dynamic_cast<nonlinear_function*>(func));
} }
// End of correction // End of correction
//*/ */
return func; return func;
} }
data* plugins_manager::get_data(const std::string& n) data* plugins_manager::get_data(const std::string& n)
......
...@@ -98,22 +98,29 @@ void vertical_segment::load(const std::string& filename, const arguments& args) ...@@ -98,22 +98,29 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
linestream >> v[i] ; linestream >> v[i] ;
/* // /*
// Correction of the data by 1/cosine(theta_L) // Correction of the data by 1/cosine(theta_L)
double factor = 1.0; double factor = 1.0;
if(args.is_defined("data-correct-cosine")) if(args.is_defined("data-correct-cosine"))
{ {
double cart[6]; double cart[6];
params::convert(&v[0], input_parametrization(), params::CARTESIAN, cart); params::convert(&v[0], input_parametrization(), params::CARTESIAN, cart);
factor = 1.0/cart[5]; if(cart[5] > 0.0 && cart[2] > 0.0)
{
factor = 1.0/cart[5]*cart[2];
}
else
{
continue;
}
} }
// End of correction // End of correction
*/ // */
for(int i=0; i<dimY(); ++i) for(int i=0; i<dimY(); ++i)
{ {
linestream >> v[dimX() + i]; linestream >> v[dimX() + i];
// v[dimX() + i] /= factor; v[dimX() + i] /= factor;
} }
// Check if the data containt a vertical segment around the mean // Check if the data containt a vertical segment around the mean
...@@ -181,6 +188,9 @@ void vertical_segment::load(const std::string& filename, const arguments& args) ...@@ -181,6 +188,9 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
} }
} }
} }
if(args.is_defined("data-correct-cosine"))
save("/tmp/data-corrected.dat");
std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ; std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ;
std::cout << "<<INFO>> data inside " << _min << " ... " << _max << std::endl ; std::cout << "<<INFO>> data inside " << _min << " ... " << _max << std::endl ;
......
...@@ -160,7 +160,7 @@ vec diffuse_function::parametersJacobian(const vec& x) const ...@@ -160,7 +160,7 @@ vec diffuse_function::parametersJacobian(const vec& x) const
void diffuse_function::bootstrap(const data* d, const arguments& args) void diffuse_function::bootstrap(const data* d, const arguments& args)
{ {/*
// Set the diffuse component // Set the diffuse component
if(params::is_cosine_weighted(d->output_parametrization())) if(params::is_cosine_weighted(d->output_parametrization()))
{ {
...@@ -185,12 +185,12 @@ void diffuse_function::bootstrap(const data* d, const arguments& args) ...@@ -185,12 +185,12 @@ void diffuse_function::bootstrap(const data* d, const arguments& args)
} }
} }
else else
*/
{ {
vec x0 = d->get(0);
for(int i=0; i<d->dimY(); ++i) for(int i=0; i<d->dimY(); ++i)
_kd[i] = x0[d->dimX() + i]; _kd[i] = std::numeric_limits<double>::max();
for(int i=1; i<d->size(); ++i) for(int i=0; i<d->size(); ++i)
{ {
vec xi = d->get(i); vec xi = d->get(i);
for(int j=0; j<d->dimY(); ++j) for(int j=0; j<d->dimY(); ++j)
......
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include <core/function.h> #include <core/function.h>
#include <core/fitter.h> #include <core/fitter.h>
#include <core/plugins_manager.h> #include <core/plugins_manager.h>
#include <core/vertical_segment.h>
#include <iostream> #include <iostream>
#include <vector> #include <vector>
...@@ -70,25 +71,61 @@ int main(int argc, char** argv) ...@@ -70,25 +71,61 @@ int main(int argc, char** argv)
// Import data // Import data
data* d_in = NULL ; data* d_in = NULL ;
d_in = plugins_manager::get_data(args["in-data"]) ; d_in = plugins_manager::get_data(args["in-data"]) ;
d_in->load(args["input"]); d_in->load(args["input"], args);
data* d_out = NULL; data* d_out = NULL;
d_out = plugins_manager::get_data(args["out-data"]) ; d_out = plugins_manager::get_data(args["out-data"]) ;
if(d_out->dimY() != d_in->dimY()) if(d_in == NULL && d_out == NULL)
{ {
std::cerr << "<<ERROR>> data types have incompatible output dimensions" << std::endl; std::cerr << "<<ERROR>> cannot import or export data" << std::endl ;
return 1;
} }
if(d_in != NULL && d_out != NULL) if(dynamic_cast<vertical_segment*>(d_out) != NULL)
{ {
params::input param = params::parse_input(args["param"]);
if(param == params::UNKNOWN_INPUT)
{
std::cerr << "<<ERROR>> unable to parse the parametrization" << std::endl;
return -1;
}
d_out->setParametrization(param);
d_out->setDimX(params::dimension(param));
d_out->setDimY(d_in->dimY());
std::cout << "<<INFO>> output DIM = " << d_out->dimX() << ", " << d_out->dimY() << std::endl;
vec temp(d_out->dimX() + d_out->dimY());
for(int i=0; i<d_in->size(); ++i)
{
// Copy the input vector
vec x = d_in->get(i);
params::convert(&x[0], d_in->parametrization(), d_out->parametrization(), &temp[0]);
for(int j=0; j<d_in->dimY(); ++j)
{
temp[d_out->dimX() + j] = x[d_in->dimX() + j];
}
d_out->set(temp);
}
}
else
{
if(d_out->dimY() != d_in->dimY())
{
std::cerr << "<<ERROR>> data types have incompatible output dimensions" << std::endl;
}
vec temp(d_in->dimX()); vec temp(d_in->dimX());
for(int i=0; i<d_out->size(); ++i) for(int i=0; i<d_out->size(); ++i)
{ {
// Copy the input vector // Copy the input vector
vec x = d_out->get(i); vec x = d_out->get(i);
params::convert(&x[0], d_out->parametrization(), d_in->parametrization(), &temp[0]); params::convert(&x[0], d_out->parametrization(), d_in->parametrization(), &temp[0]);
vec y = d_in->value(temp); vec y = d_in->value(temp);
for(int j=0; j<d_in->dimY(); ++j) for(int j=0; j<d_in->dimY(); ++j)
...@@ -98,14 +135,7 @@ int main(int argc, char** argv) ...@@ -98,14 +135,7 @@ int main(int argc, char** argv)
d_out->set(x); d_out->set(x);
} }
d_out->save(args["output"]);
return 0 ;
} }
else d_out->save(args["output"]);
{ return 0 ;
std::cerr << "<<ERROR>> cannot import or export data" << std::endl ;
return 1;
}
} }
...@@ -36,7 +36,7 @@ int main(int argc, char** argv) ...@@ -36,7 +36,7 @@ int main(int argc, char** argv)
if(args.is_defined("help")) { if(args.is_defined("help")) {
std::cout << "Usage: data2moments --input data.file --output gnuplot.file --data loader.so" << std::endl ; std::cout << "Usage: data2moments --input data.file --output gnuplot.file --data loader.so" << std::endl ;
std::cout << " -> input, output and data are mandatory parameters" << std::endl ; std::cout << "Compute the statistical moments, up to order 4, on the data file" << std::endl ;
return 0 ; return 0 ;
} }
...@@ -107,6 +107,7 @@ int main(int argc, char** argv) ...@@ -107,6 +107,7 @@ int main(int argc, char** argv)
vec k_yyyy(nY); vec k_yyyy(nY);
#ifndef OLD #ifndef OLD
// Number of elements per dimension // Number of elements per dimension
std::vector<int> dim; std::vector<int> dim;
if(args.is_defined("dim")) if(args.is_defined("dim"))
...@@ -156,7 +157,6 @@ int main(int argc, char** argv) ...@@ -156,7 +157,6 @@ int main(int argc, char** argv)
// For the global index i of the sample, compute the index in the various // For the global index i of the sample, compute the index in the various
// dimensions and then compute the position within the dimension. All the // dimensions and then compute the position within the dimension. All the
// samples are taken uniformly. // samples are taken uniformly.
// \todo Ad Metropolis integration
int indices[nX]; int indices[nX];
int global = i; int global = i;
for(int k=0; k<nX; ++k) for(int k=0; k<nX; ++k)
...@@ -228,7 +228,7 @@ int main(int argc, char** argv) ...@@ -228,7 +228,7 @@ int main(int argc, char** argv)
#ifdef NOT #ifdef NOT
k_xxxx = m_xxxx - 4*m_xxx*m_x - 3*m_xx*m_xx + 12*m_xx*m_x*m_x - 6*m_x*m_x*m_x*m_x; k_xxxx = m_xxxx - 4*m_xxx*m_x - 3*m_xx*m_xx + 12*m_xx*m_x*m_x - 6*m_x*m_x*m_x*m_x;
k_xxxy = m_xxxy - 3*m_xxy*m_x - m_xxx*m_y - 3*m_xx*m_xy + 6*m_xx*m_x*m_y + 6*m_xy*m_x*m_x - 6*m_x*m_x*m_x*m_y; k_xxxy = m_xxxy - 3*m_xxy*m_x - m_xxx*m_y - 3*m_xx*m_xy + 6*m_xx*m_x*m_y + 6*m_xy*m_x*m_x - 6*m_x*m_x*m_x*m_y;
k_xxyy = m_xxyy - 2*m_xxy*m_y - m_xyy*m_x - 2*m_xy*m_xy - m_xx*m_yy + 8*m_xy*m_x*m_y + 2*m_xx*m_y*m_y + 2*m_yy*m_x*m_x - 6*m_x*m_x*m_y*m_y; k_xxyy = m_xxyy - 2*m_xxy*m_y - 2*m_xyy*m_x - 2*m_xy*m_xy - m_xx*m_yy + 8*m_xy*m_x*m_y + 2*m_xx*m_y*m_y + 2*m_yy*m_x*m_x - 6*m_x*m_x*m_y*m_y;
k_xyyy = m_xyyy - 3*m_xyy*m_y - m_yyy*m_x - 3*m_yy*m_xy + 6*m_yy*m_x*m_y + 6*m_xy*m_y*m_y - 6*m_x*m_y*m_y*m_y; k_xyyy = m_xyyy - 3*m_xyy*m_y - m_yyy*m_x - 3*m_yy*m_xy + 6*m_yy*m_x*m_y + 6*m_xy*m_y*m_y - 6*m_x*m_y*m_y*m_y;
k_yyyy = m_yyyy - 4*m_yyy*m_y - 3*m_yy*m_yy + 12*m_yy*m_y*m_y - 6*m_y*m_y*m_y*m_y; k_yyyy = m_yyyy - 4*m_yyy*m_y - 3*m_yy*m_yy + 12*m_yy*m_y*m_y - 6*m_y*m_y*m_y*m_y;
#endif #endif
......
Markdown is supported
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