R/FedorovWynnAlgorithm.R

In PFIM: Population Fisher Information Matrix

#' Class "FedorovWynnAlgorithm"
#'
#' @description Class \code{FedorovWynnAlgorithm} represents an initial variable for ODE model.
#'
#' @name FedorovWynnAlgorithm-class
#' @aliases FedorovWynnAlgorithm
#' @docType class
#' @include OptimizationAlgorithm.R
#' @include GenericMethods.R
#' @include Design.R
#' @export
#'
#' @section Objects from the class \code{FedorovWynnAlgorithm}:
#' Objects form the class \code{FedorovWynnAlgorithm} can be created by calls of the form \code{FedorovWynnAlgorithm(...)}
#' where (...) are the parameters for the \code{FedorovWynnAlgorithm} objects.
#'
#'@section Slots for \code{FedorovWynnAlgorithm} objects:
#' \describe{
#' \item{\code{elementaryProtocols}:}{A list of vector for the initial elementary protocols.}
#' \item{\code{numberOfSubjects}:}{A vector for the number of subjects.}
#' \item{\code{proportionsOfSubjects}:}{A vector for the number of subjects.}
#' \item{\code{OptimalDesign}:}{A object Design giving the optimal Design.}
#' \item{\code{showProcess}:}{A boolean to show the process or not.}
#' \item{\code{FisherMatrix}:}{A vector giving the Fisher Information }
#' \item{\code{optimalFrequencies}:}{A vector of the optimal frequencies.}
#' \item{\code{optimalSamplingTimes}:}{A list of vectors for the optimal sampling times.}
#' \item{\code{optimalDoses}:}{A vector for the optimal doses.}
#' }

FedorovWynnAlgorithm = setClass(
  Class = "FedorovWynnAlgorithm",
  contains = "OptimizationAlgorithm",
  representation = representation(
    elementaryProtocols = "list",
    numberOfSubjects = "vector",
    proportionsOfSubjects = "vector",
    showProcess = "logical",
    optimalDesign = "Design",
    FisherMatrix = "vector",
    optimalFrequencies = "vector",
    optimalSamplingTimes = "list",
    optimalDoses = "vector")
)

#' initialize
#' @param .Object .Object
#' @param elementaryProtocols elementaryProtocols
#' @param numberOfSubjects numberOfSubjects
#' @param proportionsOfSubjects proportionsOfSubjects
#' @param showProcess showProcess
#' @return FedorovWynnAlgorithm
#' @export

setMethod(f="initialize",
          signature="FedorovWynnAlgorithm",
          definition= function (.Object, elementaryProtocols, numberOfSubjects, proportionsOfSubjects, showProcess)
          {
            if(!missing(elementaryProtocols))
              .Object@elementaryProtocols=elementaryProtocols

            if(!missing(numberOfSubjects))
              .Object@numberOfSubjects=numberOfSubjects

            if(!missing(proportionsOfSubjects))
              .Object@proportionsOfSubjects=proportionsOfSubjects

            if(!missing(showProcess))
              .Object@showProcess=showProcess

            validObject(.Object)
            return (.Object )
          }
)

##########################################################################################################
#' Fedorov-Wynn algorithm in Rcpp.
#'
#' @name FedorovWynnAlgorithm_Rcpp
#' @description Run the FedorovWynnAlgorithm in Rcpp
#' @param protocols_input parameter protocols_input
#' @param ndimen_input parameter ndimen_input
#' @param nbprot_input parameter nbprot_input
#' @param numprot_input parameter numprot_input
#' @param freq_input parameter freq_input
#' @param nbdata_input parameter nbdata_input
#' @param vectps_input parameter vectps_input
#' @param fisher_input parameter fisher_input
#' @param nok_input parameter nok_input
#' @param protdep_input parameter protdep_input
#' @param freqdep_input parameter freqdep_input
#' @return A list giving the results of the outputs of the FedorovWynn algorithm.
#' @export

