Skip to content

Commit bcba4d0

Browse files
author
Albert Akhriev
committed
Amdados application is done in general; need to test and debug thoroughly
1 parent e0c2a1e commit bcba4d0

13 files changed

Lines changed: 449 additions & 348 deletions

File tree

‎code/app/include/amdados/app/utils/amdados_utils.h‎

Lines changed: 126 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -48,91 +48,150 @@ T Bound(const T & v, const T & vmin, const T & vmax)
4848
}
4949

5050
//-------------------------------------------------------------------------------------------------
51-
// Function returns a squared value.
51+
// Function rounds the value to the nearest integer.
5252
//-------------------------------------------------------------------------------------------------
53-
template<typename T>
54-
T Square(const T & v)
53+
inline int Round(double val)
5554
{
56-
return (v * v);
55+
assert_true(std::fabs(val) < numeric_limits<int>::max());
56+
return static_cast<int>(std::floor(val + 0.5));
5757
}
5858

5959
//-------------------------------------------------------------------------------------------------
60-
// Functor converts 2D index (x,y) into a plain one.
61-
//-------------------------------------------------------------------------------------------------
62-
template<int SizeX, int SizeY>
63-
struct Sub2Ind {
64-
Sub2Ind() {
65-
static_assert(SizeX * SizeY < static_cast<int>(numeric_limits<int>::max()), "overflow");
66-
}
67-
68-
int operator()(int ix, int iy) const {
69-
assert_true((0 <= ix) && (ix < SizeX));
70-
assert_true((0 <= iy) && (iy < SizeY));
71-
//return (x + static_cast<int>(SizeX) * y);
72-
return (ix * SizeY + iy); // TODO: we already have observations based
73-
} // on this indexing: y changes faster, swap?
74-
};
75-
76-
//-------------------------------------------------------------------------------------------------
77-
// Functor converts a plane index into 2D one (x,y).
78-
//-------------------------------------------------------------------------------------------------
79-
template<int SizeX, int SizeY>
80-
struct Ind2Sub {
81-
Ind2Sub() {
82-
static_assert(SizeX * SizeY < static_cast<int>(numeric_limits<int>::max()), "overflow");
83-
}
84-
85-
void operator()(const int idx, int & x, int & y) const {
86-
assert_true((0 <= idx) && (idx < SizeX * SizeY));
87-
std::div_t divresult = std::div(idx, SizeY); // if y is the fastest
88-
x = divresult.quot;
89-
y = divresult.rem;
90-
}
91-
};
92-
93-
//-------------------------------------------------------------------------------------------------
94-
// Function reshapes a vector into 2D grid structure,
95-
// in Matlab notation: grid = reshape(vec, [SizeX, SizeY]).
60+
// Function creates a new directory or does nothing if it already exists.
61+
// TODO: this will not work on Windows, use STL "experimental" instead.
9662
//-------------------------------------------------------------------------------------------------
97-
template<int SizeX, int SizeY, typename GRID>
98-
void Reshape1Dto2D(GRID & grid, const allscale::utils::grid<double, SizeX * SizeY> & vec)
63+
inline void MakeDirectory(const char * dir)
9964
{
100-
// TODO: check grid sizes: must be SizeX by SizeY
101-
Sub2Ind<SizeX, SizeY> sub2ind;
102-
for (int i = 0; i < static_cast<int>(SizeX); i++) {
103-
for (int j = 0; j < static_cast<int>(SizeY); j++) {
104-
grid[{i,j}] = vec[{sub2ind(i,j)}];
105-
}
106-
}
65+
assert_true(dir != nullptr);
66+
std::string cmd("mkdir -p ");
67+
cmd += dir;
68+
int retval = std::system(cmd.c_str());
69+
retval = std::system("sync");
70+
(void)retval;
10771
}
10872

