R/hgm.c2wishart.R

In hgm: Holonomic Gradient Method and Gradient Descent

# $OpenXM: OpenXM/src/R/r-packages/hgm/R/hgm.c2wishart.R,v 1.5 2017/03/22 00:42:57 takayama Exp $
"hgm.tk.p2wishart" <-
function(m=3,n1=5,n2=10,beta=c(1,2,4),q0=0.3,approxdeg=-1,h=0.001,dp=-1,q=4,
         mode=c(1,1,0),method="a-rk4",err=c(-1.0,-1.0),
         automatic=1,assigned_series_error=0.00001,verbose=0,autoplot=0) { 
  x<-q; x0<-q0;
  nstrategy<-0;
  mm<-charmatch(method,c("rk4","a-rk4","a-rk4-m"));
  if (!is.na(mm)) nstrategy<- (mm-1);

  if ((m>=200) || (m<=0)) stop("Invalid size of m."); #200 is M_m_MAX in jack-n.c

  if (!is.loaded("hgm")) library.dynam("hgm",package="hgm",lib.loc=NULL);

  .C("Rmh_set_strategy",as.integer(nstrategy),as.double(err),retv=double(1),
     PACKAGE="hgm")$retv ;

  if (m<1) m=1;
  rank <- 2^m;
  rsize <- rank+1;
  if (autoplot==1) {
     dp<-floor(q/(h*100));
     mode<-c(1,1,(rank+1)*floor(q/(h*dp)+1));
  }
  if (mode[3] > 0) rsize <- rsize+mode[3]; 
  if (approxdeg < 0) approxdeg <- 6;
##argchecks
  if (class(beta)[1] != "numeric") stop("beta must be a vector.");
  if (length(beta) != m) stop("The length beta must be equal to m.");
  for (i in 1:m) {
    if (beta[i] <= 0) stop("beta[i] must be positive.");
  }
  if (m != 1) {
   for (i in 1:(m-1)) {
    for (j in  (i+1):m) {
       if (beta[i] == beta[j]) stop("beta's must be different.");
    }
   } 
  }
  if (class(mode)[1] != "numeric") stop("mode must be a vector of length 3.");
  if (length(mode) != 3) stop("mode must be a vector of length 3.");
  if ((n1 < m) || (n2 < m)) stop("n1 and n2 must be >=m.");
##end of argchecks

  pp<-2; qq<-1;
  a<-c((m+1)/2,(n1+n2)/2);
  b<-c((n1+m+1)/2);
  ef_type<-2;

  ans <-
  .C("Rmh_pFq_gen",as.integer(m),
     as.integer(pp), as.double(a),
     as.integer(qq), as.double(b),
     as.integer(ef_type),
     as.double(beta),as.double(x0),
     as.integer(approxdeg),
     as.double(h),
     as.integer(dp),as.double(x),
     as.integer(mode),as.integer(rank),
     as.integer(automatic),as.double(assigned_series_error),as.integer(verbose),
     retv=double(rsize),PACKAGE="hgm")$retv ;
  if (autoplot == 0) {
    return(ans);
  }else {
    return(matrix(ans,ncol=rank+1,byrow=1));
  }
}