FedorovWynnAlgorithm_Rcpp = function( protocols_input,  ndimen_input, nbprot_input,
                                      numprot_input, freq_input, nbdata_input,
                                      vectps_input, fisher_input, nok_input,
                                      protdep_input, freqdep_input ){
  incltxtFedorovWynnAlgorithm = '

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h> 	/* Mathematical functions */
#include <time.h>	/* Function time used to initialise the random number generator */
#include <float.h>	/* Implementation related constants */
#include <signal.h>	/* Signal handling used to detect arithmetic errors */
#include <Rcpp.h> /* Rcpp */

using namespace Rcpp;

#define FTOL 1e-6 /* Minimal change in criterion during optimisation step*/
#define FREQMIN 1e-8 /* Minimal frequency for a protocol to be kept */
#define EPSILON DBL_EPSILON
#define CRI_MAX 1e30 /* maximum value for criterium, used to test correct comp*/
#define NRANSI
#define tol 0.0001
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
#define TINY 1.0e-20
#define SMALL 1.0e-7
#define IA 16807  	/*definition des constantes pour ran1 */
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
#define FREE_ARG char*

typedef struct PROTOC /* Structure defi
ning an individual or elementary protocol */
{
	int	ntps; /* Number of sampling times */
	double	*tps; /* Vector of the sampling times */
	int	ndim; /* Number of random effects (=> size of matrices)*/
	double	*fisher; /* Fisher information matrix of the individual protocol */
} PROTOC;

typedef struct POPROT /* Structure defining a population protocol */
{
	int 	maxnp; /* Max number of individual protocols, needed for memory allocation*/
	int	np; /* Number of individual protocols WARNING np<=maxnp !*/
	int	ndim; /* Number of random effects (=> size of matrices)*/
	int	*num; /* In this program we keep trace of the index of individual protocols*/
	PROTOC	*pind; /* In NPML we would only need the individual protocols themselves*/
	double	*freq; /* Frequency of each individual protocols */
	double 	*fisher; /* Fisher information matrix of population protocol */
	double 	*finv; /* Inverse of the Fisher information matrix */
	double	det; /* Determinant of the matrix */
} POPROT;

typedef struct matrix
{
	int 	nrow;
	int	ncol;
	double	**m;
} matrix;


/******************************************************************************
   Construction and destruction of structures
*******************************************************************************/

void PROTOC_alloc( PROTOC *p, int ntps, int ndim)
{
	p->ntps = ntps;
	p->tps = (double *)calloc( p->ntps, sizeof(double));
	assert(p->tps!=NULL);
	p->ndim=ndim;
	p->fisher = (double *)calloc( p->ndim*(p->ndim+1)/2, sizeof(double));
	assert(p->fisher!=NULL);
	return;
}

POPROT *POPROT_alloc(int np,int ndim,int maxnp)
{	POPROT *p;
	p = (POPROT *)malloc(sizeof(POPROT));
	assert(p!=NULL);

	p->maxnp=maxnp;
	p->np = np;
	p->ndim=ndim;
/* Keep trace of the index of the individual protocols */
	p->num = (int *)calloc( p->maxnp, sizeof(int));
	assert(p->num!=NULL);
	p->pind = (PROTOC *)calloc( p->maxnp, sizeof(PROTOC));
	assert(p->pind!=NULL);
	p->freq = (double *)calloc( p->maxnp, sizeof(double));
	assert(p->freq!=NULL);
	p->fisher = (double *)calloc( p->ndim*(p->ndim+1)/2, sizeof(double));
	assert(p->fisher!=NULL);
	p->finv = (double *)calloc( p->ndim*(p->ndim+1)/2, sizeof(double));
	assert(p->finv!=NULL);
	p->det=1.0;
	return p;
}

void PROTOC_copy(PROTOC *p,PROTOC *p1)
{
	int i,ndat;

	for(i=0;i<p->ntps;i++) {p1->tps[i]=p->tps[i];}
	ndat=(int)(p->ndim*(p->ndim+1))/2;
	for(i=0;i<ndat;i++) {p1->fisher[i]=p->fisher[i];

		}
	return;
}

matrix *matrix_create(int nrow, int ncol)
{	int i,j;
	matrix *mat;
	mat = (matrix *)malloc(sizeof(matrix));
	assert(mat!=NULL);
	mat->nrow=nrow;
	mat->ncol=ncol;
	mat->m = (double **)calloc( nrow, sizeof(double *));
	assert(mat->m!=NULL);
	for( i=0 ; i<nrow ; i++ )
	{
	mat->m[i] = (double *)calloc( ncol, sizeof(double));
	assert(mat->m[i]!=NULL);
	for( j=0 ; j<ncol ; j++ )
		mat->m[i][j]=0.0;
	}
	return mat;

}

void PROTOC_print(PROTOC *p)
{	int i,icas=0,j;
	fprintf(stderr,"Temps du protocole :");
	for (i=0;i<p->ntps;i++) {}
	for (i=0;i<p->ndim;i++) {
		for (j=0;j<=i;j++) {
			icas++;}
	}
	return;
}

void POPROT_print(POPROT *pop,int nofish)
{	int i,icas=0,j,k;

	for(i=0;i<pop->np;i++)
	{
		for (j=0;j<pop->pind[i].ntps;j++)
		{
		}

		if(nofish==1) {

		icas=0;
		for (j=0;j<pop->pind[i].ndim;j++)
		{
			for (k=0;k<=j;k++)
			{

				icas++;
			}

		}

		}
	}
	return;
}

void POPROT_printfisher(POPROT *pop)
{	int icas,j,k;


	icas=0;
	for (j=0;j<pop->ndim;j++)
	{
		for (k=0;k<=j;k++)
		{

		icas++;
		}

	}

	icas=0;
	for (j=0;j<pop->ndim;j++)
	{
		for (k=0;k<=j;k++)
		{

		icas++;
		}

	}
	return;
}

void POPROT_calculefisher(POPROT *pop)
{	int i,iprot;

	for(i=0;i<pop->ndim*(pop->ndim+1)/2;i++)
	{
		pop->fisher[i]=0.0;
		for(iprot=0;iprot<pop->np;iprot++)
		{
	pop->fisher[i]=pop->fisher[i]+pop->pind[iprot].fisher[i]*pop->freq[iprot];
		}
	}
}

void matrix_print(matrix *mat)
{	int i,j;

	for (i=0;i<mat->nrow;i++) {
		for (j=0;j<mat->ncol;j++) {
}
}
	return;
}

void PROTOC_destroy(PROTOC *p)
{
	free(p->fisher);
	free(p->tps);
}

void POPROT_destroy(POPROT *p)
{
	free(p->finv);
	free(p->fisher);
	free(p->freq);
	free(p->pind);
	free(p->num);
}

void matrix_destroy(matrix *mat)
{
	int i;
	for( i=0 ; i<mat->nrow ; i++ )
	{
	free(mat->m[i]);
	}
	free(mat->m);
}

/******************************************************************************
   Fonctions calculatoires
*******************************************************************************/

double matrace(int iprot,PROTOC *prot,POPROT *pop);
double lik(POPROT *pop);
void lubksb(double **a, int n, int *indx, double b[]);
int ludcmp(double **a, int n, int *indx, double *d);

/* Random number generator */
double randu(double x);
double gauss1(long *idum);
float ran1(long *idum);

/******************************************************************************
   Defining global variables (see variables_npml.txt for details)
*******************************************************************************/

volatile sig_atomic_t ierr=0;	/* Parameter signaling a SIGFPE error */
				/* to be used with a signal handler */
long seed=-10;     	/* graine du generateur aleatoire */

/******************************************************************************

/******************************************************************************
   Initialisation et ajout de protocole ?l?mentaire
*******************************************************************************/

POPROT *initprot(PROTOC *allprot, IntegerVector protdep, NumericVector freqdep)

{
	int i,ij,np,nmax;

	POPROT *mypop;

	np=protdep[0];

	nmax=allprot[0].ndim*(allprot[0].ndim+1)/2+1;
	mypop=POPROT_alloc(np,allprot[0].ndim,nmax);

	for(i=0;i<np;i++) {
		mypop->freq[i]=freqdep[i];
		mypop->num[i]=protdep[i+1]-1;
		ij=mypop->num[i];
		PROTOC_alloc(&mypop->pind[i],allprot[ij].ntps,allprot[ij].ndim);
		PROTOC_copy(&allprot[ij],&mypop->pind[i]);
	}

return mypop;

}

/******************************************************************************
   Optimisation de protocole
*******************************************************************************/

void tassement(POPROT *pop)
{	int i,idec=0,np;
	np=pop->np;
	for(i=0;i<pop->np;i++)
	{
		if(pop->freq[i]<FREQMIN)
		{
			idec++;
			np--;
		}
		else
		{
			if (idec>0)
			{
				pop->freq[i-idec]=pop->freq[i];
				pop->num[i-idec]=pop->num[i];

				PROTOC_destroy(&pop->pind[i-idec]);
	PROTOC_alloc(&pop->pind[i-idec],pop->pind[i].ntps,pop->pind[i].ndim);
				PROTOC_copy(&pop->pind[i],&pop->pind[i-idec]);

			}
		}
	}

	for(i=np;i<pop->np;i++) PROTOC_destroy(&pop->pind[i]);
	pop->np=np;
	return;
}




int ajout(PROTOC *allprot,POPROT *pop,int nprot)
{	double cri,xmul,xtest,ifin=0,xdim,rcalc;
	int iprot,j,deja,qajout=-1;

	cri=pop->det;

	xtest=(double)(pop->ndim);
	xdim=xtest;

	for(iprot=0;iprot<nprot;iprot++)
	{


		deja=0;
		for(j=0;j<pop->np;j++)
		{
			if (iprot==pop->num[j]) deja=1;
		}
		if(fabs(deja)<EPSILON)
		{
			xmul=matrace(iprot,allprot,pop);
			if(xmul>=xtest)
			{
				xtest=xmul;
				qajout=iprot;
			}
		}
	}
	if(fabs(xtest-xdim)<EPSILON)
	{
		ifin=1;
	}
	else
	{
		rcalc=(xtest-xdim)/(xdim*(xtest-1));

		PROTOC_alloc(&pop->pind[pop->np],allprot[qajout].ntps,allprot[qajout].ndim);
		PROTOC_copy(&allprot[qajout],&pop->pind[pop->np]);
		pop->num[pop->np]=qajout;
		pop->freq[pop->np]=rcalc;


		for(j=0;j<pop->np;j++) pop->freq[j]=pop->freq[j]*(1-rcalc);
		pop->np++;
	}

	cri=lik(pop);
	return ifin;
}

int takeout_k(POPROT *pop,PROTOC *allprot,int *nnul)
{
	int ir=-1,j;
	double xmul;
	xmul=lik(pop); /*update Mf-1 */
	for(j=0;j<pop->np;j++)
	{
		xmul=0;
		if(pop->freq[j]<=FREQMIN)
		{
			xmul=matrace(pop->num[j],allprot,pop);
			if(xmul>(double)pop->ndim)
			{
				ir=j;
				*nnul=*nnul-1;
				return ir;
			}
		}
	}
	return ir;
}

int project_grad(POPROT *pop, PROTOC *allprot, double *gal,int ir,double
		*vmgal,int nnul)
{
/* Computes the projection of the gradient of ln(V)=ln(det(H(pop))) on the
surface sum(frequencies)=1 and freq[j]=0 for all j such as pop->freq[j]<FREQMIN
gal[j] is the jth component of the projected gradient (the direction)
gal[j]=0 if pop->freq[j]<FREQMIN and j<>ir
ir is -1 or the index removed from the active set at the previous step
vmgal is the largest abs(gal[j])
returns in, the last index j such as gal[j]<0
*/
	int j,in=-1;
	double sg,sgal;

	sg=lik(pop); /* Check that Mf-1 has been computed */
	sg=0;
	for(j=0;j<pop->np;j++)
	{
		gal[j]=0;
		if((pop->freq[j]>=FREQMIN) | (j==ir))
			gal[j]=matrace(pop->num[j],allprot,pop);
		sg=sg+gal[j];
	}
	sgal=0;*vmgal=0;
	for(j=0;j<pop->np;j++)
	{
		if((pop->freq[j]>=FREQMIN) | (j==ir))
		{
			gal[j]=gal[j]-sg/(double)(pop->np-nnul);
			if(fabs(gal[j])>*vmgal) *vmgal=(double)fabs(gal[j]);
			if(gal[j]<0) in=j;
			sgal=sgal+gal[j];
		}
	}
	return in;
}


int optim_lambda(POPROT *pop, double *gal, double *alkeep, double *crit,
	double dmax, double crit2, double vmgal)
{
/* dmax is too large a move along the direction
=> move along the direction by dmax/2**k until
-the move is too small : no new frequencies have been found (return 1)
-a better likelihood is found => update pop with the new frequencies
	-if the criteria has changed by more than FTOL, return 0
	-else return 1; the frequencies correspond to the last (smallest) move
*/
	int j;
	double ro=1,crit1;
	do
	{
		ro=ro*0.5;
		for(j=0;j<pop->np;j++) pop->freq[j]=alkeep[j]+ro*dmax*gal[j];
		crit1=lik(pop); /* Updates Mf-1 for new frequencies */
		if(crit1>*crit)
		{
			for(j=0;j<pop->np;j++) alkeep[j]=pop->freq[j];
			*crit=crit1;
			if(fabs((crit2-*crit)/crit2)>FTOL) {return 0;}
			else {return 1;}
		}
	} while((dmax*ro*vmgal)>1e-4);
	for(j=0;j<pop->np;j++) pop->freq[j]=alkeep[j];
	return 1;
}



int doptimal(POPROT *pop,PROTOC *allprot)
{
	int nnul=0, ir=-1,in=-1,nit=0;
	int j,imax,iopt;
	double crit,crit1,crit2,vmgal,dmax;
	double *gal,*alkeep;
	double sal,sal1;

	gal=(double *)calloc(pop->np,sizeof(double));
	alkeep=(double *)calloc(pop->np,sizeof(double));
	crit=lik(pop); /*Calcul de Mf-1 */

	for(nit=0;nit<500;nit++)
	{
/* alkeep stores the initial frequencies of the elementary protocols
*/
		for(j=0;j<pop->np;j++)
			alkeep[j]=pop->freq[j];
		crit2=crit;
		if(ir==(-1)) /* First time through */
		{
			nnul=0; /* Computing the number of nul frequencies */
			for(j=0;j<pop->np;j++)
			{
				if(pop->freq[j]<FREQMIN) nnul++;
			}
		}
/* Compute the projected gradient (direction)*/
		in=project_grad(pop,allprot,gal,ir,&vmgal,nnul);
		if((in<0) | (gal[in]==0))
		{
			free(gal);
			free(alkeep);
			return nnul;
		}
/* Compute the optimal change along the gradient dmax (lambda)*/
		dmax=-pop->freq[in]/gal[in];
		imax=in;
		if((ir!=(-1)) & (gal[ir]<=0))
		{
			free(gal);
			free(alkeep);
			return nnul;
		}
		ir=-1;
		for(j=0;j<pop->np;j++)
		{
			if((gal[j]<0) & (fabs(pop->freq[j]/gal[j])<dmax))
			{
				dmax=-pop->freq[j]/gal[j];
				imax=j;
			}
		}
/* Updates the frequencies as beta*(t,j)=beta(t,j)+dmax*gal(j)
Compute the likelihood for the updated vector of frequencies
*/
		sal=0;
		for(j=0;j<pop->np;j++)
		{
			pop->freq[j]=pop->freq[j]+dmax*gal[j];
			sal=sal+pop->freq[j];
		}
		sal1=sal-pop->freq[imax];
		pop->freq[imax]=0;
		crit1=lik(pop);
/* If better likelihood, update alkeep (pop->freq has been updated already)
*/
		if(crit1>crit)
		{
			for(j=0;j<pop->np;j++) alkeep[j]=pop->freq[j];
			crit=crit1;
			nnul++;
		}
		else
		{
			for(j=0;j<pop->np;j++) pop->freq[j]=alkeep[j];
/* If not, try dmax*ro for ro=1/2**k until dmax*ro too small
*/
			iopt=optim_lambda(pop,gal,alkeep,&crit,dmax,crit2,vmgal);
			if(iopt==1)
			{
				ir=takeout_k(pop,allprot,&nnul);
/* ir=-1 if for all j, matrace(q[j]) <= ndim*(ndim+1)/2, ie optimal protocol
=> in that case exit the subroutine
*/
				if(ir==(-1))
				{
					free(gal);
					free(alkeep);
					return nnul;
				}
			}
		}
	}
	free(gal);
	free(alkeep);
	return -10;
}


void index_limit(POPROT *pop)
{
	return;
}

// ******************************************************************************************************************

int main_opt(PROTOC *allprot, POPROT *pop, IntegerVector nprot)
{
  int nit=0,ifin=0,nstop,ntest;
	double cri;

	cri=lik(pop);

	for(nit=0;nit<100;nit++)
	{

  nstop=0;
do
{
  ifin=ajout(allprot,pop,nprot[0]);
  nstop++;
} while ((cri<(-1e10))& (nstop<100));

if((nstop>99)& (cri<(-1e10)))

{
  return -10;
}
if(ifin==1) break; /* Si ifin=1, on sort de la boucle sur nit */
  ntest=doptimal(pop,allprot);
  if(ntest<0)
    return -100;
  tassement(pop);

  if(pop->np>(pop->ndim*(pop->ndim+1)/2))
  {
    index_limit(pop);
    return -1;
  }
  cri=lik(pop); /* Updates matrices */
    }
return nit;
}

// ******************************************************************************************************************

/******************************************************************************
  Fonctions calculatoires
	matrace:
	lik:determinant de hh et inverse de hh
*******************************************************************************/

double matrace(int iprot,PROTOC *prot,POPROT *pop)
{/* Computes tr(Mf(iprot)*Mf-1(poprot))
WARNING pop->finv must be computed beforehand (with a call to lik)*/
	int i,j,ij;
	double xcal=0;

	for(i=0;i<(prot->ndim);i++)
	{
		for(j=0;j<i;j++)
		{

			ij=i*(i+1)/2+j;
			xcal=xcal+2*(prot[iprot].fisher[ij])*(pop->finv[ij]);
			}
		ij=i*(i+3)/2;
		xcal=xcal+(prot[iprot].fisher[ij])*(pop->finv[ij]);
		}
	return xcal;
}

double lik(POPROT *pop)
{/* Function to compute the inverse of the Fisher information matrix of the
population protocol pop. Returns the determinant of the matrix, also stored in
pop->det. The Fisher  information matrix of the population protocol is stored
in pop->fisher and the inverse is stored in pop->finv
*/
	int ndim=pop->ndim;

	int ncase=(int)(ndim*(ndim+1)/2);

	int i,j,jj,ifail;
	int *indx;
	double *col;
	matrix *xa;
	double cri;

	indx=(int *)calloc(ndim,sizeof(int));
	col=(double *)calloc(ndim,sizeof(double));
	xa=matrix_create(ndim,ndim);


	for(i=0;i<ncase;i++)
	{
		pop->fisher[i]=0;
		for(j=0;j<pop->np;j++)
		{
pop->fisher[i]=pop->fisher[i]+pop->freq[j]*pop->pind[j].fisher[i];
		}
	}
	jj=0;
	for(i=0;i<ndim;i++)
	{
		for(j=0;j<=i;j++)
		{
			xa->m[i][j]=pop->fisher[jj];
			xa->m[j][i]=pop->fisher[jj];

			jj++;
		}
	}

	ifail=ludcmp(xa->m,ndim,indx,&cri);

	if(ifail==1) return CRI_MAX;
	for(i=0;i<ndim;i++) cri=cri*xa->m[i][i];
	pop->det=cri;
	for(i=0;i<ndim;i++)
	{
		for(j=0;j<ndim;j++) col[j]=0.0;
		col[i]=1.0;
		lubksb(xa->m,ndim,indx,col);
		for(j=0;j<ndim;j++)
		{
			jj=i*(i+1)/2+j;
			pop->finv[jj]=col[j]; /*=xb[j][i] using another matrix*/
			}
		}
	/* desallocation */
	matrix_destroy(xa);
	free(indx);
	return cri;
}

void lubksb(double **a, int n, int *indx, double b[])
{
	int i,ii=-1,ip,j;
	double sum;

	for (i=0;i<n;i++) {
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		if (ii>=0)
			for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
		else if (sum) ii=i;
		b[i]=sum;
	}
	for (i=n-1;i>=0;i--) {
		sum=b[i];
		for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}

int ludcmp(double **a, int n,int *indx, double *d)
{
	int i,imax,j,k;
	double big,dum,sum,temp;
	double *vv;

	vv=(double *)calloc(n,sizeof(double));
	*d=1.0;
	for (i=0;i<n;i++) {
		big=0.0;
		for (j=0;j<n;j++)

			if ((temp=fabs(a[i][j])) > big) big=temp;

		if (big == 0.0) {

			return 1;}
		vv[i]=1.0/big;
	}
	for (j=0;j<n;j++) {
		for (i=0;i<j;i++) {
			sum=a[i][j];
			for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<n;i++) {
			sum=a[i][j];
			for (k=0;k<j;k++)
				sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big) {
				big=dum;
				imax=i;
			}
		}
		if (j != imax) {
			for (k=0;k<n;k++) {
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			*d = -(*d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != (n-1)) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<n;i++) a[i][j] *= dum;
		}
	}
	free(vv);
	return 0;
}

/* Random number generator */
double randu(double x)
{	double xrand;
	xrand=rand();
	return xrand;
}
double gauss1(long *idum)
{
/* 	float ran1(long *idum);
*/	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*ran1(idum)-1.0;
			v2=2.0*ran1(idum)-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return (double)gset;
	}
}

