Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!uunet!mcsun!hp4nl!fwi.uva.nl!smagt
From: smagt@fwi.uva.nl (Patrick van der Smagt)
Newsgroups: comp.ai.neural-nets
Subject: ART1&2 source code (WAS: Detailed description of ART-1 & 2 wanted)
Keywords: ART NN
Message-ID: <1991Feb27.153904.10455@fwi.uva.nl>
Date: 27 Feb 91 15:39:04 GMT
References: <132@blakex.RAIDERNET.COM>
Sender: news@fwi.uva.nl
Distribution: comp
Organization: FWI, University of Amsterdam
Lines: 998
Nntp-Posting-Host: wendy.fwi.uva.nl

root@blakex.RAIDERNET.COM (Blake McBride) writes:

>Am looking for a source for a detailed description of Adaptive Resonance
>Theory (ART) by Grossberg & Carpenter.  Both 1 & 2.  In addition, if 
>anyone has any source it sure would be appreciated!

Well, I do have a description of ART-1, whereas one of my students
tried to disentangle ART itself.  But one had better read Lippmann
and reconstruct ART-1 from that.  Anyway, here is some source code
for ART-1 and ART-2; I wrote neither but looked at both.  I also
have another ART-2 source code, but I will not clobber the net
by including that one, too.
						Patrick van der Smagt 
								    /\/\
                                                                    \  /
Organisation: Faculty of Mathematics & Computer Science             /  \
              University of Amsterdam, Kruislaan 403,            _  \/\/  _
              NL-1098 SJ  Amsterdam, The Netherlands            | |      | |
Phone:        +31 20  525 7524                                  | | /\/\ | |
Fax:          +31 20  525 7490                                  | | \  / | |
                                                                | | /  \ | |
email:        smagt@fwi.uva.nl                                  | | \/\/ | |
                                                                | \______/ |
                                                                 \________/

								    /\/\
``The opinions expressed herein are the author's only and do        \  /
not necessarily reflect those of the University of Amsterdam.''     /  \
                                                                    \/\/

---------art1.c-------------cut here-----------------art1.c-----------------
/* Simiulation of ART-1 Neural network (Ref: Lippman 87)
   Written by Jagath Samarabandu <bandu@cs.buffalo.edu> 
   To compile, type cc art1.c -o art1  */

#include <stdio.h>

#define N      25   /* No. of nodes in F1 */
#define M      20   /* No. of nodes in F2 */

initialize(top_down,bot_up,n,m) /* initialize the LTM traces */
     double top_down[],bot_up[]; /* top down and bot. up LTM traces */
     int n,m;                   /* No. of nodes in F1 and F2 */
{
  int i;
  double t;

  t = 1.0/(1.0+n);
  for (i=0;i<n*m;i++) {
    top_down[i] = 1.0; bot_up[i] = t;
  }
}

compute_matching_scores(mu,bot_up,input,n,m) /* calculate matching scores */
     double mu[],bot_up[],input[];            /* mu - matching score */
     int n,m;                                /* No. of F1 and F2 nodes */
{
  int i,j;

  for (j=0; j<m; j++) 
    for (i=0, mu[j]=0.0; i<n; i++)
      mu[j] += bot_up[i*m+j]*input[i];
}

double vigilance(top_down,input,jmax,n,m) /* returns |T.X|/|X| */
     double top_down[],input[];
     int n,m,jmax;
{
  int i;
  double x,t;

  for (i=0,x=0.0; i<n; i++)
    x += input[i];
  for (i=0,t=0.0; i<n; i++)
    t += top_down[i*m+jmax]*input[i];
  return(t/x);
}

int find_max(array,len) /* find the max of array and return the index */
     double array[];
     int len;
{
  int j,jmax;

  for (j=0,jmax=0; j<len; j++)
    if (array[j]>array[jmax])
      jmax = j;
  return (jmax);
}

adapt_LTM(top_down,bot_up, input,jmax,n,m) /* change top down and bot.up LTM */
     double top_down[],bot_up[],input[];
     int n,m,jmax;
{
  int i,ij;
  double sum,t;
  
  for (i=0,sum=0.5; i<n; i++)
    sum += top_down[i*m+jmax]*input[i];

  for (i=0,ij=jmax; i<n; i++,ij += m) {
    t = top_down[ij]*input[i];
    bot_up[ij] = t/sum;
    top_down[ij] = t;
  }
}

