/*
   File:        aggress.c
   Author:      Andrew W. Moore
   Created:     Wed Apr 12 20:55:15 EDT 1995
   Description: Aggressive optimization methods

   Copyright (C) 1995, Andrew W. Moore
*/

#include <stdio.h>
#include <math.h>
#include "aggress.h"
#include "simplex.h"
#include "parse.h"

void autonp_agg_params(char *m,agg_params *ap,double verb)
{
  if ( Verbosity > verb )
  {
    fprintf_agg_params(stdout,m,ap,"\n");
    wait_for_key();
  }
}

/* OE0 - Its evaluation is composed of two things.  First is its estimate
 *       of the oneshot return for the experiment.  Second is the incremental
 *	 increase in future returns given a future policy of selecting the
 *	 point with the highest predicted return (even though that's not the
 *	 future policy).  The first term is a reasonable probabilistic
 *	 computation.  The second is a gross oversimplification and
 *	 approximation.  Finally, note that this evaluation ignores the
 *	 function end_value, so it should only be used when that function is
 *	 zero for all points
 */
typedef struct oe0_workspace_struct
{
  dyv *best_sin;
  double value;
} oe0_workspace;

oe0_workspace *mk_oe0_workspace()
{
  oe0_workspace *result = AM_MALLOC(oe0_workspace);
  result->best_sin = NULL;
  result->value = 0.0;
  return result;
}

void free_oe0_workspace(char *workspace)
{
  oe0_workspace *ws = (oe0_workspace *)workspace;
  if (ws->best_sin) free_dyv(ws->best_sin);
  AM_FREE(ws,oe0_workspace);
}

/* updates the optimization method specific stuff in workspace */
void oe0_update(dyv *sin, dyv *sout, sofar *sf, OptSpec *os, aprior *apr)
{
  posterior *po;
  double tmp;
  dyv *usin,*usout,*lsout;
  oe0_workspace *ws = (oe0_workspace *)sf->workspace;

  if (!sf->workspace) /* then we'll create it */
  {
    ws = mk_oe0_workspace();
    sf->workspace = (char *)ws;
    sf->workspace_free_fn = free_oe0_workspace;
  }
/* update best point seen so far, correct way would be another simplex search
 * over the whole space, but we'll settle for looking at the last point and
 * a recomputation of the value of the current best
 */
  usin = mk_dyv(sf->at->sm->usin_size);
  sin = mk_dyv(sf->at->sin_size);
  usout = mk_dyv(1);
  lsout = mk_dyv(1);
  /* first recompute the current best's estimated value */
  if (ws->best_sin)
  {
    copy_dyv(ws->best_sin,sin);
    po=mk_posterior_from_atree(apr,sf->at,sf->breq,sin,0);
    dyv_set(lsout,0,prediction_from_posterior(po,sin,0.0,NULL,NULL));
    usout_from_sout_sm(sf->at->sm,lsout,usout);
    usin_from_sin_sm(sf->at->sm,sin,usin);
    ws->value = (os->oneshot_value)(usin,usout,os->oneshot_data);
    free_posterior(po);
  }
  /* then check the most recent point */
  copy_dyv(sf->at->darr->array[sf->at->num_points-1]->sin,sin);
  po=mk_posterior_from_atree(apr,sf->at,sf->breq,sin,0);
  dyv_set(lsout,0,prediction_from_posterior(po,sin,0.0,NULL,NULL));
  usout_from_sout_sm(sf->at->sm,lsout,usout);
  usin_from_sin_sm(sf->at->sm,sin,usin);
  tmp = (os->oneshot_value)(usin,usout,os->oneshot_data);
  if (!ws->best_sin)
  {
    ws->best_sin = mk_copy_dyv(sin);
    ws->value = tmp;
  }
  else if (tmp > ws->value)
  {
    copy_dyv(sin,ws->best_sin);
    ws->value = tmp;
  }
  free_dyv(usin); free_dyv(sin); free_dyv(usout); free_dyv(lsout);
  free_posterior(po);
}

