computation.hpp

Go to the documentation of this file.

 
#ifndef RTCTK_EXAMPLEDATATASK_COMPUTATION_HPP
#define RTCTK_EXAMPLEDATATASK_COMPUTATION_HPP
 
#include "rtctk/componentFramework/exceptions.hpp"
#include "rtctk/componentFramework/logger.hpp"
#include "rtctk/componentFramework/matrixBuffer.hpp"
#include "rtctk/exampleTopics/topics.hpp"
 
#include <chrono>
#include <numeric>
#include <string>
#include <vector>
// #include <cblas.h>
 
using namespace rtctk::componentFramework;
 
class Computation {
public:
    using TopicType = rtctk::exampleTopic::ScaoLoopTopic;
 
    Computation() {
        m_last_sample_id = 0;
        m_sample_idx = 0;
        m_iteration_idx = 0;
    }
 
    void SetStaticConfig(unsigned to_read, unsigned n_slopes, unsigned n_modes) {
        // read anything that is unlikely to change.
        m_n_samples_to_read = to_read;
        m_n_slopes = n_slopes;
        m_n_modes = n_modes;
 
        // set vectors and matrix to correct size.
        m_samples.resize(m_n_samples_to_read);
        m_avg_slopes.resize(m_n_slopes);
        m_s2m_matrix.resize(m_n_modes, m_n_slopes);
        m_avg_modes.resize(m_n_modes);
 
        // fill with zeros where needed.
        std::fill(m_avg_slopes.begin(), m_avg_slopes.end(), 0.0);
        std::fill(m_avg_modes.begin(), m_avg_modes.end(), 0.0);
 
        // some simple debug output
        LOG4CPLUS_DEBUG(GetLogger(), "m_n_slopes: " << m_n_slopes);
        LOG4CPLUS_DEBUG(GetLogger(), "m_n_modes: " << m_n_modes);
        LOG4CPLUS_DEBUG(GetLogger(), "m_n_samples_to_read: " << m_n_samples_to_read);
    }
 
    void SetDynamicConfig(MatrixBuffer<float>&& s2m_matrix) {
        m_s2m_matrix = std::move(s2m_matrix);
 
        if (m_n_modes != m_s2m_matrix.GetNrows() || m_n_slopes != m_s2m_matrix.GetNcols()) {
            std::stringstream err_text;
            err_text << "s2m wrong shape: ";
            err_text << " "
                     << "expected : " << m_n_modes << " x " << m_n_slopes;
            err_text << " "
                     << "received : " << m_s2m_matrix.GetNrows() << " x "
                     << m_s2m_matrix.GetNcols();
            using rtctk::componentFramework::RtctkException;
            CII_THROW(RtctkException, err_text.str());
        }
    }
 
    void Reset() {
        // This is invoked to clean up in case of stopped during buffering of data.
        m_last_sample_id = 0;
        m_sample_idx = 0;
        m_iteration_idx = 0;
    }
 
    void OnDataAvailable(TopicType const& sample) {
        m_last_sample_id = sample.sample_id;
        m_samples[m_sample_idx++] = sample.wfs.slopes;
    }
 
    struct Result {
        std::vector<float> const& avg_slopes;
        std::vector<float> const& avg_modes;
 
        struct {
            unsigned last_sample_id;
            unsigned iteration_idx;
            std::chrono::duration<double> elapsed;
        } stats;
    };
 
    Result Compute() {
        using rtctk::componentFramework::GetLogger;
        using rtctk::componentFramework::RtctkException;
 
        if (m_sample_idx != m_n_samples_to_read) {
            CII_THROW(RtctkException, "not all data available");
        }
 
        LOG4CPLUS_INFO(GetLogger(), "Computation::Compute() - samples: " << m_sample_idx);
 
        auto time_start = std::chrono::system_clock::now();
 
        ComputeAverageSlopes();
 
        ComputeAverageModes();
 
        auto elapsed = std::chrono::system_clock::now() - time_start;
 
        m_iteration_idx++;
        m_sample_idx = 0;
 
        return {m_avg_slopes, m_avg_modes, m_last_sample_id, m_iteration_idx, elapsed};
    }
 
private:
    void ComputeAverageSlopes() {
        // average matrix to avg_array this is overly simplistic and not optimised
        for (unsigned i = 0; i < m_n_samples_to_read; i++) {
            for (unsigned j = 0; j < m_n_slopes; j++) {
                m_avg_slopes[j] += m_samples[i][j];
            }
        }
 
        for (unsigned j = 0; j < m_n_slopes; j++) {
            m_avg_slopes[j] /= static_cast<float>(m_n_samples_to_read);
        }
    }
 
    void ComputeAverageModes() {
        // multiply avg_array with S2M matrix to project the average slopes onto modes.
        // TODO: openblas not available on Centos 7 and no point in making it a dependancy.
        //      for this example I have added a very basic matrix vector multiplication of my own
        //      It is single threaded and not optimised. It is here only for demonstration purposes.
        // cblas_sgemv(CblasRowMajor, CblasNoTrans, m_n_modes, m_n_slopes, 1.0, m_s2m_matrix.data(),
        // m_n_slopes, m_avg_slopes.data(), 1, 0.0, m_avg_modes.data(), 1);
        my_sgemv(
            m_s2m_matrix.data(), m_avg_slopes.data(), m_avg_modes.data(), m_n_modes, m_n_slopes);
    }
 
    // single threaded sgemv for row major matrix
    void
    my_sgemv(const float* mat, const float* in_vec, float* out_vec, uint32_t rows, uint32_t cols) {
        for (uint32_t i = 0; i < rows; ++i) {
            out_vec[i] = my_dot_product(mat + i * cols, in_vec, cols);
        }
    }
 
    // single threaded dot product
    float my_dot_product(const float* x, const float* y, uint32_t n) {
        double res = 0.0;
        for (uint32_t i = 0; i < n; ++i) {
            res += x[i] * y[i];
        }
        return res;
    }
 
    unsigned m_last_sample_id;
    unsigned m_iteration_idx;
    unsigned m_sample_idx;
 
    unsigned m_n_slopes;
    unsigned m_n_modes;
    unsigned m_n_samples_to_read;
 
    std::vector<float> m_avg_slopes;
    std::vector<float> m_avg_modes;
    std::vector<decltype(decltype(TopicType::wfs)::slopes)> m_samples;
    MatrixBuffer<float> m_s2m_matrix;
};
 
#endif  // RTCTK_EXAMPLEDATATASK_COMPUTATION_HPP