BNRC.h

/*
  score/BNRC.{h,c} : BNRC network score function.
  Copyright (C) 2018-2021, Yoshinori Tamada <tamada A T ytlab.jp>
  All rights reserved.
 
  THE FEATURES OF ECv CALCULATION ARE PROTECTED BY PATENT. 
  KYOTO UNIVERSITY HOLDS THE COMMERCIAL RIGHTS TO THIS SOFTWARE.
  DO NOT COPY, MODIFY OR REDISTRIBUTE WITHOUT PERMISSION.
 
  See LICENSE.txt for details of the licensing agreement.
*/
 
#ifndef __INGOR_BNRC_H
#define __INGOR_BNRC_H
 
#include "util/ytData.h"
#include "util/ytKeyValues.h"
#include "net/ytNetwork.h"
#include "net/ytNode.h"
#include "net/ytEdge.h"
 
/* In mega bytes */
//#define BNRC_DEFAULT_MAX_MEM 500
#define BNRC_DEFAULT_MAX_MEM 1000
#ifdef INGOR_PUB
#define BNRC_DEF {"BNRC",\
   BNRC_init,\
   BNRC_score,\
   BNRC_reinit,\
   Score_default_cache,\
   BNRC_edgeProp,\
   BNRC_nodeProp,\
   BNRC_finalize,\
   BNRC_partialResidual,\
   NULL,\
   BNRC_setEdgeProp,\
   BNRC_setNodeProp,\
   BNRC_status,\
   BNRC_setParam,\
  }
#else
#define BNRC_DEF {"BNRC",\
   BNRC_init,\
   BNRC_score,\
   BNRC_reinit,\
   Score_default_cache,\
   BNRC_edgeProp,\
   BNRC_nodeProp,\
   BNRC_finalize,\
   BNRC_partialResidual,\
   BNRC_calcEdgeContrib,\
   BNRC_setEdgeProp,\
   BNRC_setNodeProp,\
   BNRC_status,\
   BNRC_setParam,\
  }
#endif /* INGOR_PUB */
 
typedef struct {
  int p; /* cols (variables) of data matrices X and Y */
  int N; /* Column length (stride width) of X, Y, Xm, Ym */
  double * X; /* Column-major (n, p) input (explanatory) matrix */
  double * Y; /* Column-major (n, p) output (target) matrix */
  int * Yn; /* The number of non-missing samples in Y. */
  char * Xm; /* Column-major (n, p) input missing value matrix.
               A value of 1 represents non-missing, and 0 missing values. */
  char * Ym; /* Column-major (n, p) output matrix */
  double * xl; /* left-most minimum value vector of X */
  double * xr; /* right-most maximum value vector of X */
  double * B; /* p x (4 * n) row-major B-spline design matrix.
                 B + j * 4 * n corresponds to the design matrix for the j-th
                 variable */
  int * Bindex; /* p x n idndex vector for design matrices.
                  'Bindex + j * m' corresponds to the index vector for the
                  j-th variable */
} BNRC_Data;
 
typedef struct {
  int M; /* The number of B-splines */
  int H; /* The number of hyperparamter (length of 'beta') */
 
  double * beta; /* B-spline hyperparameters to optimize. */
  double * log_beta; /* log of 'beta' */
 
  double * K; /* M x M matrix to store K (= Dk' Dk). */
  double log_det_K_plus; /* |K|_+ : product of non zero eigenvalues. */
 
  int * pa; /* dummy parents */
 
  /* level 2 only: */ 
  /*** LambdaInv[j] points to 'h' of 'M x M' matrices. */
  /***  (LambdaInv = (B'B + n beta K)^-1 */
  double ** LambdaInv;
 
  /* level 2 & 3; log|Lambda| (= B'B + n beta K) */
  double ** logDetLambda; /* log|Lambda_jk| : p * H */
 
  /* level 3 */
  /*** LambdaInvBt = (B'B + n beta K)^-1 B' (M x n)*/
  /*** LambdaInvBt[j] points to h of (M x n) matrices. */
  double ** LambdaInvBt;
 
  int maxLoops; /* max loops for optimization */
  int stop; /* stops if parameter estimation is not converged. */
  double outer; /* outer region of the input value range */
} BNRC_Global;
 
