Source

is-solver / siman_solve / nm_solve_smmddw.c

Full commit
/*
 * nm_solve_smmddw.c
 *
 * Use Nelder-Mead simplex to repeatedly invoke Stephen's solver.
 *
 * 
 * This is for the special case of running SMM with DDW
 * 
 * Usage:
 *
 *	nm_solve --rho=DDD --sigma_v=DDD --a=DDD --gamma=DDD --s=DDD \
 *        --lambda_1=DDD --lambda_2=DDD --delta=DDD --epsilon_P=DDD  \
 *	  --theta=DDD  [--prefix=PPPP]                               \
 *
 *	--prefix is optional prefix for generated files.
 */

#include <math.h>
#include <tgmath.h>
#include <string.h>
#include <stdio.h>
#include <gsl/gsl_multimin.h>

/* Command to run */
#define SOLVE_COMMAND		"solve_isC-icpc"

/* Solve parameters */
#define MODEL_VALUE_COUNT	8

/*
 * Model value ranges -
 *
 * 	These values are used to bound the stepping during the
 *	simulated annealuing processing. They are also used to
 *	validate the input parameters.
 */

/* STEPHEN - FILL IN */
double ModelValueMin[MODEL_VALUE_COUNT] =
{
  -10.0,	/* rho  	*/
  -10.0,	/* sigma_v  	*/
  -10.0,	/* a	       	*/
  -10.0,	/* gamma  	*/
  -10.0,	/* s 		*/
  -10.0,	/* lambda_1 	*/
  -10.0,	/* lambda_2 	*/
  /* -10.0, */	/* delta 	*/
  /* -10.0,	/\* epsilon_P	*\/ */
  -10.0		/* theta	*/
};

/* STEPHEN - FILL IN */
double ModelValueMax[MODEL_VALUE_COUNT] =
{
  10.0,		/* rho  	*/
  10.0,		/* sigma_v 	*/
  10.0,		/* a 	 	*/
  10.0,		/* gamma     	*/
  10.0,		/* s		*/
  10.0,		/* lambda_1 	*/
  10.0,		/* lambda_2	*/
  /* 10.0,	 */	/* delta 	*/
  /* 10.0,		/\* epsilon_P	*\/ */
  10.0		/* theta	*/
};

char *ModelValueNames[MODEL_VALUE_COUNT] = 
{
  "rho",	
  "sigma_v",	
  "a",		
  "gamma",	
  "s",		
  "lambda_1",	
  "lambda_2",	
  /* "delta",	 */
  /* "epsilon_P", */
  "theta"
};

/* Global variables */
static char	*Prefix = "Pref";
static char     *Momentlist = "Momentlist";
static char     *Datamoments = "Datamoments";
static char     *Gmmweightmatrix = "Gmmweightmatrix";
static char     *DoSMM = "--do_smm_for_ddw --add_ddw='1' --delta='0.15' --epsilon_P='-4.0' --log_sigma_v ";
int		 Iter;

double	SolveEnergy(const gsl_vector *, void *);

/* Utility functions */
void PrintModel(char *, gsl_vector *);

