/*
   m__o
   Any point-->Orbit by iterations
*/

#include "common.h"
#include "corefunc.h"

#include "stagen.h"

#include "_util.h"
#include "_linalg.h"

#include "m__o.h"
#include "_iter.h"	/* iterator's (i.e. generator's) data */


/* Numbers of visible name classes;
   the order must be the same as in /names section */
#define CL_TIME		1
#define CL_PHASE	2
#define CL_PAR		3
#define CL_UTEST	4


Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(m__o_DataPtr) staData=NULL;	/* global copy of our own data area */
Local(StdLinAlgDataPtr) staFunc=NULL;	/* global copy of pointer to standard linea algebra functions */
Local(UtilitiesDataPtr) staUtil=NULL;	/* global copy of pointer to utility functions */
Local(Int2) npar,nphase;		/* dims of par and phase spaces */
Local(Vector) x1=NULL;			/* ptr to vector - the first point produced by starter */
Local(FloatHiPtr) P=NULL;		/* current set of params used by Decoder */

Local(Vector) x0=NULL;			

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(m__o_DataPtr)stagenptr->StaPtr;
  /* Use standard linea algebra */
  staFunc=(StdLinAlgDataPtr)stagenptr->StdLinAlg;
  staUtil=(UtilitiesDataPtr)stagenptr->Utilities;
  /* Get dimensions of Phase and Par spaces */
  npar=*stagenptr->ClassDim[CL_PAR];
  nphase=*stagenptr->ClassDim[CL_PHASE];
  /* Call udf SetInitPoint function, if any */
  if (!stagenData->ContDataGenLen && FuncParActive(staData->SetInitPoint)) {	
                                                        /* there is UDF to set InitPoint */
    Int2 Index_X,Index_P;
    Index_X=stagenptr->ParIndex(staData,X);
    Index_P=stagenptr->ParIndex(staData,P);
    FuncParCall(staData->SetInitPoint)(staData);
    stagenptr->UpdatePar(Index_X);
    stagenptr->UpdatePar(Index_P);
  }
  /* Create a copy parameters */
  if (!P) {
    P=MemNew(sizeof(FloatHi)*npar);
    MemCopy(P,staData->P,sizeof(FloatHi)*npar);
    /* Allocate memory for values of rhs used by Decoder/Define/Tests */
    x0=staFunc->CreateVector(nphase+1);
    x1=staFunc->CreateVector(nphase+1);
  }
  stagenData->BifDataOutLen=0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  return 0;
}

Local(void) Term (Boolean OK) {
    if (OK&&P&&(stagenData->Reason==RF_CLEAN)) stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
    if (P) {
      P=MemFree(P);
      x0=staFunc->DestroyVector(x0);
      x1=staFunc->DestroyVector(x1);
    }
}


/* Any point  -> Orbit Starter */

Entry(Int2) o_o_Starter(StagenDataPtr stagenptr) {
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  if (Init(stagenptr)) return 1;	/* initialize local variables */
  
  /* Create the first point for generator */
  x1[0]=*stagenptr->ClassDim[CL_TIME] ? *staData->T : *staData->Tdef;
  MemCopy(x1+1,staData->X,sizeof(FloatHi)*nphase);
  ((ContDataPtr)stagenptr->GenPtr)->X0=x1;
  /* That's all */
  return 0;
}

/***********/
/* Decoder */
Entry(Int2) o_Decoder(VoidPtr vpoint, FloatHiPtr PNTR indirect, Int2 oper) {
  switch (oper) {
    case DECODER_INIT:	/* point is a pointer to StagenData structure */
    case DECODER_TERM:
      break;
    case DECODER_DIM:	/* returns number of M-points in given G-point */
      return 1;
    default:		/* extract next M-point from G-point pointed by vpoint */
      /* vpoint is a Vector */
      indirect[CL_TIME]=(VAR *)vpoint;
      indirect[CL_PHASE]=(VAR *)vpoint+1;
      indirect[CL_PAR]=P;
  }
  return 0;
}

#define _itmap staData->itmap  			/* number of superpositions */
/**********************/
/* Defining functions */
Entry(Int2) o_DefRHS (Vector x, Vector y) {
Int2 i;
staFunc->CopyVector(x,x0);
  for (i=1; i<=_itmap; i++) {
    stagenData->Rhs(x0+1,P,&x0[0],y+1);    /* 0-component is for time */
    staFunc->CopyVector(y,x0);    
  }
  return 0;
}

/*********************/
/* Default processor */

Entry(Int2) o_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {
  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  MemCopy(x1,x,sizeof(VAR)*(nphase+1));     /* for the last point */
  return 0;
}