float ran1(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0 || !iy) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}

/******************************************************************************
 FedorovInit_R
*******************************************************************************/
// [[Rcpp::export]]

List FedorovWynnAlgorithm_Rcpp( List protocols,
                                IntegerVector ndimen,
                                IntegerVector nbprot,
                                IntegerVector numprot,
                                NumericVector freq,
                                IntegerVector nbdata,
                                NumericVector vectps,
                                NumericVector fisher,
                                IntegerVector error,
                                IntegerVector protdep,
                                NumericVector freqdep)
{

IntegerVector nb_protocols = as<IntegerVector>(protocols[0]);
IntegerVector nb_times = as<IntegerVector>(protocols[1]);
IntegerVector nb_dimensions = as<IntegerVector>(protocols[2]);
IntegerVector total_cost = as<IntegerVector>(ndimen);
NumericMatrix samplingTimes = as<NumericMatrix >(protocols[4]);
NumericMatrix fisher_matrices = as<NumericMatrix>(protocols[5]);


// ***************************************************************************************

// ***************
// init parameters
// ***************

IntegerVector nprot = as<IntegerVector>(nb_protocols);
IntegerVector ndim = as<IntegerVector>(nb_dimensions);
int Ntot = total_cost[2];

// matrix for the optimal sampling times
NumericMatrix optimal_sampling_times(nprot[0],nb_times[0]);

int i,j,ij,ik,iprot,icas;

// *******************
//protocols
// *******************

// Vector of individual protocols
PROTOC *allprot;
allprot = (PROTOC *)calloc(nprot[0], sizeof(PROTOC));

// Population protocol
POPROT *pop;

// *********************************************
// sampling times
// allprot[i].tps = sampling time i
// allprot[i].tps[j] = element j of sampling time i
// *********************************************
// fisher matrices
// allprot[iprot].fisher[k]
// *********************************************

for(iprot=0;iprot<nprot[0];iprot++){

  PROTOC_alloc(&allprot[iprot],nb_times[0],ndim[0]);

  NumericVector row_matrix_sampling_times = samplingTimes.row(iprot);

  NumericVector fisher_matrix = fisher_matrices.row(iprot);

  for ( j=0; j< allprot[iprot].ntps;j++){
  allprot[iprot].tps[j] = row_matrix_sampling_times[j];
  }

  icas=0;

  for (ij=0;ij<ndim[0];ij++) {
    for (ik=0;ik<=ij;ik++) {
      allprot[iprot].fisher[icas]=fisher_matrix[icas]*Ntot/nb_times[0];
      icas++;
    }
  }
}

pop = initprot(allprot, protdep, freqdep);

// *********************************************
// optimisation des protocoles
// nok
// *********************************************

int nok = 0;
double sumf = 0;

nok = main_opt(allprot,pop,nprot);

  if(nok>=0) {

    sumf=0;

    for(i=0;i<pop->np;i++) {

      numprot[i]=pop->num[i]+1;
      nbdata[i]=pop->pind[i].ntps;

      freq[i]=pop->freq[i]*Ntot/(nb_times[0]);

      sumf=sumf+freq[i];

      freqdep[i]=pop->freq[i];
    }

    for(i=0;i<pop->np;i++) {
      freq[i]=freq[i]/sumf;

      for(j=0;j<pop->pind[i].ntps;j++) {
        vectps[j]=pop->pind[i].tps[j];

optimal_sampling_times(i,j) = pop->pind[i].tps[j];

      }
    }
    for(i=0;i<(pop->ndim*(pop->ndim+1)/2);i++) {
      fisher[i]=pop->fisher[i];
    }
  }

  return Rcpp::List::create(
    Rcpp::Named("freq") = freq,
    Rcpp::Named("optimal_sampling_times") = optimal_sampling_times,
    Rcpp::Named("fisher") = fisher,
    Rcpp::Named("numprot") = numprot);


} // end FedorovInit_R
'
FedorovWynnAlgorithm_Rcpp = inline::cxxfunction( signature( protocols_input = "list",
                                                            ndimen_input = "integer",
                                                            nbprot_input = "integer",
                                                            numprot_input = "integer",
                                                            freq_input = "numeric",
                                                            nbdata_input = "integer",
                                                            vectps_input = "numeric",
                                                            fisher_input = "numeric",
                                                            nok_input = "integer",
                                                            protdep_input = "integer",
                                                            freqdep_input = "numeric"),
                                                 plugin = "Rcpp",
                                                 incl = incltxtFedorovWynnAlgorithm,
                                                 body = '
          List protocols = Rcpp::as<List>(protocols_input);
          IntegerVector ndimen  = Rcpp::as<IntegerVector>(ndimen_input);
          IntegerVector nbprot  = Rcpp::as<IntegerVector>(nbprot_input);
          IntegerVector numprot = Rcpp::as<IntegerVector>(numprot_input);
          NumericVector freq = Rcpp::as<NumericVector>(freq_input);
          IntegerVector nbdata = Rcpp::as<IntegerVector>(nbdata_input);
          NumericVector vectps = Rcpp::as<NumericVector>(vectps_input);
          NumericVector fisher = Rcpp::as<NumericVector>(fisher_input);
          IntegerVector nok = Rcpp::as<IntegerVector>(nok_input);
          IntegerVector protdep = Rcpp::as<IntegerVector>(protdep_input);
          NumericVector freqdep = Rcpp::as<NumericVector>(freqdep_input);

          return Rcpp::wrap( FedorovWynnAlgorithm_Rcpp(  protocols, ndimen, nbprot, numprot, freq, nbdata, vectps, fisher, nok, protdep, freqdep ) );')

output = FedorovWynnAlgorithm_Rcpp( protocols_input,  ndimen_input, nbprot_input,
                                    numprot_input, freq_input, nbdata_input,
                                    vectps_input, fisher_input, nok_input,
                                    protdep_input, freqdep_input )

return( output )

}