10973
//-------------------------------------------------------------------------------------------------
110-
// Function unrolls 2D grid structure into a vector, in Matlab notation: vec = grid(:).
74+
// Function creates an output directory inside the current one, which is supposed to be the
75+
// project root folder. If the directory is already exist all its content will be deleted.
76+
// TODO: this will not work on Windows, use STL "experimental" instead.
11177
//-------------------------------------------------------------------------------------------------
112-
template<int SizeX, int SizeY, typename GRID>
113-
void Reshape2Dto1D(allscale::utils::grid<double, SizeX * SizeY> & vec, const GRID & grid)
78+
void CreateAndCleanOutputDir(const std::string & dir)
11479
{
115-
// TODO: check grid sizes: must be SizeX by SizeY
116-
Sub2Ind<SizeX, SizeY> sub2ind;
117-
for (int i = 0; i < static_cast<int>(SizeX); i++) {
118-
for (int j = 0; j < static_cast<int>(SizeY); j++) {
119-
vec[{sub2ind(i,j)}] = grid[{i,j}];
120-
}
121-
}
80+
assert_true(!dir.empty());
81+
string cmd("mkdir -p ");
82+
cmd += dir;
83+
int retval = std::system(cmd.c_str());
84+
retval = std::system("sync");
85+
retval = std::system((string("/bin/rm -fr ") + dir + "/*.png").c_str());
86+
retval = std::system((string("/bin/rm -fr ") + dir + "/*.pgm").c_str());
87+
retval = std::system((string("/bin/rm -fr ") + dir + "/*.jpg").c_str());
88+
retval = std::system((string("/bin/rm -fr ") + dir + "/*.avi").c_str());
89+
retval = std::system("sync");
90+
(void)retval;
12291
}
12392

