[Rd] want to call R from aplatform written i strict ANSI C

[Ralf Östermark]

Hi,
I would like to connect my GHA-system to R and to deliver the
library to users of R.

could you give me a comprehensive example containing the step by step
procedure to:
1. create a callable library with R binaries (not having a standalone
version of R with its own main()). If possible, I would like to run the
system on unix,alpha,linux computers as well as on IBM's parallel supercomputers.
2. code the interface between a main-program written in C and a
user-defined script written in R. Consider the following case:
The algorithm - written in strict ANSI C (explained below) - needs to call a user-written
R-program occasionally in order to boost computations, say, in an
ill-conditioned computational problem where classical optimization tools
can be used only locally. For this, the C-program needs somehow to
transfer data, an initialized vector of parameters and the name of the
user-defined R program to some interface function. The R program in turn
can use all the resources available in the R environment. The results of
the computations (e.g., the updated parameter vector, matrices of
forecasts and forecast errors) need to be returned to the C-program from
R. The corresponding allocations need to be done before R is invoked and freed afterwards. 

Using R would increase the flexibility of the C-platform considerably
and hopefully bring new features to R.

I would be very happy if I could do this both sequentially on
alpha/unix/linux platforms and in parallel on supercomputers.

I have developed an algorithmic platform (GHA = Genetic Hybrid
Algorithm) over the past years (referee articles are listed under
www.abo.fi/~rosterma). The system is compiled and runs without warnings
under maximum warning level in sequential mode on alpha, unix and linux
mainframes and in parallel mode on IBMs supercomputers and Cray T3E (the
latter is a more or less outdated supercomputer today).

The system was built in strict ANSI C and it can be connected to
computational algorithms through an accelerator function. We can solve
problems using pure classical optimization tools on the one extreme
(global, constrained (linear,nonlinear, milp,minlp etc)) and pure
genetic search on the other. We can also use hybrid approaches where
various classical/artificial intelligence approaches are combined.

GHA may be used as a subsystem to other programs
or as the main system using other systems as subsystems (libraries). As
I am not an expert in R, in the first stage I would like to
invoke some simple user-written R routines from within GHA, using R as a
support library for GHA. If successful, I would be in a position to
write instructions for other users and to submit my package to the R
community. For the moment I have no external users of my algorithm. I
have studied web-material concerning R extensions. The file R-exts.pdf
(chapter 5), for example, describes various entry points to R from C.
Also chapter 7 discusses this issue.

However, I have not found a comprehensive example that would give me a
step by step procedure. Consider the following example where I invoke my
GHA-system as a standalone program to solve a minlp problem:

#include <stdio.h>  
#include <stdlib.h> 
#include <string.h> 
#include <time.h>   
#include <ctype.h>  
#include <math.h>   
#include <float.h>
#include "supergha.h"
#include "rpa_proj.h"
MINLP_ptr *SHAREX_PROB;
extern void ghasystem(int argc,char *argv[]);
int main(int argc,char **argv) {
 int LOAD_SYSTEM = TRUE;
 INITIALIZE_RAN1
 SHAREX_PROB  = get_data(argc,argv,LOAD_SYSTEM);
 ghasystem(argc,argv);
 exit(0);
}   /* end of main() */

The function ghasystem() controls all GHA-processes. The key user-level
functions are defined in project.c. In the below example I have used
gha() to solve a minlp-problem using classical optimization tools only.

/* project.c contains the following problem specific functions:
   (1) evaluator()
   (2) accelerator()
   (3) gradient()
   (4) hessian()
   (5) pre_processor()
   (6) post_processor()
*/
/* Global Definitions */
#include "project.h"
#include "rpa_proj.h"
#include "rpa_proj.c"
extern RPA_ptr   *SHAREX_ROOT;
extern MINLP_ptr *SHAREX_PROB;
IVECTOR SUPER_FLAG        = &(INT_STATUS[16]);
IVECTOR ALLOCATION_SWITCH = &(INT_PROTOCOL[30]);
extern void compare_NVAR(int n1,int n2,int n3,char *stage);
/* end of global definitions */
void evaluator(DVECTOR w,double Penalty,double *gf,double *F,double *Dev) {
 int i,n,n_i,n_x,n_d,m_f,ACCELERATE;
 ACCELERATE=INT_STATUS[0];
 REAL_STATUS[4] += 1.0;    /* number of evaluator() calls */
 n   = min(SHAREX_PROB->n,*(TASK.NVAR));
 n_i = min(SHAREX_PROB->n_i,*(TASK.NVAR));
 n_x = min(SHAREX_PROB->n_x,*(TASK.NVAR));
 n_d = min(SHAREX_PROB->n_d,*(TASK.NVAR));
 m_f = SHAREX_PROB->m_f;
 /* NOTE: w[0] is F */
 *F   = 0.0;
 *Dev = 0.0;
 if(m_f EQ 0) {for(i=0;i<n;i++) *F   += w[i]*SHAREX_PROB->c[i];}
 else {
  for(i=0;i<m_f;i++) *F +=
f_function(SHAREX_PROB,SHAREX_PROB->x,SHAREX_PROB->du
al_x,i);
 }
 *Dev = SHAREX_PROB->Dev;
 *gf  = *F+Penalty*(*Dev);
}  /* end of evaluator() */