typedef struct {
  double ** tmp_LambdaInvBt; /* Map to LambdaInvBt */
  double ** tmp_logDetLambda; /* Map to logDetLambda or logDetLambdaWork. */
  double * T1; /* Working memory 1 ; Dk : (M - 2) * M, B : (M + 1) + n * (deg + 1),
                             BtB : M * M, S : n * n */
  double * T2; /* Working memory 2 : M * M */
  int * S1; /* Integer working memory for a pivot table */
  double * dsytrf_work;
  int dsytrf_work_size;
 
  double * RS; /* n-residual vector */
  double * beta_gKg; /* mp-vector for beta gKg */
  double * B_gamma; /* B_gamma + k * n : B_jk gamma : n * mp */
  double * tmp_gamma; /* gamma : M * 1 */
  double * tmp_B_gamma; /* B_gamma : n * 1 */
  int * beta_jk; /* selection of beta[h] : max_parents */
 
  double * LambdaInvBtWork;  /* mp x H x (M x n) matrices. */
  double * logDetLambdaWork; /* mp x H */
 
  int * pa; /* copy of parents */
  int q; /* number of parents */
} BNRC_Work;
 
typedef struct {
  double * gamma; /* M vector of gamma for q_j parents */
  double * xl; /* q_j min range */
  double * xr; /* q_j max range */
  double mean;
  double var;
} BNRC_Model;
 
void * BNRC_init(ytData * data, ytKeyValues * args);
 
void BNRC_reinit(void * buff, ytData * data);
double BNRC_score(void * buff, const int j, const int * parents, const int q);
void BNRC_edgeProp(void * buff, const int j, const int * parents, const int q,
                     const int k, ytEdge * edge);
void BNRC_nodeProp(void * buff, int j, ytNode * node);
void BNRC_finalize(void * buff);
void BNRC_partialResidual(void * buff, const int j, const int * parents, const int q, 
                            FILE ** fp, double * ll);
void BNRC_setEdgeProp(void * buff, const int j, const int * parents, const int q,
                        const int k, const ytEdge * edge);
void BNRC_setNodeProp(void * buff, const int j, const ytNode * node);
void BNRC_calcRelContrib(void * buff, int j,
                           const int * parents, int q, double * ar);
#ifndef INGOR_PUB
void BNRC_calcEdgeContrib(void * buff, ytNetwork * network, int j, 
                          const int * parents, int q,
                          int mode, double * ar);
#endif /* INGOR_PUB */
void BNRC_status(void * buff);
int BNRC_debug(int argc, char * argv[]);
 
 
void BNRC_prepare_beta(double hyper_bg, double hyper_inc, int hyper_num,
                       double * beta, double * log_beta, int verbose);
 
size_t BNRC_Data_alloc(BNRC_Data * D, int p, int n, 
                       int allocXY, int dynamic, int level, int dry);
void BNRC_Data_finalize(BNRC_Data * D, int allocXY);
 
size_t BNRC_Global_alloc(BNRC_Global * G, int p, int n, int mp, int level,
                         int dry);
void BNRC_Global_finalize(BNRC_Global * G);
 
size_t BNRC_Work_alloc(BNRC_Work * W, int p, int n, int mp, int H,
                       int M, int level, int T1_size, int dry);
void BNRC_Work_finalize(BNRC_Work * W);
 
size_t BNRC_Model_alloc(BNRC_Model * P, int mp, int M, int dry);
void BNRC_Model_finalize(BNRC_Model * P);
 
void BNRC_range(BNRC_Data * D, 
                const ytDoubleArray * xlar, const ytDoubleArray * xrar,
                const char * type, double outer);
void BNRC_proc_mv(BNRC_Data * D, int level);
void BNRC_calc_B(BNRC_Data * D, int ofs, const char * type,
                 int n, int M, BNRC_Work * W);
void BNRC_calc_Lambda(BNRC_Data * D, BNRC_Global * G, int p, int n,
                      const char * type, BNRC_Work * W, int level);
 
double BNRC_calc(BNRC_Data * D, int ofs, int j, const int * parents, int q,
                 int n, int nn, int level, BNRC_Work * W, BNRC_Model * P, BNRC_Global * G);
 
size_t BNRC_get_T1_size(int n, int M);
 
void BNRC_setParam(void * buff, const char * key, ytObject * value);
 
#endif /* __INGOR_BNRC_H */