/* used by the simplex search to evaluate possible experiments */
double oe0_eval(posterior *po, agg_params *ap, sofar *sf, OptSpec *os, dyv *x)
{
  dist *response = mk_predict_dist(po,x);
  dist *pred = mk_prediction_dist(response,po->noise_dist,10);
  nonpar_data *d = (nonpar_data *)pred->data;
  double IntSize = ((dyv_ref(pred->xmax,0)-dyv_ref(pred->xmin,0))/
		    ((double)(d->pdf_vals->size-1)));
  double oneshotval,futureval=0.0;
  double result = 0.0;
  dyv *usin,*usout,*sout;
  oe0_workspace *ws = (oe0_workspace *)sf->workspace;
  int i;

/* compute probabilistic oneshot and future values of this experiment */
  sout = mk_dyv(1);
  usout = mk_dyv(1);
  usin = mk_usin_from_sin_sm(sf->at->sm,x);
  dyv_set(sout,0,dyv_ref(pred->xmin,0));
  for (i=0;i<d->pdf_vals->size;i++,sout->farr[0]+=IntSize)
  {
    usout_from_sout_sm(sf->at->sm,sout,usout);
    oneshotval = (os->oneshot_value)(usin,usout,os->oneshot_data);
    if (os->duration_type == FINITE_KNOWN)
    {
      if (oneshotval > ws->value)
	futureval = ((oneshotval - ws->value)*
		     (double)(os->duration_steps - sf->at->num_points));
      else futureval = 0.0;
    }
    else my_error("oe0: does not support this duration_type");
    result += dyv_ref(d->pdf_vals,i)*IntSize*(oneshotval+futureval);
  }

  free_dyv(sout); free_dyv(usout); free_dyv(usin);
  free_dist(response); free_dist(pred);
  return result;
}

/* OE1: gitten's indices ------------------------------------------- */

void oe1_update(dyv *sin, dyv *sout, sofar *sf, OptSpec *os, aprior *apr)
  {
  return;
  }

/* this "magic" function returns the advantage of pulling arm1 first then
 * behaving optimally over always pulling arm2, where
 *      arm 1 -> with specified prior, noisless rewards,
 *      arm 2 -> gauranteed/known reward=guess_index
 * value is >0 if guess_index is too low, <0 if too high. Is used to do 
 * a search (probably binary) for the true gittens's index. 
 *
 * This function was mostly modified from Jeff Schneider's oe0_eval().
 */
double eval_gittens_magic(nonpar_data *d, dist *response, dist *pred,
		sofar *sf, OptSpec *os, dyv *x, double IntSize,
		double total_weight, double future_weight, double guess_index)
  {
  double result;	/* will hold reward for pulling on arm 1 first,
			   then behaving optimally after that. */
  int i;
  dyv *sout,*usout,*usin;
  double oneshotval, futureval;

  /* An approximation: Assume we can ignore all indices < 0. (i.e. that
     all indices are  >= 0 */
  if (guess_index < 0)
	return -guess_index;	/* >0 indicates guess_index was too low */

  /* compute probabilistic oneshot and future values of this experiment */
  sout = mk_dyv(1);
  usout = mk_dyv(1);
  usin = mk_usin_from_sin_sm(sf->at->sm,x);
  dyv_set(sout,0,dyv_ref(pred->xmin,0));

  result = 0.0;
  for (i=0;i<d->pdf_vals->size;i++,sout->farr[0]+=IntSize)
    {
    usout_from_sout_sm(sf->at->sm,sout,usout);
    oneshotval = (os->oneshot_value)(usin,usout,os->oneshot_data);

    if (oneshotval > guess_index)
	futureval = oneshotval * future_weight;
    else 
	futureval = guess_index * future_weight;
    result += dyv_ref(d->pdf_vals,i)*IntSize*(oneshotval+futureval);
    }

  free_dyv(sout); 
  free_dyv(usout); 
  free_dyv(usin);

  return result - guess_index * total_weight;
  }

