/*
 *  File : lib.c
 *
 *   Andrew H. Fagg
 *   5/10/91
 *
 *    Useful library routines for work in NSL.
 *    Includes :
 *       Matrix, Vector manipulation.
 *       Random weight generators.
 *       Distributed representation manipulation (creation, interpretation,
 *              addition with, etc.)
 *       
 *
 */

#include "nsl_include.h"
#include "alib.h"

#define nl << "\n"
#define sp << " "

/*
 *  nsl_matrix vec_mult_vec_trans(nsl_vector& vecA, nsl_vector& vecB)
 *
 *                  T
 *    Return = A * B
 *
 */

nsl_matrix vec_mult_vec_trans(NslFloat1& vecA, NslFloat1& vecB)

{
  int An, Bn;
  int i, j;

  An = vecA.get_xn();
  Bn = vecB.get_xn();

  nsl_matrix out(Bn, An);

  for(i = 0; i < An; ++i)
    for(j = 0; j < Bn; ++j)
      out[j][i] = vecA[i] * vecB[j];
  return out;
}

nsl_matrix vec_mult_vec_trans(NslFloat1& vecA, NslInt1& vecB)

{
  int An, Bn;
  int i, j;

  An = vecA.get_xn();
  Bn = vecB.get_xn();

  nsl_matrix out(Bn, An);

  for(i = 0; i < An; ++i)
    for(j = 0; j < Bn; ++j)
      out[j][i] = vecA[i] * vecB[j];
  return out;
}

/*
 *  nsl_vector mat_mult_col_vec(nsl_matrix& mat, nsl_vector& vec)
 *
 *    Computes C = M * V
 *
 *      Where C and V are both column vectors and M is a matrix.
 *
 */

nsl_vector mat_mult_col_vec(NslFloat2& mat, NslFloat1& vec)
{
  int xn;
  int yn;
  int vxn;
  int i, j;

  xn = mat.get_xn();
  yn = mat.get_yn();
  vxn = vec.get_xn();

  if(xn != vxn)
    {
      cout << "mat_mult_col_vec : dimension mismatch.\n";
      cout << xn sp << yn sp << vxn nl;
      exit(0);
    }

  nsl_vector out(yn);
  for(i = 0; i < yn; ++i)
    {
      out.elem(i) = 0.0;
      for(j = 0; j < xn; ++j)
	out.elem(i) += vec.elem(j) * mat.elem(j,i);
    }
      
  return out;
}
nsl_vector mat_mult_col_vec(NslFloat2& mat, NslInt1& vec)
{
  int xn;
  int yn;
  int vxn;
  int i, j;

  xn = mat.get_xn();
  yn = mat.get_yn();
  vxn = vec.get_xn();

  if(xn != vxn)
    {
      cout << "mat_mult_col_vec : dimension mismatch.\n";
      cout << xn sp << yn sp << vxn nl;
      exit(0);
    }

  nsl_vector out(yn);
  for(i = 0; i < yn; ++i)
    {
      out.elem(i) = 0.0;
      for(j = 0; j < xn; ++j)
	out.elem(i) += vec.elem(j) * mat.elem(j,i);
    }
      
  return out;
}


/*
 *  nsl_vector mat_trans_mult_col_vec(nsl_matrix& mat, nsl_vector& vec)
 *
 *                  T
 *    Computes C = M * V
 *
 *    Where C and V are column vectors.
 *
 */

nsl_vector mat_trans_mult_col_vec(NslFloat2& mat, NslInt1& vec)

{
  int xn;
  int yn;
  int vxn;
  int i, j;

  xn = mat.get_xn();
  yn = mat.get_yn();
  vxn = vec.get_xn();

  if(yn != vxn)
    {
      cout << "mat_mult_col_vec : dimension mismatch.\n";
      cout << xn sp << yn sp << vxn nl;
      exit(0);
    }

  nsl_vector out(xn);
  for(i = 0; i < xn; ++i)
    {
      out.elem(i) = 0.0;
      for(j = 0; j < yn; ++j)
	out.elem(i) += vec.elem(j) * mat.elem(i,j);
    }
      
  return out;
}



/*
 *  void extract_vector(nsl_matrix& mat, int col, nsl_vector& vec)
 *
 *   Copy the <col>th column of <mat> into <vec>.
 *
 */

void extract_col_vector(NslInt2& mat, int col, NslInt1& vec)
{
  int row;
  int xn, yn;

  xn = vec.get_xn();
  yn = mat.get_yn();
  if(xn != yn)
    {
      cout << "Error\n";
      exit(0);
    }
  for(row = 0; row < xn; ++row)
    {
      vec.elem(row) = mat.elem(row, col);
    }
}


