BNRC.h

/*
  score/BNRC.{h,c} : BNRC network score function.
  Copyright (C) 2018-2021, Yoshinori Tamada <tamada A T ytlab.jp>
  All rights reserved.
 
  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that the following conditions
  are met:
 
  * Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.
 
  * Redistributions in binary form must reproduce the above copyright
    notice, this list of conditions and the following disclaimer in 
    the documentation and/or other materials provided with the 
    distribution.
 
  * Neither the name of Hirosaki University, Kyoto University, The 
    University of Tokyo, nor the names of its contributors may be 
    used to endorse or promote products derived from this software 
    without specific prior written permission.
 
  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
  COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  POSSIBILITY OF SUCH DAMAGE.
*/
 
#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
#define BNRC_DEF {"BNRC",\
   BNRC_init,\
   BNRC_score,\
   BNRC_reinit,\
   Score_default_cache,\
   BNRC_edgeProp,\
   BNRC_nodeProp,\
   BNRC_finalize,\
   BNRC_partialResidual,\
   BNRC_setEdgeProp,\
   BNRC_setNodeProp,\
   BNRC_status,\
   BNRC_setParam,\
  }
 
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' */
  /*** 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 */
  double * T2; /* Working memory 2 */
  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_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 */