/*
   File:        twod.c
   Author:      Andrew W. Moore
   Created:     Sat Feb 25 17:05:16 EST 1995
   Description: Testing auton with two-d functions

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "twod.h"
#include "plant.h"
#include "shootpl.h"
#include "tarmpl.h"
#include "kw.h"
#include "minteg.h"
#include "transform.h"
#include "opt.h"

double default_twod_fn(char *data,dyv *sin,char *exp_code)
{
  dyv *target = mk_constant_dyv(dyv_size(sin),0.8);
  double z;
  int i;
  double result;

  dyv_set(target,0,0.2);
  dyv_subtract(target,sin,target);
  for ( i = 0 ; i < dyv_size(sin) ; i++ )
    dyv_set(target,i,(((i % 2)==0) ? 1.0 : 2.0) * dyv_ref(target,i));
  z = dyv_magnitude(target);
  z = z / (1.0 + z);

  result = 1.0 - z;
  
  if ( char_is_in(exp_code,'p') )
    printf("default_twod_fn(%g,%g) = %g\n",dyv_ref(sin,0),dyv_ref(sin,1),
           result
          );

  free_dyv(target);
  return(result);
}

void default_free_twod_data_fn(char *data)
{
  /* Do absolutely nowt! */
}

void maybe_compute_twodexp_zones(twodexp *td)
{
  if ( td -> auto_compute_zones )
  {
    td -> good_zone = good_zone_thresh(td -> twod_fn,td -> twod_data,
                                       td->high_value_proportion,td->dims
                                      );
    td -> bad_zone = 
       real_min(real_max(td->twod_fn(td->twod_data,td->start_center,"")/2.0,
                         td->twod_fn(td->twod_data,td->start_center,"") - 0.1
                        ),
                0.75
               );
  }
  else
    printf("I notice that td->auto_compute_zones is FALSE. I won't compute the zones\n");
}

twodexp *mk_default_twodexp(int dims, arequest *areq)
{
  twodexp *td = AM_MALLOC(twodexp);

  td -> opt_type     = "rsm";
  td -> dims         = dims;
  td -> steps        = 60;
  td -> high_value_proportion = 0.05;
  td -> start_center = mk_constant_dyv(dims,0.2);
  dyv_set(td->start_center,0,0.8);
  td -> cp           = mk_default_choice_params(dims);
  td -> apr          = mk_default_aprior(max_term_index(dims,areq->regdeg,
							areq->timemode));
  td -> radius       = 0.1;
  td -> dradius      = td->radius/2.0;
  td -> minradius    = td->radius/10.0;
  td -> fixed_dradius= TRUE;
  td -> timeweight   = 0.0;
  td -> autoweight   = 0;

  td -> rsm_final_excite = TRUE;
  td -> use_steepest = FALSE;

  td -> twod_fn      = default_twod_fn;
  td -> twod_data    = (char *) NULL;
  td -> free_twod_data_fn = default_free_twod_data_fn;
  td -> noise        = 0.05;
  td -> explain      = FALSE;
  td -> auto_compute_zones = TRUE;
  td -> bad_zone     = 0.1;
  td -> good_zone    = 0.9;

  default_agg_params(td->ap);
  td -> os           = mk_default_optspec();
  td -> os -> duration_steps = td -> steps;
  td -> auto_kw_steps = 0; /* meaning don't change kw */
  td -> auto_kw_pct = 0.0; 
  td -> contour_display_steps = 5;
  return(td);
}

double npfun_twod_fn(char *data,dyv *sin,char *exp_code)
{
  npfun *npf = (npfun *) data;
  double result = eval_npfun(npf,sin);
  return(result);
}

void npfun_free_twod_data_fn(char *data)
{
  npfun *npf = (npfun *) data;
  free_npfun(npf);
}

double plant_twod_fn(char *data,dyv *sin,char *exp_code)
{
  plant *p = (plant *) data;
  dyv *sout = mk_dyv(1);
  double result;
  if ( char_is_in(exp_code,'a') ) ag_on("");
  eval_plant(p,sin,exp_code,sout);
  if ( char_is_in(exp_code,'a') ) 
  {
    wait_for_key();
    ag_off();
  }
  result = dyv_ref(sout,0);
  free_dyv(sout);
  return(result);
}

void plant_free_twod_data_fn(char *data)
{
  plant *p = (plant *) data;
  free_plant(p);
}