load_data(d,max,n,fname) /* load data from file */
     int d[],max,n;
     char *fname[];
{
  FILE *fp;
  int n_pat,n_var,i,j;
  
  if (!(fp=fopen(fname,"r"))) exit(perror(fname));
  fscanf(fp,"%d %d\n",&n_pat,&n_var);
  if (n_pat>max) {
    printf("Warning! only %2d patterns out of %d are read\n",max,n_pat);
    n_pat = max;
  }
  if (n_var!=n)
    exit(printf("wrong pattern size: should be %2d. was %2d\n",n,n_var)); 
    
  for (i=0;j<n_pat;i++)
    for (j=0;j<n_var;j++)
      fscanf(fp,"%d",&d[i*n+j]);
  fclose(fp);
}

display(in,top_down,x,y,m) /* display input and top down weights */
     double in[],top_down[];
     int x,y,m;
{
  int i,ix,iy,j;

  for (iy=0,i=0; iy<y; iy++){
    for (ix=0,i=iy*y; ix<x; ix++,i++)
      printf("%c",(in[i]>0.5)?'#':' ');
    printf(" | ");
    for (j=0; j<m; j++) {
      for (ix=0,i=iy*y; ix<x; ix++,i++)
	printf("%c",(top_down[i*m+j]>0.5)?'#':' ');
      printf(" ");
    }
    printf("\n");
  }
  printf("\n");
}
/*****************  main routine starts here  *******************************/

int data[20][N]={
#include "art1.data"
};

main(argc,argv)
     int argc;
     char *argv[];
{
char ch;
  double t[N][M],b[N][M],x[N],mu[M],rho;
  int max_j,i,j,n_pat,ok,seq[M];

  if (argc>1)
    n_pat = load_data(data,20,N,argv[1]);
  else n_pat=20;
  initialize(t,b,N,M);
  printf("Vigilance threshold: "); scanf("%lf",&rho);
  printf("\nSimulation of ART1 network with vigilance Threshold = %3.3lf\n\n",rho);

  do {
    for (i=0; i<n_pat; i++) {
      for (j=0; j<N; j++)
	x[j] = (double)data[i][j];
      compute_matching_scores(mu,b,x,N,M);
      bzero((char *)seq,M*sizeof(int)); j=1;
      do {
	max_j = find_max(mu,M); seq[max_j] = j++;
	if (vigilance(t,x,max_j,N,M)>rho) {
	  adapt_LTM(t,b,x,max_j,N,M);
	  seq[max_j] = -1;
	  break;
	}
	else 
	  mu[max_j] = 0.0;
      } while (1);
      printf("IN:%2d    ",i);
      for (j=0;j<M; j++) {
	if (seq[j]>0)
	  printf("%1d     ",seq[j]);
	else if (seq[j]==0)
	  printf("      ");
	else {
	  printf("R\n"); break;
	}
      }
      display(x,t,5,5,M);
    }
    printf("Another round? (1-yes): "); scanf("%d",&ok);
  } while (ok);
  ok =1;

do {
   
   printf("\n[A-T] or ESC to quit :"); scanf("%s",&ch);
   printf("\n");
   if ((ch == 27) || ((ch < 'A') || (ch > 'T'))) ok =0;
   
   else
      {
      for (j=0; j<N; j++)
	x[j] = (double)data[ch -'A'][j];
      compute_matching_scores(mu,b,x,N,M);
      bzero((char *)seq,M*sizeof(int)); j=1;
      do {
	max_j = find_max(mu,M); seq[max_j] = j++;
	if (vigilance(t,x,max_j,N,M)>rho) {
	  adapt_LTM(t,b,x,max_j,N,M);
	  seq[max_j] = -1;
	  break;
	}
	else 
	  mu[max_j] = 0.0;
      } while (1);
      printf("IN:%2d    ",i++);
      for (j=0;j<M; j++) {
	if (seq[j]>0)
	  printf("%1d     ",seq[j]);
	else if (seq[j]==0)
	  printf("      ");
	else {
	  printf("R\n"); break;
	}
      }
      display(x,t,5,5,M);

 }
 } while (ok);
}
----------------------end-------------end------------------------------------
------art1 data file--------cut here-------------art1 data file--------------
/* art1 data file - 20x25 */

{1,1,1,1,1,
 1,0,0,0,1,
 1,1,1,1,1,
 1,0,0,0,1,
 1,0,0,0,1,},

{1,1,1,1,0,
 1,0,0,0,1,
 1,1,1,1,0,
 1,0,0,0,1,
 1,1,1,1,0,},

