/*sequence.c - Generate a random sequence of stimuli on alternating sides, so
  that stimuli on the same side will be far enough apart to disambiguate
  responses.
  Copyright (c) 1996 Matthew Belmonte

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  If you find this program useful, please send mail to Matthew Belmonte.
  <mkb4@Cornell.edu>.  If you base a publication on data processed by this
  program, please notify Matthew Belmonte and include the following citation
  in your publication:

	Matthew Belmonte, `A Software System for Analysis of
	Steady-State Evoked Potentials', Association for Computing
	Machinery SIGBIO Newsletter 17:1:9-14 (April 1997).
*/

#include <stdio.h>
#include <errno.h>
#include <math.h>
#include "fastshft.h"

#define H ((ST_SLOWEST_RESPONSE-ST_FASTEST_RESPONSE+ST_MSECS_PER_TICK/2)/ST_MSECS_PER_TICK-1)
#define K (ST_LONGEST_ISI/ST_MSECS_PER_TICK)

/*Find probabilities and cumulative expectations on successive iterations 'i' of
  all elements 't(i)' in a distribution {p[0], p[1], ..., p[H-1], p[H]/(K-H+1),
  p[H]/(K-H+1), ..., p[H]/(K-H+1)} over [0, K] conditioned to t(i+1)-t(i) >= H.

  Since the current selection t(i+1) is restricted to the interval
  [max(0, H-t(i)), K], t(i+1) has no chance of being 0 unless t(i) >= H.  Since
  the die is fair, the probability that t(i) >= H is simply
  (K-H+1)*(p[H]/(K-H+1)), that is, the number of possible values of t(i) that
  are greater than or equal to H, times the probability, for each of these
  values X, that t(i)=X.  Cancelling the K-H+1, we're left with p[H].

  If, based on the foregoing computation, t(i+1) is allowed to be 0, then 0 is
  one of K+1 possible values that t(i+1) can assume.  So the _conditional_
  probability that t(i+1)=0 is 1/(K+1) (since the die is fair), and that makes
  the probability that t(i+1)=0 p[H]/(K+1).  This is exactly the computation
  implemented below for p[0].

  Whenever t(i+1) is allowed to be 0, it's also allowed to be 1.  So the
  probability that t(i+1)=1 must be at least p[H]/(K+1).  However, t(i+1) is
  also allowed to be 1 when t(i)=H-1.  The probability of this is p[H-1].  If it
  is the case that t(i)=H-1, then t(i+1) is restricted to [1, K].  So the
  _conditional_ probability that t(i+1)=1 is 1/K, since there are now only K
  (instead of K+1 as above) values from which to choose (and, again, since the
  die is fair).  This makes the probability that t(i+1)=1 p[H]/(K+1)+p[H-1]/K
  (since the two events t(i) >= H and t(i)=H-1 are independent).  Note that the
  first term in this sum is exactly what was just computed as the new value of
  p[0].

  In general, the new value of p[n] is the sum of the new value of p[n-1], and
  1/(K+1-n) times the old value of p[H-n].  What this is implementing, in
  shorthand, is the iteration of a Markov chain with the following (K+1)*(K+1)
  transition probability matrix (illustrated here for the case H=3, K=4):

  0		0		0		1/(K-H+1)	1/(K-H+1)
  0		0		1/(K-(H-1)+1)	1/(K-(H-1)+1)	1/(K-(H-1)+1)
  0		1/(K-(H-2)+1)	1/(K-(H-2)+1)	1/(K-(H-2)+1)	1/(K-(H-2)+1)
  1/(K+1)	1/(K+1)		1/(K+1)		1/(K+1)		1/(K+1)
  1/(K+1)	1/(K+1)		1/(K+1)		1/(K+1)		1/(K+1)

  The initial probability vector is uniform, that is,
  [1/(K+1)	1/(K+1)		1/(K+1)		1/(K+1)		1/(K+1)]

  Thanks to Bruce Fischer for explaining conditioned distributions, and to
  Anthony Gamst for explaining Markov chains.
*/
int iter_probs(desired_duration_of_trial, smallest_expectation)
int desired_duration_of_trial;
double *smallest_expectation;
  {
  register int elt, iter;
  double p[H+1], prev_p[H+1], e[H+1], expected_duration;
/*initialise for a fair die*/
  printf("ITERATION 0\n");
  for(elt = 0; elt <= H; elt++)
    {
    p[elt] = 1.0/(K+1);
    e[elt] = p[elt];
    printf("%lf %lf\n", p[elt], e[elt]);
    }
  p[H] *= K-H+1;
  *smallest_expectation = e[0];
  expected_duration = (K+1)/2.0;
  iter = 1;
  while(expected_duration < desired_duration_of_trial)
  /*inv: p[0..H-1] are the probabilities of obtaining the values 0, 1, ..., H-1,
    respectively, on the 'iter'th selection, and p[H] is K-H+1 times the
    probability of obtaining H (which is equal to the probability of obtaining
    H+1, and to that of obtaining H+2, and so on up to K) on the 'iter'th
    selection.  (In practice, this computation converges sufficiently within
    about 10 iterations.)  e[0..H] are the expected counts of each value 0, 1,
    ..., H, respectively, after 'iter' selections have been made.  (e[H] is
    actually the expectation for _all_ values in the interval [H, K].)
    '*smallest_expectation' is the least element of e[], and 'expected_duration'
    is the expected duration of the trial so far, that is,

       H-1
      -----
      \
    ( /    (e[elt]*elt) ) + (K-H+1)*e[H]*H + (iter+1)/2
      -----
      elt=0

    The final term in the expression above is necessary in order to account for
    the i+1 target stimuli that occur at the beginnings and ends of the i
    sequences of nontarget stimuli that are being generated here.  Each target
    is half a time period in duration.
*/
    {
    for(elt = 0; elt <= H; elt++)
      prev_p[elt] = p[elt];
    p[0] = prev_p[H]/(K+1);
    e[0] += p[0];
    *smallest_expectation = e[0];
    printf("\nITERATION %d\n%lf %lf\n", iter++, p[0], e[0]);
    expected_duration = (iter+1)/2.0;
    for(elt = 1; elt <= H; elt++)
      {
      p[elt] = p[elt-1]+prev_p[H-elt]/(K+1-elt);
      e[elt] += p[elt];
      if(e[elt] < *smallest_expectation)
	*smallest_expectation = e[elt];
      expected_duration += e[elt]*elt;
      printf("%lf %lf\n", p[elt], e[elt]);
      }
    p[H] *= K-H+1;
    expected_duration += e[H]*(K-H);
    }
/*'iter' selections, with initially uniform probabilities, produce an expected
  count of at least '*smallest_expectation' for each possible outcome.  The
  expected duration of an experimental trial based on this sequence is not much
  more than 'desired_duration_of_trial'.
*/
  return(iter);
  }