93+
//@{
12494
//-------------------------------------------------------------------------------------------------
125-
// Function creates a new directory or does nothing if it already exists.
95+
// Functions for global reduction across all the subdomains.
96+
// A T T E N T I O N: these functions must be used ONLY for testing, debugging or visualization.
12697
//-------------------------------------------------------------------------------------------------
127-
inline void MakeDirectory(const char * dir)
98+
double ReduceMean(const ::allscale::api::user::data::Grid<double,2> & grid)
12899
{
129-
assert_true(dir != nullptr);
130-
std::string cmd("mkdir -p "); // TODO: not portable, use STL "experimental" instead
131-
cmd += dir;
132-
int retval = std::system(cmd.c_str()); // TODO: mutex
133-
retval = std::system("sync");
134-
(void)retval;
100+
double sum = 0.0;
101+
for (int x = 0; x < SubDomGridSize[_X_]; ++x) {
102+
for (int y = 0; y < SubDomGridSize[_Y_]; ++y) { sum += grid[{x,y}]; }}
103+
return (sum / static_cast<double>(SubDomGridSize[_X_] * SubDomGridSize[_Y_]));
104+
}
105+
double ReduceAbsMin(const ::allscale::api::user::data::Grid<double,2> & grid)
106+
{
107+
double v = std::fabs(grid[{0,0}]);
108+
for (int x = 0; x < SubDomGridSize[_X_]; ++x) {
109+
for (int y = 0; y < SubDomGridSize[_Y_]; ++y) { v = std::min(v, std::fabs(grid[{x,y}])); }}
110+
return v;
111+
}
112+
double ReduceAbsMax(const ::allscale::api::user::data::Grid<double,2> & grid)
113+
{
114+
double v = std::fabs(grid[{0,0}]);
115+
for (int x = 0; x < SubDomGridSize[_X_]; ++x) {
116+
for (int y = 0; y < SubDomGridSize[_Y_]; ++y) { v = std::max(v, std::fabs(grid[{x,y}])); }}
117+
return v;
135118
}
119+
//@}
120+
121+
/*//-------------------------------------------------------------------------------------------------*/
122+
/*// Function returns a squared value.*/
123+
/*//-------------------------------------------------------------------------------------------------*/
124+
/*template<typename T>*/
125+
/*T Square(const T & v)*/
126+
/*{*/
127+
/*return (v * v);*/
128+
/*}*/
129+
130+
/*//-------------------------------------------------------------------------------------------------*/
131+
/*// Functor converts 2D index (x,y) into a plain one.*/
132+
/*//-------------------------------------------------------------------------------------------------*/
133+
/*template<int SizeX, int SizeY>*/
134+
/*struct Sub2Ind {*/
135+
/*Sub2Ind() {*/
136+
/*static_assert(SizeX * SizeY < static_cast<int>(numeric_limits<int>::max()), "overflow");*/
137+
/*}*/
138+
139+
/*int operator()(int ix, int iy) const {*/
140+
/*assert_true((0 <= ix) && (ix < SizeX));*/
141+
/*assert_true((0 <= iy) && (iy < SizeY));*/
142+
/*//return (x + static_cast<int>(SizeX) * y);*/
143+
/*return (ix * SizeY + iy); // TODO: we already have observations based*/
144+
/*} // on this indexing: y changes faster, swap?*/
145+
/*};*/
146+
147+
/*//-------------------------------------------------------------------------------------------------*/
148+
/*// Functor converts a plane index into 2D one (x,y).*/
149+
/*//-------------------------------------------------------------------------------------------------*/
150+
/*template<int SizeX, int SizeY>*/
151+
/*struct Ind2Sub {*/
152+
/*Ind2Sub() {*/
153+
/*static_assert(SizeX * SizeY < static_cast<int>(numeric_limits<int>::max()), "overflow");*/
154+
/*}*/
155+
156+
/*void operator()(const int idx, int & x, int & y) const {*/
157+
/*assert_true((0 <= idx) && (idx < SizeX * SizeY));*/
158+
/*std::div_t divresult = std::div(idx, SizeY); // if y is the fastest*/
159+
/*x = divresult.quot;*/
160+
/*y = divresult.rem;*/
161+
/*}*/
162+
/*};*/
163+
164+
/*//-------------------------------------------------------------------------------------------------*/
165+
/*// Function reshapes a vector into 2D grid structure,*/
166+
/*// in Matlab notation: grid = reshape(vec, [SizeX, SizeY]).*/
167+
/*//-------------------------------------------------------------------------------------------------*/
168+
/*template<int SizeX, int SizeY, typename GRID>*/
169+
/*void Reshape1Dto2D(GRID & grid, const allscale::utils::grid<double, SizeX * SizeY> & vec)*/
170+
/*{*/
171+
/*// TODO: check grid sizes: must be SizeX by SizeY*/
172+
/*Sub2Ind<SizeX, SizeY> sub2ind;*/
173+
/*for (int i = 0; i < static_cast<int>(SizeX); i++) {*/
174+
/*for (int j = 0; j < static_cast<int>(SizeY); j++) {*/
175+
/*grid[{i,j}] = vec[{sub2ind(i,j)}];*/
176+
/*}*/
177+
/*}*/
178+
/*}*/
179+
180+
/*//-------------------------------------------------------------------------------------------------*/
181+
/*// Function unrolls 2D grid structure into a vector, in Matlab notation: vec = grid(:).*/
182+
/*//-------------------------------------------------------------------------------------------------*/
183+
/*template<int SizeX, int SizeY, typename GRID>*/
184+
/*void Reshape2Dto1D(allscale::utils::grid<double, SizeX * SizeY> & vec, const GRID & grid)*/
185+
/*{*/
186+
/*// TODO: check grid sizes: must be SizeX by SizeY*/
187+
/*Sub2Ind<SizeX, SizeY> sub2ind;*/
188+
/*for (int i = 0; i < static_cast<int>(SizeX); i++) {*/
189+
/*for (int j = 0; j < static_cast<int>(SizeY); j++) {*/
190+
/*vec[{sub2ind(i,j)}] = grid[{i,j}];*/
191+
/*}*/
192+
/*}*/
193+
/*}*/
194+
136195

137196
} // end namespace utils
138197
} // end namespace app

‎code/app/include/amdados/app/utils/common.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,3 +30,5 @@ using ::std::string;
3030
} // end namespace app
3131
} // end namespace amdados
3232

33+
// Useful message for experimentation: "!!!!!!! B E W A R E: TEMPORARY CODE !!!!!!!"
34+

