Search code examples
c++boost

What datatype would is expected for the arrays in c++?


I've been trying to assign a datatype to the three-dimensional arrays, e.g., double, but keep getting the error

                   error: request for member 'size' in 'msd_x', which is of non-class 
                          type 'const sample_type' {aka 'const long unsigned int'}
                   -> size_t N = msd_x.size();
class mean_square_displacement
{
public:
    typedef boost::multi_array_types::size_type sample_type;
    //typedef fixed_vector<double, 3> result_type;
    typedef double result_type;

    result_type operator() (sample_type const& msd_x, sample_type const& msd_y, sample_type const& msd_z) const
    {
        size_t N = msd_x.size();
        result_type msd = 0;
        for (unsigned int i = 0; i < N; ++i) {
            for (unsigned int j = 0; j < N; ++j) {
                auto dr = (msd_x[i+1][j] - msd_x[i][j])  + (msd_y[i+1][j+1] - msd_y[i][j+1]) + (msd_z[i+1][j+2] - msd_z[i][j+2]); 
            }
            // accumulate squared displacements
            msd += dr * dr;
        }
        return msd / N;

    }
}

Along with that it throws one more error:

                   error: invalid types 'const sample_type {aka const long  
                          unsigned int}[unsigned int]' for array subscript

I guess the way code accessing the arrays is fine. But it still throws the error. The data code reads is of double datatype, shown in the code below:

int main(int argc, char const* argv[])
{
    correlator::multi_tau_correlator<double> corr(    // TODO replace sample type
        nsamples * sampling_interval / 30            trajectory length
      , sampling_interval                        // time resolution at lowest level
      , 10   
    );
    // define time correlation functions
    auto msd = make_correlation(correlator::mean_square_displacement(), corr);
    corr.add_correlation(msd);
    // main loop for all coordinates
    for (size_t i = 1; i < nsamples; ++i) {
        auto position_array = first_rows[i]; 
        // append data to the correlator, which possibly computes some time correlations
        corr.sample(position_array);
    }
    corr.finalise();

    return 0;
}

Does anyone point out what is wrong with the datatypes?


Solution

  • I'm guessing here based on the many questions leading up to this.

    It is pretty obvious that you do not want sample type to be a scalar, but a 2-dimensional array.

    I will sketch a generic short-cut that would allow you to write mean_square_displacement::operator() to accept those, potentially even without knowing the concrete type(s) of the arguments.

    Caution

    However, I want to first caution that there is no way to know whether that will work with the multi_tau_correlator framework that you've never described. It seems like you are trying to apply the correlator in a parallel fashion.

    Firstly it's unknown whether the correlator framework allows you to do that (technically) and secondly it is likely not going to be fast, as this kind of parallelism is the realm of highly optimized libraries (that use actually vectorized operations like SIMD, AVX, OpenGL/Cuda).

    All I see here is a raw quadratic loop, which does not bode well for efficiency.

    The blind fix

    Since you don't want sample_type to be what you define it to be, simply define it as something else, preferrably, the thing you need it to be!

    typedef array_2d_t sample_type;
    

    There's a good chance that this is not what you actually are passing either (I suspect a sub_array view or (optionally strided) multi_array_view slice of a 3D multi array). But you can let the compiler figure it out based on the fact that you know you can index the arguments like this arg[int][int]:

    template <typename Sample>
    result_type operator()(Sample const& msd_x, Sample const& msd_y,
                           Sample const& msd_z) const
    

    This lets the compiler deduce the concrete type, assuming all arguments have the same statical type. Otherwise, you can go even more open:

    template <typename Sx, typename Sy, typename Sz>
    result_type operator()(Sx const& msd_x, Sy const& msd_y,
                           Sz const& msd_z) const
    

    More Caution

    I can see you're out of your depth. E.g.

    • in your mean_square_displacement we see auto dr being used outside its scope. It will not compile. For the above I have GUESSED that you meant to place the accumulation into the inner loop.
    • you're also explicitly addressing msd_x out of bounds (index N doesn't exist, since N == msd_x.size()).
    • Likely msd_y/msd_z have the same extents, so they too have this issue, but even with +2
    • If course the naive fix to loop till N-2 risks UB again unless you make sure that N is never <2

    Here's the minimum fixes that make this code look like it might technically work, though:

    struct mean_square_displacement {
        //typedef fixed_vector<double, 3> result_type;
        typedef double result_type;
    
        template <typename Sample>
        result_type operator()(Sample const& msd_x, Sample const& msd_y,
                               Sample const& msd_z) const
        {
            size_t N = msd_x.size();
            assert(N>=2);
            result_type msd = 0;
            for (unsigned int i = 0; i < N-1; ++i) {
                for (unsigned int j = 0; j < N-2; ++j) {
                    auto dr = //
                        (msd_x[i + 1][j + 0] - msd_x[i][j + 0]) +
                        (msd_y[i + 1][j + 1] - msd_y[i][j + 1]) +
                        (msd_z[i + 1][j + 2] - msd_z[i][j + 2]);
                    // accumulate squared displacements
                    msd += dr * dr;
                }
            }
            return msd / N;
    
        }
    };