/*Produce a pseudorandom integer from a uniform distribution on [lb, ub].*/
int choose(lb, ub)
int lb, ub;
  {
  return(lb+(ub-lb+1)*(random()&((1L<<26)-1))/(1L<<26));
  }

void main()
  {
  register int tick, isi, side, stimulus, trial;
  int num_stimuli, num_trials;
  double instances_per_trial;
  FILE *stim_file;
  fprintf(stderr, "Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  num_stimuli = iter_probs((ST_MSECS_PER_TRIAL-2*ST_LONGEST_ISI)/ST_MSECS_PER_TICK, &instances_per_trial);
  num_trials = (int)((AV_BIN_SIZE+0.5*instances_per_trial)/instances_per_trial);
  printf("Generating %d trials each containing %d targets.\n", num_trials, num_stimuli);
  if((stim_file = fopen("stimuli.dat", "w")) == NULL)
    {
    perror((char *)0);
    exit(errno);
    }
  for(trial = 0; trial != num_trials; trial++)
  /*inv: 'trial' trials have been output*/
    {
  /*begin with an instance of the longest possible sequence of nontargets*/
    for(tick = 0; tick != K; tick++)
      fprintf(stim_file, "1\n3\n");
  /*flip a coin in order to decide on which side the first target will appear*/
    side = 1+(random()&2);
  /*if first target is on the right then insert an extra nontarget on the left*/
    if(side == 3)
      fprintf(stim_file, "1\n");
  /*initialise 'isi' s.t. the interval for the initial choice is unrestricted*/
    isi = K;
    for(stimulus = 0; stimulus != num_stimuli; stimulus++)
    /*inv: 'stimulus' target stimuli have been generated, and each of these is
      followed immediately by t(i) pairs of nontarget stimuli, where t(i) is
      selected uniformly over over [0, K] conditioned to t(i+1)-t(i) >= H.*/
      {
    /*issue target*/
      fprintf(stim_file, "%d\n", side+1);
    /*choose inter-target interval*/
      isi = choose((H-isi>0? H-isi: 0), K);
    /*toggle side*/
      side = (side+2)%4;
    /*issue nontargets*/
      for(tick = 0; tick != isi; tick++)
	fprintf(stim_file, "%d\n%d\n", side, (side+2)%4);
      }
  /*issue the final target*/
    fprintf(stim_file, "%d\n", side+1);
  /*end with an instance of the longest possible sequence of nontargets*/
    for(tick = 0; tick != K; tick++)
      fprintf(stim_file, "%d\n%d\n", (side+2)%4, side);
  /*issue an end-of-trial code*/
    fprintf(stim_file, "254\n");
    }
  fclose(stim_file);
  exit(0);
  }