{1,1,1,1,1,
 1,0,0,0,0,
 1,0,0,0,0,
 1,0,0,0,0,
 1,1,1,1,1,},

{1,1,1,1,0,
 1,0,0,0,1,
 1,0,0,0,1,
 1,0,0,0,1,
 1,1,1,1,0,},

{1,1,1,1,1,
 1,0,0,0,0,
 1,1,1,1,1,
 1,0,0,0,0,
 1,1,1,1,1,},

{1,1,1,1,1,
 1,0,0,0,0,
 1,1,1,1,1,
 1,0,0,0,0,
 1,0,0,0,0,},

{1,1,1,1,1,
 1,0,0,0,0,
 1,0,1,1,1,
 1,0,0,0,1,
 1,1,1,1,1,},

{1,0,0,0,1,
 1,0,0,0,1,
 1,1,1,1,1,
 1,0,0,0,1,
 1,0,0,0,1,},

{1,1,1,1,1,
 0,0,1,0,0,
 0,0,1,0,0,
 0,0,1,0,0,
 1,1,1,1,1,},

{1,1,1,1,1,
 0,0,1,0,0,
 0,0,1,0,0,
 0,0,1,0,0,
 1,1,1,0,0,},

{1,0,0,0,1,
 1,0,0,1,0,
 1,1,1,0,0,
 1,0,0,1,0,
 1,0,0,0,1,},

{1,0,0,0,0,
 1,0,0,0,0,
 1,0,0,0,0,
 1,0,0,0,0,
 1,1,1,1,1,},

{1,0,0,0,1,
 1,1,0,1,1,
 1,0,1,0,1,
 1,0,0,0,1,
 1,0,0,0,1,},

{1,0,0,0,1,
 1,1,0,0,1,
 1,0,1,0,1,
 1,0,0,1,1,
 1,0,0,0,1,},

{1,1,1,1,1,
 1,0,0,0,1,
 1,0,0,0,1,
 1,0,0,0,1,
 1,1,1,1,1,},

{1,1,1,1,1,
 1,0,0,0,1,
 1,1,1,1,1,
 1,0,0,0,0,
 1,0,0,0,0,},

{1,1,1,1,1,
 1,0,0,0,1,
 1,0,1,0,1,
 1,0,0,1,1,
 1,1,1,1,1,},

{1,1,1,1,1,
 1,0,0,0,1,
 1,1,1,1,1,
 1,0,0,1,0,
 1,0,0,0,1,},

{1,1,1,1,1,
 1,0,0,0,0,
 1,1,1,1,1,
 0,0,0,0,1,
 1,1,1,1,1,},

{1,1,1,1,1,
 0,0,1,0,0,
 0,0,1,0,0,
 0,0,1,0,0,
 0,0,1,0,0,},
--------------------end------------------end---------------------------------
------------art2.c-----------cut here-----------------art2.c-----------------
#include <stdio.h>
#include <math.h>
#include <strings.h>			/* dans addition */

#define M 10
#define N 8
#define STREQ(s1,s2) (strcmp((s1),(s2)) == 0)
  
/*****************************************************************/
/**	art2 algorithm						**/
/** simulates the art2 algorithm for clustering vectors		**/
/**	Steve Lehar Sept 89					**/
/*****************************************************************/

/**** global variables & functions ****/
int reset;
float w[M],p[M],oldp[M],x[M],q[M],v[M],u[M],r[M],topdown[M];
float w1[M],p1[M],x1[M],q1[M],v1[M],u1[M];
float rho,normp,normu,normv,normw,normr;
float normq1,normp1,normu1,normv1,normw1,normr1;
float zup[M][N],zdn[N][M],oldzdn[M];
float I[M], F1[M], F2[N],F2val[N];
int F2winner,F2reset[N];
float pchange=1.0,zchange=1.0,diff();
float a=5.0, b=5.0, c=0.225, e=0.000001, theta = 0.3,del_t = 0.1;
char pattern[80];
FILE *fp;
int iter,learning,nreset=0;
float f(),floatval();