#define OE1_EPSILON (2.0e-3)	/* approximate allowed absolute error in 
				   gitten's index. (BUT: IntSize will 
			           probably be a much larger source of 
                                   error than this.) */

double oe1_eval(posterior *po, agg_params *ap, sofar *sf, OptSpec *os, dyv *x)
  {
  dist *response = mk_predict_dist(po,x);
  dist *pred = mk_prediction_dist(response,po->noise_dist,10);
  nonpar_data *d = (nonpar_data *)pred->data;
  double IntSize = ((dyv_ref(pred->xmax,0)-dyv_ref(pred->xmin,0))/
		    ((double)(d->pdf_vals->size-1)));

  double a, b, c;	/* for doing a binary search for gitten's index */
  double diff, eval;
  double total_weight, future_weight;
  double beta;		/* discount rate */ 
  static int warned=0;

  if (os->duration_type != INFINITE && !warned)
	{
	fprintf(stderr,"WARNING: OE1 (gitten's indices) being used "
		        " with non-infinite horizon.\n"
                       "WARNING: (Optimality proofs don't hold!)\n");
	warned=1;
	}

  if (os->duration_type == FINITE_KNOWN)
	{
	future_weight = (double)(os->duration_steps - sf->at->num_points);
	total_weight = (double)(os->duration_steps - sf->at->num_points) + 1;
	}
  else if (os->duration_type == INFINITE)
	{
	beta = os->discount_rate;
	future_weight = beta/(1.0-beta);
	total_weight = 1.0/(1.0-beta);
	}
  else
	{
	my_error("oe1: Supports only FINITE_KNOWN and INFINITE "
		  "horizon types. Sorry.");
	}

  /* Start doing a binary search for the true gittens index. 
     Our bracket will always be (a,b). First, establish the bracket */
  a = -0.1;
  while (eval_gittens_magic(d, response, pred, sf, os, x, IntSize, 
		total_weight, future_weight, a) < 0)
  	a *= 1.5;

  diff=0.2;
  b = a + diff;
  while (eval_gittens_magic(d, response, pred, sf, os, x, IntSize, 
		total_weight, future_weight, b) > 0)
	{
	diff *= 1.5;
  	b = a+diff;
	}

  printf("{(%f,%f)}", (float)a, (float)b);

  /* Commence binary search */
  while (fabs(b-a) > OE1_EPSILON)
	{
	c = (a+b)/2.0;
	eval= eval_gittens_magic(d, response, pred, sf, os, x, IntSize, 
			total_weight, future_weight, c);
	if (eval > 0)
		a = c;
	else
		b = c;
	}

  free_dist(response); free_dist(pred);

  printf("(%.2f,%.2f->%.5f)", (float)dyv_ref(x,0), (float)dyv_ref(x,1),
			(float)(a+b)/2.0);

  return (a+b)/2.0;
  }

/* ----------------------------------------------------------------- */