/*
 *  void extract_row_vector(nsl_matrix& mat, int row, nsl_vector& vec)
 *
 *   Copy the <row>th row of <mat> into <vec>.
 *
 */

void extract_row_vector(NslInt2& mat, int row, NslInt1& vec)
{
  int col;
  int xn, yn;

  xn = vec.get_xn();
  yn = mat.get_yn();
  if(xn != yn)
    {
      cout << "Error\n";
      exit(0);
    }
  for(col = 0; col < xn; ++col)
    {
      vec.elem(col) = mat.elem(row, col);
    }
}



/*
 *  void insert_vector(nsl_matrix& mat, int col, nsl_vector& vec)
 *
 *   Copy <vec> into the <col>th column of <mat>.
 *
 */

void insert_col_vector(NslInt2& mat, int col, NslInt1& vec)
{
  int row;
  int xn, yn;

  xn = vec.get_xn();
  yn = mat.get_yn();
  if(xn != yn)
    {
      cout << "Error\n";
      exit(0);
    }
  for(row = 0; row < xn; ++row)
    {
      mat.elem(row, col) = vec.elem(row);
    }
}


/*
 *  void insert_row_vector(nsl_matrix& mat, int row, nsl_vector& vec)
 *
 *   Copy <vec> into the <row>th row of <mat>.
 *
 */

void insert_row_vector(NslInt2& mat, int row, NslInt1& vec)
{
  int col;
  int xn, yn;

  xn = vec.get_xn();
  yn = mat.get_yn();
  if(xn != yn)
    {
      cout << "Error\n";
      exit(0);
    }
  for(col = 0; col < xn; ++col)
    {
      mat.elem(row, col) = vec.elem(col);
    }
}



/*
 * float random_value()
 *
 *  Return = a random number between [-1.0, 1.0].
 *
 */

float random_value()
{
  return(((float)(rand() % 30001)/ 15000.0) - 1.0);
}


/*
 * float random()
 *
 *  Return = a random number between [0.0, 1.0).
 *
 */

float randomv()
{
  return((float)(rand() % 30000)/ 30000.0);
}


/*
 *  nsl_vector distributed_sum(nsl_vector& a, nsl_vector& b)
 *
 *   Given distributed representations of values stored in vectors
 *  <a> and <b>, compute the distributed sum of the two.  The resulting
 *  vector length is the sum of the lengths of both less 1, with the
 *  lowest index (represented by the 0th position), being the sum
 *  of the lowest indices of each.
 *
 *
 */

nsl_vector distributed_sum(nsl_vector& a, nsl_vector& b)

{
  int alen;
  int blen;
  int i, j;

  alen = a.get_xn();
  blen = b.get_xn();

  nsl_vector result(alen + blen - 1);

  result = 0.0;
  for(i = 0; i < alen; ++i)
    {
      for(j = 0; j < blen; ++j)
	{
	  result.elem(i+j) = result.elem(i+j) + a.elem(i) * b.elem(j);
	}
    }
  return(result);
}




/*
 *  nsl_vector biased_distributed_sum(nsl_vector& a, nsl_vector& b)
 *
 *   Given distributed representations of values stored in vectors
 *  <a> and <b>, compute the distributed sum of the two.  The resulting
 *  vector length is the sum of the lengths of both less 1, with the
 *  lowest index (represented by the 0th position), being the sum
 *  of the lowest indices of each.
 *
 *
 */

nsl_vector biased_distributed_sum(nsl_vector& a, nsl_vector& b,
				  nsl_matrix& bias, int size)

{
  int alen;
  int blen;
  int i, j;

  alen = a.get_xn();
  blen = b.get_xn();

           /*  Check dimension of bias vector.  */

  if(alen != bias.get_xn() || blen != bias.get_yn())
    {
      cout << "biased_distributed_sum : bias vector size does not match input, "
	<< bias.get_name() sp << alen sp << blen sp << bias.get_xn() sp
	  << bias.get_yn() nl;
      exit(0);
    }

  nsl_vector result(size);


  result = 0.0;
  for(i = 0; i < alen; ++i)
    {
      for(j = 0; j < blen; ++j)
	{
           /*  Check bias to make sure it is in range of output.  */

	  if((int)bias.elem(i,j) < 0 || (int)bias.elem(i,j) >= size)
	    {
	      cout << "biased_distributed_sum : bias vector element error, "
		<< bias.get_name() sp << (int)bias.elem(i,j) nl;
	    }
	  result.elem((int)bias.elem(i,j)) =
	    result.elem((int)bias.elem(i,j)) + a.elem(i) * b.elem(j);
	}
    }
  return(result);
}