#' Resize the fisher Matrix from a vector to a matrix.
#'
#' @name resizeFisherMatrix
#' @param nbOfDimensions : a numeric for the dimensions of the fisher matrix.
#' @param fisherMatrix : a vector that contain the low triangular Fisher matrix + its main diagonal.
#' @return The Fisher matrix of size nbOfDimensions*nbOfDimensions
#' @export

setGeneric("resizeFisherMatrix",
           function( nbOfDimensions, fisherMatrix )
           {
             standardGeneric("resizeFisherMatrix")
           }
)

#' @rdname resizeFisherMatrix
#' @export

setMethod( f="resizeFisherMatrix",
           definition = function( nbOfDimensions, fisherMatrix )
           {
             M = matrix(0, nrow = nbOfDimensions, ncol = nbOfDimensions)
             k=1
             for (i in 1:nbOfDimensions){
               for (j in 1:i) {
                 M[i,j] = fisherMatrix[k]
                 k=k+1
               }}
             M = (M+t(M))
             diag(M) = diag(M)/2

             return( M )
           })

# ======================================================================================================

#' @rdname setParameters
#' @export

setMethod("setParameters",
          "FedorovWynnAlgorithm",
          function( object, parameters ) {
            object@name = "FedorovWynnAlgorithm"
            object@elementaryProtocols = parameters$elementaryProtocols
            object@numberOfSubjects = parameters$numberOfSubjects
            object@proportionsOfSubjects = parameters$proportionsOfSubjects
            object@showProcess = parameters$showProcess
            return( object )
          })