main(argc,argv)
int argc;
char *argv[];
{
  int i,j,k;
  char key;
  
  /**** initialize ****/
  init_weights();
  get_commandline_args(argc,argv);

  /**** learn training patterns ****/
  for(i=0;i<10;i++){
    fscanf(fp,"%s ",pattern);
    learning = 1;
/****
    if(STREQ(pattern,"E") || STREQ(pattern,"F")) learning = 0;
    else learning = 1;
****/
    read_input();
    init_reset();
    print_input();
  restart:;
    init_nodes();

    /**** equilibrate f0 ****/
    for(iter=0;iter<10;iter++){
      update_f0();
    }

    /**** equilibrate f1 ****/
    for(iter=0;iter<10;iter++){
      update_f1();
    }

    /**** equilibrate f1 & f2 ****/
    for(iter=0;iter<10;iter++){
      update_f1();
      update_f2();
      update_r();
    }

    /**** check for reset ****/
    if(reset){
      printf("RESET\n");
      reset = 0;
      nreset++;
      if(nreset==N)goto nextpattern;
      if(learning) goto restart;
    }

    if(!learning) goto skip;

    /**** learn resonant pattern ****/
    for(k=0;k<1000;k++){
      update_zup();
      update_zdn();
      update_f1();
    }


skip:;
    print_vals();
    if(is_arg("-v",argc,argv)) print_f0();
    print_f1(); 
    print_f2();
    print_z(F2winner);

}
    printf("\n");
nextpattern:;

  fclose(fp);
}


/****
    zchange = diff(zdn[F2winner],oldzdn,N);
    pchange = diff(p,oldp,M);
    print_f2(); 
    print_all();
****/

/*****************************************************************/
/**	function f()						**/
/** returns 0 if x<theta, else returns theta			**/
/*****************************************************************/
float f(x)
float x;
{
  if(x<theta) return(0.0); else return(x);
}

/*****************************************************************/
/**	diff							**/
/*****************************************************************/
float diff(v1,v2,n)
float v1[],v2[];
int n;
{
  int i;
  float d;

  d = 0.0;
  for(i=0;i<n;i++){
    d += (v1[i] - v2[i])*(v1[i] - v2[i]);
    v2[i] = v1[i];
  }
  return(d);
}

/*****************************************************************/
/**	initialize weights					**/
/*****************************************************************/
init_weights()
{
  int i,j;

  /**** initialize weights ****/
  for(j=0;j<N;j++){
    for(i=0;i<M;i++){
      zup[i][j] = 1.0/(0.2*sqrt((float)M));
      zdn[j][i] = 0.0;
    }
  }
  
}

/*****************************************************************/
/**	init_reset						**/
/*****************************************************************/
init_reset()
{
  int j;

  for(j=0;j<N;j++){
    F2reset[j] = 0;
  }
  nreset = 0;
}

/*****************************************************************/
/**	initialize nodes					**/
/*****************************************************************/
init_nodes()
{
  int i,j;

  /**** initialize nodes ****/
  for(i=0;i<M;i++){
    F1[i] = 0.0;
    w1[i] = 0.0;
    x1[i] = 0.0;
    v1[i] = 0.0;
    u1[i] = 0.0;
    p1[i] = 0.0;
    q1[i] = 0.0;
    w[i] = 0.0;
    x[i] = 0.0;
    v[i] = 0.0;
    u[i] = 0.0;
    p[i] = 0.0;
    q[i] = 0.0;
  }
  for(j=0;j<M;j++){
    F2[j] = 0.0;
  }
  
}

/*****************************************************************/
/**	get_commandline_args					**/
/*****************************************************************/
get_commandline_args(argc,argv)
int argc;
char *argv[];
{

  /**** open data file & get rho ****/
  /* if (!is_arg("-i",argc,argv)){
    printf("ERROR no input file given [-i file]\n");
    exit(1); 
  }
  if (!is_arg("-r",argc,argv)){
    printf("ERROR no rho value given [-r rho]\n");
    exit(1);
  }
   fp = fopen(strval("-i",argc,argv),"r"); */
  fp = fopen("data","r");
  /* rho = floatval("-r",argc,argv); */
  rho = 0.7; 
}  

/*****************************************************************/
/**	read_input						**/
/*****************************************************************/
read_input()
{
  int i,j;

  /**** read input vector ****/
  for(i=0;i<M;i++){
    fscanf(fp,"%f",&I[i]);
  }
}
    