void update_twodexp_from_args(twodexp *td,arequest *areq,int argc,char *argv[])
{
  char *osv_string;
  char *osv_fname;
  dyv *old_start_center = td->start_center;
  int i;

  td -> steps = int_from_args("steps",argc,argv,td -> steps);
  td -> os -> duration_steps = td->steps;
  td -> high_value_proportion = 
    double_from_args("high_value_proportion",
                    argc,argv,td -> high_value_proportion
                   );

  td -> start_center = 
    mk_dyv_from_args("start_center",argc,argv,td -> start_center);

  autonp_dyv("old_start_center",old_start_center,1.0);
  autonp_dyv("new_start_center",td -> start_center,1.0);

  free_dyv(old_start_center);

  update_choice_params_from_args(td->cp,argc,argv);

/* efficiency problem: this will always remove the old apr and put in another
 * even if the command line contains no modification to the default
 */
  if (td->apr) free_aprior(td->apr);
  td->apr = mk_aprior_from_args(max_term_index(td->dims,areq->regdeg,
					       areq->timemode),
				argc,argv);

  td -> radius = double_from_args("radius",argc,argv,td -> radius);

  td -> dradius = double_from_args("dradius",argc,argv,td -> dradius);
  td -> minradius = double_from_args("minradius",argc,argv,td -> minradius);
  td -> fixed_dradius = bool_from_args("fixed_dradius",argc,argv,td -> fixed_dradius);

  td -> noise = double_from_args("noise",argc,argv,td -> noise);
  td -> rsm_final_excite = 
    bool_from_args("rsm_final_excite",argc,argv,td->rsm_final_excite);
  td -> use_steepest = 
    bool_from_args("use_steepest",argc,argv,td->use_steepest);
  td -> auto_compute_zones = 
    bool_from_args("auto_compute_zones",argc,argv,td->auto_compute_zones);

  td -> good_zone = double_from_args("good_zone",argc,argv,td -> good_zone);
  td -> bad_zone = double_from_args("bad_zone",argc,argv,td -> bad_zone);
  td -> explain = index_of_arg("explain",argc,argv) > 0;

  if ( index_of_arg("shoot",argc,argv) >= 0 )
  {
    td -> twod_fn = plant_twod_fn;
    td -> free_twod_data_fn = plant_free_twod_data_fn;
    td -> twod_data = (char *) mk_shoot_plant();
  }
  if ( index_of_arg("tarm",argc,argv) >= 0 )
  {
    td -> twod_fn = plant_twod_fn;
    td -> free_twod_data_fn = plant_free_twod_data_fn;
    td -> twod_data = (char *) mk_tarm_plant(td->dims);
    td -> good_zone = 0.5;
    td -> bad_zone = -0.7;
    td -> auto_compute_zones = FALSE;
  }
  else
  {
    int npfun_seed = int_from_args("npfun",argc,argv,-1);
    int npfun_modes = int_from_args("npfunmodes",argc,argv,1);
    if ( npfun_seed >= 0 )
    {
      npfun *npf;
      npf = mk_random_seeded_npfun(td->dims,npfun_seed,npfun_modes);

      td -> twod_fn = npfun_twod_fn;
      td -> free_twod_data_fn = npfun_free_twod_data_fn;
      td -> twod_data = (char *) npf;
/*
      compute_twodexp_zones(td);
*/
    }
  }

  td->opt_type = string_from_args("opt_type",argc,argv,"rsm");

  update_agg_params_from_args(td->ap,argc,argv);

/* the string is right on the command line */
  osv_string = string_from_args("1shotvals",argc,argv,NULL);
/* its given in the file */
  osv_fname = string_from_args("1shotvalf",argc,argv,NULL);
  if (osv_fname && !osv_string)
  {
    int nread,done,max,oldsize;
    char *tmp = NULL;
    FILE *fp = safe_fopen(osv_fname,"r");

    max = 10000; oldsize = 0; done = 0; osv_string = NULL;
    while(!done)
    {
      tmp = osv_string;
      osv_string = (char *)am_malloc(max * sizeof(char));
      if (oldsize) 
      {
	for (i=0;i<oldsize;i++) osv_string[i]=tmp[i];
	am_free(tmp,oldsize);
      }
      nread = fread(osv_string+oldsize,sizeof(char),max-oldsize,fp);
      if (nread < (max-oldsize)) 
      {
	done = 1;
	osv_string[oldsize+nread] = '\0';
      }
      else
      {
	oldsize = max;
	max *= 2;
      }
    }
    set_oneshot_value(td->os,osv_string);
    am_free(osv_string,max);
  }
  else if (osv_string) set_oneshot_value(td->os,osv_string);
  /* and don't free osv_string because it came from the command line */
  td->auto_kw_steps=int_from_args("auto_kw_steps",argc,argv,td->auto_kw_steps);
  td->auto_kw_pct=double_from_args("auto_kw_pct",argc,argv,td->auto_kw_pct);
  td->contour_display_steps=int_from_args("contour_display_steps",argc,argv,td->contour_display_steps);
}