# ======================================================================================================

#' Get the optimal frequencies
#' @name getOptimalFrequencies
#' @param object An object from the class \linkS4class{FedorovWynnAlgorithm}.
#' @return A vector giving the optimal frequencies
#' @export

setGeneric("getOptimalFrequencies",
           function(object )
           {
             standardGeneric("getOptimalFrequencies")
           })

#' @rdname getOptimalFrequencies
#' @export

setMethod(f="getOptimalFrequencies",
          signature="FedorovWynnAlgorithm",
          definition = function( object )
          {
            return( object@optimalFrequencies )
          }
)

# ======================================================================================================

#' @rdname optimize
#' @export

setMethod(f = "optimize",
          signature = "FedorovWynnAlgorithm",
          definition = function( object, optimizerParameters, optimizationObject )
          {
            # get fim type
            fim = getFim( optimizationObject )

            # generate Fims from constraints
            fims = generateFimsFromConstraints( optimizationObject )

            # generate initial elementary protocols
            elementaryProtocolsFW = getElementaryProtocols( optimizationObject, fims )

            # parameters for FedorovWynn algorithm in Rcpp
            numberOfprotocols = elementaryProtocolsFW$numberOfprotocols
            totalCost = elementaryProtocolsFW$totalCost
            numberOfTimes = elementaryProtocolsFW$numberOfTimes
            nbOfDimensions = elementaryProtocolsFW$nbOfDimensions
            samplingTimes = elementaryProtocolsFW$samplingTimes
            fisherMatrices = elementaryProtocolsFW$fisherMatrices

            ndim = nbOfDimensions
            ndimen = c( numberOfprotocols, nbOfDimensions, totalCost )
            npInit = numberOfprotocols
            numprot = rep(0,ndim*2)
            freq = rep(0,ndim*2)
            nbdata = rep(0,ndim*2)
            vectps = rep(0,length(numprot))
            fisher = rep(0,ndim*(ndim + 1)/2)
            nok = 0

            # indices for initial oversampling in matrix of initialElementaryProtocols
            optimizerParameters = getOptimizerParameters( optimizationObject )
            initialElementaryProtocols = unlist( optimizerParameters$elementaryProtocols )

            indexElementaryProtocols = rowSums(samplingTimes == initialElementaryProtocols[col(samplingTimes)]) == ncol(samplingTimes)
            indexElementaryProtocols = which( indexElementaryProtocols == TRUE )
            numberOfElementaryProtocol = length( indexElementaryProtocols )
            proportionsOfSubjects = optimizerParameters$proportionsOfSubjects

            # initial protocols, number and indices
            zeprot = rep(0,ndim*2)
            zeprot = c( numberOfElementaryProtocol, indexElementaryProtocols )

            # initials frequencies
            zefreq = rep(0,ndim*2)
            proportionsOfSubjects = getProportionsOfSubjects( optimizationObject )
            zefreq[1:length(proportionsOfSubjects)] = proportionsOfSubjects

            # run the FedorovWynn algorithm
            output = FedorovWynnAlgorithm_Rcpp( elementaryProtocolsFW, ndimen, npInit, numprot, freq, nbdata, vectps,
                                                fisher, nok, zeprot, zefreq )

            # optimalSamplingTimes : remove rows of 0 and convert matrix to a list of vector
            optimalSamplingTimes = output$optimal_sampling_times
            indexOptimalSamplingTimes = which(rowSums( optimalSamplingTimes )>0)

            # check if FW has converged
            if ( length( indexOptimalSamplingTimes ) == 0 )
            {
              print ( " ==================================================================================================== ")
              print ( paste0( " The algorithm has not converged " ) )
              print ( " ==================================================================================================== ")
              stop()
            }

            optimalSamplingTimes = optimalSamplingTimes[indexOptimalSamplingTimes,]

            # optimal frequencies
            optimalFrequencies = output$freq[output$freq>0]

            # indices optimalSamplingTimes in sampling Times
            indexOptimalSamplingTime = output$numprot
            indexOptimalSamplingTime = indexOptimalSamplingTime[indexOptimalSamplingTime>0]

            # optimal design
            optimalDesign = Design( name = c( "Design optimized" ) )

            # optimal sampling times
            designs = getDesigns( optimizationObject )
            design = designs[[1]]
            arms = getArms( design )
            arm = arms[[1]]

            samplingTimesConstraints = getSamplingTimesConstraints( arm )

            numberOfsamplingsOptimisable = lapply( samplingTimesConstraints, function(x) getNumberOfsamplingsOptimisable (x) )
            numberOfsamplingsOptimisable = unlist( numberOfsamplingsOptimisable )

            outcomes = unlist( lapply( samplingTimesConstraints, function(x) getOutcome( x ) ) )

            splitOptimalSamplingTimes = list()

            if ( length( indexOptimalSamplingTimes ) > 1 )
            {
              for ( i in 1:length( indexOptimalSamplingTimes ) )
              {
                splitOptimalSamplingTimes[[i]] = split(optimalSamplingTimes[i,],
                                                       rep( 1:length(numberOfsamplingsOptimisable ), numberOfsamplingsOptimisable ) )

                names( splitOptimalSamplingTimes[[i]] ) = outcomes
              }
            }else {
              splitOptimalSamplingTimes[[1]] = split( optimalSamplingTimes,
                                                      rep( 1:length(numberOfsamplingsOptimisable ), numberOfsamplingsOptimisable ) )

              names( splitOptimalSamplingTimes[[1]] ) = outcomes
            }

            # optimal doses
            optimalDoses = unlist( fims$designArmDose )
            optimalDoses = optimalDoses[indexOptimalSamplingTime]

            # case : populationFIM

            # number of individuals and freq
            numberOfIndividuals = optimizerParameters$numberOfSubjects * optimalFrequencies

            #  optimal sampling times and responses
            listArms = list()
            k = 1

            for( ind in indexOptimalSamplingTime )
            {
              # set size
              arm = fims$listArms[[ind]]
              arm = setSize( arm, numberOfIndividuals[k] )
              armName = paste0("Arm",ind)
              arm = setName( arm, armName )

              # set administration
              administration = getAdministrations( arm )
              administration = setDose( administration[[1]], optimalDoses[k] )
              arm = setAdministrations( arm, list( administration ) )

              # set sampling times
              samplingTimesArms = list()

              for ( outcome in outcomes )
              {
                # get samplingTimes and samplings
                samplings =  splitOptimalSamplingTimes[[k]][[outcome]]

                samplingTime = getSamplingTime( arm, outcome )
                samplingTime = setSamplings( samplingTime, samplings )
                arm = setSamplingTime( arm, samplingTime )
              }
              listArms = append( listArms, arm )
              k=k+1
            }

            optimalDesign = setArms( optimalDesign, listArms )

            # Fisher Matrix : from vector to matrix
            output$fisher = resizeFisherMatrix(  nbOfDimensions, output$fisher )

            if ( class ( fim ) %in% c( "IndividualFim", "BayesianFim" ) )
            {
              listArms = list()

              # case : individualFIM & bayesianFIM
              indexOptimalSamplingTime = indexOptimalSamplingTime[1]
              optimalDoses = optimalDoses[ which.max( indexOptimalSamplingTime ) ]
              numberOfIndividuals = optimizerParameters$numberOfSubjects

              for ( ind in indexOptimalSamplingTime )
              {
                # set size and name
                arm = fims$listArms[[ind]]
                arm = setSize( arm, 1 )
                armName = paste0("Arm",ind)
                arm = setName( arm, armName )

                # set administration
                administration = getAdministrations( arm )
                administration = setDose( administration[[1]], optimalDoses[1] )
                arm = setAdministrations( arm, list( administration ) )

                for ( outcome in outcomes )
                {
                  # get samplingTimes and samplings
                  samplings = splitOptimalSamplingTimes[[1]][[outcome]]

                  samplingTimes = getSamplingTime( arm, outcome )
                  samplingTimes = setSamplings( samplingTimes, samplings )

                  arm = setSamplingTime( arm, samplingTimes )
                }
                listArms = append( listArms, arm )
              }

              optimalDesign = setArms( optimalDesign, listArms )

              # Fisher Matrix : from vector to matrix
              output$fisher = elementaryProtocolsFW$fisherMatrices[indexOptimalSamplingTime,]
              output$fisher = resizeFisherMatrix( ndim,output$fisher )
              output$freq = c(1.0)
            }

            # =====================================
            # outputs
            # =====================================

            object = setOptimalDesign( object, optimalDesign )
            object@optimalDoses = optimalDoses
            object@FisherMatrix = output$fisher
            object@optimalFrequencies = list( listArms = listArms, optimalFrequencies = optimalFrequencies )
            object@optimalSamplingTimes = splitOptimalSamplingTimes

            return( object )
          })