/*****************************************************************/
/**	update_f0						**/
/*****************************************************************/
update_f0()
{
    int i,j;

    /**** update w nodes ****/
    normw1 = 0.0;
    for(i=0;i<M;i++){
      w1[i] = I[i]  + a*u1[i];
      normw1 += w1[i] * w1[i];
    }
    normw1 = sqrt(normw1);

    /**** update p nodes ****/
    normp1 = 0.0;
    for(i=0;i<M;i++){
      p1[i] = u1[i];
      normp1 += p1[i] * p1[i];
    }
    normp1 = sqrt(normp1);

    /**** update x nodes ****/
    for(i=0;i<M;i++){
      x1[i] = w1[i] / (e+normw1);
    }

    /**** update q nodes ****/
    normq1 = 0.0;
    for(i=0;i<M;i++){
      q1[i] = p1[i] / (e+normp1);
      normq1 += q1[i] * q1[i];
    }
    normq1 = sqrt(normq1);

    /**** update v nodes ****/
    normv1 = 0.0;
    for(i=0;i<M;i++){
      v1[i] = f(x1[i]) + b*f(q1[i]);
      normv1 += v1[i] * v1[i];
    }
    normv1 = sqrt(normv1);

    /**** update u nodes ****/
    for(i=0;i<M;i++){
      u1[i] = v1[i] / (e+normv1);
    }

}
    
/*****************************************************************/
/**	update_f1						**/
/*****************************************************************/
update_f1()
{
    int i,j;

    /**** update w nodes ****/
    normw = 0.0;
    for(i=0;i<M;i++){
      w[i] = q1[i]  + a*u[i];
      normw += w[i] * w[i];
    }
    normw = sqrt(normw);

    /**** update top-down value ****/
    for(i=0;i<M;i++){
      topdown[i] = 0.0;
      for(j=0;j<N;j++){
        topdown[i] += F2[j] * zdn[j][i];
      }
    }

    /**** update p nodes ****/
    normp = 0.0;
    for(i=0;i<M;i++){
      p[i] = u[i] + topdown[i];
      normp += p[i] * p[i];
    }
    normp = sqrt(normp);

    /**** update x nodes ****/
    for(i=0;i<M;i++){
      x[i] = w[i] / (e+normw);
    }

    /**** update q nodes ****/
    for(i=0;i<M;i++){
      q[i] = p[i] / (e+normp);
    }

    /**** update v nodes ****/
    normv = 0.0;
    for(i=0;i<M;i++){
      v[i] = f(x[i]) + b*f(q[i]);
      normv += v[i] * v[i];
    }
    normv = sqrt(normv);

    /**** update u nodes ****/
    normu = 0.0;
    for(i=0;i<M;i++){
      u[i] = v[i] / (e+normv);
      normu += u[i] * u[i];
    }
    normu = sqrt(normu);

}  

/*****************************************************************/
/**	update_f2						**/
/*****************************************************************/
update_f2()
{
   int i,j;
   float max;

  /**** compute F2 ****/
  for(j=0;j<N;j++){
    F2[j] = 0.0;
    for(i=0;i<M;i++){
      F2[j]+= p[i] * zup[i][j];
    }
    F2val[j] = F2[j];
  }

  /**** compute F2winner ****/
  max = 0.0;
  for(j=0;j<N;j++){
    if(F2[j] > max && !F2reset[j]){
      max = F2[j];
      F2winner = j;
    }
    F2[j] = 0.0;
  }
  F2[F2winner] = 0.8;
}

/*****************************************************************/
/**	update_r						**/
/*****************************************************************/
update_r()
{
  int i,j;
  float cnormp;

  /**** compute cnormp ****/
  cnormp = 0.0;
  for(i=0;i<M;i++){
    cnormp += c * p[i] * p[i];
  }
  cnormp = sqrt(cnormp);


  /**** compute reset ****/
  normr = 0.0;
  for(i=0;i<M;i++){
    r[i] = (q1[i] + c*p[i]) / (e + normq1 + c*normp);
    normr += r[i] * r[i];
  }
  normr = sqrt(normr);

  if(rho/(e+normr) > 1.0) reset = 1; else reset = 0;
  if(reset){
    F2reset[F2winner] = 1.0;
  }
    
}
/*****************************************************************/
/**	update_zup						**/
/*****************************************************************/
update_zup()
{
  int i,j;

  for(i=0;i<M;i++){
    zup[i][F2winner] += del_t * (p[i] - zup[i][F2winner]);
  }

}

/*****************************************************************/
/**	update_zdn						**/
/*****************************************************************/
update_zdn()
{
   int i,j;

  /**** set top-down weights ****/
  for(i=0;i<M;i++){
    zdn[F2winner][i] += del_t * (p[i] - zdn[F2winner][i]);
  }

}

/*****************************************************************/
/** 	print_f0						**/
/*****************************************************************/
print_f0()
{
  int i,j;

  /**** print values ****/
  printf("**** F 0  iteration %2d ****\n",iter);
  for(i=0;i<M;i++){
    printf("%s%2d: I=%5.2f  w1=%5.2f  x1=%5.2f  v1=%5.2f  ",
                   pattern, i, I[i], w1[i], x1[i], v1[i]);
    printf("u1=%5.2f  p1=%5.2f  q1=%5.2f\n",
            u1[i], p1[i], q1[i]);
  }
}

