text: /* file: linear_model_train_normeq_kernel.h */
  /*******************************************************************************
  * Copyright 2014-2021 Intel Corporation
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
  * You may obtain a copy of the License at
  *
  *     http://www.apache.org/licenses/LICENSE-2.0
  *
  * Unless required by applicable law or agreed to in writing, software
  * distributed under the License is distributed on an "AS IS" BASIS,
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *******************************************************************************/
  
  /*
  //++
  //  Declaration of common base classes for normal equations model training.
  //--
  */
  
  #ifndef __LINEAR_MODEL_TRAIN_NORMEQ_KERNEL_H__
  #define __LINEAR_MODEL_TRAIN_NORMEQ_KERNEL_H__
  
  #include "services/env_detect.h"
  #include "data_management/data/numeric_table.h"
  #include "src/data_management/service_numeric_table.h"
  
  namespace daal
  {
  namespace algorithms
  {
  namespace linear_model
  {
  namespace normal_equations
  {
  namespace training
  {
  namespace internal
  {
  using namespace daal::services;
  using namespace daal::data_management;
  using namespace daal::internal;
  using namespace daal::services::internal;
  
  /**
   * Abstract class that defines interface for the helper function that computes the regression coefficients.
   */
  template <typename algorithmFPType, CpuType cpu>
  class KernelHelperIface
  {
  public:
      /**
       * Computes regression coefficients by solving the system of linear equations
       * \param[in] p         Size of the system of linear equations
       * \param[in] a         Matrix of size P x P with semifinished left hand side of the system
       * \param[in,out] aCopy Auxiliary matrix of size P x P with a copy of the matrix a
       * \param[in] ny        Number of right hand sides of the system
       * \param[in,out] b     Matrix of size Ny x P.
       *                      On input, the right hand sides of the system of linear equations
       *                      On output, the regression coefficients
       * \param[in] interceptFlag Flag. If true, then it is required to compute an intercept term
       * \return Status of the computations
       */
      virtual Status computeBetasImpl(DAAL_INT p, const algorithmFPType * a, algorithmFPType * aCopy, DAAL_INT ny, algorithmFPType * b,
                                      bool inteceptFlag) const = 0;
  };
  
  /**
   * Implements the common part of the regression coefficients computation from partial result
   */
  template <typename algorithmFPType, CpuType cpu>
  class FinalizeKernel
  {
      typedef WriteRows<algorithmFPType, cpu> WriteRowsType;
      typedef WriteOnlyRows<algorithmFPType, cpu> WriteOnlyRowsType;
      typedef ReadRows<algorithmFPType, cpu> ReadRowsType;
  
  public:
      /**
       * Computes regression coefficients by solving the symmetric system of linear equations
       *      - X' - matrix of size N x P' that contains input data set of size N x P
       *             and optionally a column of 1's.
       *             Column of 1's is added when it is required to compute an intercept term
       *      - P' - number of columns in X'.
       *             P' = P + 1, when it is required to compute an intercept term;
       *             P' = P, otherwise
       * \param[in]  xtx      Input matrix \f$X'^T \times X'\f$ of size P' x P'
       * \param[in]  xty      Input matrix \f$X'^T \times Y\f$ of size Ny x P'
       * \param[out] xtxFinal Resulting matrix \f$X'^T \times X'\f$ of size P' x P'
       * \param[out] xtyFinal Resulting matrix \f$X'^T \times Y\f$ of size Ny x P'
       * \param[out] beta     Matrix with regression coefficients of size Ny x (P + 1)
       * \param[in]  interceptFlag    Flag. True if intercept term is not zero, false otherwise
       * \param[in]  helper   Object that implements the differences in the regression
       *                      coefficients computation
       * \return Status of the computations
       */
      static Status compute(const NumericTable & xtx, const NumericTable & xty, NumericTable & xtxFinal, NumericTable & xtyFinal, NumericTable & beta,
                            bool interceptFlag, const KernelHelperIface<algorithmFPType, cpu> & helper);
  
      static Status copyDataToTable(const algorithmFPType * data, size_t dataSizeInBytes, NumericTable & table);
  
      /**
       * Solves the symmetric system of linear equations
       * \param[in] p             Size of the system of linear equations
       * \param[in] a             Matrix of size P x P with the left hand side of the system
       * \param[in] ny            Number of right hand sides of the system
       * \param[in,out] b         Matrix of size Ny x P.
       *                          On input, the right hand sides of the system of linear equations
       *                          On output, the regression coefficients
       * \param[in] internalError Error code that have to be returned in case incorrect parameters
       *                          are passed into lapack routines
       * \return Status of the computations
       */
      static Status solveSystem(DAAL_INT p, algorithmFPType * a, DAAL_INT ny, algorithmFPType * b, const ErrorID & internalError);
  };
  
  /**
   * Thread local storage used on the partial results update stage
   */
  template <typename algorithmFPType, CpuType cpu>
  class ThreadingTask
  {
      typedef ReadRows<algorithmFPType, cpu> ReadRowsType;
  
  public:
      DAAL_NEW_DELETE();
  
      /**
       * Creates thread local storage of the requested size
       * \param[in] nBetasIntercept   Number of colums in the partial result
       * \param[in] nResponses        Number of responses
       * \return Pointer on the thread local storage object if the object was created successfully, NULL otherwise
       */
      static ThreadingTask<algorithmFPType, cpu> * create(size_t nBetasIntercept, size_t nResponses);
      virtual ~ThreadingTask();
  
      /**
       * Updates local partial result with the new block of data
       * \param[in] startRow  Index of the starting row of the block
       * \param[in] nRows     Number of rows in the block of data
       * \param[in] xTable    Input data set of size N x P
       * \param[in] yTable    Input array of responses of size N x Ny
       * \return Status of the computations
       */
      Status update(DAAL_INT startRow, DAAL_INT nRows, const NumericTable & xTable, const NumericTable & yTable);
  
      /**
       * Reduces thread local partial results into global partial result
       * \param[out] xtx Global partial result of size P' x P'
       * \param[out] xty Global partial result of size Ny x P'
       */
      void reduce(algorithmFPType * xtx, algorithmFPType * xty);
  
  protected:
      /**
       * Construct thread local storage of the requested size
       * \param[in]  nBetasIntercept  Number of colums in the partial result
       * \param[in]  nResponses       Number of responses
       * \param[out] st               Status of the object construction
       */
      ThreadingTask(size_t nBetasIntercept, size_t nResponses, Status & st);
      algorithmFPType * _xtx;    /*!< Partial result of size P' x P' */
      algorithmFPType * _xty;    /*!< Partial result of size Ny x P' */
      ReadRowsType _xBlock;      /*!< Object that manages memory block of the input data set */
      ReadRowsType _yBlock;      /*!< Object that manages memory block of the input array of responses */
      DAAL_INT _nBetasIntercept; /*!< P' - number of columns in the partial result */
      DAAL_INT _nResponses;      /*!< Ny - number of responses */
  };
  
  /**
   * Implements the common part of the partial results update with new block of input data
   */
  template <typename algorithmFPType, CpuType cpu>
  class UpdateKernel
  {
      typedef WriteRows<algorithmFPType, cpu> WriteRowsType;
      typedef ReadRows<algorithmFPType, cpu> ReadRowsType;
      typedef ThreadingTask<algorithmFPType, cpu> ThreadingTaskType;
  
  public:
      /**
       * Updates normal equations model with the new block of data
       * \param[in]  x        Input data set of size N x P
       * \param[in]  y        Input responses of size N x Ny
       * \param[out] xtx      Matrix \f$X'^T \times X'\f$ of size P' x P'
       * \param[out] xty      Matrix \f$X'^T \times Y\f$ of size Ny x P'
       * \param[in]  initializeResult Flag. True if results initialization is required, false otherwise
       * \param[in]  interceptFlag    Flag.
       *                              - True if it is required to compute an intercept term and P' = P + 1
       *                              - False otherwis, P' = P
       * \return Status of the computations
       */
      static Status compute(const NumericTable & x, const NumericTable & y, NumericTable & xtx, NumericTable & xty, bool initializeResult,
                            bool interceptFlag);
  };
  
  /**
   * Implements the common part of computations that merges together several partial results
   */
  template <typename algorithmFPType, CpuType cpu>
  class MergeKernel
  {
      typedef WriteOnlyRows<algorithmFPType, cpu> WriteOnlyRowsType;
      typedef ReadRows<algorithmFPType, cpu> ReadRowsType;
  
  public:
      /**
       * Merges an array of partial results into one partial result
       * \param[in] n          Number of partial resuts in the input array
       * \param[in] partialxtx Array of n numeric tables of size P x P
       * \param[in] partialxtx Array of n numeric tables of size Ny x P
       * \param[out] xtx       Numeric table of size P x P
       * \param[out] xty       Numeric table of size Ny x P
       * \return Status of the computations
       */
      static Status compute(size_t n, NumericTable ** partialxtx, NumericTable ** partialxty, NumericTable & xtx, NumericTable & xty);
  
  protected:
      /**
       * Adds input numeric table to the partial result
       * \param[in]  partialTable       Numeric table with partial sums
       * \param[out] result             Resulting array with full sums
       * \param[in]  threadingCondition Flag. If true, then the operation is performed in parallel
       * \return Status of the computations
       */
      static Status merge(const NumericTable & partialTable, algorithmFPType * result, bool threadingCondition);
  };
  
  } // namespace internal
  } // namespace training
  } // namespace normal_equations
  } // namespace linear_model
  } // namespace algorithms
  } // namespace daal
  #endif
  