void fprintf_twodexp(FILE *s,char *m1,twodexp *td,char *m2)
{
  char buff[100];
  fprintf(s,"%s -> opt_type = %s%s",m1,td->opt_type,m2);
  fprintf(s,"%s -> dims = %d%s",m1,td->dims,m2);
  fprintf(s,"%s -> steps = %d%s",m1,td->steps,m2);
  fprintf(s,"%s -> high_value_proportion = %g%s",m1,
          td->high_value_proportion,m2
         );
  sprintf(buff,"%s -> start_center",m1);
  fprintf_dyv(s,buff,td->start_center,m2);
  sprintf(buff,"%s -> cp",m1);
  fprintf_choice_params(s,buff,td->cp,m2);
  sprintf(buff,"%s -> pr",m1);
  fprintf_aprior(s,buff,td->apr,m2);
  fprintf(s,"%s -> radius = %g%s",m1,td->radius,m2);
  fprintf(s,"%s -> noise = %g%s",m1,td->noise,m2);
  fprintf(s,"%s -> rsm_final_excite = %d%s",m1,td->rsm_final_excite,m2);
  fprintf(s,"%s -> use_steepest = %d%s",m1,td->use_steepest,m2);
  fprintf(s,"%s -> auto_compute_zones = %d%s",m1,td->auto_compute_zones,m2);
  fprintf(s,"%s -> good_zone = %g%s",m1,td->good_zone,m2);
  fprintf(s,"%s -> bad_zone = %g%s",m1,td->bad_zone,m2);
  fprintf(s,"%s -> explain = %d%s",m1,td->explain,m2);
  sprintf(buff,"%s -> ap",m1);
  fprintf_agg_params(s,buff,td->ap,m2);
}

void free_sofar(sofar *sf)
{
  free_atree(sf->at);
  free_arequest(sf->areq);
  free_brequest(sf->breq);
  free_dyv(sf->base);
  if (sf->lastdir) free_dyv(sf->lastdir);
  if (sf->workspace_free_fn && sf->workspace) 
    (sf->workspace_free_fn)(sf->workspace);
  AM_FREE(sf,sofar);
}

void free_twodexp(twodexp *td)
{
  td->free_twod_data_fn(td->twod_data);
  free_dyv(td->start_center);
  free_choice_params(td->cp);
  free_aprior(td->apr);
  free_optspec(td->os);
  AM_FREE(td,twodexp);
}

void random_dyv_in_unit_cube(dyv *x)
{
  int i;
  for ( i = 0 ; i < dyv_size(x) ; i++ )
    dyv_set(x,i,range_random(0.0,1.0));
}

double good_zone_thresh(
    double (*twod_fn)(char *data,dyv *sin,char *exp_code),
    char *twod_data,
    double high_value_proportion,
    int dims
  )
{
  dyv *sin = mk_dyv(dims); 
  int size = 5000;
  int i;
  dyv *scores = mk_dyv(size);
  double result;

  for ( i = 0 ; i < size ; i++ )
  {
    double v;
    random_dyv_in_unit_cube(sin);
    v = twod_fn(twod_data,sin,"");
    dyv_set(scores,i,v);
  }

  free_dyv(sin);

  dyv_sort(scores,scores);
  result = 
    dyv_ref(scores,(int) (dyv_size(scores) * (1.0 - high_value_proportion)));

  free_dyv(scores);
  return(real_min(result,0.975));
}

void init_twodres(twodres *tr)
{
  tr -> steps = 0;
}

void twodres_observe(
    twodres *tr,
    dyv *chosen_sin,
    double noiseless_sout,
    double sout,
    double value
  )
{
  if ( tr -> steps >= MAX_TWODRES_STEPS )
    my_error("MAX_TWODRES_STEPS too small");

  tr->e_choice_response[tr->steps] = noiseless_sout;
  tr->actual_response[tr->steps] = sout;
  tr->actual_value[tr->steps] = value;
  tr->steps += 1;
}

void fprintf_twodres(FILE *s,char *m1,twodres *tr,char *m2)
{
  char buff[100];
  sprintf(buff,"%s -> e_choice_response",m1);
  fprintf_realnums(s,buff,tr->e_choice_response,tr->steps,m2);
  sprintf(buff,"%s -> actual_response",m1);
  fprintf_realnums(s,buff,tr->actual_response,tr->steps,m2);
}

