1. stevejb
  2. is-solver

Commits

stevejb  committed 3dbf9ee

Adding an additional estimator

  • Participants
  • Parent commits 08ffa59
  • Branches default

Comments (0)

Files changed (1)

File siman_solve/nm_solve_smmddw_loggamma.c

View file
  • Ignore whitespace
+/*
+ * 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_gamma --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;
+}