double aggress_eval(posterior *po,agg_params *ap,sofar *sf,OptSpec *os,dyv *x)
{
  dist *ds = mk_predict_dist(po,x);
  double result;

  autonp_agg_params("aggress-eval-ap",ap,1.987);
  autonp_dyv("aggress-eval-x",x,1.987);
  if ( Verbosity > 1.987 )
  {
    printf("dist_univariate_mean(ds) = %g\n",
            dist_univariate_mean(ds)
           );
    printf("invert_univariate_cdf(ds,ap->opt_conf) = %g\n",
            invert_univariate_cdf(ds,ap->opt_conf)
           );
    printf("dist_univariate_cdf(ds,ap->min_resp) = %g\n",
            dist_univariate_cdf(ds,ap->min_resp)
           );
    wait_for_key();
  }

  if ( ap->agg_type == PREDICTED_MAX )
    result = dist_univariate_mean(ds);
  else if ( ap->agg_type == OPTIMISTIC_MAX )
    result = invert_univariate_cdf(ds,ap->opt_conf);
  else if ( ap->agg_type == CAUTIOUS_OPTIMISTIC_MAX )
  {
    double opt_case = invert_univariate_cdf(ds,ap->opt_conf);
    double pess_case = invert_univariate_cdf(ds,ap->pess_conf);
    
    result = opt_case - 100.0 * real_max(0.0,ap->min_resp - pess_case);
  }
  else if ( ap->agg_type == OE0)
    result = oe0_eval(po,ap,sf,os,x);
  else if ( ap->agg_type == OE1)
    result = oe1_eval(po,ap,sf,os,x);
  else 
  {
    result = 0.0;
    my_error("unknown ap type");
  }

  autonp_realnum("aggess-eval-result",result,1.987);
  free_dist(ds);
  return(result);
}

void default_agg_params(agg_params *ap)
{
  ap -> agg_type = PREDICTED_MAX;
  ap -> opt_conf = 0.95;
  ap -> min_resp = 0.2;
  ap -> pess_conf = 0.075;
}

void update_agg_params_from_args(agg_params *ap,int argc,char *argv[])
{
  bool pmax  = index_of_arg("pmax",argc,argv)  >= 0;
  bool iemax = index_of_arg("iemax",argc,argv) >= 0;
  bool comax = index_of_arg("comax",argc,argv) >= 0;
  bool oe0   = index_of_arg("oe0",argc,argv)   >= 0;
  bool oe1   = index_of_arg("oe1",argc,argv)   >= 0;
  int s = ((pmax)  ? 1 : 0) + ((iemax) ? 1 : 0) + 
          ((comax) ? 1 : 0) + ((oe0)   ? 1 : 0)   +
          ((oe1)   ? 1 : 0);		/*AYN*/

  if ( s > 1 ) my_error("more than one of pmax, iemax and comax on command");

  ap -> agg_type = (comax) ? CAUTIOUS_OPTIMISTIC_MAX :
                   (iemax) ? OPTIMISTIC_MAX : 
                   (pmax) ? PREDICTED_MAX :
		   (oe0) ? OE0 :
		   (oe1) ? OE1 : 		/*AYN*/
                   ap -> agg_type;

  ap -> opt_conf = double_from_args("opt_conf",argc,argv,ap->opt_conf);
  ap -> min_resp = double_from_args("min_resp",argc,argv,ap->min_resp);
  ap -> pess_conf = double_from_args("pess_conf",argc,argv,ap->pess_conf);
}

char *string_from_agg_type(int at)
{
  char *result =
          (at == PREDICTED_MAX) ? "PREDICTED_MAX" :
          (at == OPTIMISTIC_MAX) ? "OPTIMISTIC_MAX" :
          (at == CAUTIOUS_OPTIMISTIC_MAX) ? "CAUTIOUS_OPTIMISTIC_MAX" :
	  (at == OE0) ? "OE0" :
	  (at == OE1) ? "OE1" : 		/* AYN */
          "<UnKnOwN AgGtYpE>";
  return(result);
}

void fprintf_agg_params(FILE *s,char *m1,agg_params *ap,char *m2)
{
  fprintf(s,"%s -> agg_type = %s%s",
          m1,string_from_agg_type(ap->agg_type),m2
         );
  fprintf(s,"%s -> opt_conf = %g%s",m1,ap->opt_conf,m2);
  fprintf(s,"%s -> min_resp = %g%s",m1,ap->min_resp,m2);
  fprintf(s,"%s -> pess_conf = %g%s",m1,ap->pess_conf,m2);
}

double eval_aggress_exp(aprior *apr,agg_params *ap, sofar *sf,OptSpec *os,
			dyv *query)
{
  posterior *po;
  double result;

  po = mk_posterior_from_atree(apr,sf->at,sf->breq,query,0);
  result = aggress_eval(po,ap,sf,os,query);
  free_posterior(po);
  return(result);
}
  