/* Main program */
int main(int argc, char *argv[])
{
  int		i;
  int		j;
  double	Model[MODEL_VALUE_COUNT];

  /* Initialize model values */
  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    Model[i] = INFINITY;
  }

  /* Capture model values arguments */
  for (i = 1; i < argc; i++)
  {
    if (strncmp(argv[i], "--prefix=", 9) == 0)
    {
      Prefix = &argv[i][9];
      printf("Prefix set to %s\n", Prefix);
      continue;
    }

    if (strncmp(argv[i], "--momentlist=", 13) == 0) {
	Momentlist = &argv[i][13];
	printf("Momentlist set to %s\n", Momentlist);
	continue;
    }

    if (strncmp(argv[i], "--datamoments=", 14) == 0) {
	Datamoments = &argv[i][14];
	printf("Datamoments set to %s\n", Datamoments);
	continue;
    }

    if (strncmp(argv[i], "--gmm_weight_matrix=", 20) == 0) {
	Gmmweightmatrix = &argv[i][20];
	printf("Gmmweightmatrix set to %s\n", Gmmweightmatrix);
	continue;
    }


    j = ModelNameIndex(argv[i] + 2);
    
    if (j != -1)
    {
      Model[j] = atof(argv[i] + 3 + strlen(ModelValueNames[j]));
      continue;
    }

    printf("Usage: %s --rho=d.dd ... --epsilon_P=d.dd [--prefix=PPPP]\n", argv[0]);
    exit(1);
  }

  /* Validate model values */
  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    if (Model[i] == INFINITY)
    {
      printf("Error: No value supplied for parameter '%s'\n", ModelValueNames[i]);
      exit(2);
    }

    if (Model[i] < ModelValueMin[i])
    {
      printf("Error: Value supplied for parameter '%s' less than min.\n", ModelValueNames[i]);
      exit(2);
    }

    if (Model[i] > ModelValueMax[i])
    {
      printf("Error: Value supplied for parameter '%s' more than max.\n", ModelValueNames[i]);
      exit(3);
    }
  }
  
  /* // test the argument parsing */
  /* exit(1); //  */

  /* Set IEEE floating point mode from environment */
  gsl_ieee_env_setup();

  /* Get set up for Neller Mead */
  const gsl_multimin_fminimizer_type	*T = gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer       	*S = NULL;
  gsl_vector                    	*StartVec;
  gsl_vector                    	*StepVec;
  gsl_multimin_function          	 MinexFunc;
  int 					 Status;
  double				 Size;

  /* Create starting point */
  StartVec = gsl_vector_alloc(MODEL_VALUE_COUNT);
  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    gsl_vector_set(StartVec, i, Model[i]);
  }
  
  /* Set initial step size to 0.1 */
  StepVec = gsl_vector_alloc(MODEL_VALUE_COUNT);
  gsl_vector_set_all(StepVec, 0.1);

  /* Initialize method */
  MinexFunc.n      = MODEL_VALUE_COUNT;
  MinexFunc.f      = &SolveEnergy;
  MinexFunc.params = NULL;

  S = gsl_multimin_fminimizer_alloc(T, MODEL_VALUE_COUNT);
  gsl_multimin_fminimizer_set(S, &MinexFunc, StartVec, StepVec);

  do
  {
    Iter++;
    Status = gsl_multimin_fminimizer_iterate(S);
           
    if (Status) 
      break;
     
    Size   = gsl_multimin_fminimizer_size(S);
    Status = gsl_multimin_test_size(Size, 1e-2);
     
    if (Status == GSL_SUCCESS)
    {
      printf ("converged to minimum at\n");
    }
     
    printf ("%5d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e   %10.3e f() = %7.3f size = %.3f\n", 
	    Iter,
	    gsl_vector_get(S->x, 0), 
	    gsl_vector_get(S->x, 1), 
	    gsl_vector_get(S->x, 2), 
	    gsl_vector_get(S->x, 3), 
	    gsl_vector_get(S->x, 4), 
	    gsl_vector_get(S->x, 5), 
	    gsl_vector_get(S->x, 6), 
	    gsl_vector_get(S->x, 7), 
	    S->fval, Size);
  }
  while (Status == GSL_CONTINUE && Iter < 1000);

  PrintModel("Solved", S->x);

  gsl_vector_free(StartVec);
  gsl_vector_free(StepVec);
  gsl_multimin_fminimizer_free(S);

}

/*
 * PrintModel -
 *
 *	Display the model with the specified label. No return value.
 */

void PrintModel(char *Label, gsl_vector *Model)
{
  int	i;

  printf("%s:\n", Label);
  printf("  ");

  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    printf("%s=%f", ModelValueNames[i], gsl_vector_get(Model, i));
    if (i != (MODEL_VALUE_COUNT - 1))
    {
      printf(", ");
    }
  }
  printf("\n");
}

/*
 * SolveEnergy -
 */

double SolveEnergy(const gsl_vector *v, void *params)
{
  int		 i;
  char		 Command[2048];
  char		 CommandArg[1024];
  char		 CommandOut[2048];
  char           MlCommand[1024];
  char           DMCommand[1024];
  char           GMWCommand[1024];
  char		*EnPtr;
  FILE		*FP;
  double 	 Energy;

  /* Create basic command line */
  strcpy(Command, "");
  strcat(Command, SOLVE_COMMAND);
  strcat(Command, " ");

  /* Add iteration */
  sprintf(CommandArg, "--uuid_i='%s%d'", Prefix, Iter);
  strcat(Command, CommandArg);
  strcat(Command, " ");

  /* Add model values */
  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    sprintf(CommandArg, "--%s='%f'", ModelValueNames[i], gsl_vector_get(v, i));
    strcat(Command, CommandArg);
    strcat(Command, " ");
  }

  /* Add smm flag */
  strcat(Command, DoSMM);
  strcat(Command, " ");

  /* Add momentlist */
  sprintf(MlCommand, "--momentlist='%s'", Momentlist);
  strcat(Command, MlCommand);
  strcat(Command, " ");

  /* Add datamoments */
  sprintf(DMCommand, "--datamoments='%s'", Datamoments);
  strcat(Command, DMCommand);
  strcat(Command, " ");

  /* Add gmm_weight_matrix */
  sprintf(GMWCommand, "--gmm_weight_matrix='%s'", Gmmweightmatrix);
  strcat(Command, GMWCommand);
  strcat(Command, " ");


  /* Check to see if this worked */
  /* printf("Running command plus ML + DM: %s\n", Command); */
  printf("Running command: %s\n", Command);


  Energy = INFINITY;
  FP = popen(Command, "r");
  while (fgets(CommandOut, sizeof CommandOut, FP))
  {
    if (strstr(CommandOut, "Hey Dad") != NULL)
    {
      EnPtr = strstr(CommandOut, "Energy=");
      if (EnPtr != NULL)
      {
	if (sscanf(EnPtr + 7, "%lf", &Energy) == 1)
	{
	  printf("Found energy: %f\n", Energy);
	  pclose(FP);
	  break;
	}
      }
    }
  }
  return Energy;
}

int ModelNameIndex(char *ModelValueName)
{
  int i;
  for (i = 0; i < MODEL_VALUE_COUNT; i++)
  {
    if (strncmp(ModelValueName, ModelValueNames[i], strlen(ModelValueNames[i])) == 0)
    {
      return i;
    }
  }
  return -1;
}