‎code/app/include/amdados/app/utils/image_writer.h‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
//-----------------------------------------------------------------------------
2+
// Author : Albert Akhriev, albert_akhriev@ie.ibm.com
3+
// Copyright : IBM Research Ireland, 2017
4+
//-----------------------------------------------------------------------------
5+
6+
#pragma once
7+
8+
namespace amdados {
9+
namespace app {
10+
namespace utils {
11+
12+
//=================================================================================================
13+
// Class (member function) converts a state field into image and saves it into "pgm" file.
14+
//=================================================================================================
15+
class ImageWriter
16+
{
17+
private:
18+
std::string m_output_dir; // output directory
19+
std::vector<unsigned char> m_image_buffer; // temporary image buffer
20+
21+
public:
22+
//-------------------------------------------------------------------------------------------------
23+
// Constructor.
24+
//-------------------------------------------------------------------------------------------------
25+
explicit ImageWriter(const std::string & output_dir) : m_output_dir(output_dir)
26+
{
27+
if (m_output_dir.empty()) m_output_dir.append(".");
28+
}
29+
30+
//-------------------------------------------------------------------------------------------------
31+
// Function writes the whole property field into a file in binary, gray-scaled PGM format,
32+
// where all the values are adjusted to [0..255] interval.
33+
// A T T E N T I O N: the function must be used outside pfor(..) loop.
34+
// \param title file title (extension will be "pgm").
35+
// \param field global grid of subdomains to be saved as an image.
36+
// \param time_index discrete timestamp.
37+
// \param write_image true, if saving into file was intended.
38+
// \return reference to image that can be plotted, for example, by Gnuplot.
39+
//-------------------------------------------------------------------------------------------------
40+
const std::vector<unsigned char> & Write(const char * title,
41+
const ::allscale::api::user::data::Grid<
42+
Matrix<NELEMS_X,NELEMS_Y>,2> & field,
43+
int time_index,
44+
bool write_image)
45+
{
46+
using ::allscale::api::user::pfor;
47+
48+
const double TINY = numeric_limits<double>::min() /
49+
std::pow(numeric_limits<double>::epsilon(),3);
50+
const bool flipY = false;
51+
52+
const int ImageSize = GLOBAL_NELEMS_X * GLOBAL_NELEMS_Y;
53+
if (m_image_buffer.capacity() < ImageSize) { m_image_buffer.reserve(ImageSize); }
54+
m_image_buffer.resize(ImageSize);
55+
56+
// Compute the minimum and maximum values of the property field.
57+
double lower = 0.0, upper = 0.0;
58+
{
59+
::allscale::api::user::data::Grid<double,2> minvalues(SubDomGridSize),
60+
maxvalues(SubDomGridSize);
61+
pfor(Origin, SubDomGridSize, [&](const point2d_t & idx) {
62+
minvalues[idx] = *(std::min_element(field[idx].begin(), field[idx].end()));
63+
maxvalues[idx] = *(std::max_element(field[idx].begin(), field[idx].end()));
64+
});
65+
lower = ReduceAbsMin(minvalues);
66+
upper = ReduceAbsMax(maxvalues);
67+
//assert_true(upper - lower > TINY);
68+
}
69+
70+
// Convert the property field into one-byte-per-pixel representation.
71+
pfor(Origin, SubDomGridSize, [&](const point2d_t & idx) {
72+
const auto & subfield = field[idx];
73+
for (int y = 0; y < NELEMS_Y; ++y) { int yg = Sub2GloY(idx, y);
74+
for (int x = 0; x < NELEMS_X; ++x) { int xg = Sub2GloX(idx, x);
75+
int pos = xg + GLOBAL_NELEMS_X * (flipY ? yg : (GLOBAL_NELEMS_Y - 1 - yg));
76+
if (!(static_cast<unsigned>(pos) < static_cast<unsigned>(ImageSize))) assert_true(0);
77+
m_image_buffer[pos] = static_cast<unsigned char>(
78+
Round(255.0 * (subfield(x,y) - lower) / std::max(upper - lower, TINY)) );
79+
}}
80+
});
81+
82+
// Write a file in binary, gray-scaled PGM format.
83+
if (write_image) {
84+
std::stringstream filename;
85+
filename << m_output_dir << "/" << title
86+
<< std::setfill('0') << std::setw(5) << time_index << ".pgm";
87+
std::ofstream f(filename.str(),
88+
std::ios_base::out | std::ios_base::binary | std::ios_base::trunc);
89+
assert_true(f.good()) << "failed to open file for writing: " << filename.str() << endl;
90+
f << "P5\n" << GLOBAL_NELEMS_X << " " << GLOBAL_NELEMS_Y << "\n" << int(255) << "\n";
91+
f.write(reinterpret_cast<const char*>(m_image_buffer.data()), m_image_buffer.size());
92+
f << std::flush;
93+
}
94+
return m_image_buffer;
95+
}
96+
97+
};
98+
99+
} // end namespace utils
100+
} // end namespace app
101+
} // end namespace allscale
102+

‎code/app/include/amdados/app/utils/kalman_filter.h‎