typedef struct ba_exp_struct
{
  aprior *apr;
  agg_params *ap;
  sofar *sf;
  dproj *dp;
  int tnum;
  OptSpec *os;
} ba_exp;

extern dproj draw_auton_state_global_dp[1];

double eval_ba_exp(char *data,dyv *x)
{
  ba_exp *ba = (ba_exp *) data;
  double result = eval_aggress_exp(ba->apr,ba->ap,ba->sf,ba->os,x);
  char buff[10];

  if ( Verbosity > 45.0 )
  {
    fprintf_dyv(stdout,"eval_ba_exp",x,"\n");
    printf("eval = %g\n",result);
    if ( ba->dp != NULL )
      dproj_dyv_dot(ba->dp,x,result);
    wait_for_key();
  }
  if ( Verbosity >= 0.0 )
  {
    if (ba->tnum >= 0)
    {
      sprintf(buff,"%d",ba->tnum);
      printf("%s",buff);
      dproj_dyv_string(draw_auton_state_global_dp,x,0.0,buff);
      fflush(stdout);
    }
  }
  return(result);
}

double eval_ba_exp_2(char *data, double x, double y)
{
  double result;
  dyv *d = mk_dyv(2);
  dyv_set(d,0,x);
  dyv_set(d,1,y);
  result = eval_ba_exp(data,d);
  free_dyv(d);
  if (result < 0.0) return 0.0; /* hack just to make the graphics nice */
  else return result;
}

double best_aggress_exp(aprior *apr, agg_params *ap, sofar *sf, OptSpec *os, 
			double radius,dyv *best_exp)
