Commit f8843e56 authored by Ludovic Courtès's avatar Ludovic Courtès

core: vertical_segment: Store data in an 'Eigen::MatrixXd'.

parent d172fa6d
...@@ -36,6 +36,8 @@ ...@@ -36,6 +36,8 @@
#include <Eigen/Core> #include <Eigen/Core>
typedef Eigen::VectorXd vec; typedef Eigen::VectorXd vec;
typedef Eigen::Ref<vec> vecref;
typedef Eigen::Ref<const vec> const_vecref;
// Convenience functions. // Convenience functions.
static inline double norm(const vec& v) static inline double norm(const vec& v)
......
...@@ -21,31 +21,43 @@ ...@@ -21,31 +21,43 @@
# include <endian.h> # include <endian.h>
#endif #endif
//using namespace alta; using namespace Eigen;
static bool cosine_correction(vec& v, unsigned int dimX, unsigned int dimY, // A deleter for arrays, to work around the lack of array support in C++11's
// 'shared_ptr':
// <http://stackoverflow.com/questions/8947579/why-isnt-there-a-stdshared-ptrt-specialisation#8947700>.
static void delete_array(double *thing)
{
delete[] thing;
}
static bool cosine_correction(vecref v, unsigned int dimX, unsigned int dimY,
alta::params::input in_param) alta::params::input in_param)
{ {
using namespace alta; using namespace alta;
assert(v.size() == dimX + dimY);
// 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;
double cart[6]; double cart[6];
params::convert(&v[0], in_param, params::CARTESIAN, cart); params::convert(&v(0), in_param, params::CARTESIAN, cart);
if(cart[5] > 0.0 && cart[2] > 0.0) if(cart[5] > 0.0 && cart[2] > 0.0)
{ {
factor = 1.0/cart[5]*cart[2]; factor = 1.0/cart[5]*cart[2];
for(unsigned int i = 0; i < dimY; ++i) for(unsigned int i = 0; i < dimY; ++i)
{ {
v[i + dimX] /= factor; std::cout << "i = " << i << " dimx = " << dimX
<< " + = " << i + dimX << " | " << v.size()
<< std::endl;
v(i + dimX) /= factor;
} }
return true; return true;
} }
else return false; else return false;
} }
typedef Eigen::Ref<const vec> vecref;
// Return true if V is between MIN and MAX. // Return true if V is between MIN and MAX.
static bool within_bounds(const vecref v, static bool within_bounds(const vecref v,
const vec& min, const vec& max) const vec& min, const vec& max)
...@@ -64,11 +76,13 @@ enum ci_kind ...@@ -64,11 +76,13 @@ enum ci_kind
// Read a confidence interval on the output parameters from INPUT into V. // Read a confidence interval on the output parameters from INPUT into V.
static void read_confidence_interval(std::istream& input, static void read_confidence_interval(std::istream& input,
vec& v, ci_kind kind, vecref v, ci_kind kind,
unsigned int dimX, unsigned int dimX,
unsigned int dimY, unsigned int dimY,
const alta::arguments& args) const alta::arguments& args)
{ {
assert(v.size() == dimX + 3 * dimY);
for(int i=0; i<dimY; ++i) for(int i=0; i<dimY; ++i)
{ {
double min_dt = 0.0, max_dt = 0.0; double min_dt = 0.0, max_dt = 0.0;
...@@ -77,8 +91,8 @@ static void read_confidence_interval(std::istream& input, ...@@ -77,8 +91,8 @@ static void read_confidence_interval(std::istream& input,
{ {
input >> min_dt ; input >> min_dt ;
input >> max_dt ; input >> max_dt ;
min_dt = min_dt-v[dimX+i]; min_dt = min_dt-v(dimX + i);
max_dt = max_dt-v[dimX+i]; max_dt = max_dt-v(dimX + i);
} }
else if(i == 0 && kind == SYMMETRICAL_CONFIDENCE_INTERVAL) else if(i == 0 && kind == SYMMETRICAL_CONFIDENCE_INTERVAL)
{ {
...@@ -97,18 +111,18 @@ static void read_confidence_interval(std::istream& input, ...@@ -97,18 +111,18 @@ static void read_confidence_interval(std::istream& input,
if(args.is_defined("dt-relative")) if(args.is_defined("dt-relative"))
{ {
v[dimX + dimY+i] = v[dimX + i] * (1.0 + min_dt) ; v(dimX + dimY+i) = v(dimX + i) * (1.0 + min_dt) ;
v[dimX + 2*dimY+i] = v[dimX + i] * (1.0 + max_dt) ; v(dimX + 2*dimY+i) = v(dimX + i) * (1.0 + max_dt) ;
} }
else if(args.is_defined("dt-max")) else if(args.is_defined("dt-max"))
{ {
v[dimX + dimY+i] = v[dimX + i] + std::max(v[dimX + i] * min_dt, min_dt); v(dimX + dimY+i) = v(dimX + i) + std::max(v(dimX + i) * min_dt, min_dt);
v[dimX + 2*dimY+i] = v[dimX + i] + std::max(v[dimX + i] * max_dt, max_dt); v(dimX + 2*dimY+i) = v(dimX + i) + std::max(v(dimX + i) * max_dt, max_dt);
} }
else else
{ {
v[dimX + dimY+i] = v[dimX + i] + min_dt ; v(dimX + dimY+i) = v(dimX + i) + min_dt ;
v[dimX + 2*dimY+i] = v[dimX + i] + max_dt ; v(dimX + 2*dimY+i) = v(dimX + i) + max_dt ;
} }
// You can enforce the vertical segment to stay in the positive // You can enforce the vertical segment to stay in the positive
...@@ -116,9 +130,9 @@ static void read_confidence_interval(std::istream& input, ...@@ -116,9 +130,9 @@ static void read_confidence_interval(std::istream& input,
// that the data point is also clamped to zero if negative. // that the data point is also clamped to zero if negative.
if(args.is_defined("dt-positive")) if(args.is_defined("dt-positive"))
{ {
v[dimX + i] = std::max(v[dimX + i], 0.0); v(dimX + i) = std::max(v(dimX + i), 0.0);
v[dimX + dimY+i] = std::max(v[dimX + dimY+i], 0.0); v(dimX + dimY+i) = std::max(v(dimX + dimY+i), 0.0);
v[dimX + 2*dimY+i] = std::max(v[dimX + 2*dimY+i], 0.0); v(dimX + 2*dimY+i) = std::max(v(dimX + 2*dimY+i), 0.0);
} }
#ifdef DEBUG #ifdef DEBUG
...@@ -147,6 +161,7 @@ alta::data* alta::load_data_from_text(std::istream& input, ...@@ -147,6 +161,7 @@ alta::data* alta::load_data_from_text(std::istream& input,
params::output out_param = params::parse_output(header.get_string("PARAM_OUT", "UNKNOWN_OUTPUT")); params::output out_param = params::parse_output(header.get_string("PARAM_OUT", "UNKNOWN_OUTPUT"));
std::pair<int, int> dim = header.get_pair<int>("DIM"); std::pair<int, int> dim = header.get_pair<int>("DIM");
size_t row_count = dim.first + 3 * dim.second;
min = args.get_vec("min", dim.first, -std::numeric_limits<float>::max()) ; min = args.get_vec("min", dim.first, -std::numeric_limits<float>::max()) ;
max = args.get_vec("max", dim.first, std::numeric_limits<float>::max()) ; max = args.get_vec("max", dim.first, std::numeric_limits<float>::max()) ;
...@@ -165,8 +180,7 @@ alta::data* alta::load_data_from_text(std::istream& input, ...@@ -165,8 +180,7 @@ alta::data* alta::load_data_from_text(std::istream& input,
vs_value == 2 ? ASYMMETRICAL_CONFIDENCE_INTERVAL vs_value == 2 ? ASYMMETRICAL_CONFIDENCE_INTERVAL
: (vs_value == 1 ? SYMMETRICAL_CONFIDENCE_INTERVAL : (vs_value == 1 ? SYMMETRICAL_CONFIDENCE_INTERVAL
: NO_CONFIDENCE_INTERVAL); : NO_CONFIDENCE_INTERVAL);
std::vector<double> content;
std::vector<vec> content;
// Now read the body. // Now read the body.
while(input.good()) while(input.good())
...@@ -182,14 +196,24 @@ alta::data* alta::load_data_from_text(std::istream& input, ...@@ -182,14 +196,24 @@ alta::data* alta::load_data_from_text(std::istream& input,
} }
else else
{ {
auto start = content.size();
// Read the data point x and y coordinates // Read the data point x and y coordinates
vec v = vec::Zero(dim.first + 3*dim.second) ; for(int i=0; i < dim.first + dim.second; ++i)
for(int i=0; i<dim.first+dim.second; ++i)
{ {
linestream >> v[i] ; double item;
linestream >> item;
content.push_back(item);
} }
// Save room for the confidence interval.
for (int i = 0; i < 2 * dim.second; i++)
content.push_back(0.);
// If data is not in the interval of fit // If data is not in the interval of fit
// std::cout << " start = " << start << " rows = " << row_count
// << " size = " << content.size() << "\n";
Map<VectorXd> v(&content[start], row_count);
if (!(within_bounds(v.segment(0, dim.first), min, max) if (!(within_bounds(v.segment(0, dim.first), min, max)
&& within_bounds(v.segment(dim.first, dim.second), && within_bounds(v.segment(dim.first, dim.second),
ymin, ymax))) ymin, ymax)))
...@@ -197,15 +221,14 @@ alta::data* alta::load_data_from_text(std::istream& input, ...@@ -197,15 +221,14 @@ alta::data* alta::load_data_from_text(std::istream& input,
if(args.is_defined("data-correct-cosine")) if(args.is_defined("data-correct-cosine"))
{ {
if (!cosine_correction(v, dim.first, dim.second, in_param)) if (!cosine_correction(v.segment(0, dim.first + dim.second),
dim.first, dim.second, in_param))
continue; continue;
} }
// Read the confidence interval data if available. // Read the confidence interval data if available.
read_confidence_interval(linestream, v, kind, read_confidence_interval(linestream, v, kind,
dim.first, dim.second, args); dim.first, dim.second, args);
content.push_back(std::move(v));
} }
} }
...@@ -215,8 +238,15 @@ alta::data* alta::load_data_from_text(std::istream& input, ...@@ -215,8 +238,15 @@ alta::data* alta::load_data_from_text(std::istream& input,
<< " -> R^" << dim.second << std::endl ; << " -> R^" << dim.second << std::endl ;
std::cout << "<<INFO>> " << content.size() << " elements to fit" << std::endl ; std::cout << "<<INFO>> " << content.size() << " elements to fit" << std::endl ;
double *raw_content = new double[content.size()];
memcpy(raw_content, content.data(), content.size() * sizeof(double));
parameters param(dim.first, dim.second, in_param, out_param); parameters param(dim.first, dim.second, in_param, out_param);
data* result = new vertical_segment(param, std::move(content)); size_t element_count = content.size() / row_count;
data* result = new vertical_segment(param, element_count,
std::shared_ptr<double>(raw_content,
delete_array));
if(args.is_defined("data-correct-cosine")) if(args.is_defined("data-correct-cosine"))
result->save("/tmp/data-corrected.dat"); result->save("/tmp/data-corrected.dat");
...@@ -276,6 +306,7 @@ void alta::save_data_as_binary(std::ostream &out, const alta::data& data) ...@@ -276,6 +306,7 @@ void alta::save_data_as_binary(std::ostream &out, const alta::data& data)
for(int i=0; i < data.size(); ++i) for(int i=0; i < data.size(); ++i)
{ {
// FIXME: Currently we are not saving the confidence interval here.
vec sample = data.get(i); vec sample = data.get(i);
const double *numbers = sample.data(); const double *numbers = sample.data();
...@@ -311,26 +342,33 @@ alta::data* alta::load_data_from_binary(std::istream& in, const alta::arguments& ...@@ -311,26 +342,33 @@ alta::data* alta::load_data_from_binary(std::istream& in, const alta::arguments&
} }
// TODO: Arrange to use mmap and make it zero-alloc and zero-copy. // TODO: Arrange to use mmap and make it zero-alloc and zero-copy.
std::vector<vec> content(sample_count); size_t rows = dim.first + 3 * dim.second;
for (int i = 0; i < sample_count; i++) double *content = new double[sample_count * rows];
for (std::streamsize i = 0; i < sample_count; i++)
{ {
vec row = vec::Zero(dim.first + dim.second); double* start = &content[i * rows];
std::streamsize expected = row.size() * sizeof(double); std::streamsize byte_count =
(dim.first + dim.second) * sizeof(*content);
for (std::streamsize total = 0; for (std::streamsize total = 0;
total < expected && !in.eof(); total < byte_count && !in.eof();
total += in.gcount()) total += in.gcount())
{ {
char* ptr = (char*)row.data(); in.read((char *) start + total, byte_count - total);
in.read(ptr + total, expected - total); }
}
// XXX: Fill in the confidence interval values with zeros since we
content[i] = row; // don't have those values.
content[i * rows + dim.first + dim.second] = 0.;
content[i * rows + dim.first + dim.second + 1] = 0.;
} }
parameters param(dim.first, dim.second, parameters param(dim.first, dim.second,
params::parse_input(header["PARAM_IN"]), params::parse_input(header["PARAM_IN"]),
params::parse_output(header["PARAM_OUT"])); params::parse_output(header["PARAM_OUT"]));
return new alta::vertical_segment(param, std::move(content)); return new alta::vertical_segment(param, sample_count,
std::shared_ptr<double>(content,
delete_array));
} }
...@@ -18,131 +18,138 @@ ...@@ -18,131 +18,138 @@
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <cassert> #include <cassert>
#include <cstring>
using namespace alta; using namespace alta;
using namespace Eigen;
//#define RELATIVE_ERROR //#define RELATIVE_ERROR
// Note: For some reason returning an rvalue here doesn't work, so let's hope // Note: For some reason returning an rvalue here doesn't work, so let's hope
// RVO comes into play. // RVO comes into play.
static vec data_min(unsigned int size, const std::vector<vec>& data) static vec data_min(Ref<MatrixXd> data)
{ {
vec min(size); vec min(data.rows());
for(int k = 0; k < size; ++k) for(int k = 0; k < data.rows(); ++k)
{ {
min[k] = std::numeric_limits<double>::max(); min[k] = std::numeric_limits<double>::max();
} }
for (auto&& item: data) for(int i = 0; i < data.cols(); i++)
{ {
for (int k = 0; k < size; ++k) auto col = data.col(i);
for (int k = 0; k < data.rows(); ++k)
{ {
min[k] = std::min(min[k], item[k]); min[k] = std::min(min[k], col[k]);
} }
} }
return min; return min;
} }
static vec data_max(unsigned int size, const std::vector<vec>& data) static vec data_max(Ref<MatrixXd> data)
{ {
vec max(size); vec max(data.rows());
for(int k = 0; k < size; ++k) for(int k = 0; k < data.rows(); ++k)
{ {
max[k] = -std::numeric_limits<double>::max(); max[k] = -std::numeric_limits<double>::max();
} }
for (auto&& item: data) for(int i = 0; i < data.cols(); i++)
{ {
for (int k = 0; k < size; ++k) auto col = data.col(i);
for (int k = 0; k < data.rows(); ++k)
{ {
max[k] = std::max(max[k], item[k]); max[k] = std::max(max[k], col[k]);
} }
} }
return max; return max;
} }
// Return a matrix view of PTR showing only the 'dimX' first rows of that
// matrix.
static Ref<MatrixXd>
x_view(double *ptr, size_t cols, const parameters& params)
{
return Map<MatrixXd, 0, OuterStride<> >(ptr, params.dimX(), cols,
OuterStride<>(params.dimX()
+ 3 * params.dimY()));
}
vertical_segment::vertical_segment(const parameters& params, vertical_segment::vertical_segment(const parameters& params,
std::vector<vec>&& input_data) size_t size,
: data(params, input_data.size(), std::shared_ptr<double> input_data)
data_min(params.dimX(), input_data), : data(params, size,
data_max(params.dimX(), input_data)), data_min(x_view(input_data.get(), size, params)),
data_max(x_view(input_data.get(), size, params))),
_data(input_data), _is_absolute(true), _dt(0.1) _data(input_data), _is_absolute(true), _dt(0.1)
{ {
} }
vertical_segment::vertical_segment(const parameters& params, unsigned int size): vertical_segment::vertical_segment(const parameters& params, unsigned int rows):
data(params, size), _is_absolute(true), _dt(0.1) data(params, rows), _is_absolute(true), _dt(0.1)
{ {
initializeToZero(size); _data =
std::shared_ptr<double>(new double[rows * column_number()]);
memset(_data.get(), '\0', sizeof(double) * rows * column_number());
} }
void
vertical_segment::initializeToZero( unsigned int number_of_data_elements )
{
_data.clear();
unsigned int const size_of_one_element = _parameters.dimX() + 3*_parameters.dimY();
for( unsigned int i=0; i < number_of_data_elements; i++)
{
vec initial_element = vec::Zero( size_of_one_element );
_data.push_back( initial_element );
}
}
void vertical_segment::get(int i, vec& x, vec& yl, vec& yu) const void vertical_segment::get(int i, vec& x, vec& yl, vec& yu) const
{ {
auto matrix = matrix_view();
#ifdef DEBUG #ifdef DEBUG
assert(i >= 0 && i < _data.size()); assert(i >= 0 && i < matrix.size());
#endif #endif
x.resize(_parameters.dimX()); yl.resize(_parameters.dimY()) ; yu.resize(_parameters.dimY()) ; x.resize(_parameters.dimX()); yl.resize(_parameters.dimY()) ; yu.resize(_parameters.dimY()) ;
for(int j=0; j<_parameters.dimX(); ++j) for(int j=0; j<_parameters.dimX(); ++j)
{ {
x[j] = _data[i][j]; x[j] = matrix(i, j);
} }
for(int j=0; j<_parameters.dimY(); ++j) for(int j=0; j<_parameters.dimY(); ++j)
{ {
yl[j] = _data[i][_parameters.dimX() + 1*_parameters.dimY() + j]; yl[j] = matrix(i, _parameters.dimX() + 1*_parameters.dimY() + j);
yu[j] = _data[i][_parameters.dimX() + 2*_parameters.dimY() + j]; yu[j] = matrix(i, _parameters.dimX() + 2*_parameters.dimY() + j);
} }
} }
void vertical_segment::get(int i, vec& yl, vec& yu) const void vertical_segment::get(int i, vec& yl, vec& yu) const
{ {
yl.resize(_parameters.dimY()) ; yu.resize(_parameters.dimY()) ; auto matrix = matrix_view();
for(int j=0; j<_parameters.dimY(); ++j)
{ yl.resize(_parameters.dimY()) ; yu.resize(_parameters.dimY()) ;
yl[j] = _data[i][_parameters.dimX() + _parameters.dimY() + j] ; for(int j=0; j<_parameters.dimY(); ++j)
yu[j] = _data[i][_parameters.dimX() + 2*_parameters.dimY() + j] ; {
} yl[j] = matrix(i, _parameters.dimX() + _parameters.dimY() + j);
yu[j] = matrix(i, _parameters.dimX() + 2*_parameters.dimY() + j);
}
} }
vec vertical_segment::get(int i) const vec vertical_segment::get(int i) const
{ {
//SLOW !!!!! and useless // Omit the confidence interval data.
// call _data[i] auto relevant_rows = _parameters.dimX() + _parameters.dimY();
const int n = _parameters.dimX() + _parameters.dimY(); auto matrix =
vec res = _data[i].head(n); Map<MatrixXd, 0, OuterStride<> >(_data.get(), relevant_rows, size(),
OuterStride<>(column_number()));
return res ;
return matrix.col(i);
} }
void vertical_segment::set(int i, const vec& x) void vertical_segment::set(int i, const vec& x)
{ {
// Check if the input data 'x' has the size of a vertical segment (i.e. dimX+3*dimY), // Check if the input data 'x' has the size of a vertical segment (i.e. dimX+3*dimY),
// if it only has the size of a single value, then create the segment. // if it only has the size of a single value, then create the segment.
if(x.size() == _parameters.dimX() + 3*_parameters.dimY()) { if(x.size() == column_number()
_data[i] = x; || x.size() == _parameters.dimX() + _parameters.dimY())
{
} else if(x.size() == _parameters.dimX() + _parameters.dimY()) { auto matrix = matrix_view();
_data[i] = vs(x); auto row = matrix.row(i);
row.head(x.size()) = x;
} else { } else {
std::cerr << "<<ERROR>> Passing an incorrect element to vertical_segment::set" << std::endl; std::cerr << "<<ERROR>> Passing an incorrect element to vertical_segment::set" << std::endl;
throw; throw;
......
...@@ -81,7 +81,8 @@ class vertical_segment : public data ...@@ -81,7 +81,8 @@ class vertical_segment : public data
public: // methods public: // methods
vertical_segment(const parameters& params, vertical_segment(const parameters& params,
std::vector<vec>&& data); size_t size,