void get_twod_summary(
    twodres *tr,
    double good_zone,
    double bad_zone,
    twod_summary *tsum
  )
{
  double total_return = 0.0;
  double sum_actual_regret = 0.0;
  double sum_expected_regret = 0.0;
  double dis_counter = 0;
  double n = (double) int_max(1,tr->steps);
  int steps_to_good_zone;
  int i;
  double sum_final_ten = 0.0;

  for ( i = 0 ; i < tr -> steps ; i++ )
  {
    total_return += tr->actual_value[i];
    sum_actual_regret += 1 - tr->actual_response[i];
    sum_expected_regret += 1 - tr->e_choice_response[i];
    if ( i >= tr->steps - 10 )
      sum_final_ten += 1.0 - tr->e_choice_response[i];
    if ( tr->e_choice_response[i] < bad_zone ) dis_counter += 1;
  }

  for ( i = tr->steps-1 ; 
        i >= 0 && tr->e_choice_response[i] >= good_zone ; 
        i-- 
      )
  {
    /* skip */
  }

  steps_to_good_zone = i + 1;

  for ( i = 0 ; i < tr->steps && tr->e_choice_response[i] < good_zone ; i++ )
  {
    /* skip */
  }

  tsum -> total_return         = total_return;
  tsum -> avg_return           = total_return / n;
  tsum -> mean_actual_regret   = sum_actual_regret / n;
  tsum -> mean_expected_regret = sum_expected_regret / n;
  tsum -> final_ten_regret     = sum_final_ten / 
                                 (double) int_max(1,int_min(tr->steps,10));
  tsum -> num_disasters        = dis_counter;
  tsum -> steps_to_perm_goal   = ( steps_to_good_zone >= tr->steps ) ?
                                 tr->steps :
                                 steps_to_good_zone;
  tsum -> steps_before_goal    = i;
  tsum -> steps                = tr -> steps;
}
  
void explain_twodres(
    FILE *s,
    char *message,
    twodres *tr,
    double good_zone,
    double bad_zone
  )
{
  twod_summary ts[1];
  get_twod_summary(tr,good_zone,bad_zone,ts);

  fprintf(s,"After taking %d step%s, here's how the optimizer is doing:\n",
          ts->steps,(ts->steps==1)?"":"s"
         );
  fprintf(s,"Total actual return                     = %g\n",ts->total_return);
  fprintf(s,"Average actual return                   = %g\n",ts->avg_return);
  fprintf(s,"Mean actual regret                      = %g\n",
          ts->mean_actual_regret
         );
  fprintf(s,"Mean expected regret                    = %g\n",
          ts->mean_expected_regret
         );
  fprintf(s,"Mean regret of final ten experiments    = %g\n",
          ts->final_ten_regret
         );
  fprintf(s,"percentage of disasters                 = %4.1f%%\n",
          100.0 * ts->num_disasters / (double) int_max(1,ts->steps)
         );
  fprintf(s,"steps til base permanently in goal zone = ");
  if ( ts->steps_to_perm_goal >= ts->steps )
    fprintf(s,"not-happened-yet\n");
  else
    fprintf(s,"%d\n",ts -> steps_to_perm_goal);

  fprintf(s,"steps til base first entered goal zone = ");
  if ( ts->steps_before_goal >= ts->steps )
    fprintf(s,"not-happened-yet\n");
  else
    fprintf(s,"%d\n",ts -> steps_before_goal);
}
  
void explain_oneline_twodres(
    FILE *s,
    twodres *tr,
    twodexp *td,
    char *meth_name,
    int seed
  )
{
  twod_summary ts[1];

  get_twod_summary(tr,td->good_zone,td->bad_zone,ts);

  fprintf(s,"mth %18s seed %5d ",meth_name,seed);
  fprintf(s,"mxr %7.5f pd %4.1f stpg %3d sbg %3d ftr %6.3f\n",
          ts->mean_expected_regret,
          ((double)ts->num_disasters) * 100.0 / (double) int_max(1,ts->steps),
          ts->steps_to_perm_goal,
          ts->steps_before_goal,
          ts->final_ten_regret
         );
}

sofar *mk_simple_sofar(twodexp *td)
{
  sofar *sf = AM_MALLOC(sofar);
  sf -> at = NULL;
  sf -> base = mk_copy_dyv(td->start_center);
  sf -> radius = td -> radius;
  sf -> timecutoff = 0;
  sf -> lastdir = NULL;
  sf -> workspace = NULL;
  sf -> workspace_free_fn = NULL;
  return(sf);
}

sofar *mk_sofar(twodexp *td)
{
  sofar *sf;

  if ( eq_string(td->opt_type,"rsm") )
    sf = mk_simple_sofar(td);
  else if ( eq_string(td->opt_type,"pmax") )
    sf = mk_simple_sofar(td);
  else if ( eq_string(td->opt_type,"iemax") )
    sf = mk_simple_sofar(td);
  else if ( eq_string(td->opt_type,"comax") )
    sf = mk_simple_sofar(td);
  else if ( eq_string(td->opt_type,"kw"))
  {
    sf = mk_simple_sofar(td);
    sf->lastdir = mk_dyv(td->dims);
  }
  else if ( eq_string(td->opt_type,"oe0"))
    sf = mk_simple_sofar(td);
  else
  {
    sf = NULL;
    printf("mk_sofar: td->opt_type = %s\n",td->opt_type);
    my_error("unknown opt_type");
  }

  return(sf);
}