# ======================================================================================================

#' @rdname getDataFrameResults
#' @export

setMethod(f="getDataFrameResults",
          signature = "FedorovWynnAlgorithm",
          definition = function( object )
          {
            optimalFrequenciesAndArms = getOptimalFrequencies( object )

            optimalFrequencies = optimalFrequenciesAndArms$optimalFrequencies
            optimalFrequencies = round( optimalFrequencies, 2 )

            arms = optimalFrequenciesAndArms$listArms

            armNames = unlist( lapply( arms, function(x) getName( x ) ) )
            armSizes = unlist( lapply( arms, function(x) getSize( x ) ) )
            armSizes = round( armSizes, 2 )

            armNames = unique( armNames )
            armSizes = unique( armSizes )

            armsAndOptimalFrequencies = data.frame( armNames = armNames,
                                                    armSizes = armSizes,
                                                    optimalFrequencies = optimalFrequencies )

            colnames( armsAndOptimalFrequencies ) = c("Arm", "Size", "Frequency" )
            rownames( armsAndOptimalFrequencies ) = NULL

            return( armsAndOptimalFrequencies )
          })

# ======================================================================================================

#' @rdname plotFrequencies
#' @export

setMethod(f="plotFrequencies",
          signature = "FedorovWynnAlgorithm",
          definition = function( object )
          {
            data = getDataFrameResults( object )
            data = data[order(data$Frequency), ]
            data$Arm = factor(data$Arm, levels = data$Arm)

             plotData = ggplot( data, aes( x = data[,1], y = data[,3] ) ) +

              geom_bar( stat = "identity", width = 0.5, position = "identity" ) +

              scale_y_continuous( name = "\n Frequencies\n", limits = c(0, 1.05),
                                  breaks = scales::pretty_breaks(n = 10), expand = c(0, 0) ) +

              scale_x_discrete( name = "Arms \n" ) +

              theme(
                legend.position = "none",
                axis.title.x.top = element_text(color = "red", vjust = 2.0),
                axis.text.x.top = element_text(angle = 90, hjust = 0, color = "red"),
                plot.title = element_text(size = 16, hjust = 0.5),
                axis.title.x = element_text(size = 16),
                axis.title.y = element_text(size = 16),
                axis.text.x = element_text(size = 16, angle = 0, vjust = 0.5, color = "black"),
                axis.text.y = element_text(size = 16, angle = 0, vjust = 0.5, hjust = 0.5, color = "black"),
                strip.text.x = element_text(size = 16)
              ) +

              coord_flip()

            return( plotData )
          })

