/* IO(1) */

load x9_10;

/*
 * Reference: J. K. Wolf, A. M. Viterbi, and G. S. Dixon,
 * Finding the Best Set of K Paths Through a Trellis
 * with Application to Multitarget Tracking,
 * IEEE AES-25, No. 2, March 1989, pp. 287-296
 */

void get_x( int _K)
{
/* data setup for simple x/y/velocity targets */

  int i, j;
  vector jj;
  matrix ppp;

  if( narg == 0) {

    LOADED=9;

    N = 21;	/* N data samples */
    M = 4;	/* M symbol values */
    K = 1;	/* find K best completely unmerged paths */

    MM = M;	/* save original M, used when K > 1 */

    if( undef(example_num))
      example_num = 1;

/*** note: dt is global and is used in get_cost() below ***/

    dt = 4/(N-1);	/* time will go from 0 to 4 */
     t = [0:N-1]*dt;

  if( example_num == 1) {

    xy1 = [ t; 2*(t.*(1-t/4))]';
     v1 = [ ones(1,N); 2-t]';

    xy2 = [ t; 4-2*(t.*(1-t/4))]';
     v2 = [ ones(1,N); t-2]';

    /* some random clutter samples, uniform over the same ranges as xy1, v1, xy2, v2 */


if(0) /*** testing with clutter as HUGE_VAL's ***/

  xyr1 = vr1 = xyr2 = vr2 = HUGE_VAL * ones( rows(xy1), cols(xy1));

else {
    xyr1 = 4 * rand( rows(xy1), cols(xy1));
     vr1 = 4 * rand( rows(v1), cols(v1)) - 2;

    xyr2 = 4 * rand( rows(xy1), cols(xy1));
     vr2 = 4 * rand( rows(v1), cols(v1)) - 2;
}

  } else {

    print( "bad example_num", example_num);
    return;

  }

/*
 * xv = [ xy1(1) v1(1); xy2(1) v2(1);  r  r  r  r;  r  r  r  r;	-- t = 1
 *	  xy1(2) v1(2); xy2(2) v2(2);  r  r  r  r;  r  r  r  r;	-- t = 2
 *	  ...
 *	  xy1(N) v1(N); xy2(N) v2(N);  r  r  r  r;  r  r  r  r;	-- t = N
 *	]
 */
    xv = reshape( [xy1 v1 xy2 v2 xyr1 vr1 xyr2 vr2], M*N, 4);

/* add some noise */

/* 0.01 stdev was used for the Asilomar summary */

   Kn = 0.0001;	/* measurement noise variance */

   xv += sqrt(Kn) * rand( rows(xv), cols(xv), "gauss");

/* just add noise to the x/y positions, not velocities */
/*
   xv(:,[1 2]) = xv(:,[1 2]) + 0.2 * rand( rows(xv), 2, "gauss");
*/

/* randomize the observation rows */

/*
    ppp = permutations(1:M);

    for( i = 1; i <= N; ++i) {
      j = 1 + (i-1) * M;
      jj = [j:j+M-1];
      xv(jj) = xv( jj( ppp( rand(1,1,"int",M))));
    }
*/

    zstate = 0;

    z = zeros(1,N);	/* dummy */

    return;
  }

  /* narg > 0 */

  if( LOADED != 9) {
    print( "get_x: call get_x with no args first.");
    return;
  }

  if( _K > MM) {
    printf( "get_x: K > M not allowed!");
    return;
  }

  K = _K;
  combos = combinations( 1:MM, K);
  perms = permutations( 1:K);
  M = rows(combos);
  MP = rows(perms);

  zstate = 0;
}