void fprintf_sofar(FILE *s,char *m1,sofar *sf,char *m2)
{
  char buff[100];
  sprintf(buff,"%s -> atree",m1);
  fprintf_atree(s,buff,sf->at,m2);
  sprintf(buff,"%s -> base",m1);
  fprintf_dyv(s,buff,sf->base,m2);
  fprintf(s,"%s -> radius = %g%s",m1,sf->radius,m2);
  fprintf(s,"%s -> workspace = %x%s",m1,(unsigned int)sf->workspace,m2);
  fprintf(s,"%s -> workspace_free_fn = %x%s",m1,(unsigned int)sf->workspace_free_fn,m2);
}

dyv *twod_rsm_choose(twodexp *td,sofar *sf)
/*
   sf gets updated.
  (i)  Chooses an experiment
  (ii) Performs it, and updates sf accordingly.

  This is only ever called when there is at least one result in sofar.
*/
{
  dyv *new_choice;
  posterior *po;
  bool move_center;
  norm_query *nq;
  norm_query_request *nqr = mk_norm_query_request();
  apost *apo;

/* we could've used mk_posterior_from_atree here (it does exactly these
 * operations), but we need the xtx for the call to mk_choice
 */
  nq = mk_norm_query_from_sin(sf->at,sf->breq,sf->base,nqr);
  apo = mk_apost(td->apr,nq->xtx,nq->xty,nq->yty);
  po = mk_posterior(apo,0,sf->at->sin_size,
		    sf->at->regdeg,sf->at->timemode);

  new_choice = mk_choice(po,
			 sf->at,
			 td->cp,
			 sf->base,
			 nq->xtx,
			 sf->breq,
			 sf->radius,
			 td->rsm_final_excite,
			 td->use_steepest,
			 td->explain,
			 &move_center
			 );
  
  if (Verbosity > 9.0) fprintf_posterior(stdout,"po",po,"");

  if ( move_center )
    copy_dyv(new_choice,sf->base);

  free_norm_query(nq);
  free_norm_query_request(nqr);
  free_apost_and_contents(apo);
  free_posterior(po);

  return(new_choice);
}

dyv *twod_aggress_choose(twodexp *td,sofar *sf)
/*
   sf gets updated.
  (i)  Chooses an experiment
  (ii) Performs it, and updates sf accordingly.

  This is only ever called when there is at least one result in sofar.
*/
{
  dyv *new_choice = mk_dyv(td->dims);
  double best_val = best_aggress_exp(td->apr,td->ap,sf,td->os,td->radius,
				     new_choice);
  if ( Verbosity >= 0.0 )
    printf("Using the %s type aggress method, the best value is %g\n",
           string_from_agg_type(td->ap->agg_type),
           best_val
          );

  copy_dyv(new_choice,sf->base);
  return(new_choice);
}

dyv *twod_choose(twodexp *td,sofar *sf)
/*
  Chooses an experiment. makes it and returns it.

  It is fine for sofar params to be updated by the choosing method here
  if necessary.

  This is only ever called when there is at least one result in sofar.
*/
{
  dyv *result = NULL;

  if ( eq_string(td->opt_type,"rsm") )
    result = twod_rsm_choose(td,sf);
  else if ( eq_string(td->opt_type,"pmax") ||
           eq_string(td->opt_type,"iemax") ||
           eq_string(td->opt_type,"comax") ||
	   eq_string(td->opt_type,"oe0")
         )
    result = twod_aggress_choose(td,sf);
  else if (eq_string(td->opt_type,"kw"))
  {
    result = mk_dyv(td->dims);
    keifWolfAS(sf->at,td->minradius,&sf->radius,&td->dradius,
	       td->fixed_dradius,sf->lastdir,0,result);
  }
  else
  {
    printf("td->opt_type = %s\n",td->opt_type);
    my_error("twod_cycle: unknown opt_type");
  }

  return(result);
}

void twod_rsm_observe(twodexp *td,dyv *sin,double sout,sofar *sf)
{
  dyv *inlo,*inhi,*outlo,*outhi,*local_sout;

  local_sout = mk_constant_dyv(1,sout);
  if (!sf->at)  /* make an atree with this point */
  {
    dym *local_dym_sin = mk_dym(1,sin->size);
    dym *local_dym_sout = mk_constant_dym(1,1,sout);
    copy_dyv_to_dym_row(sin,local_dym_sin,1);
    if (!sf->areq) my_error("twod_rsm_observe: no atree or arequest given\n");
    inlo = mk_constant_dyv(sin->size,0.0);
    inhi = mk_constant_dyv(sin->size,1.0);
    outlo = mk_constant_dyv(local_sout->size,0.0);
    outhi = mk_constant_dyv(local_sout->size,1.0);
 /* mk_atree wants unscaled parms, but its ok since we have identity scaling */
    sf->at = mk_atree(sf->areq,local_dym_sin,local_dym_sout,
		      inlo,inhi,outlo,outhi);
    free_dyv(inlo); free_dyv(inhi); free_dyv(outlo); free_dyv(outhi);
    free_dym(local_dym_sin); free_dym(local_dym_sout);
  }
  else /* just add the point */
    add_sin_and_sout_to_atree(sf->at,sin,local_sout);

  free_dyv(local_sout); 
}