void accelerator(int FIX_Flip,int ROW,double Penalty,DVECTOR LOWER,
              DVECTOR UPPER,DVECTOR w_IN,DVECTOR w_OUT) {
 double MILP_CALLS;
 if(FIX_Flip > 0) return;
 SHAREX_PROB->INITIALIZED = FALSE;
 SHAREX_PROB->TREE_COUNTER = 0.0;
 SHAREX_PROB->ACTIVE_TREE  = 0.0;
 REAL_STATUS[5] += 1.0;    /* number of accelerator() calls */

 MILP_CALLS = milp_caller(SHAREX_PROB);

  /*
   NOTE: if(ni_g>0 OR nx_g>0) then w must be restructured.
         if(n_change>0 and we want to change TASK.NVAR in full GHA-search, 
         then POP,HISTORY,w_in,START,IMPROVED,BEST must be restructured.
         GHA_SYS.SUMDIM and GHA_SYS.MAXDIM must be adjusted as well
         as TASK.DERIVABLES etc.
  */
 compare_NVAR(SHAREX_PROB->n,*(TASK.NVAR),*(TASK.NVAR),"accelerator() 1");
 memcpy(w_OUT,SHAREX_PROB->x,*(TASK.NVAR)*sizeof(double));
}      /* end of accelerator() */

void gradient(DVECTOR w,int n_w,double Penalty,DVECTOR d,int n_d,IVECTOR
pos) {
 /* NOTE: FIXED_h=1,VARIABLE_h=2 */
 int h_TYPE=2;
 num_gradient(h_TYPE,w,n_w,Penalty,d,n_d,pos);
}   /* end of gradient */

void hessian(DVECTOR w,int n_w,double Penalty,DMATRIX G,int n_G,IVECTOR
pos) {
 /* NOTE: FIXED_h=1,VARIABLE_h=2 */
 int h_TYPE=2;
 num_hessian(h_TYPE,w,n_w,Penalty,G,n_G,pos);
}        /* end of hessian() */

void pre_processor (struct SOLUTION *START) {
 int PHASE = INT_STATUS[17];
 if((PHASE EQ 1) AND (*(TASK.ACTIVE_SYSTEM) EQ 0) AND
(*ALLOCATION_SWITCH EQ FALSE)) {
  *ALLOCATION_SWITCH = TRUE;
  initlp_();
 }   /* end of if() */
}    /* end of pre_processor() */

void post_processor(struct SOLUTION *BEST) {
 double Penalty;
 int PHASE = INT_STATUS[17];
 Penalty = REAL_STATUS[0];
 if((*SUPER_FLAG EQ TRUE) AND (PHASE EQ GHA_SYS.SYSTEM_CALLS)){
  deallocate_MINLP(SHAREX_PROB);
 }   /* end of if() */
}    /* end of post_processor() */

The key connection between R and GHA - from the viewpoint of GHA - is
the accelerator()-function. This function is activated by GHA in
critical stages of the main iterations (a flag indicates where precisely
the function is invoked). If I could call a user-written R-script from
the accelerator() at a critical stage, then I could - for example -
split the solution process into the following sequence:

1. given the current local environment of the problem as defined by
gha(), call the R script
2. the R script returns a local solution to the problem.
3. the system returns to step 1 until convergence.

Since gha works both sequentially and in parallel, I can boost the
search process in difficult problems using the parallel programming
facilities of gha (tested on the parallel supercomputers of the Centre
of Scientific Comuting in Helsinki).

I am grateful for any help on this matter.
regards
Ralf Östermark
Professor of Accounting and Optimization Systems