Lines changed: 19 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -16,15 +16,15 @@ template<int PROBLEM_SIZE, int NUM_OBSERVATIONS>
1616
class KalmanFilter
1717
{
1818
public:
19-
using matrix_t = Matrix<PROBLEM_SIZE, PROBLEM_SIZE>;
20-
using matrix_OxN_t = Matrix<NUM_OBSERVATIONS, PROBLEM_SIZE>;
21-
using matrix_NxO_t = Matrix<PROBLEM_SIZE, NUM_OBSERVATIONS>;
22-
using matrix_OxO_t = Matrix<NUM_OBSERVATIONS, NUM_OBSERVATIONS>;
23-
using vector_t = Vector<PROBLEM_SIZE>;
24-
using vector_obs_t = Vector<NUM_OBSERVATIONS>;
19+
typedef Matrix<PROBLEM_SIZE, PROBLEM_SIZE> matrix_t;
20+
typedef Matrix<NUM_OBSERVATIONS, PROBLEM_SIZE> matrix_OxN_t;
21+
typedef Matrix<PROBLEM_SIZE, NUM_OBSERVATIONS> matrix_NxO_t;
22+
typedef Matrix<NUM_OBSERVATIONS, NUM_OBSERVATIONS> matrix_OxO_t;
23+
typedef Vector<PROBLEM_SIZE> vector_t;
24+
typedef Vector<NUM_OBSERVATIONS> vector_obs_t;
2525

2626
private:
27-
Cholesky<PROBLEM_SIZE> m_chol; // Cholesky decomposition solver
27+
Cholesky<NUM_OBSERVATIONS> m_chol; // Cholesky decomposition solver
2828
LUdecomposition<PROBLEM_SIZE> m_lu; // LU decomposition solver
2929

3030
vector_t m_x_prior; // placeholder for the vector x_{k|k-1} = A * x
@@ -99,14 +99,19 @@ void IterateInverse(
9999
// x_prior = A * x, P_prior = A * P * A^t, where A is avaliable via its inversion B = A^{-1}.
100100
m_lu.Init(B); // decompose: B = L * U
101101
m_lu.Solve(m_x_prior, x); // x_prior = B^{-1} * x_{t}
102-
m_lu.SolveBatch(m_P_prior, P); // P_prior = B^{-1} * P (P symmetric!)
102+
#if 1
103+
m_lu.BatchSolve(m_P_prior, P); // P_prior = B^{-1} * P (P symmetric!)
103104
P = m_P_prior; // use P as a temporary matrix
104-
m_lu.SolveBatchTr(m_P_prior, P); // P_prior = B^{-1} * (B^{-1} * P)^t = A * P * A^t
105+
m_lu.BatchSolveTr(m_P_prior, P); // P_prior = B^{-1} * (B^{-1} * P)^t = A * P * A^t
105106
AddMatrices(m_P_prior, m_P_prior, Q); // P_prior = A * P * A^t + Q
106107
Symmetrize(m_P_prior); // correct loss of symmetry due to round-off errors
107108

108109
// Estimate posterior state and covariance.
109110
PosteriorEstimation(H, R, z, x, P);
111+
#else
112+
#pragma message("!!!!!!! B E W A R E: TEMPORARY CODE !!!!!!! NO KALMAN FILTERING")
113+
x = m_x_prior;
114+
#endif
110115
}
111116

112117
private:
@@ -118,8 +123,11 @@ void IterateInverse(
118123
// \param x out: new state.
119124
// \param P out: new covariance.
120125
//-------------------------------------------------------------------------------------------------
121-
void PosteriorEstimation(const matrix_OxN_t & H, const matrix_OxO_t & R, const vector_obs_t & z,
122-
vector_t & x, matrix_t & P)
126+
void PosteriorEstimation(const Matrix<NUM_OBSERVATIONS, PROBLEM_SIZE> & H,
127+
const Matrix<NUM_OBSERVATIONS, NUM_OBSERVATIONS> & R,
128+
const VectorView<NUM_OBSERVATIONS> & z,
129+
VectorView<PROBLEM_SIZE> & x,
130+
Matrix<PROBLEM_SIZE, PROBLEM_SIZE> & P)
123131
{
124132
// y = z - H * x_prior
125133
MatVecMult(m_y, H, m_x_prior);

0 commit comments

Comments
 (0)