void twod_observe(twodexp *td,dyv *sin,double sout,sofar *sf)
/*
   Sofar gets updated. Possibly in an opt_type-specific way.
*/
{
  if ( eq_string(td->opt_type,"rsm") )
    twod_rsm_observe(td,sin,sout,sf);
  else if ( eq_string(td->opt_type,"pmax") ||
           eq_string(td->opt_type,"iemax") ||
           eq_string(td->opt_type,"comax")
         )
    twod_rsm_observe(td,sin,sout,sf);  
        /* Note aggress just does the same as rsm */
  else if (eq_string(td->opt_type,"kw"))
    twod_rsm_observe(td,sin,sout,sf);
  else if (eq_string(td->opt_type,"oe0"))
  {
    dyv *lsout = mk_constant_dyv(1,sout);
    twod_rsm_observe(td,sin,sout,sf);
    oe0_update(sin,lsout,sf,td->os,td->apr);
    free_dyv(lsout);
  }
  else
  {
    printf("opt_type %s ?\n",td->opt_type);
    my_error("twod_observe: unknown opt_type");
  }
}

void twod_perform_and_record(
   twodexp *td,
   sofar *sf,
   dyv *sin,
   char *show_code,
   char *run_exp_code,
   twodres *tr
  )
{
  double noiseless_sout;
  double sout = sout;
  dyv *usin = mk_dyv(sin->size);
  dyv *usout = mk_dyv(1);
  dyv *lsout = mk_dyv(1);

  bool show = ((sf->at->num_points == 0) || 
	       (td->contour_display_steps < 2) ||
	       ((sf->at->num_points%td->contour_display_steps) == 1));

  if ( show && !char_is_in(show_code,'q') )
  {
    char buff[100];
    sprintf(buff,"are/%d.ps",tr->steps);
    ag_on(buff);
  }

  if ( Verbosity > 0.00001 )
  {
    fprintf_dyv(stdout,"base                ",sf->base,"\n");
    fprintf_dyv(stdout,"next_sin_to_be_tried",sin,"\n");
  }

  draw_auton_state(td->apr,
		   sf->at,
                   sf->base,
                   sf->radius,
                   sin,
                   TRUE,
                   show ? show_code : "",
		   sf->breq,
                   td -> twod_fn,
                   td -> twod_data,
                   td -> good_zone
                  );

  ag_off();
  if ( !char_is_in(show_code,'q') )
    wait_for_key();

  noiseless_sout = td->twod_fn(td->twod_data,sin,run_exp_code);
  sout = noiseless_sout + td->noise * gen_gauss();

  dyv_set(lsout,0,sout);
  usout_from_sout_sm(sf->at->sm,lsout,usout);
  usin_from_sin_sm(sf->at->sm,sin,usin);
  twod_observe(td,sin,sout,sf);
  twodres_observe(tr,sin,noiseless_sout,sout,
		  (td->os->oneshot_value)(usin,usout,td->os->oneshot_data));
  free_dyv(usout); free_dyv(usin); free_dyv(lsout);

  if ( Verbosity >= 0.0 )
  {
    int k;
    printf("Observe: sin ( ");
    for ( k = 0 ; k < dyv_size(sin) ; k++ )
      printf("%4.2f ",dyv_ref(sin,k));
    printf(") --> sout %6g [ noiseless = %6g ]\n",sout,noiseless_sout);
  }

  if ( !char_is_in(show_code,'q') )
  {
    explain_twodres(stdout,"current-twodres",tr,td->good_zone,td->bad_zone);
  }
  else
  {
    char c = 'a'+(char) int_max(0,int_min(25,(int)floor(26 * noiseless_sout)));
    if ( noiseless_sout < td -> bad_zone || noiseless_sout > td->good_zone )
      c = c - 'a' + 'A';
    printf("%c",c);
    fflush(stdout);
  }
}

