Is it possible to add/subtract arrays read from a data file?-CodePudding

I want to perform some mathematical operations such as addition/subtraction, mean etc. on the data read from the file. The data is two dimensional, looks something like this: So now for example I want to subtract first row from the 3rd or 4th, it is very simple in python using indexing, but I have no idea how to do in c . Here is the code I'm using.

#include <boost/multi_array.hpp>
#include <boost/timer/timer.hpp>
#include <boost/range/irange.hpp>
#include <h5xx/h5xx.hpp>
#include <iostream>
#include <vector>
#include <algorithm>
#include <iterator>
#include <string>

using array_2d_t = boost::multi_array<float, 2>;

template <typename T> 
void print_array(T const& array)
{
    for (auto const& row : array) 
        { for (auto v : row)
            printf("f ", v);
        printf("\n"); //prints a new line similar t0 \n
    }
    std::cout << "\n End of file" << std::endl;
}

array_2d_t read_frame(std::string const& filename, unsigned frame_no) {
    //h5xx::file xaa("../../data/xaa.h5", h5xx::file::mode::in);
    h5xx::file xaa(filename, h5xx::file::mode::in);

    h5xx::group   g(xaa, "particles/lipids/box/positions");
    h5xx::dataset ds(g, "value");

    // determine dataset shape: frames, particle count, space dimension
    auto ds_shape = h5xx::dataspace(ds).extents<3>();
    array_2d_t arr(boost::extents[ds_shape[1]][ds_shape[2]]);

    std::vector<hsize_t> offsets{frame_no, 0, 0};
    std::vector<hsize_t> counts{1, arr.shape()[0], arr.shape()[1]};
    //std::vector<hsize_t> strid{1, arr.shape()[0], arr.shape()[1]};
    // slice ohne stride:
    h5xx::slice slice(offsets, counts);
    h5xx::read_dataset(ds, arr, slice);
    return arr;
}

int main(int argc, char const* argv[])
{
    if (argc < 2) {
        std::cout << "Usage: " << argv[0] << " input.h5" << std::endl;
        return -1;
    }
    std::string filename(argv[1]);
    print_array(read_frame(filename, 1));
    return 0; 
}

read_frame function takes the argument filename and frame_no. The file is divided into frames each containing (11214) rows, but it is not important to know here. So read_frame returns the two dimensional data for every frame. In the main() function I can print the data for every frame, but now suppose I want to subtract the first frame from the second, how can I do that? Does anyone know ho

CodePudding user response：

So, what the data looks like in text form is very irrelevant. What's relevant is that you have to arrays of the same dimensions and want to apply element-wise arithmetic operations to it.

You can use a proper Matrix library (like Boost UBlas, Eigen, Armadillo etc.).

Or you can write them for Boost MultiArray: How to execute mathematical operations between two boost::multi_arrays (C )?

Here's a sample:

#include <algorithm>
#include <boost/multi_array.hpp>
#include <boost/range/irange.hpp>
#include <boost/timer/timer.hpp>
#include <h5xx/h5xx.hpp>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>

namespace ArrayOperators {
    template <typename Op> struct ArrayOp {
        Op op;
        constexpr explicit ArrayOp(Op op) : op(op) {}

        template <typename T, typename Scalar, size_t Dim> 
        constexpr inline auto operator()(
            boost::multi_array<T, Dim> const& l,
            Scalar const& v) const
        {
            std::array<int, Dim> shape;
            std::copy_n(l.shape(), Dim, shape.data());

            using R = boost::multi_array<decltype(op(T{}, v)), Dim>;
            R result(shape);
            
            std::transform(
               l.data(), l.data() l.num_elements(),
               result.data(),
               [&op=op,v](auto const& el) { return op(el, v); });

            return result;
        }

        template <typename T, typename U, size_t Dim> 
        constexpr inline auto operator()(
            boost::multi_array<T, Dim> const& l,
            boost::multi_array<U, Dim> const& r) const
        {
            std::array<int, Dim> shape;
            std::copy_n(l.shape(), Dim, shape.data());
            assert(std::equal(shape.begin(), shape.end(), r.shape()));

            using R = boost::multi_array<decltype(op(T{}, U{})), Dim>;
            R result(shape);
            
            std::transform(
               l.data(), l.data() l.num_elements(),
               r.data(), result.data(),
               [&op=op](auto const& v1, auto const& v2) { return op(v1, v2); });

            return result;
        }
    };

    template <typename L, typename R>
    static constexpr inline auto operator (L const& l, R const& r) {
        return ArrayOp {std::plus<>{}} (l, r); }

    template <typename L, typename R>
    static constexpr inline auto operator-(L const& l, R const& r) {
        return ArrayOp {std::minus<>{}} (l, r); }

    template <typename L, typename R>
    static constexpr inline auto operator/(L const& l, R const& r) {
        return ArrayOp {std::divides<>{}} (l, r); }

    template <typename L, typename R>
    static constexpr inline auto operator*(L const& l, R const& r) {
        return ArrayOp {std::multiplies<>{}} (l, r); }
}

using array_2d_t = boost::multi_array<float, 2>;

template <typename T> 
void print_array(T const& array)
{
    for (auto const& row : array) 
        { for (auto v : row)
            printf("f ", v);
        printf("\n"); //prints a new line similar t0 \n
    }
    std::cout << "\n End of file" << std::endl;
}

array_2d_t read_frame(std::string const& filename, unsigned frame_no) {
    //h5xx::file xaa("../../data/xaa.h5", h5xx::file::mode::in);
    h5xx::file xaa(filename, h5xx::file::mode::in);

    h5xx::group   g(xaa, "particles/lipids/box/positions");
    h5xx::dataset ds(g, "value");

    // determine dataset shape: frames, particle count, space dimension
    auto ds_shape = h5xx::dataspace(ds).extents<3>();
    array_2d_t arr(boost::extents[ds_shape[1]][ds_shape[2]]);

    std::vector<hsize_t> offsets{frame_no, 0, 0};
    std::vector<hsize_t> counts{1, arr.shape()[0], arr.shape()[1]};
    //std::vector<hsize_t> strid{1, arr.shape()[0], arr.shape()[1]};
    // slice ohne stride:
    h5xx::slice slice(offsets, counts);
    h5xx::read_dataset(ds, arr, slice);
    return arr;
}

int main(int argc, char const* argv[])
{
    if (argc < 2) {
        std::cout << "Usage: " << argv[0] << " input.h5" << std::endl;
        return -1;
    }
    std::string filename(argv[1]);

    using namespace ArrayOperators;
    auto avg4 = //
        (read_frame(filename, 0)   read_frame(filename, 1)  
         read_frame(filename, 2)   read_frame(filename, 3)) /
        4.0;

    print_array(read_frame(filename, 4) - avg4);
    return 0; 
}

Prints:

 -0.610001   0.410004  -0.150000 
  0.372501  -0.395000  -0.037500 
 -0.062500   0.027500  -0.122500 
  0.307503   0.184998  -0.212500 
  0.150000   0.009995   0.362500 
 -0.497500   0.197500   0.217500 
  0.405003  -0.185000   0.115000 
 -0.072500   0.247498   0.237500 
 -0.480000   0.034998   0.107500 
 -0.130000  -0.162499  -0.087500 
(...)
 -0.340000  -0.280001  -0.040000 
 -0.572502  -0.747505  -0.005000 
  0.392494   0.012500   0.045000 
  0.500000  -0.040000  -0.162500 
  0.355000   0.700000   0.065000 

 End of file