Commit 25aefd05 authored by Jan Oliver Oelerich's avatar Jan Oliver Oelerich

Intermediate state of something. I will now try to use cereal for serialization...

parent 6666faea
The MIT License (MIT)
Copyright (c) 2015 Mohamed-Yassine MADDOURI
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
[![Build Status](https://travis-ci.org/maddouri/hyper_array.svg?branch=master)](https://travis-ci.org/maddouri/hyper_array)
# hyper_array
`hyper_array::array` is a simple, templated, multi-dimensional array class that is inspired by [`orca_array`](https://github.com/astrobiology/orca_array). It uses modern C++ techniques in order to achieve good performance, code clarity and platform independence.
`hyper_array` is header-only, contained in a single file -- [`hyper_array.hpp`](include/hyper_array/hyper_array.hpp) -- and doesn't depend on external libraries. Its only requirement is a C++11-compliant compiler.
<!--
grep -e '^#\+' README.md | sed 's/#\{3\}/ */' | sed 's/#\{2\}/ */' | sed 's/#\{1\}/*/' | sed 's/^\([^A-Za-z]\+\)\(.\+\)$/\1[\2](#\2)/'
-->
* [Basics](#basics)
* [Main Class](#main-class)
* [Array Order](#array-order)
* [Construction](#construction)
* [Assignment](#assignment)
* [Element Access](#element-access)
* [Standard Library Compatibility](#standard-library-compatibility)
* [Development](#development)
## Basics
### Main Class
The class template `hyper_array::array<ValueType, Dimensions, Order>` represents a `Dimensions`-dimension array of `ValueType` elements. Therefore, the type and number of dimensions are specified at compile-time. The length along each dimension can be set at run-time.
### Array Order
The third template argument --`Order`-- designates the [array order](https://en.wikipedia.org/wiki/Row-major_order) which can be either **row-major** _(i.e. C convention)_ or **column-major** _(i.e. Fortran convention)_ (cf. `hyper_array::array_order` enum).
By default, if `Order` is not specified, the order is row-major.
### Construction
A new array can be instantiated using one of the following constructors:
```c++
#include "hyper_array/hyper_array.hpp"
using hyper_array::array;
/// the usual way of constructing hyper arrays
array(DimensionLengths... dimensionLengths);
// usage example
array<double, 3> my3DArray{32, 64, 128};
/// copy constructor
array(const array_type& other);
// usage example
auto arrayCopy = my3DArray;
/// move constructor
array(array_type&& other);
// usage example
auto consumer = std::move(my3DArray);
/// create a new hyper array from "raw data"
array(::std::array<size_type, Dimensions> lengths, value_type* rawData);
// usage example
double* rawData = new double[262144];
array<double, 3> dataWrapper{{32, 64, 128}, rawData};
/// create and initialize a hyper array
/// given the dimension lengths and the array elements
array(::std::array<size_type, Dimensions> lengths, // dimension lenths
std::initializer_list<value_type> values, // array elements (you can provide less than size() elements)
const value_type& defaultValue = {}); // default initialization value (in case values.size() < size())
// usage example
const array<double, 2> constantArray{
{2, 3}, // dimension lengths
{11, 12, 13,
21, 22, 23} // array elements
};
```
### Assignment
`hyper_array` can be assigned for performing both deep copying and moving:
```c++
array<double, 4> hyperArray{42, 42, 42, 42};
/// copy assignment
array_type& operator=(const array_type& other);
// usage example
array<double, 4> someArray{32, 64, 128, 256};
// ...
hyperArray = someArray;
/// move assignment
array_type& operator=(array_type&& other);
// usage example
array<double, 4> temporaryArray{32, 64, 128, 256};
// ...
hyperArray = std::move(temporaryArray);
```
### Element Access
Single elements can be accessed for reading and assignment using various methods:
```c++
/// access using an index tuple
value_type& operator()(Indices... indices);
const value_type& operator()(Indices... indices) const;
// usage example
array<double, 3> arr{4, 5, 6};
arr(3, 1, 4) = 3.14;
std::cout << "arr(3, 1, 4) == " << arr(3, 1, 4) << std::endl; // arr(3, 1, 4) = 3.14
/// access using an index tuple, with index bounds checking
value_type& at(Indices... indices);
const value_type& at(Indices... indices) const;
// usage example
array<double, 3> arr{4, 5, 6};
arr.at(3, 1, 4) = 3.14;
std::cout << "arr.at(3, 1, 4) == " << arr.at(3, 1, 4) << std::endl; // arr.at(3, 1, 4) = 3.14
/// access using the index of the element in the underlying data array
value_type& operator[](const index_type idx);
const value_type& operator[](const index_type idx) const;
// usage example
array<int, 3> arr{4, 5, 6};
arr[100] = 314;
cout << "arr[100] == arr(3, 1, 4): " << std::boolalpha << (arr[100] == arr(3, 1, 4)) << endl; // arr[100] == arr(3, 1, 4): true
```
### Standard Library Compatibility
Currently, `hyper_array::array` implements the same iterators as `std::array`, which makes it compatible with most of the C++ Standard Library's algorithms and containers, as well as the range-based for loop syntax introduced in C++11. In addition, `operator<<()` is overloaded in order to provide an easy way to visualize array information.
```c++
// range-based for loop
array<double, 3> aa{4, 5, 6};
for (auto& x : aa) {
x = 1; // initialize the array
}
// algorithms
std::iota(aa.begin(), aa.end(), 1); // fill aa with a sequence of consecutive numbers
array<double, 3> bb{aa.lengths()}; // array::lengths() returns the lengths along each direction
std::copy(aa.begin(), aa.end(), bb.rbegin()); // reverse-copy
array<double, 3> cc{aa.lengths()};
// cc = aa + bb
std::transform(aa.begin(), aa.end(),
bb.begin(),
cc.begin(),
[](double a, double b) { return a + b; });
// stream operator
cout << "cc: " << cc << endl;
// cc: [dimensions: 3 ][lengths: 4 5 6 ][coeffs: 30 6 1 ][size: 120 ][data: 121 121 121 ...]
```
## Development
`hyper_array` is in constant development and new features will be added when appropriate. The goal is to keep it simple but useful, and efficient but maintainable.
This diff is collapsed.
/*
Copyright (c) 2015, 2016 Florian Dang
Permission is hereby granted, free of charge, to any person obtaining a copy of this
software and associated documentation files (the "Software"), to deal in the Software
without restriction, including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or
substantial portions of the Software.
*/
#ifndef MULTIDIM_GRID_HPP_
#define MULTIDIM_GRID_HPP_
#include <ostream>
#include <type_traits>
#include <utility>
#include <iterator>
#include <iostream>
#include <array>
template<typename T, size_t N>
using array = std::array<T,N>;
using std::next;
namespace multidim
{
namespace _impl_multi_grid
{
// meta functions
template<typename T>
constexpr T meta_prod(T x) { return x; }
template<typename T, typename... Ts>
constexpr T meta_prod(T x, Ts... xs) { return x * meta_prod(xs...); }
template<typename Iter, size_t D0, size_t... DIMS> struct flatten_t;
template <typename Iter, size_t D0>
constexpr size_t meta_flatten(Iter first, std::index_sequence<D0> seq)
{
return flatten_t<Iter, D0>::compute(first, seq);
}
template <typename Iter, size_t D0, size_t... DIMS>
constexpr size_t meta_flatten(Iter first, Iter last, std::index_sequence<D0, DIMS...> seq)
{
return flatten_t<Iter, D0, DIMS...>::compute(first, last, seq);
}
template<typename Iter, size_t D0>
struct flatten_t<Iter, D0>
{
static constexpr size_t compute(Iter first, Iter, std::index_sequence<D0>)
{
return *first;
}
};
template <typename Iter, size_t D0, size_t... DIMS>
struct flatten_t
{
static constexpr size_t compute(Iter first, Iter last, std::index_sequence<D0, DIMS...>)
{
return *first * meta_prod(DIMS...) +
meta_flatten(next(first), last, std::index_sequence<DIMS...>{} );
}
};
template<typename T, size_t... DIMS>
constexpr size_t flatten_to_index(const array<size_t, sizeof...(DIMS)>& coord)
{
// The problem here is that begin and end are constexpr since C++17 only...
return meta_flatten(coord.cbegin(), coord.cend(), std::index_sequence<DIMS...>{});
}
template<typename T, size_t... DIMS>
constexpr auto unflatten_to_coordinates(size_t idx) -> array<size_t, sizeof...(DIMS)>
{
const size_t num_dims = sizeof...(DIMS);
constexpr auto size_per_dim = array<size_t, num_dims>{{DIMS...}};
auto coord = array<size_t, num_dims>{};
size_t prod_dims = meta_prod(DIMS...);
size_t r = 0;
for(size_t i = 0; i <= num_dims - 2; i ++) {
// prod_dims /= size_per_dim[num_dims - i - 1];
// coord[num_dims-i-1] = idx / prod_dims;
prod_dims /= size_per_dim[i];
coord[i] = idx / prod_dims;
r = idx % prod_dims;
idx = r;
}
// coord[0] = r;
coord[num_dims - 1] = r;
return coord;
}
template<typename T, size_t... DIMS>
class Grid
{
public:
static constexpr size_t num_dims = sizeof...(DIMS);
static constexpr size_t prod_dims = meta_prod(DIMS...);
static_assert(num_dims > 0, "Grid dimension needs to be > 0");
static_assert(prod_dims > 0, "All dimension size must be > 0");
using ArrayValues = array<T,prod_dims>;
using ArrayCoord = array<size_t,num_dims>;
using value_type = typename ArrayValues::value_type; // T
using reference = typename ArrayValues::reference; // T&
using const_reference = typename ArrayValues::const_reference; // const T&
using size_type = typename ArrayValues::size_type; // size_t
using iterator = typename ArrayValues::iterator; // random access iterator
using const_iterator = typename ArrayValues::const_iterator;
using difference_type = typename ArrayValues::difference_type;
Grid() : Grid(ArrayValues{}) {}
Grid(const ArrayValues& values) : values_(values) {}
Grid(const Grid&) = default;
Grid& operator=(const Grid&) noexcept = default;
Grid(Grid&&) = default;
Grid& operator=(Grid&&) noexcept = default;
~Grid() = default;
bool operator==(const Grid& other) { return (values_ == other.values_); }
bool operator!=(const Grid& other) { return (values_ != other.values_); }
size_type size() { return prod_dims; }
size_type max_size() { return values_.max_size(); }
void swap(const Grid& other) noexcept { values_.swap(other.values_); }
constexpr iterator begin() { return values_.begin(); }
constexpr const_iterator begin() const { return values_.begin(); }
constexpr const_iterator cbegin() const { return values_.cbegin(); }
constexpr iterator end() { return values_.end(); }
constexpr const_iterator end() const { return values_.end(); }
constexpr const_iterator cend() const { return values_.cend(); }
constexpr reference operator[] (size_type idx) { return values_[idx]; };
constexpr const_reference operator[] (size_type idx) const { return values_[idx]; };
constexpr reference operator[] (const ArrayCoord& coord)
{
return values_[flatten_to_index<T,DIMS...>(coord)];
}
constexpr const_reference operator[] (const ArrayCoord& coord) const
{
return values_[flatten_to_index<T,DIMS...>(coord)];
};
// reference operator()(std::initializer_list<size_type> l) {
// return (*this)[ArrayCoord(l)];
// }
// const_reference operator()(std::initializer_list<size_type> l) const {
// return const_cast<reference>(static_cast<const Grid&>(*this)(l));
// };
constexpr auto unflatten(size_type idx) const -> ArrayCoord
{
return unflatten_to_coordinates<T,DIMS...>(idx);
}
constexpr size_type flatten(const ArrayCoord& coord) const
{
return flatten_to_index<T,DIMS...>(coord);
}
private:
ArrayValues values_;
};
} // namespace _impl_multi_grid
using _impl_multi_grid::Grid;
using _impl_multi_grid::flatten_to_index;
using _impl_multi_grid::unflatten_to_coordinates;
} // namespace multidim
#endif // MULTIDIM_GRID_HPP_
......@@ -34,9 +34,10 @@ set(INTERNAL_FILES
classes/Slice.cpp classes/Slice.hpp
libatomic/Atom.hpp
libatomic/Element.hpp
libatomic/Cell.hpp
libatomic/elements_json.hpp
classes/SimulationState.hpp
utilities/Wave.hpp
libatomic/elements_json.hpp
utilities/mpi.hpp
utilities/memory.hpp
utilities/output.hpp
......
......@@ -236,7 +236,7 @@ void IO::initBuffers(const shared_ptr<GridManager> &gridman) {
}
void IO::writeParams(const shared_ptr<GridManager> &gridman, const string & conf_file_string) {
void IO::writeParams(const shared_ptr<GridManager> &gridman, const string &conf_file_string) {
Params &p = Params::getInstance();
auto &mpi_env = mpi::Environment::getInstance();
......@@ -808,8 +808,8 @@ void IO::dumpWave(string filename, const Wave &wave) {
fclose(f);
}
void IO::initCrystalFromNCFile(const shared_ptr<Crystal> &crystal) {
Params &p = Params::getInstance();
std::shared_ptr<stemsalabim::atomic::Cell> IO::initCrystalFromXYZFile(const std::string &filename) {
atomic::ElementProvider elp = atomic::ElementProvider::getInstance();
atomic::Scattering &sf = atomic::Scattering::getInstance();
......@@ -825,7 +825,20 @@ void IO::initCrystalFromNCFile(const shared_ptr<Crystal> &crystal) {
vector<size_t> a_len;
if(mpi_env.isMaster()) {
NCFile f(p.paramsFileName(), false, true);
ifstream input(filename);
string line;
vector<string> tokens;
double x, y, z, msd;
size_t c = 0;
while(std::getline(input, line)) {
if(c == 0) {
// first line, number of atoms
algorithms::trim(line);
int n_atoms = std::stoi(line);
}
}
auto g = f.group("AMBER");
auto lv = g.var("lattice_coordinates");
......@@ -865,25 +878,199 @@ void IO::initCrystalFromNCFile(const shared_ptr<Crystal> &crystal) {
}
}
for(const auto & at: atom_types) {
crystal->addElement(elp.elementBySymbol(at));
for(const auto &at: atom_types) {
cell->addElement(elp.elementBySymbol(at));
sf.initPotential(elp.elementBySymbol(at));
}
for(size_t i_atom = 0; i_atom < e_len[1]; i_atom++) {
crystal->addAtom(make_shared<atomic::Atom>(lattice_coords[3 * i_atom + 0],
lattice_coords[3 * i_atom + 1],
lattice_coords[3 * i_atom + 2],
msds[i_atom],
elp.elementBySymbol(atom_types[element_ids[i_atom]]),
i_atom));
cell->addAtom(make_shared<atomic::Atom>(lattice_coords[3 * i_atom + 0],
lattice_coords[3 * i_atom + 1],
lattice_coords[3 * i_atom + 2],
msds[i_atom],
elp.elementBySymbol(atom_types[element_ids[i_atom]]),
i_atom));
}
cell->setLengths(system_lengths[0], system_lengths[1], system_lengths[2]);
// this order is very important!
atomic::ElementProvider p = atomic::ElementProvider::getInstance();
// some useful variable declarations
ifstream input(filename);
string line;
vector<string> tokens;
double x, y, z, msd;
unsigned int slice = 0;
// the first line is only the number of atoms.
getline(iss, line);
// from the second line, extract the system size
// and convert from nm to angstrom
getline(iss, line);
//determine, if we have extended xyz or normal xyz format. In normal xyz format, this
// line contains only 3 numbers, corresponding to cell size in x y and z direction. In
// extended xyz, it is like this: Lattice="lx 0.0 0.0 0.0 ly 0.0 0.0 0.0 lz"
if(line.find("Lattice=") != string::npos) {
regex e("Lattice=\"([0-9\\seE+-\\.]+)\"", regex_constants::icase);
smatch m;
if(regex_search(line, m, e)) {
tokens = algorithms::split(m[1], ' ');
algorithms::trim(tokens[0]);
algorithms::trim(tokens[4]);
algorithms::trim(tokens[8]);
_size_x = std::stod(tokens[0]);
_size_y = std::stod(tokens[4]);
_size_z = std::stod(tokens[8]);
} else {
output::error("Wrong syntax of extended XYZ file line 2.");
}
} else {
tokens = algorithms::split(line, ' ');
algorithms::trim(tokens[0]);
algorithms::trim(tokens[1]);
algorithms::trim(tokens[2]);
_size_x = std::stod(tokens[0]);
_size_y = std::stod(tokens[1]);
_size_z = std::stod(tokens[2]);
}
atomic::Scattering &sf = atomic::Scattering::getInstance();
// the rest of the lines is atoms with positions
size_t atom_id = 0;
while(getline(iss, line)) {
if(line.length() == 0)
continue;
tokens = algorithms::split(line, ' ');
if(tokens.size() != 5 && tokens.size() != 6) {
output::error("Init file needs to have 5 or 6 columns!\n");
}
bool with_slices = (tokens.size() == 6);
algorithms::trim(tokens[0]);
algorithms::trim(tokens[1]);
algorithms::trim(tokens[2]);
algorithms::trim(tokens[3]);
algorithms::trim(tokens[4]);
if(with_slices) {
algorithms::trim(tokens[5]);
slice = std::stoi(tokens[5]);
}
// x y z are columns 1,2,3 and the mean square displacement is column 4. Convert to angstrom
x = std::stod(tokens[1]);
y = std::stod(tokens[2]);
z = std::stod(tokens[3]);
msd = std::stod(tokens[4]);
// set the site's occupant
if(with_slices) {
_atoms.push_back(std::make_shared<atomic::Atom>(x,
y,
z,
msd,
p.elementBySymbol(tokens[0]),
atom_id,
slice));
} else {
_atoms.push_back(std::make_shared<atomic::Atom>(x, y, z, msd, p.elementBySymbol(tokens[0]), atom_id));
}
_elements.insert(p.elementBySymbol(tokens[0]));
sf.initPotential(p.elementBySymbol(tokens[0]));
atom_id++;
}
}
std::shared_ptr<stemsalabim::atomic::Cell> IO::initCrystalFromNCFile(const std::string &filename) {
atomic::ElementProvider elp = atomic::ElementProvider::getInstance();
atomic::Scattering &sf = atomic::Scattering::getInstance();
auto &mpi_env = mpi::Environment::getInstance();
vector<short> element_ids;
vector<float> lattice_coords;
vector<float> msds;
vector<float> system_lengths;
vector<string> atom_types;
vector<size_t> e_len;
vector<size_t> a_len;
if(mpi_env.isMaster()) {
NCFile f(filename, false, true);
auto g = f.group("AMBER");
auto lv = g.var("lattice_coordinates");
auto ev = g.var("element");
auto av = g.var("atom_types");
auto mv = g.var("msd");
auto sv = g.var("slice");
e_len = ev.len();
a_len = av.len();
atom_types.resize(a_len[0]);
char buff[100];
for(size_t i = 0; i < a_len[0]; ++i) {
av.get(vs({i, 0}), vs({1, a_len[1]}), buff);
atom_types[i] = buff;
}
element_ids = ev.get<short>(vs({0, 0}), vs({1, e_len[1]}));
lattice_coords = lv.get<float>(vs({0, 0, 0}), vs({1, e_len[1], 3}));
msds = mv.get<float>(vs({0, 0}), vs({1, e_len[1]}));
system_lengths = g.var("system_lengths").get<float>();
}
if(mpi_env.isMpi()) {
mpi_env.broadcast(a_len, 0);
mpi_env.broadcast(e_len, 0);
mpi_env.broadcast(element_ids, 0);
mpi_env.broadcast(lattice_coords, 0);
mpi_env.broadcast(msds, 0);
mpi_env.broadcast(system_lengths, 0);
if(mpi_env.isSlave())
atom_types.resize(a_len[0]);
for(size_t i = 0; i < a_len[0]; ++i) {
mpi_env.broadcast(atom_types[i], 0);
}
}
for(const auto &at: atom_types) {
cell->addElement(elp.elementBySymbol(at));
sf.initPotential(elp.elementBySymbol(at));
}
for(size_t i_atom = 0; i_atom < e_len[1]; i_atom++) {
cell->addAtom(make_shared<atomic::Atom>(lattice_coords[3 * i_atom + 0],
lattice_coords[3 * i_atom + 1],
lattice_coords[3 * i_atom + 2],
msds[i_atom],
elp.elementBySymbol(atom_types[element_ids[i_atom]]),
i_atom));
}
crystal->setLengths(system_lengths[0], system_lengths[1], system_lengths[2]);
crystal->generateSlices();
cell->setLengths(system_lengths[0], system_lengths[1], system_lengths[2]);
}
void IO::initDisplacementsFromNCFile(const shared_ptr<FPConfManager> &fpman) {
......@@ -969,7 +1156,7 @@ vector<vector<float>> IO::readStoredPotential(const std::shared_ptr<FPConfManage
for(size_t i = 0; i < len[2]; ++i) {
for(size_t j = 0; j < len[3]; ++j) {
pot[i][j] = flattened_potential[i*len[3]+j];
pot[i][j] = flattened_potential[i * len[3] + j];
}
}
......
......@@ -33,6 +33,7 @@
#include "../utilities/algorithms.hpp"
#include "../utilities/netcdf_utils.hpp"
#include "../utilities/Wave.hpp"
#include "../libatomic/Cell.hpp"
namespace stemsalabim {
......@@ -152,7 +153,8 @@ namespace stemsalabim {
*/
static void dumpWave(std::string filename, const Wave &wave);
void initCrystalFromNCFile(const std::shared_ptr<Crystal> &crystal);
static std::shared_ptr<atomic::Cell> initCrystalFromNCFile(const std::string & filename);
static std::shared_ptr<atomic::Cell> initCrystalFromXYZFile(const std::string & filename);
void initDisplacementsFromNCFile(const std::shared_ptr<FPConfManager> &fpman);
......
......@@ -616,8 +616,6 @@ void Simulation::initBuffers() {
Params &prms = Params::getInstance();
typedef typename memory::buffer::number_buffer<float> buf_type;
auto size_type_size = sizeof(typename vector<float>::size_type);
// The serialization buffer is for sending work packages around. It should be larger than the
// theoretical maximal size of any work package filled with data.
unsigned long number_intensities_per_pixel = _gridman->adfDetectorGrid().size() * _gridman->adfSliceCoords().size();
......@@ -625,12 +623,6 @@ void Simulation::initBuffers() {
_gridman->storedCbedSizeY() *
_gridman->cbedSliceCoords().size()) : 0;
unsigned long max_size_work_package = sizeof(ScanPoint) +
2 * size_type_size +
number_intensities_per_pixel * sizeof(float) +
number_cbed_per_pixel * sizeof(float);
_serialization_buffer.resize(prms.workPackageSize() * max_size_work_package + /*Add 10MB*/10 * 1024 * 1024);
_adf_intensity_buffer = std::make_shared<buf_type>(number_intensities_per_pixel, 2 * prms.workPackageSize());
_cbed_intensity_buffer = std::make_shared<buf_type>(number_cbed_per_pixel, 2 * prms.workPackageSize());
}
......
......@@ -107,7 +107,6 @@ namespace stemsalabim {
/// memory fragmentation errors.
std::shared_ptr<memory::buffer::number_buffer<float>> _adf_intensity_buffer;
std::shared_ptr<memory::buffer::number_buffer<float>> _cbed_intensity_buffer;
std::vector<char> _serialization_buffer;