/*
 *  void gauss_to_value(nsl_vector& a, float& total, float& div)
 *
 *    Convert a distributed representation into the weighted sum
 *  (<total>) and activity (<div>) values.  The coded value is
 *  then just <total>/<div>. 
 *  Computes the location relative to the array indicies.
 *
 */

void gauss_to_value(nsl_vector& a, float& total, float& div)
{
  int alen;

  int i;

  alen = a.get_xn();

  total = div = 0.0;

  for(i = 0; i < alen; ++i)
    {
      total += i * a.elem(i);
      div += a.elem(i);
    }
}


/*
 *  void gauss_to_value(float min, float max, int pad,
 *                  nsl_vector& vec, float& total, float& div)
 *
 *    Convert a distributed representation into the weighted sum
 *  (<total>) and activity (<div>) values.  The coded value is
 *  then just <total>/<div>.  The value is computed relative
 *  to the <min> and <max> values - <min> corresponds to a
 *  a gaussian centered at array position <pad>, and <max>
 *  corresponds to an array position <array_size> - <pad> - 1.
 *  
 *
 */

void gauss_to_value(float min, float max, int pad,
		      nsl_vector& vec, float& total, float& div)

{
  float step;
  int i;
  float current;

  int size = vec.get_xn();
  step = (max - min) / ((float) (size - pad * 2 - 1.0));
  for(total = div = 0.0, current = min - (float) pad * step, i = 0; i < size;
      ++i, current += step)
    {
      total += current * vec.elem(i);
      div += vec.elem(i);
    }
}



/*
 *  float gauss(float mean, float position, float std)
 *
 *    Return the value of a gaussian at <position>.  The
 *  gaussian has a mean of <mean> and a standard deviation of <std>.
 *
 */

float gauss(float mean, float position, float std)

{
  return((1.0/(NEWCONST * std)) * Exp(-0.5 * Square((position - mean) / std)));
}


/*
 *  nsl_vector value_to_gauss(float val, float min, float max, 
 *	       int pad, float std, nsl_vector& vec) *
 *
 *    Compute a distributed representation for <val>, in the
 *  shape of a gaussian centered about <val>.  <val> must
 *  fall between <min> and <max>.  At <val> = <min>, the
 *  output is a gaussian centered about the <pad>th
 *  cell in the array.  When <val> = <max>, the gaussian
 *  is centered about element number <array_size> - <pad> - 1.
 *
 */

nsl_vector value_to_gauss(float val, float min, float max, 
	       int pad, float std, nsl_vector& vec)

{
  int size;
  int i;
  float step;
  float current;
  float gain;

  size = vec.get_xn();
  
  gain = 1.0 / gauss(0.0, 0.0, std);
  step = (max - min) / ((float) (size - pad * 2 - 1.0));
  for(current = min - (float) pad * step, i = 0; i < size;
      ++i, current += step)
    vec.elem(i) = gauss(val, current, std) * gain;
  return(vec);
}




/*
 *   Other usefull routines.
 *
 */

nsl_vector square(nsl_vector& a)
{
  int i;
  int j;

  i = a.get_xn();

  nsl_vector out(i);

  for(j = 0; j < i; ++j)
    {
      out.elem(j) = a.elem(j) * a.elem(j);
    }
  return(out);
}


nsl_vector sqrt(nsl_vector& a)
{
  int i;
  int j;

  i = a.get_xn();

  nsl_vector out(i);

  for(j = 0; j < i; ++j)
    {
      out.elem(j) = sqrt(a.elem(j));
    }
  return(out);
}


nsl_vector abs(nsl_vector& a)
{
  int i;
  int len;

  len = a.get_xn();

  nsl_vector out(len);

  for(i = 0; i < len; ++i)
    out.elem(i) = fabs(a.elem(i));
  return(out);
}


/*
 *  void randomize(nsl_vector& vec)
 *
 *    Set the values of <vec> to random numbers in the range [-1.0, 1.0].
 *
 */

void randomize(nsl_vector& vec)

{
  int i;
  int init;
  int limit;

  init = vec.get_x0();
  limit = vec.get_x1();
  for(i = init; i <= limit; ++i)
    {
      vec.elem(i) = random_value();
    }
}


/*
 *  overload void randomize(nsl_matrix& mat)
 *
 *    Set the values of <mat> to random numbers in the range [-1.0, 1.0].
 *
 */

void randomize(nsl_matrix& mat)

{
  int i, j;
  int x0, x1;
  int y0, y1;

  x0 = mat.get_x0();
  x1 = mat.get_x1();
  y0 = mat.get_y0();
  y1 = mat.get_y1();

  for(i = x0; i <= x1; ++i)
    for(j = y0; j <= y1; ++j)
      {
	mat.elem(i,j) = random_value();
      }
}