void run_twod_experiment(twodexp *td, char *show_code, char *exp_code,
			 twodres *tr, sofar *sf)
{
  int iter,i;
  dyv *initial_chosen_sin = sf->base;

  if ((td->os->duration_type == FINITE_KNOWN) && 
      (td->os->duration_steps != td->steps))
    my_error("run_twod_experiment: os->duration->steps != td->steps\n");
  maybe_compute_twodexp_zones(td);
  td -> ap->min_resp = td->bad_zone;
  init_twodres(tr);

/* noiseless_sout and sout are same cause we don't know the noiseless value */
  if (sf->at)
  {
    dyv *usin = mk_dyv(sf->at->sin_size);
    dyv *usout = mk_dyv(sf->at->out_size);
    for (i=0;i<sf->at->num_points;i++)
    {
      usout_from_sout_sm(sf->at->sm,sf->at->darr->array[i]->sout,usout);
      usin_from_sin_sm(sf->at->sm,sf->at->darr->array[i]->sin,usin);
      twodres_observe(tr,sf->at->darr->array[i]->sin,
		      dyv_ref(sf->at->darr->array[i]->sout,0),
		      dyv_ref(sf->at->darr->array[i]->sout,0),
		      (td->os->oneshot_value)(usin,usout,td->os->oneshot_data));
    }
    free_dyv(usin); free_dyv(usout);
  }

  if ( Verbosity >= 0.01 )
  {
    fprintf_twodexp(stdout,"td",td,"\n");
    wait_for_key();
  }

  twod_perform_and_record(td,sf,initial_chosen_sin,show_code,exp_code,tr);

  for ( iter = 0 ; iter < td->steps ; iter++ )
  {
    dyv *chosen_sin = twod_choose(td,sf);
    twod_perform_and_record(td,sf,chosen_sin,show_code,exp_code,tr);
    free_dyv(chosen_sin);

    if ((td->auto_kw_steps>0)&&!((iter+1)%td->auto_kw_steps))
    {
      double kw = sf->breq->wght->kwidth;
      /* note: when output clips are added to rsm, the NULL below must be
	       replaced with the output_clips dym
      find_best_kwidth(sf->at,sf->breq,td->apr,NULL,SQUARED_LOO,
		       td->auto_kw_pct,&kw);
       */
      find_best_kwidth(sf->at,sf->breq,td->apr,SQUARED_LOO,td->auto_kw_pct,&kw);
      printf("Using new kwidth %8.3f\n",kw);
      sf->breq->wght->kwidth = kw;
      sf->radius = kw; /* for "rsm" : reasonable, but may overide users request
			* for others: does nothing, but modify the graphics */
    }

  }

  printf("\n");
  explain_twodres(stdout,"final",tr,td->good_zone,td->bad_zone);

  if ( !char_is_in(show_code,'q') )
  {
    int i;
    char expname[100];
    char fname[100];
    printf("\n\nI will now save some data and graphs in a file for you.\n");
    printf("I will save several files, with the same base name, but\n");
    printf("differeing suffixes.\n");
    (void) input_string("What base name shall I use? > ",expname,99);

    sprintf(fname,"%s-true.ps",expname);
    ag_on(fname);
    draw_auton_state(td->apr,
		     sf->at,
                     sf->base,
                     sf->radius,
                     sf->base,
                     TRUE,
                     "t",
		     sf->breq,
                     td -> twod_fn,
                     td -> twod_data,
                     td->good_zone
                    );
    ag_off();
    printf("Saved this diagram in %s\n",fname);
    wait_for_key();

    sprintf(fname,"%s-predict.ps",expname);
    ag_on(fname);

    draw_auton_state(td->apr,
		     sf->at,
                     sf->base,
                     sf->radius,
                     sf->base,
                     TRUE,
                     "f",
		     sf->breq,
                     td -> twod_fn,
                     td -> twod_data,
                     td->good_zone
                    );
    ag_off();
    printf("Saved this diagram in %s\n",fname);
    wait_for_key();
 
    sprintf(fname,"%s-optimistic.ps",expname);
    ag_on(fname);
    draw_auton_state(td->apr,
		     sf->at,
                     sf->base,
                     sf->radius,
                     sf->base,
                     TRUE,
                     "o",
		     sf->breq,
                     td -> twod_fn,
                     td -> twod_data,
                     td->good_zone
                    );
    ag_off();
    printf("Saved this diagram in %s\n",fname);
    wait_for_key();

    sprintf(fname,"%s-pessimistic.ps",expname);
    ag_on(fname);
    draw_auton_state(td->apr,
		     sf->at,
                     sf->base,
                     sf->radius,
                     sf->base,
                     TRUE,
                     "p",
		     sf->breq,
                     td -> twod_fn,
                     td -> twod_data,
                     td->good_zone
                    );
    ag_off();
    printf("Saved this diagram in %s\n",fname);
    wait_for_key();

    sprintf(fname,"%s-data.txt",expname);

    for ( i = 0 ; i < 2 ; i++ )
    { 
      FILE *s = (i == 0) ? stdout : safe_fopen(fname,"w");
      fprintf_twodexp(s,"td",td,"\n");
      fprintf(s,"\n\n");
      explain_twodres(s,expname,tr,td->good_zone,td->bad_zone);

      if ( i != 0 ) fclose(s);
    }

    printf("Saved this data in %s\n",fname);
    wait_for_key();
  }
}