/*
   Result stored in best_exp. Choose the best experiment from a grid of
   experiment points.

   Uses simplex hill-climber.

   If this is too depressingly slow damn well use
   a more sensible optimization method.
*/
{
  double best_value = 0.0;
  int trials = 3;
  ba_exp ba[1];
  dproj dp[1];
  int tnum;
  dyv *delta;

  ba -> apr = apr;
  ba -> ap = ap;
  ba -> sf = sf;
  ba -> os = os;

  draw_2d_border(dp,0.0,0.0,1.0,1.0);
  ba -> dp = dp;


  ba->tnum = -1; /* just makes these graphics come out ok */
  ag_on("");
  draw_auton_state(NULL,sf->at,sf->base,sf->radius,sf->base,FALSE,"",
		   NULL,NULL,NULL,0.0);
  draw_2d_function(eval_ba_exp_2,(char *)ba,12,15,"Experiment Value",
		   0.0,0.0,1.0,1.0,NULL);
  wait_for_key();


  for ( tnum = 0 ; tnum < trials ; tnum++ )
  {
    dyv *start_x;
    dyv *exp_x = mk_dyv(sf->at->sin_size);
    double value;
    double ftol = 0.01;
    double step = 0.3;

    if ( tnum == 0 && sf->at->num_points > 0 )  /* start at the last point */  
    {                                         
      start_x = mk_copy_dyv(sf->at->darr->array[sf->at->num_points-1]->sin);
    }
    else if (tnum == 1 && ap->agg_type == OE0 && sf->workspace)
    {                                      /* start at the best we've seen */
      oe0_workspace *ws = (oe0_workspace *)sf->workspace;
      start_x = mk_copy_dyv(ws->best_sin);
    }
    else                               /* start at a random previous point */
    {
      start_x=mk_copy_dyv(sf->at->darr->array[int_random(sf->at->num_points)]->sin);
    }

    if ( Verbosity > 8.0 )
    {
      fprintf_dyv(stdout,"start simplex search ",start_x,"\n");
      dproj_dyv_small_circle(dp,start_x,value);
      wait_for_key();
    }

    ba -> tnum = tnum;

    amoeba_search(start_x,step,ftol,eval_ba_exp,(char *)ba,
                  FALSE,TRUE,exp_x,&value
                 );

    if ( Verbosity > 8.0 )
    {
      printf("\n*\n*\n");
      fprintf_dyv(stdout,"exp from simplex search ",exp_x,"\n");
      printf("that was trial number %d with result %g\n",tnum,value);
      dproj_dyv_small_cross(dp,exp_x,value);
      wait_for_key();
    }

    printf("\nbest_aggress_exp()[inside]: exp_x: %f,%f-->%f \n", 
			dyv_ref(exp_x,0), dyv_ref(exp_x,1), value);

    if ( tnum == 0 || value > best_value )
    {
      best_value = value;
      copy_dyv(exp_x,best_exp);
      printf("\nNEW!! best_aggress_exp(): best_exp: %f,%f-->%f \n", 
			dyv_ref(best_exp,0), dyv_ref(best_exp,1), best_value);
    }

    free_dyv(exp_x);
    free_dyv(start_x);
  }

  if ( Verbosity >= 0.0 ) printf("\n");

  printf("\nEnd of 3 trials: best_aggress_exp(): best_sxp: %f,%f-->%f \n", 
			dyv_ref(best_exp,0), dyv_ref(best_exp,1), best_value);

  delta = mk_random_unit_dyv(sf->at->sin_size);
  printf("Original Delta = (%f,%f)\n", dyv_ref(delta,0), dyv_ref(delta,1));

  printf("Check 0: (b:%f,%f  d:%f,%f) \n", 
	best_exp->farr[0], best_exp->farr[1], delta->farr[0], delta->farr[1]);

  printf("radius=%f\n", radius);
  dyv_scalar_mult(delta, 0.2*radius, delta);
  printf("radius=%f\n", radius);

  printf("Check 1: (b:%f,%f  d:%f,%f) \n", 
	best_exp->farr[0], best_exp->farr[1], delta->farr[0], delta->farr[1]);

  dyv_plus(delta,best_exp,best_exp);

  printf("Check 2: (b:%f,%f  d:%f,%f) \n", 
	best_exp->farr[0], best_exp->farr[1], delta->farr[0], delta->farr[1]);

  printf("New Delta = (%f,%f)\n", dyv_ref(delta,0), dyv_ref(delta,1));
  printf("radius = %f\n");

  free_dyv(delta);

  printf("\nbest_aggress_exp(): %f,%f-->%f \n", 
			dyv_ref(best_exp,0), dyv_ref(best_exp,1), best_value);

  
  dyv_clip_in_unit(best_exp);

  return(best_value);
}

typedef struct PLO_data_struct
{
  dyv *in;
  dyv *out;
} PLO_data;

double piecewise_linear_oneshot_value(dyv *usin, dyv *usout, char *data)
{
  int i;
  PLO_data *plo = (PLO_data *)data;

  if (plo->in->size == 0) return 0.0;  /* no function */
  if (plo->in->size == 1) return dyv_ref(plo->out,0);  /* constant function */
  for (i=0;((i<plo->in->size)&&(dyv_ref(usout,0)>dyv_ref(plo->in,i)));i++);
  if (i==0) i++;
  else if (i == plo->in->size) i--;
  if (dyv_ref(plo->in,i) == dyv_ref(plo->in,i-1))
    return dyv_ref(plo->in,i);     /* use righthand value on steps */
/* otherwise do linear interpolation/extrapolation */
  return dyv_ref(plo->out,i-1)+((dyv_ref(usout,0) - dyv_ref(plo->in,i-1)) *
				(dyv_ref(plo->out,i) - dyv_ref(plo->out,i-1))/
				(dyv_ref(plo->in,i) - dyv_ref(plo->in,i-1)));
}

void piecewise_linear_data_free_fn(char *data)
{
  PLO_data *plo = (PLO_data *)data;
  if (plo->in) free_dyv(plo->in);
  if (plo->out) free_dyv(plo->out);
  AM_FREE(data,PLO_data);
}