/*****************************************************************/
/**	print_f1						**/
/*****************************************************************/
print_f1()
{
  int i,j;

  /**** print F1 values ****/
  printf("w:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",w[i]);
  }
  printf("\n");
  printf("x:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",x[i]);
  }
  printf("\n");
  printf("v:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",v[i]);
  }
  printf("\n");
  printf("u:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",u[i]);
  }
  printf("\n");
  printf("p:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",p[i]);
  }
  printf("\n");
  printf("q:   ");
  for(i=0;i<M;i++){
    printf("%5.2f ",q[i]);
  }
  printf("\n");
}

/*****************************************************************/
/**	print_f2						**/
/*****************************************************************/
print_f2()
{
  int i,j;

  /**** print f2 values ****/
  printf("F2:  ");
  for(j=0;j<N;j++){
    printf("%5.2f ",F2val[j]);
  }
  printf("\n");

  /**** underline winner ****/
  printf("     ");
  for(j=0;j<N;j++){
    if(j==F2winner) printf("^WINS^");
    else printf("      ");
  }
  printf("\n");
}

/*****************************************************************/
/**	print_z							**/
/*****************************************************************/
print_z(node)
int node;
{
  int i,j;

  printf("zup%1d:",node);
  for(i=0;i<M;i++){
    printf("%5.2f ",zup[i][node]);
  }
  printf("\n");
  printf("zdn%1d:",node);
  for(i=0;i<M;i++){
    printf("%5.2f ",zdn[node][i]);
  }
  printf("\n");
}

/*****************************************************************/
/**	print_input						**/
/*****************************************************************/
print_input()
{
  int i,j;

  /**** print input ****/
  printf("\n %s:  ",pattern);
  for(i=0;i<M;i++){
    printf("%5.2f ",I[i]);
  }
  printf("\n");
}

/*****************************************************************/
/** 	print_vals						**/
/*****************************************************************/
print_vals()
{
  printf("rho=%5.2f normr=%5.2f\twinner=%1d",rho,normr,F2winner);
  if(reset)printf(" RESET");
  if(learning) printf(" learning");
  printf("\n");
}

/*****************************************************************/
/** 	print_all						**/
/*****************************************************************/
print_all(iter)
int iter;
{
  int i,j;

  /**** print values ****/
  printf("**** iteration %2d ****\n",iter);
  for(i=0;i<3;i++){
    printf("%s I=%5.2f  v1=%5.2f  q1=%5.2f  v=%5.2f  ",
                   pattern,I[i], v1[i], q1[i], v[i]);
    printf("p=%5.2f  r=%5.2f \n",
            p[i], r[i]);
  }
  printf("normr=%5.2f ",normr);
  if(reset) printf("RESET\n");
  printf("\n");

}


/*----------------------------------------dans additions --------*/

int is_arg(str,arg_c,arg_v)
char str[],*arg_v[];
int arg_c;
{
   int result;
   result = 0;
   while (arg_c > 0) 
      {
      result = (strcmp(str,arg_v[--arg_c]) == 0) ? 1 : 0;
      }
   return(result);
}
      
float floatval(str,arg_c,arg_v)
char str[],*arg_v[];
int arg_c;
{
   float result;
   result = 0.0;
   while (arg_c > 0)
      {  
      if (strcmp(str,arg_v[--arg_c]) == 0) break;
      }  
   if (arg_c > 0)
         {
         if (sscanf(arg_v[arg_c+1],"%lf",result) < 1)
         printf("\n\7ERR: syntax error - expecting a float, got %s\n",
                 arg_v[arg_c+1]);
         }
   return(result);
}

char strval(str,arg_c,arg_v)
char str[],*arg_v[];
int arg_c;
{
   char  result[128];
   int ll;
   for (ll = 0; ll < 128; result[ll++] = 'A'); /* clear out result str */
   result[127] = '\n';				/* null terminate  */

   while (arg_c > 0)
      {  
      if (strcmp(str,arg_v[--arg_c]) == 0) break;
      }  
   if (arg_c > 0)
         {
         if (sscanf(arg_v[arg_c+1],"%s",result) < 1)
           printf("\n\7ERR: syntax error - expecting a string, got nothing \n");
         }
   return(result);
}