# ======================================================================================================

#' @title show
#' @rdname show
#' @param object object
#' @export

setMethod(f="show",
          signature = "FedorovWynnAlgorithm",
          definition = function( object )
          {
            armsAndOptimalFrequencies = getDataFrameResults( object )

            armsAndOptimalFrequencies = armsAndOptimalFrequencies[ order( armsAndOptimalFrequencies$Frequency, decreasing = TRUE ), ]

            colnames( armsAndOptimalFrequencies ) = c("Arm","Size", "Optimal frequency")

            cat( " ************************************************* ")
            cat("\n")
            cat( " Arm, size and optimal frequency")
            cat("\n")
            cat( " ************************************************* ")
            cat("\n\n")

            print( armsAndOptimalFrequencies )

          })

# ======================================================================================================
# generateReportOptimization
# ======================================================================================================

#' @rdname generateReportOptimization
#' @export

setMethod( "generateReportOptimization",
           signature = "FedorovWynnAlgorithm",
           definition = function( object, optimizationObject, outputPath, outputFile, plotOptions )
           {
             # projectName and outputs tables
             projectName = getName( optimizationObject )

             evaluationFIMResults = getEvaluationFIMResults( optimizationObject )
             fimType = is( getFim( evaluationFIMResults ) )[1]

             evaluationFIMIntialDesignResults = getEvaluationInitialDesignResults( optimizationObject )

             tablesEvaluationFIMIntialDesignResults = generateTables( evaluationFIMIntialDesignResults, plotOptions )

             tablesOptimizationObject = generateTables( optimizationObject, plotOptions )

             # markdown template
             path = system.file(package = "PFIM")
             path = paste0( path, "/rmarkdown/templates/skeleton/" )
             nameInputFile = paste0( path, "template_FedorovAlgorithm.rmd" )

             rmarkdown::render( input = nameInputFile,
                                output_file = outputFile,
                                output_dir = outputPath,
                                params = list(
                                  object = "object",
                                  plotOptions = "plotOptions",
                                  projectName = "projectName",
                                  fimType = "fimType",
                                  tablesEvaluationFIMIntialDesignResults = "tablesEvaluationFIMIntialDesignResults",
                                  tablesOptimizationObject = "tablesOptimizationObject" ) )

           })

##########################################################################################################
# END Class "FedorovWynnAlgorithm"
##########################################################################################################