void set_piecewise_linear_oneshot_value(OptSpec *os,char *fname)
{
  dym *ins = NULL;
  dym *outs = NULL;
  PLO_data *plo;

  os->oneshot_value = piecewise_linear_oneshot_value;
  plo = AM_MALLOC(PLO_data);
  read_io_dyms(fname,NULL,&ins,&outs);
  plo->in = mk_dyv_from_dym_col(ins,0);
  plo->out = mk_dyv_from_dym_col(outs,0);
  os->oneshot_data = (char *)plo;
  os->oneshot_data_free_fn = piecewise_linear_data_free_fn;
  free_dym(ins); free_dym(outs);
}

typedef struct OSV_data_struct
{
  parse_node *pn;
  bindings *bs;
} OSV_data;

double oneshot_value_fn(dyv *usin, dyv *usout, char *data)
{
  dym *in = mk_col_dym_from_dyv(usin);
  dym *out = mk_col_dym_from_dyv(usout);
  dym *result;
  double res;
  OSV_data *d = (OSV_data *)data;

  add_binding(d->bs,"in",in);
  add_binding(d->bs,"out",out);
  result = eval_parse_node(d->pn,d->bs);
  res = dym_ref(result,0,0);
  free_dym(in); free_dym(out); free_dym(result);
  return res;
}

void OSV_data_free_fn(char *data)
{
  OSV_data *d = (OSV_data *)data;
  free_parse_node(d->pn);
  free_bindings(d->bs);
  AM_FREE(data,OSV_data);
}

void set_oneshot_value(OptSpec *os, char *s)
{
  OSV_data *d;

  os->oneshot_value = oneshot_value_fn;
  d = AM_MALLOC(OSV_data);
  d->pn = top_level_parse(s);
  if (!d->pn) my_error("set_oneshot_value: parse failed\n");
  d->bs = make_empty_bindings();
  os->oneshot_data = (char *)d;
  os->oneshot_data_free_fn = OSV_data_free_fn;
}

double default_oneshot_value(dyv *usin, dyv *usout, char *data)
{
  return dyv_ref(usout,0);
}

double default_end_value(dym *usins, dym *usouts, char *data)
{
  return 0.0;
}

void default_print_summary_fn(OptSpec *os, dym *usins, dym *usouts)
{
  double value,total_value = 0.0;
  int i;
  dyv *usin,*usout;

  usin = mk_dyv(usins->cols);
  usout = mk_dyv(usouts->cols);
  for (i=0;i<usins->rows;i++)
  {
    copy_dym_row_to_dyv(usins,usin,i);
    copy_dym_row_to_dyv(usouts,usout,i);
    value = (os->oneshot_value)(usin,usout,os->oneshot_data);
    printf("Experiment %4d has value %8.3f\n",i,value);
    total_value += value;
  }
  value = (os->end_value)(usins,usouts,NULL);
  total_value += value;
  printf("End value of experiment set %8.3f\n",value);
  printf("Total value of experiments  %8.3f\n",total_value);
  free_dyv(usin); free_dyv(usout);
}

OptSpec *mk_default_optspec()
{
  OptSpec *os = AM_MALLOC(OptSpec);
  os->duration_type = FINITE_KNOWN;
  os->duration_steps = 0;
  os->discount_rate = 1.0;
  os->oneshot_value = default_oneshot_value;
  os->oneshot_data = NULL;
  os->oneshot_data_free_fn = NULL;
  os->end_value = default_end_value;
  os->print_summary_fn = default_print_summary_fn;
  return os;
}

void free_optspec(OptSpec *os)
{
  if (os->oneshot_data)
  {
    if (!os->oneshot_data_free_fn) 
      my_error("free_optspec: oneshot_data exists, but free_fn does not\n");
    else (os->oneshot_data_free_fn)(os->oneshot_data);
  }
  AM_FREE(os,OptSpec);
}