dym *mk_scaled_dym_from_unscaled(dym *us,dyv *los, dyv *his)
{
  int i,j;
  dym *result = mk_dym(us->rows,us->cols);

  for (i=0;i<us->rows;i++)
    for (j=0;j<us->cols;j++)
      dym_set(result,i,j,(dym_ref(us,i,j) - dyv_ref(los,j))/
	                 (dyv_ref(his,j) - dyv_ref(los,j)));
  return result;
}

/* Optionally read a file of past points, then choose and run experiments */
void rsm_main(int argc, char *argv[])
{
  atree *at=NULL;
  dym *usins=NULL, *usouts=NULL;
  dyv *in_hi,*in_lo,*out_hi,*out_lo;
  dyv *start_center=NULL;
  char *fname = string_from_args("fname",argc,argv,NULL);
  char *format = string_from_args("fo",argc,argv,NULL);
  facode *fc=NULL;
  arequest *areq=NULL;
  brequest *breq;
  int dims;
  twodexp *td;
  sofar *sf;
  char *show_code = string_from_args("show_code",argc,argv,"f");
  char *exp_code = string_from_args("exp_code",argc,argv,"ad");
  twodres tr[1];

  if ( char_is_in(show_code,'q') && index_of_arg("verbose",argc,argv) < 0 )
    Verbosity = -1.0;

  if (fname)
  {
    read_io_dyms(fname,format,&usins,&usouts);
    if (usins && dym_cols(usins))
    {
      make_robust_lo_hi(usins,&in_lo,&in_hi);
      make_robust_lo_hi(usouts,&out_lo,&out_hi);
      dyv_set(in_lo,0,0.0); dyv_set(in_lo,1,0.0);
      dyv_set(in_hi,0,1.0); dyv_set(in_hi,1,1.0);
      dyv_set(out_lo,0,0.0); dyv_set(out_hi,0,1.0);
      fc = mk_facode_from_args(argc,argv,dym_cols(usins));
  /* for rsm we want ellipses and no neighbors to be the default local model */
      if (index_of_arg("gm",argc,argv) < 0) 
      {
	fc->regdeg = ELLIPSES_REGDEG;
	fc->neighs_code = 0;
      }
      areq = mk_arequest_from_facode(fc);
      at = mk_atree(areq,usins,usouts,in_lo,in_hi,out_lo,out_hi);
      start_center = mk_copy_dyv(at->darr->array[at->num_points-1]->sin);
    }
  }

  if (at) 
  {
    dims = at->sin_size;
    if (dims != int_from_args("dims",argc,argv,dims))
      my_error("rsm_main: command line dims not equal dims in given file");
  }
  else /* prepare an fc and arequest and make an empty tree */
  {
    dims = int_from_args("dims",argc,argv,2);
    fc = mk_facode_from_args(argc,argv,dims);
  /* for rsm we want ellipses and no neighbors to be the default local model */
    if (index_of_arg("gm",argc,argv) < 0) 
    {
      fc->regdeg = ELLIPSES_REGDEG;
      fc->neighs_code = 0;
    }
    areq = mk_arequest_from_facode(fc);
    usins = mk_dym(0,dims); usouts = mk_dym(0,1);
    in_lo = mk_constant_dyv(dims,0.0); in_hi = mk_constant_dyv(dims,1.0);
    out_lo = mk_constant_dyv(1,0.0); out_hi = mk_constant_dyv(1,1.0);
    at = mk_atree(areq,usins,usouts,in_lo,in_hi,out_lo,out_hi);
  }

  breq = mk_brequest_from_facode(fc);
/* make it easy to specify a new kwidth superceding the facode */
  breq->wght->kwidth = double_from_args("kwidth",argc,argv,breq->wght->kwidth);
  td = mk_default_twodexp(dims,areq);
  if (start_center) td->start_center = mk_copy_dyv(start_center);

  update_twodexp_from_args(td,areq,argc,argv);

  if (!eq_string(td->opt_type,"rsm")) /* just makes the graphics nicer */
    td->radius = breq->wght->kwidth;

  sf = mk_sofar(td);
  sf->at = at;
  sf->areq = areq;
  sf->breq = breq;

  run_twod_experiment(td,show_code,exp_code,tr,sf);

  free_twodexp(td);

  if (usins) free_dym(usins); if (usouts) free_dym(usouts);
  free_dyv(in_lo); free_dyv(in_hi); free_dyv(out_lo); free_dyv(out_hi);
  if (start_center) free_dyv(start_center);
  free_facode(fc);
/* don't free atree,arequest,brequest free_sofar will handle them all */
  free_sofar(sf);

  close_statistics();

  if ( !char_is_in(show_code,'q') )
  {
    am_malloc_report();
    dym_malloc_report();
  }

  return;
}