/***********************
 * Simulation of TODAM *
 ***********************
 *
 * October 1992
 *
 * Peter Nobel & David Huber
 * Indiana University
 * Bloomington, IN  47405
 *
 * pnobel@ucs.indiana.edu & dhuber@ucs.indiana.edu
 *
 * This program simulates the following test-conditions:
 *      single item recognition
 *      associative recognition: intact - rearranged
 *      associative recognition: intact - mixed
 *      associative recognition: intact - new
 *      cued recall
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <sys/types.h>
#include <sys/timeb.h>
#include <math.h>
#include <malloc.h>
#include <string.h>

#define MAX_ENTRIES 10000

typedef struct {
    char *name;
    char *val;
} entry;

char *makeword(char *line, char stop);
char *fmakeword(FILE *f, char stop, int *len);
char x2c(char *what);
void unescape_url(char *url);
void plustospace(char *str);


#define RAND_MAX 2147483648.0         /* the maximum number returned by rand */
#define NSIM 500                           /* maximum number of simulations */
#define NMAX 32                        /* maximum number of pairs on a list */
#define FMAX 63                               /* maximum number of features */
#define SMAXW 3                              /* maximum strength weak items */
#define SMAXS 5                            /* maximum strength strong items */
#define TRUE 1
#define SIZE1 7                                   /* number of items menu 1 */
#define SIZE2 6                                   /* number of items menu 2 */
#define SIZE3 2                                   /* number of items menu 3 */
#define SIZE4 4                                   /* number of items menu 4 */
#define SIZE5 8                                   /* number of items menu 5 */
#define SIZE6 3                                   /* number of items menu 6 */

#define integer int                             /* used in d-prime function */
#define real float
static integer c__9 = 9;
static integer c__1 = 1;
static integer c__4 = 4;
float hrate;
float farate;
float dp;

FILE *fp;                                    /* file pointer to output file */
char fout[100];                                         /* output file name */
FILE *fptemp;

/***********************
 * function prototypes *
 ***********************/

void menu(void);
void making_matrix(void);
void study(void);
void recognition_single(void);
void recognition_ass_rearr(void);
void recognition_ass_mixed(void);
void recognition_ass_new(void);
void cued_recall(void);
void ranset(void);
void perm(int *, int);
void simulate(void);
void results(void);
void settings(void);
void screen(void);
void action1(int);
void action2(int);
void action3(int);
void action4(int);
void action5(int);
void action6(int);
float rand1(void);
float gau(double, double);
int randw(int);
int cor_win(int);
float single_fam(int);
float double_fam(int, int);
int decision(float);
void convolve(int, int);
void correlate(int);


/********************
 * menu definitions *
 ********************/
char *menu1[] =
			{ "Select type of test",
				"Set length of list",
				"Set strength of list",
				"Set parameters",
				"Simulation",
				"Write results",
				"Exit",                   };

char *menu2[] =
			{ "Single item recognition",
				"Assoc. recog. - rearranged",
				"Assoc. recog. - mixed",
				"Assoc. recog. - new",
				"Cued recall",
				"Return to Main Menu",    };

char *menu3[] =
			{ "Set length of list (number of pairs)",
				"Return to Main Menu",    };

char *menu4[] =
			{ "Pure weak",
				"Pure strong",
				"Mixed",
				"Return to Main Menu",    };

char *menu5[] =
			{ "Number of Features",
				"Power",
				"Forgetting parameter",
				"Weights",
				"Criterion (recog.)",
				"Decision system",
				"Criterion (recall)",
				"Return to Main Menu",    };

char *menu6[] =
			{ "Number of simulations",
				"Simulate",
				"Return to Main Menu",    };

/******************************
 * variables for menu control *
 ******************************/
int curpos;                 /* keeps track of the cursor in the menu (rows) */
int code;                   /* reads in the (extended) code of key pressed  */
int back;                   /* switch between menus: 0 = stay, 1 = return   */

/***************************************
 * global variables for the simulation *
 ***************************************/
int nsim=50;                                       /* number of simulations */
int cond=1;                                                 /* type of test */
int n=8;             /* the length of the list (# of pairs): multiples of 8 */
int length;                               /* actual number of presentations */
int strength=1; /*strength of the list: 1-pure weak, 2-pure strong, 3-mixed */
int s1=1;                                         /* strength of weak items */
int s2=2;                                       /* strength of strong items */

/********************************** 
 * global variables for the model *
 **********************************/


float	mu; /* mean of noise */
float	sigma; /* std dev of noise */
float	new_memory[FMAX][FMAX];
float	matrix1[FMAX][FMAX];
float	crit;
float	recall_crit;
int	features;
float	learn;

float items[NMAX*4][FMAX];   /* max no. of different items of FMAX features */

float memory[FMAX];                                        /* memory vector */

float convolved[FMAX];                               /* two convolved items */
float correlated[FMAX];               /* item correlated with memory vector */

int order[NMAX*SMAXS];                    /* order of presentation of pairs */

int count,                         /* counter for the number of simulations */
		count_cw,                   /* counter for weak corrects in cued recall */
		count_cs,                 /* counter for strong corrects in cued recall */
		count_hw,                       /* counter for weak hits in recognition */
		count_hs,                     /* counter for strong hits in recognition */
		count_f,                     /* counter for false alarms in recognition */
		count_fs;                      /* counter for strong false alarms (ass) */

float famtw,                         /* total sum familiarity: weak targets */
			famtw2,                /* total sum familiarity squared: weak targets */
			famts,                       /* total sum familiarity: strong targets */
			famts2,              /* total sum familiarity squared: strong targets */
			famd,                           /* total sum familiarity: distractors */
			famd2,                  /* total sum familiarity squared: distractors */
			famds,                   /* total sum familiarity: strong distractors */
			famds2;          /* total sum familiarity squared: strong distractors */

float VARTW[NSIM],          /* array for variance familiarity: weak targets */
			VARTS[NSIM],        /* array for variance familiarity: strong targets */
			VARD[NSIM],            /* array for variance familiarity: distractors */
			VARDS[NSIM],           /* array for variance familiarity: distractors */
			DPRPW[NSIM],         /* array for d-prime (weak) based on proportions */
			DPRPS[NSIM],       /* array for d-prime (strong) based on proportions */
			DPRDW[NSIM],       /* array for d-prime (weak) based on distributions */
			DPRDS[NSIM];     /* array for d-prime (strong) based on distributions */


/****************************
 * parameters for the model *
 ****************************/
int features=15;                               /* number of features chosen */
float power=1.0;                                         /* power of vector */
float alpha=0.9;                                    /* forgetting parameter */
float gamma1=1.0;                                          /* weight item 1 */
float gamma2=1.0;                                          /* weight item 2 */
float gamma3=1.0;                         /* weight associative information */
float upper=0.1;                             /* upper criterion recognition */
float lower=0.1;                             /* lower criterion recognition */
float noise=0.0;                            /* variance for decision system */
float rate=0.9;                                         /* convergence rate */
float giveup=0.3;                      /* cut-off for correlation in recall */

/************************
 * variables for output *
 ************************/

/* cued recall */
float pcw;                               /* final result recall: weak items */
float pcs;                             /* final result recall: strong items */
float corw;      /* output variable of average correlation for weak targets */
float cors;    /* output variable of average correlation for strong targets */
float mycorw;      /* output variable of average correlation for weak targets */
float mycors;    /* output variable of average correlation for strong targets */

/* recognition: variance and d' for both averaged and total */
float hw;                            /* final result recognition: weak hits */
float mtw;                                 /* mean familiarity weak targets */
float vartw;                    /* variance familiarity weak targets (avg.) */
float tvartw;                  /* variance familiarity weak targets (total) */
float hs;                          /* final result recognition: strong hits */
float mts;                               /* mean familiarity strong targets */
float varts;                  /* variance familiarity strong targets (avg.) */
float tvarts;                /* variance familiarity strong targets (total) */
float f;                          /* final result recognition: false alarms */
float fs;                  /* final result recognition: strong false alarms */
float md;                                   /* mean familiarity distractors */
float mds;                           /* mean familiarity strong distractors */
float vard;                      /* variance familiarity distractors (avg.) */
float vards;              /* variance familiarity strong distractors (avg.) */
float tvard;                    /* variance familiarity distractors (total) */
float tvards;            /* variance familiarity strong distractors (total) */
float dprpw;              /* dprime based on proportions: weak items (avg.) */
float tdprpw;            /* dprime based on proportions: weak items (total) */
float dprps;            /* dprime based on proportions: strong items (avg.) */
float tdprps;          /* dprime based on proportions: strong items (total) */
float dprdw;            /* dprime based on distributions: weak items (avg.) */
float tdprdw;          /* dprime based on distributions: weak items (total) */
float dprds;          /* dprime based on distributions: strong items (avg.) */
float tdprds;        /* dprime based on distributions: strong items (total) */

/********
 * MAIN *
 ********/
main()
{
    entry entries[MAX_ENTRIES];
    register int x,m=0;
    int cl;

    printf("Content-type: text/html%c%c",10,10);

    if(strcmp(getenv("REQUEST_METHOD"),"POST")) {
        printf("This script should be referenced with a METHOD of POST.\n");
        printf("If you don't understand this, see this ");
        printf("<A HREF=\"http://www.ncsa.uiuc.edu/SDG/Software/Mosaic/Docs/fill-out-forms/overview.htm
l\">forms overview</A>.%c",10);
        exit(1);
    }
    if(strcmp(getenv("CONTENT_TYPE"),"application/x-www-form-urlencoded")) {
        printf("This script can only be used to decode form results. \n");
        exit(1);
    }
    cl = atoi(getenv("CONTENT_LENGTH"));

    for(x=0;cl && (!feof(stdin));x++) {
        m=x;
        entries[x].val = fmakeword(stdin,'&',&cl);
        plustospace(entries[x].val);
        unescape_url(entries[x].val);
        entries[x].name = makeword(entries[x].val,'=');

    if (strcmp(entries[x].name, "TypeOfTest") == 0){
        if (strcmp(entries[x].val, "Single item recognition") == 0)
            cond = 1;

        else if (strcmp(entries[x].val, "Associative recognition - rearranged") == 0)
            cond = 2;
        else if (strcmp(entries[x].val, "Associative recognition - mixed") == 0)
            cond = 3;
        else if (strcmp(entries[x].val, "Associative recognition - new") == 0)
            cond = 4;
        else if (strcmp(entries[x].val, "Cued recall") == 0)
            cond = 5;
        }
    else if (strcmp(entries[x].name, "ListLength") == 0){
        sscanf(entries[x].val, "%d", &n);
        if (n > NMAX) {
            printf("Maximum number of pairs is %d!", NMAX);
            n = NMAX;
            }
        if (n % 8 != 0) {                           /* if not a multiple of 8 */
            printf("Number of pairs should be a multiple of 8! Setting list length to 8.");
            n=8;
            }
        }
    else if (strcmp(entries[x].name, "ListStrength") == 0){
        if (strcmp(entries[x].val, "Pure Weak") == 0)
            strength = 1;
        else if (strcmp(entries[x].val, "Pure Strong") == 0)
            strength = 2;
        else if (strcmp(entries[x].val, "Mixed") == 0)
            strength = 3;
        }
    else if (strcmp(entries[x].name, "WeakStrength") == 0){
        sscanf(entries[x].val, "%d", &s1);
        if (s1 > SMAXW) {
            printf("Maximum strength for weak items is %d!", SMAXW);
            s1 = SMAXW;
            }
           else if (s1 < 1) {
            printf("Minimum strength for weak items is 1!");
            s1 = 1;
            }
        }
    else if (strcmp(entries[x].name, "StrongStrength") == 0) {
        sscanf(entries[x].val, "%d", &s2);
        if (s2 > SMAXS) {
            printf("Maximum strength for strong items is %d!", SMAXS);
            s2 = SMAXS;
            }
        else if (s2 < 2) {
            printf("Minimum strength for strong items is 2!");
            s2 = 2;
            }
        }
    else if (strcmp(entries[x].name, "NumberOfFeatures") == 0) {
        sscanf(entries[x].val, "%d", &features);
        if (features > 32) {
            printf("Maximum number of features is 32!");
            features = 32;
            }
        else if (features < 4) {
            printf("Minimum number of features is 4!");
            features = 4;
            }
        }
    else if (strcmp(entries[x].name, "power") == 0)
        sscanf(entries[x].val, "%f", &power);
    else if (strcmp(entries[x].name, "alpha") == 0)
        sscanf(entries[x].val, "%f", &alpha);
    else if (strcmp(entries[x].name, "gamma1") == 0)
        sscanf(entries[x].val, "%f", &gamma1);
    else if (strcmp(entries[x].name, "gamma2") == 0)
        sscanf(entries[x].val, "%f", &gamma2);
    else if (strcmp(entries[x].name, "gamma3") == 0)
        sscanf(entries[x].val, "%f", &gamma3);
    else if (strcmp(entries[x].name, "upper") == 0)
        sscanf(entries[x].val, "%f", &upper);
    else if (strcmp(entries[x].name, "lower") == 0)
        sscanf(entries[x].val, "%f", &lower);
    else if (strcmp(entries[x].name, "noise") == 0)
        sscanf(entries[x].val, "%f", &noise);
    else if (strcmp(entries[x].name, "rate") == 0)
        sscanf(entries[x].val, "%f", &rate);
    else if (strcmp(entries[x].name, "giveup") == 0)
        sscanf(entries[x].val, "%f", &giveup);
    else if (strcmp(entries[x].name, "NumberOfSimulations") == 0){
        sscanf(entries[x].val, "%d", &nsim);
        if (nsim > NSIM) {
            printf("Maximum number simulatoins is %d!", NSIM);
            nsim = NSIM;
            }
        }
    }

    printf("<H1>Query Results</H1>");
    if (s2 <= s1)
        printf("Strong items should be stronger than weak ones!<BR>\n");

    printf("<PRE>\n");
    ranset();                                 /* init the random no generator */
    simulate();
    screen();
    printf("</PRE>\n");
}


/*********************
 * FUNCTION SIMULATE *
 *********************
 *
 * main simulation routine
 *
 */
void simulate(void)
{
	int i;

	count_cw  = 0;                                 /* initialize the counters */
	count_cs  = 0;
	count_hw  = 0;
	count_hs  = 0;
	count_f   = 0;
	count_fs  = 0;
	count     = 0;

	famtw     = 0.0;
	famtw2    = 0.0;
	famts     = 0.0;
	famts2    = 0.0;
	famd      = 0.0;
	famd2     = 0.0;
	famds     = 0.0;
	famds2    = 0.0;

	corw      = 0.0;
	cors      = 0.0;
																					 /* initialize the storage arrays */
	for (i=0; i<NSIM; i++) {
		VARTW[i]  = 0.0;
		VARTS[i]  = 0.0;
		VARD[i]   = 0.0;
		VARDS[i]  = 0.0;
		DPRPW[i]  = 0.0;
		DPRPS[i]  = 0.0;
		DPRDW[i]  = 0.0;
		DPRDS[i]  = 0.0;
	}

	if (strength==1)
		length=s1*n;
	else if (strength==2)
		length=s2*n;
	else
		length=(s1*n/2)+(s2*n/2);
	for (i=1; i<=nsim; i++) {
		making_matrix();                               /* initialize the matrix */
		study();
		switch(cond) {
			case 1:
				recognition_single();
				break;
			case 2:
				recognition_ass_rearr();
				break;
			case 3:
				recognition_ass_mixed();
				break;
			case 4:
				recognition_ass_new();
				break;
			case 5:
				cued_recall();
				break;
		}
	}
}

/**************************
 * FUNCTION MAKING_MATRIX *
 **************************
 *
 * This function sets up a new lexicon
 *
 */
void making_matrix(void)
{
	int i, j, k;
	int check;
	float temp, choose;
	float pwer, totpower;

												/* creates the features of each item in the lexicon */
	for (i=0; i<(4*n); i++) {

		L1:

		for (k=0; k<features; k++) {
			items[i][k] = gau((double)0,(double)(power/features));
		}
											 /* makes sure that each pair of items item is unique */
		for (j=0; j<i; j++) {
			check=0;
			for (k=0; k<features; k++) {
				if (items[i][k] != items[j][k])
					check=1;
			}
			if (check == 0)
				goto L1;
		}
	}
}

/******************
 * FUNCTION STUDY *
 ******************
 *
 * creation of memory due to list presentation
 *
 */
void study(void)
{
	int i, j, k;
	int str;

	for (i=0; i<features; i++)
		memory[i] = 0.0;

	if (strength == 1)
		str=s1;
	if (strength == 2)
		str=s2;

	if (strength == 1 || strength == 2) {                        /* pure list */
		i=0;
		for (j=0; j<n; j++) {                                    /* pair number */
			for (k=0; k<str; k++) {                      /* number of repetitions */
				order[i] = j;
				i++;
			}
		}
	}
	else {                                                      /* mixed list */
		i=0;
		for (j=0; j<n/2; j++) {
			for (k=0; k<s1; k++) {
				order[i] = j;
				i++;
			}
		}
		for (j=n/2; j<n; j++) {
			for (k=0; k<s2; k++) {
				order[i] = j;
				i++;
			}
		}
	}
	perm(order, length);

	for (i=0; i<length; i++) {
		for (j=0; j<features; j++) {            /* apply alpha to memory vector */
			memory[j] *= alpha;
		}
		for (j=0; j<features; j++) {                /* add first member of pair */
			memory[j] += (gamma1*items[2*order[i]][j]);
		}
		for (j=0; j<features; j++) {               /* add second member of pair */
			memory[j] += (gamma2*items[(2*order[i])+1][j]);
		}
		convolve(2*order[i], (2*order[i])+1);
		for (j=0; j<features; j++) {                         /* add convolution */
			memory[j] += (gamma3*convolved[j]);
		}
	}
}                                                     /* end function STUDY */

/*******************************
 * FUNCTION RECOGNITION_SINGLE *
 *******************************
 *
 * single item recognition simulation routine
 *
 */
void recognition_single(void)
{
	static real hrate, farate, dp;
	extern int dprime_();
	int i;
	int win;
	float t;
	float fam;                                    /* familiarity of test-item */
	int chw  = 0,
			chs  = 0,
			cf   = 0;
	float ftw   = 0.0,          /* variables to compute means and variances */
				ftw2  = 0.0,          /* for the target and distractor distibutions */
				fts   = 0.0,
				fts2  = 0.0,
				fd    = 0.0,
				fd2   = 0.0,
				vtw   = 0.0,
				vts   = 0.0,
				vdw   = 0.0;

	for (i=0; i<n; i++) {                                          /* targets */
		fam=single_fam(2*i);
		if ((strength == 1) || (strength == 3 && i < n/2)) {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}
		if ((strength == 2) || (strength == 3 && i >= n/2)) {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}
		win = decision(fam);
		if (win == 1) {
			if ((strength == 1) || (strength == 3 && i < n/2)) {
				count_hw++;
				chw++;
			}
			if ((strength == 2) || (strength == 3 && i >= n/2)) {
				count_hs++;
				chs++;
			}
		}
	}
	for (i=2*n; i<(2*n)+n; i++) {                              /* distractors */
		fam = single_fam(i);

		famd += fam;
		famd2 += fam * fam;
		fd += fam;
		fd2 += fam * fam;

		win = decision(fam);
		if (win == 1) {
			count_f++;
			cf++;
		}
	}
																		 /* compute results and store in arrays */
	if (strength == 3)
		t = (float)n/2;
	else
		t = (float)n;

	vtw = ((t*ftw2)-(ftw*ftw))/(t*(t-1));
	if (vtw < 0.0001)
		vtw = 0.0001;
	VARTW[count]  = vtw;

	vts = ((t*fts2)-(fts*fts))/(t*(t-1));
	if (vts < 0.0001)
		vts = 0.0001;
	VARTS[count]  = vts;

	vdw = (((float)n * fd2)-(fd*fd))/((float)n*((float)n-1));
	if (vdw < 0.0001)
		vdw = 0.0001;
	VARD[count]   = vdw;

	hrate = (float)chw/t;
	farate = (float)cf/(float)n;
	dprime_(&hrate, &farate, &dp);
	DPRPW[count]  = dp;

	hrate = (float)chs/t;
	farate = (float)cf/(float)n;
	dprime_(&hrate, &farate, &dp);
	DPRPS[count]  = dp;

	DPRDW[count]  = ((ftw/t)-(fd/(float)n))/(sqrt(VARD[count]));
	DPRDS[count]  = ((fts/t)-(fd/(float)n))/(sqrt(VARD[count]));

	count++;
}

/**********************************
 * FUNCTION RECOGNITION_ASS_REARR *
 **********************************
 *
 * associative recognition simulation routine
 * targets    : A-B (intact)
 * distractors: A-D (rearranged)
 *
 */
void recognition_ass_rearr(void)
{
	static real hrate, farate, dp;
	extern int dprime_();
	int i, j;
	int win;
	float t;
	float fam;                                    /* familiarity of test-item */
	int chw   = 0,
			chs   = 0,
			cf    = 0,
			cfs   = 0;
	float ftw   = 0.0,          /* variables to compute means and variances */
				ftw2  = 0.0,          /* for the target and distractor distibutions */
				fts   = 0.0,
				fts2  = 0.0,
				fd    = 0.0,
				fd2   = 0.0,
				fds   = 0.0,
				fds2  = 0.0,
				vtw   = 0.0,
				vts   = 0.0,
				vdw   = 0.0,
				vds   = 0.0;

	for (i=0; i<n/4; i++) {                                      /* targets 1 */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 2) {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}
		else {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 2) {
				count_hs++;
				chs++;
			}
			else {
				count_hw++;
				chw++;
			}
		}
	}

	for (i=n/4; i<n/2; i++) {                                /* distractors 1 */
		j=(2*i)+3;
		if (i==((n/2)-1))
			j=(2*n/4)+1;
		fam = double_fam((2*i),j);
		fam += single_fam(2*i);
		fam += single_fam(j);

		famd += fam;
		famd2 += fam * fam;
		fd += fam;
		fd2 += fam * fam;

		win = decision(fam);
		if (win == 1) {
			count_f++;
			cf++;
		}
	}

	for (i=n/2; i<3*n/4; i++) {                                  /* targets 2 */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 1) {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}
		else {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 1) {
				count_hw++;
				chw++;
			}
			else {
				count_hs++;
				chs++;
			}
		}
	}

	for (i=3*n/4; i<n; i++) {                                /* distractors 2 */
		j=(2*i)+3;
		if (i==(n-1))
			j=(2*3*n/4)+1;
		fam = double_fam((2*i),j);
		fam += single_fam(2*i);
		fam += single_fam(j);

		if (strength == 3) {
			famds += fam;
			famds2 += fam * fam;
			fds += fam;
			fds2 += fam * fam;
		}
		else {
			famd += fam;
			famd2 += fam * fam;
			fd += fam;
			fd2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 3) {
				count_fs++;
				cfs++;
			}
			else {
				count_f++;
				cf++;
			}
		}
	}
																		 /* compute results and store in arrays */
	if (strength == 3)
		t = (float)n/4;
	else
		t = (float)n/2;

	vtw = ((t*ftw2)-(ftw*ftw))/(t*(t-1));
	if (vtw < 0.0001)
		vtw = 0.0001;
	VARTW[count]  = vtw;

	vts = ((t*fts2)-(fts*fts))/(t*(t-1));
	if (vts < 0.0001)
		vts = 0.0001;
	VARTS[count]  = vts;

	vdw = ((t*fd2)-(fd*fd))/(t*(t-1));
	if (vdw < 0.0001)
		vdw = 0.0001;
	VARD[count]   = vdw;

	if (strength == 3) {
		vds = ((t*fds2)-(fds*fds))/(t*(t-1));
		if (vds < 0.0001)
			vds = 0.0001;
		VARDS[count]  = vds;
	}

	hrate = (float)chw/t;
	farate = (float)cf/t;
	dprime_(&hrate, &farate, &dp);
	DPRPW[count]  = dp;

	hrate = (float)chs/t;
	if (strength == 3 )
		farate = (float)cfs/t;
	else
		farate = (float)cf/t;
	dprime_(&hrate, &farate, &dp);
	DPRPS[count]  = dp;

	DPRDW[count]  = ((ftw/t)-(fd/t))/(sqrt(VARD[count]));
	if (strength == 3 )
		DPRDS[count]  = ((fts/t)-(fds/t))/(sqrt(VARDS[count]));
	else
		DPRDS[count]  = ((fts/t)-(fd/t))/(sqrt(VARD[count]));

	count++;
}

/**********************************
 * FUNCTION RECOGNITION_ASS_MIXED *
 **********************************
 *
 * associative recognition simulation routine
 * targets    : A-B (intact)
 * distractors: A-X (mixed)
 *
 */
void recognition_ass_mixed(void)
{
	static real hrate, farate, dp;
	extern int dprime_();
	int i, j;
	int win;
	float t;
	float fam;                                    /* familiarity of test-item */
	int chw   = 0,
			chs   = 0,
			cf    = 0,
			cfs   = 0;
	float ftw   = 0.0,          /* variables to compute means and variances */
				ftw2  = 0.0,          /* for the target and distractor distibutions */
				fts   = 0.0,
				fts2  = 0.0,
				fd    = 0.0,
				fd2   = 0.0,
				fds   = 0.0,
				fds2  = 0.0,
				vtw   = 0.0,
				vts   = 0.0,
				vdw   = 0.0,
				vds   = 0.0;

	for (i=0; i<n/4; i++) {                                      /* targets 1 */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 2) {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}
		else {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 2) {
				count_hs++;
				chs++;
			}
			else {
				count_hw++;
				chw++;
			}
		}
	}

	for (i=n/4; i<n/2; i++) {                                /* distractors 1 */
		j=(2*n)+i;
		fam = double_fam(2*i,j);
		fam += single_fam(2*i);
		fam += single_fam(j);

		famd += fam;
		famd2 += fam * fam;
		fd += fam;
		fd2 += fam * fam;

		win = decision(fam);
		if (win == 1) {
			count_f++;
			cf++;
		}
	}

	for (i=n/2; i<3*n/4; i++) {                                  /* targets 2 */
		fam = double_fam(2*i,(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 1) {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}
		else {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 1) {
				count_hw++;
				chw++;
			}
			else {
				count_hs++;
				chs++;
			}
		}
	}

	for (i=3*n/4; i<n; i++) {                                /* distractors 2 */
		j=(2*n)+i;
		fam = double_fam(2*i,j);
		fam += single_fam(2*i);
		fam += single_fam(j);

		if (strength == 3) {
			famds += fam;
			famds2 += fam * fam;
			fds += fam;
			fds2 += fam * fam;
		}
		else {
			famd += fam;
			famd2 += fam * fam;
			fd += fam;
			fd2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 3) {
				count_fs++;
				cfs++;
			}
			else {
				count_f++;
				cf++;
			}
		}
	}
																		 /* compute results and store in arrays */
	if (strength == 3)
		t = (float)n/4;
	else
		t = (float)n/2;

	vtw = ((t*ftw2)-(ftw*ftw))/(t*(t-1));
	if (vtw < 0.0001)
		vtw = 0.0001;
	VARTW[count]  = vtw;

	vts = ((t*fts2)-(fts*fts))/(t*(t-1));
	if (vts < 0.0001)
		vts = 0.0001;
	VARTS[count]  = vts;

	vdw = ((t*fd2)-(fd*fd))/(t*(t-1));
	if (vdw < 0.0001)
		vdw = 0.0001;
	VARD[count]   = vdw;

	if (strength == 3) {
		vds = ((t*fds2)-(fds*fds))/(t*(t-1));
		if (vds < 0.0001)
			vds = 0.0001;
		VARDS[count]  = vds;
	}

	hrate = (float)chw/t;
	farate = (float)cf/t;
	dprime_(&hrate, &farate, &dp);
	DPRPW[count]  = dp;

	hrate = (float)chs/t;
	if (strength == 3 )
		farate = (float)cfs/t;
	else
		farate = (float)cf/t;
	dprime_(&hrate, &farate, &dp);
	DPRPS[count]  = dp;

	DPRDW[count]  = ((ftw/t)-(fd/t))/(sqrt(VARD[count]));
	if (strength == 3 )
		DPRDS[count]  = ((fts/t)-(fds/t))/(sqrt(VARDS[count]));
	else
		DPRDS[count]  = ((fts/t)-(fd/t))/(sqrt(VARD[count]));

	count++;
}

/********************************
 * FUNCTION RECOGNITION_ASS_NEW *
 ********************************
 *
 * associative recognition simulation routine
 * targets    : A-B
 * distractors: X-Y
 *
 */
void recognition_ass_new(void)
{
	static real hrate, farate, dp;
	extern int dprime_();
	int i, j;
	int win;
	float t;
	float fam;                                    /* familiarity of test-item */
	int chw   = 0,
			chs   = 0,
			cf    = 0;
	float ftw   = 0.0,          /* variables to compute means and variances */
				ftw2  = 0.0,          /* for the target and distractor distibutions */
				fts   = 0.0,
				fts2  = 0.0,
				fd    = 0.0,
				fd2   = 0.0,
				vtw   = 0.0,
				vts   = 0.0,
				vdw   = 0.0;

	for (i=0; i<n/2; i++) {                                      /* targets 1 */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 2) {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}
		else {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 2) {
				count_hs++;
				chs++;
			}
			else {
				count_hw++;
				chw++;
			}
		}
	}

	for (i=n/2; i<n; i++) {                                      /* targets 2 */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		if (strength == 1) {
			famtw += fam;
			famtw2 += fam * fam;
			ftw += fam;
			ftw2 += fam * fam;
		}
		else {
			famts += fam;
			famts2 += fam * fam;
			fts += fam;
			fts2 += fam * fam;
		}

		win = decision(fam);
		if (win == 1) {
			if (strength == 1) {
				count_hw++;
				chw++;
			}
			else {
				count_hs++;
				chs++;
			}
		}
	}

	for (i=n; i<2*n; i++) {                                    /* distractors */
		fam = double_fam((2*i),(2*i)+1);
		fam += single_fam(2*i);
		fam += single_fam((2*i)+1);

		famd += fam;
		famd2 += fam * fam;
		fd += fam;
		fd2 += fam * fam;

		win = decision(fam);
		if (win == 1) {
			count_f++;
			cf++;
		}
	}
																		 /* compute results and store in arrays */
	if (strength == 3)
		t = (float)n/2;
	else
		t = (float)n;

	vtw = ((t*ftw2)-(ftw*ftw))/(t*(t-1));
	if (vtw < 0.0001)
		vtw = 0.0001;
	VARTW[count]  = vtw;

	vts = ((t*fts2)-(fts*fts))/(t*(t-1));
	if (vts < 0.0001)
		vts = 0.0001;
	VARTS[count]  = vts;

	vdw = (((float)n * fd2)-(fd*fd))/((float)n*((float)n-1));
	if (vdw < 0.0001)
		vdw = 0.0001;
	VARD[count]   = vdw;

	hrate = (float)chw/t;
	farate = (float)cf/(float)n;
	dprime_(&hrate, &farate, &dp);
	DPRPW[count]  = dp;

	hrate = (float)chs/t;
	farate = (float)cf/(float)n;
	dprime_(&hrate, &farate, &dp);
	DPRPS[count]  = dp;

	DPRDW[count]  = ((ftw/t)-(fd/(float)n))/(sqrt(VARD[count]));
	DPRDS[count]  = ((fts/t)-(fd/(float)n))/(sqrt(VARD[count]));

	count++;
}


/************************
 * FUNCTION CUED_RECALL *
 ************************
 *
 * the cued recall simulation routine
 */
void cued_recall(void)
{
	int i;
	int win;

	for (i=0; i<n; i++) {
		correlate(2*i);
		win = cor_win(i);

		if (win == ((2*i)+1)) {
			if (strength == 1)
				count_cw++;
			if (strength == 2)
				count_cs++;
			if (strength == 3 && i<n/2)
				count_cw++;
			if (strength == 3 && i>=n/2)
				count_cs++;
		}
	}
}                                                           /* end function */

/*********************
 * function CONVOLVE *
 *********************
 *
 * convolves two vectors, truncates to the number of features,
 * and centers the result
 *
 */
void convolve(int word1, int word2)
{
	int i, j, k;
	int shift;                             /* shift for truncated convolution */
	float temp[(FMAX*2)-1];         /* temporary storage for full convolution */

	for (i=0; i<(2*features)-1; i++)                            /* initialize */
		temp[i] = 0.0;
	for (i=0; i<features; i++)
		convolved[i] = 0.0;

	for (i=0; i<(2*features)-1; i++) {                            /* convolve */
		for (j=0; j<features; j++) {
			for (k=0; k<features; k++) {
				if (j+k == i) {
					temp[i] += items[word1][j]*items[word2][k];
				}
			}
		}
	}

	shift = (int)features/2;                                      /* truncate */
	for (i=0; i<features; i++)
		convolved[i] = temp[i+shift];
}                                                  /* end function CONVOLVE */

/**********************
 * function CORRELATE *
 **********************
 *
 * Correlates a vector with the memory vector, truncates to the
 * number of features, and centers the result.
 *
 * Invert the word first and than use convolution
 *
 */
void correlate(int word)
{
	int i, j, k;
	int shift;                             /* shift for truncated correlation */
	float temp[(FMAX*2)-1];         /* temporary storage for full correlation */
	float inv[FMAX];                                  /* inverted word vector */

	for (i=0; i<(2*features)-1; i++)                            /* initialize */
		temp[i] = 0.0;
	for (i=0; i<features; i++) {
		correlated[i] = 0.0;
		inv[i] = 0.0;
	}

	for (i=0; i<features; i++)                               /* invert vector */
		inv[i] = items[word][features-i-1];

	for (i=0; i<(2*features)-1; i++) {                            /* convolve */
		for (j=0; j<features; j++) {
			for (k=0; k<features; k++) {
				if (j+k == i) {
					temp[i] += (memory[j] * inv[k]);
				}
			}
		}
	}

	shift = (int)features/2;                                      /* truncate */
	for (i=0; i<features; i++)
		correlated[i] = temp[i+shift];
}                                                 /* end function CORRELATE */

/***********************
 * FUNCTION SINGLE_FAM *
 ***********************
 *
 * for a given single item probe, a familiarity is computed
 *
 */
float single_fam(int i)
{
	float fam;                      /* familiarity of test-item (dot-product) */
	int k;

	fam=0.0;
														/* familiarity values from first items in pairs */
	for (k=0; k<features; k++) {
		fam += (items[i][k] * memory[k]);
	}
	return fam;
}

/***********************
 * FUNCTION DOUBLE_FAM *
 ***********************
 *
 * for a pair probe, a familiarity is computed
 *
 */
float double_fam(int first, int second)
{
	float fam;                                    /* familiarity of test-item */
	int j,k;

	fam=0.0;
	convolve(first, second);
	for (k=0; k<features; k++) {
		fam += (convolved[k] * memory[k]);
	}
	return fam;
}

/*********************
 * FUNCTION DECISION *
 *********************
 *
 * decision system for recognition
 *
 */
int decision(float fam)
{
	int win;
	int done;
	float hi, lo;
	float temp;
	float adjust;

	hi = upper;
	lo = lower;

	if (noise != 0)
		temp = fam + gau((double)0,(double)noise);                 /* add noise */
	else
		temp = fam;

	if (temp >= hi)                                         /* "yes" response */
		win = 1;
	else if (temp < lo)                                      /* "no" response */
		win = 0;
	else {                                                            /* wait */
		done=0;
		do {
			temp = fam + gau((double)0,(double)noise);
			adjust = ((hi-lo) - ((hi-lo)*rate))/2;
			hi -= adjust;
			lo += adjust;
			if (temp > hi) {
				win = 1;
				done = 1;
			}
			if (temp < lo) {
				win = 0;
				done = 1;
			}
		} while (done == 0);
	}
	return win;
}                                                  /* end function DECISION */

/********************
 * FUNCTION COR_WIN *
 ********************
 *
 * The item in memory with the highest correlation (dot-product) to the
 * retrieved vector (correlated) is found
 *
 */
int cor_win(int wordnr)
{
	int i, k, t;
	float cor, wincor, total;
	int win;

	total=0.0;
	win=64;
	wincor=-10.0;
	for (i=0; i<2*n; i++) {
		cor=0.0;
		for (k=0; k<features; k++)
			cor += (correlated[k] * items[i][k]);
		total+=cor;
		if (cor>wincor) {
			win=i;
			wincor=cor;
		}
	}

	if (strength == 1)
		corw += total/(float)n;
	else if (strength == 2)
		cors += total/(float)n;
	else {
		if (wordnr < (n/2))
			corw += total/(float)n*(float)2.0;
		else
			cors += total/(float)n*(float)2.0;
	}
	if (wincor<giveup)
		win=64;
	return win;
}

/*****************
 * FUNCTION MENU *
 *****************
 *
 * August 1992
 * Peter Nobel
 *
 */
void menu(void)
{
	curpos=0;
	while (TRUE) {
		settings();
	}
}

/********************
 * FUNCTION ACTION1 *
 ********************
 *
 * performs action for MENU1 based on cursor position
 * Main menu
 */
void action1(int pos)
{
	curpos = 0;
	switch(pos) {
		case 0:
			back = 0;
			do {
				settings();
			} while (back == 0);
			break;
		case 1:
			back = 0;
			do {
				settings();
			} while (back == 0);
			break;
		case 2:
			back = 0;
			do {
				settings();
			} while (back == 0);
			break;
		case 3:
			back =0;
			do {
				settings();
			} while (back == 0);
			break;
		case 4:
			back =0;
			do {
				settings();
			} while (back == 0);
			break;
		case 5:
			results();
			break;
		case 6:
			exit(0);
	}
}

/********************
 * FUNCTION ACTION2 *
 ********************
 *
 * performs action for MENU2 based on cursor position
 * Type of test
 */
void action2(int pos)
{
	switch(pos) {
		case 0:
			cond = 1;
			features=15;
			power=1.0;
			alpha=0.9;
			gamma1=1.0;
			gamma2=1.0;
			gamma3=1.0;
			upper=0.1;
			lower=0.1;
			noise=0.0;
			rate=0.9;
			break;
		case 1:
			cond = 2;
			features=15;
			power=1.0;
			alpha=0.9;
			gamma1=1.0;
			gamma2=1.0;
			gamma3=1.0;
			upper=0.1;
			lower=0.1;
			noise=0.0;
			rate=0.9;
			break;
		case 2:
			cond = 3;
			features=15;
			power=1.0;
			alpha=0.9;
			gamma1=1.0;
			gamma2=1.0;
			gamma3=1.0;
			upper=0.1;
			lower=0.1;
			noise=0.0;
			rate=0.9;
			break;
		case 3:
			cond = 4;
			features=15;
			power=1.0;
			alpha=0.9;
			gamma1=1.0;
			gamma2=1.0;
			gamma3=1.0;
			upper=0.1;
			lower=0.1;
			noise=0.0;
			rate=0.9;
			break;
		case 4:
			cond = 5;
			features=15;
			power=1.0;
			alpha=0.9;
			gamma1=1.0;
			gamma2=1.0;
			gamma3=1.0;
			giveup=0.3;
			break;
		case 5:
			back = 1;
			curpos = 0;
			break;
	}
}

/********************
 * FUNCTION ACTION3 *
 ********************
 *
 * performs action for MENU3 based on cursor position
 * Length of list
 */
void action3(int pos)
{
	switch(pos) {
		case 0:
			printf("Number of pairs = ");
			scanf("%d", &n);
			if (n > NMAX) {
				printf("Maximum number of pairs is %d!", NMAX);
				n = NMAX;
			}
			if (n % 8 != 0) {                           /* if not a multiple of 8 */
				printf("Number of pairs should be a multiple of 8!");
			}
			break;
		case 1:
			back = 1;
			curpos = 0;
			break;
	}
}

/********************
 * FUNCTION ACTION4 *
 ********************
 *
 * performs action for MENU4 based on cursor position
 * Strength of list
 */
void action4(int pos)
{
	int done, i;

	switch(pos) {
		case 0:
			strength = 1;
			printf("Strength selected: %s.  Value = ", menu4[pos]);
			scanf("%d", &s1);
			if (s1 > SMAXW) {
				printf("Maximum strength for weak items is %d!", SMAXW);
				s1 = SMAXW;
			}
			else if (s1 < 1) {
				printf("Minimum strength for weak items is 1!");
				s1 = 1;
			}
			break;
		case 1:
			strength = 2;
			printf("Strength selected: %s.  Value = ", menu4[pos]);
			scanf("%d", &s2);
			if (s2 > SMAXS) {
				printf("Maximum strength for strong items is %d!", SMAXS);
				s2 = SMAXS;
			}
			else if (s2 < 2) {
				printf("Minimum strength for strong items is 2!");
				s2 = 2;
			}
			break;
		case 2:
			done=0;
			strength = 3;
			while (done==0) {
				for (i=0; i<1; i++) {
					printf("Strength selected: %s.  Weak item strength = ",
								menu4[pos]);
					scanf("%d", &s1);
					done=1;
					if (s1 > SMAXW) {
						printf("Maximum strength for weak items is %d!", SMAXW);
						done=0;
						break;
					}
					else if (s1 < 1) {
						printf("Minimum strength for weak items is 1!");
						done=0;
						break;
					}
					printf("Strength selected: %s.  Strong item strength = ",
								menu4[pos]);
					scanf("%d", &s2);
					if (s2 > SMAXS) {
						printf("Maximum strength for strong items is %d!", SMAXS);
						done=0;
						break;
					}
					else if (s2 <= s1) {
						printf("Strong items must be stronger than weak ones!");
						done=0;
						break;
					}
				}
			}
			break;
		case 3:
			back = 1;
			curpos = 0;
			break;
	}
}

/********************
 * FUNCTION ACTION5 *
 ********************
 *
 * performs action for MENU5 based on cursor position
 * Parameters of the model
 */
void action5(int pos)
{
	int i, done;
	char c;

	switch(pos) {
		case 0:
			printf("%s = ",menu5[0]);
			scanf("%d", &features);
			if (features % 2 == 0) {                                   /* if even */
				printf("Number of features should be odd!");
			}
			else {
				if (features > FMAX) {
					printf("Maximum number of features is %d!", FMAX);
					features = FMAX;
				}
				else if (features < 3) {
					printf("Minimum number of features is 3!");
					features = 3;
				}
			}
			break;
		case 1:
			printf("%s = ",menu5[1]);
			scanf("%f", &power);
			if (power < 0) {
				printf("Power can't be negative!");
				power=1.0;
			}
			break;
		case 2:
			printf("Alpha = ");
			scanf("%f", &alpha);
			if (alpha <= 0 || alpha > 1) {
				printf("Alpha is between 0 and 1!");
				alpha=1.0;
			}
			break;
		case 3:
			done=0;
			while (done==0) {
				for (i=0; i<1; i++) {
					printf("Weight item 1 = ");
					scanf("%f", &gamma1);
					done=1;
					if (gamma1 < 0 || gamma1 > 1) {
						printf("Weight should be between 0 and 1!");
						done=0;
						break;
					}
					printf("Weight item 2 = ");
					scanf("%f", &gamma2);
					done=1;
					if (gamma2 < 0 || gamma2 > 1) {
						printf("Weight should be between 0 and 1!");
						done=0;
						break;
					}
					printf("Weight associative information = ");
					scanf("%f", &gamma3);
					if (gamma3 < 0 || gamma3 > 1) {
						printf("Weight should be between 0 and 1!");
						done=0;
						break;
					}
				}
			}
			break;
		case 4:
			done=0;
			while (done==0) {
				for (i=0; i<1; i++) {
					printf("Lower criterion = ");
					scanf("%f", &lower);
					done=1;
					printf("Upper criterion = ");
					scanf("%f", &upper);
					 if (upper < lower) {
						printf("Upper criterion must be higher than lower!");
						done=0;
						break;
					}
				}
			}
			break;
		case 5:
			done=0;
			while (done==0) {
				for (i=0; i<1; i++) {
					printf("Noise variance = ");
					scanf("%f", &noise);
					done=1;
					if (noise < 0) {
						printf("Variance cannot be negative!");
						done=0;
						break;
					}
					printf("Convergence rate = ");
					scanf("%f", &rate);
					 if (rate <= 0 || rate >= 1) {
						printf("Rate must be between 0 and 1!");
						done=0;
						break;
					}
				}
			}
			break;
		case 6:
			printf("%s = ",menu5[6]);
			scanf("%f", &giveup);
			break;
		case 7:
			back = 1;
			curpos = 0;
			break;
	}
}

/********************
 * FUNCTION ACTION6 *
 ********************
 *
 * performs action for MENU6 based on cursor position
 * Simulation
 */
void action6(int pos)
{
	switch(pos) {
		case 0:
			if (nsim > 500) {
				printf("Maximum number of simulations is 500!");
				nsim = 500;
			}
			break;
		case 1:
			simulate();
			screen();
			back = 1;
			curpos = 0;
			break;
		case 2:
			back = 1;
			curpos = 0;
			break;
	}
}

/*********************
 * FUNCTION SETTINGS *
 *********************
 *
 * displays values for variables under the menus
 *
 */
void settings(void)
{
	printf("Current settings:\n");
	printf("\tTest:\t\t\t\t%s (%d sim.)\n", menu2[cond-1], nsim);
	printf("\tLength:\t\t\t\t%d\n", n);
	if (strength == 1)
		printf("\tStrength:\t\t\t%s; value = %d\n", menu4[strength-1], s1);
	if (strength == 2)
		printf("\tStrength:\t\t\t%s; value = %d\n", menu4[strength-1], s2);
	if (strength == 3)
		printf("\tStrength:\t\t\t%s; weak = %d, strong = %d\n",
							menu4[strength-1], s1, s2);
	printf("\t%s:\t\t%d\n", menu5[0], features);
	printf("\t%s:\t\t\t\t%2.2f\n", menu5[1], power);
	printf("\t%s:\t\tAlpha = %2.2f\n", menu5[2], alpha);
	printf("\t%s:\t\t\tGamma1= %2.2f, Gamma2= %2.2f, Gamma3= %2.2f\n",
		menu5[3], gamma1, gamma2, gamma3);
	if (cond == 5) {
		printf("\t%s:\t\tCrit. = %2.3f\n", menu5[6], giveup);
	}
	else {
		printf("\t%s:\t\tLower = %2.3f, Upper = %2.3f\n",
				menu5[4], lower, upper);
		printf("\t%s:\t\tNoise var. = %2.3f, Conv. rate = %2.3f\n",
				menu5[5], noise, rate);
	}
}

/*******************
 * FUNCTION SCREEN *
 *******************
 *
 * put settings() and results of simulation on screen
 *
 */
void screen(void)
{
	char junk;
	int i;
	float t, N, S;
	static real hrate, farate, dp;
	extern int dprime_();

																								 /* intialize the variables */
	pcw = pcs = 0.0;
	hw = mtw = vartw = tvartw = 0.0;
	hs = mts = varts = tvarts = 0.0;
	f  = md  = vard  = tvard  = 0.0;
	fs = mds = vards = tvards = 0.0;
	dprpw = tdprpw = dprps = tdprps = 0.0;
	dprdw = tdprdw = dprds = tdprds = 0.0;


	settings();
	N = (float)nsim;
	if (cond == 2 || cond == 3) {
		if (strength == 3)
			t = (float)n/4.0;
		else
			t = (float)n/2.0;
	}
	else {
		if (strength == 3)
			t = (float)n/2.0;
		else
			t = (float)n;
	}

	S = t*N;

	if (cond == 5) {                                   /* results cued recall */
		pcw=(float)count_cw/(N*t);
		pcs=(float)count_cs/(N*t);
		if (strength == 3) {
			printf("\n\nWeak items:\tP(C) = %.3f\n", pcw);
			printf("Average correlation = %1.3f\n\n",corw/(float)nsim);
			printf("Strong items:\tP(C) = %.3f\n", pcs);
			printf("Average correlation = %1.3f\n",cors/(float)nsim);
		}
		else if (strength == 2) {
			printf("\n\nStrong items:\tP(C) = %.3f\n", pcs);
			printf("Average correlation = %1.3f\n",cors/(float)nsim);
		}
		else {
			printf("\n\nWeak items:\tP(C) = %.3f\n", pcw);
			printf("Average correlation = %1.3f\n",corw/(float)nsim);
		}
	}

	else {                                             /* results recognition */
		for (i=0; i<N; i++) {
			vartw += VARTW[i];
			varts += VARTS[i];
			vard  += VARD[i];
			vards += VARDS[i];
			dprpw += DPRPW[i];
			dprps += DPRPS[i];
			dprdw += DPRDW[i];
			dprds += DPRDS[i];
		}
		vartw   /= N;
		varts   /= N;
		vard    /= N;
		vards   /= N;
		dprpw   /= N;
		dprps   /= N;
		dprdw   /= N;
		dprds   /= N;

		hw      = (float)count_hw/S;
		mtw     = famtw/S;
		tvartw  = ((S*famtw2)-(famtw*famtw))/(S*(S-1));

		hs      = (float)count_hs/S;
		mts     = famts/S;
		tvarts  = ((S*famts2)-(famts*famts))/(S*(S-1));

		if (cond == 1 || cond == 4) {
			t     = (float)n;
			S     = t*N;
		}

		f       = (float)count_f/S;
		md      = famd/S;
		tvard   = ((S*famd2)-(famd*famd))/(S*(S-1));

		if ((cond == 2 || cond == 3) && strength == 3) {
			fs     = (float)count_fs/S;
			mds    = famds/S;
			tvards = ((S*famds2)-(famds*famds))/(S*(S-1));
		}

		hrate   = hw;
		farate  = f;
		dprime_(&hrate, &farate, &dp);
		tdprpw  = dp;

		hrate   = hs;
		if ((cond == 2 || cond == 3) && strength == 3)
			farate  = fs;
		else
			farate  = f;
		dprime_(&hrate, &farate, &dp);
		tdprps  = dp;

		tdprdw  = (mtw-md)/(sqrt(tvard));
		if ((cond == 2 || cond == 3) && strength == 3)
			tdprds  = (mts-mds)/(sqrt(tvards));
		else
			tdprds  = (mts-md)/(sqrt(tvard));

		printf("\n\nResults averaged over simulations, and for total simulation"
						 " (in brackets)\n");
		printf("\t\tTargets\t\t\t\tDistractors\n");
		if (strength == 1) {
			printf("\nproportion\t%.3f\t\t\t\t%.3f\n", hw,f);
			printf("mean fam.\t%.3f\t\t\t\t%.3f\n", mtw, md);
			printf("variance\t%.3f (%.3f)\t\t\t%.3f (%.3f)\n\n",
						vartw, tvartw, vard, tvard);
			printf("d' (prop.)\t%.3f (%.3f)\n\n", dprpw, tdprpw);
			printf("d' (distr.)\t%.3f (%.3f)\n", dprdw, tdprdw);
		}
		if (strength == 2) {
			printf("\nproportion\t%.3f\t\t\t\t%.3f\n", hs,f);
			printf("mean fam.\t%.3f\t\t\t\t%.3f\n", mts, md);
			printf("variance\t%.3f (%.3f)\t\t\t%.3f (%.3f)\n\n",
						varts, tvarts, vard, tvard);
			printf("d' (prop.)\t%.3f (%.3f)\n\n", dprps, tdprps);
			printf("d' (distr.)\t%.3f (%.3f)\n", dprds, tdprds);
		}
		if (strength == 3) {
			if (cond == 2 || cond == 3)
				printf("\t\tweak\t\tstrong\t\tweak\t\tstrong\n");
			else
				printf("\t\tweak\t\tstrong\n");

			if (cond == 2 || cond == 3) {
				printf("proportion\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", hw, hs, f, fs);
				printf("mean fam.\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",
							mtw, mts, md, mds);
				printf("variance\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f (%.3f)"
							"\n\n",vartw, tvartw,varts,tvarts,vard,tvard,vards,tvards);
			}
			else {
				printf("proportion\t%.3f\t\t%.3f\t\t%.3f\n", hw, hs, f);
				printf("mean fam.\t%.3f\t\t%.3f\t\t%.3f\n", mtw, mts, md);
				printf("variance\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f (%.3f)\n\n",
							vartw, tvartw, varts, tvarts, vard, tvard);
			}

			printf("d' (prop.)\t%.3f (%.3f)\t%.3f (%.3f)\n\n",
						dprpw, tdprpw, dprps, tdprps);
			printf("d' (distr.)\t%.3f (%.3f)\t%.3f (%.3f)\n",
						dprdw, tdprdw, dprds, tdprds);
		}
	}

}

/********************
 * FUNCTION RESULTS *
 ********************
 *
 * write the results from the simulation to a file
 *
 */
void results(void)
{
	int i;
	char junk;

	for (i=0; i<100; i++) {                        /* zero out filename array */
		fout[i] = 0;
	}

	printf("Enter name of output file: ");
	scanf("%s", fout);
	fp=fopen(fout,"a");                               /* write or append file */
	if (fp == NULL) {
		printf("\n\nERROR opening output file %s\n", fout);
	}

	fprintf(fp,"\nCurrent settings:\n");
	fprintf(fp,"\tTest:\t\t\t\t%s (%d sim.)\n", menu2[cond-1], nsim);
	fprintf(fp,"\tLength:\t\t\t\t%d\n", n);
	if (strength == 1)
		fprintf(fp,"\tStrength:\t\t\t%s; value = %d\n", menu4[strength-1], s1);
	if (strength == 2)
		fprintf(fp,"\tStrength:\t\t\t%s; value = %d\n", menu4[strength-1], s2);
	if (strength == 3)
		fprintf(fp,"\tStrength:\t\t\t%s; weak = %d, strong = %d\n",
							menu4[strength-1], s1, s2);
	fprintf(fp,"\t%s:\t\t%d\n", menu5[0], features);
	fprintf(fp,"\t%s:\t\t\t\t%2.2f\n", menu5[1], power);
	fprintf(fp,"\t%s:\t\tAlpha = %2.2f\n", menu5[2], alpha);
	fprintf(fp,"\t%s:\t\t\tGamma1=%2.2f, Gamma2=%2.2f, Gamma3=%2.2f\n",
				menu5[3], gamma1, gamma2, gamma3);

	if (cond == 5)
		fprintf(fp,"\t%s:\t\tCrit. = %2.3f\n", menu5[6], giveup);
	else {
		fprintf(fp,"\t%s:\t\tLower = %2.3f, Lower = %2.3f\n",
				menu5[4], lower, upper);
		fprintf(fp,"\t%s:\t\tNoise var. = %2.3f, Conv. rate = %2.3f\n",
				menu5[5], noise, rate);
	}

	if (cond == 5) {
		if (strength == 3) {
			fprintf(fp,"\n\nWeak items:\tP(C) = %.3f\n", pcw);
			fprintf(fp,"Average correlation = %1.3f\n\n",corw/(float)nsim);
			fprintf(fp,"Strong items:\tP(C) = %.3f\n", pcs);
			fprintf(fp,"Average correlation = %1.3f\n",cors/(float)nsim);
		}
		else if (strength == 2) {
			fprintf(fp,"\n\nStrong items:\tP(C) = %.3f\n", pcs);
			fprintf(fp,"Average correlation = %1.3f\n\n",cors/(float)nsim);
		}
		else {
			fprintf(fp,"\n\nWeak items:\tP(C) = %.3f\n", pcw);
			fprintf(fp,"Average correlation = %1.3f\n",corw/(float)nsim);
		}
	}
	else {
		fprintf(fp,"\nResults averaged over simulations, and for"
						 " total simulation (in brackets)\n");
		fprintf(fp,"\t\tTargets\t\t\t\tDistractors\n");

		if (strength == 1) {
			fprintf(fp,"\nproportion\t%.3f\t\t\t\t%.3f\n", hw,f);
			fprintf(fp,"mean fam.\t%.3f\t\t\t\t%.3f\n", mtw, md);
			fprintf(fp,"variance\t%.3f (%.3f)\t\t\t%.3f (%.3f)\n\n",
						vartw, tvartw, vard, tvard);
			fprintf(fp,"d' (prop.)\t%.3f (%.3f)\n\n", dprpw, tdprpw);
			fprintf(fp,"d' (distr.)\t%.3f (%.3f)\n", dprdw, tdprdw);
		}

		if (strength == 2) {
			fprintf(fp,"\nproportion\t%.3f\t\t\t\t%.3f\n", hs,f);
			fprintf(fp,"mean fam.\t%.3f\t\t\t\t%.3f\n", mts, md);
			fprintf(fp,"variance\t%.3f (%.3f)\t\t\t%.3f (%.3f)\n\n",
						varts, tvarts, vard, tvard);
			fprintf(fp,"d' (prop.)\t%.3f (%.3f)\n\n", dprps, tdprps);
			fprintf(fp,"d' (distr.)\t%.3f (%.3f)\n", dprds, tdprds);
		}

		if (strength == 3) {
			if (cond == 2 || cond == 3) {
				fprintf(fp,"\t\tweak\t\tstrong\t\tweak\t\tstrong\n");
				fprintf(fp,"proportion\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",
							hw, hs, f, fs);
				fprintf(fp,"mean fam.\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",
							mtw, mts, md, mds);
				fprintf(fp,"variance\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f"
							" (%.3f)\n\n",
							vartw, tvartw,varts,tvarts,vard,tvard,vards,tvards);
			}
			else {
				fprintf(fp,"\t\tweak\t\tstrong\n");
				fprintf(fp,"proportion\t%.3f\t\t%.3f\t\t%.3f\n", hw, hs, f);
				fprintf(fp,"mean fam.\t%.3f\t\t%.3f\t\t%.3f\n", mtw, mts, md);
				fprintf(fp,"variance\t%.3f (%.3f)\t%.3f (%.3f)\t%.3f (%.3f)\n\n",
							vartw, tvartw, varts, tvarts, vard, tvard);
			}
			fprintf(fp,"d' (prop.)\t%.3f (%.3f)\t%.3f (%.3f)\n\n",
						dprpw, tdprpw, dprps, tdprps);
			fprintf(fp,"d' (distr.)\t%.3f (%.3f)\t%.3f (%.3f)\n",
						dprdw, tdprdw, dprds, tdprds);
		}
	}

	fclose(fp);
}

/*******************
 * FUNCTION RANSET *
 *******************
 *
 * set up the random seed
 */
void ranset(void)
{
	struct timeb timebuffer;

	ftime(&timebuffer);
	srand(timebuffer.millitm);

}

/*****************
 * FUNCTION RAND *
 *****************
 *
 * random number generator between 0 and 1
 *
 */
float rand1(void)
{
	return((float)rand()/RAND_MAX);
}

/****************
 * FUNCTION GAU *
 ****************
 *
 * gaussian distribution, truncated at 0.0
 *
 */
float gau(double mu, double sigma)
{
	float random_array[12];
	int i;
	float total, result;

	total =0.0;
	for (i=0; i<=11; i++) {
		random_array[i] = rand1();
		total += random_array[i];
	}

	result = (sqrt(sigma) * ( total - 6.0)) + mu;

	return(result);
}                                                       /* end function GAU */

/*****************
 * function PERM *
 *****************
 *
 * Takes and array (list) of length (len) and permutes the order of its
 * elements.
 *
 * The length (len) from PERM() becomes the range delimiter (n) in RANDW()
 */
void perm(int *list, int len)
{
	int i, j, temp;

	for (i = 0; i < len; i++) {
		j = randw(len);
		temp = list[j];
		list[j] = list[i];
		list[i] = temp;
	}

}                                                      /* end function PERM */

/******************
 * function RANDW *
 ******************
 *
 * Calls the random number generator RAND() and uses the modulus (%)
 * operator to limit the random numbers to a specific range
 */
int randw(int n)
{
	return(rand()%n);

}                                                     /* end function RANDW */

/*******************
 * FUNCTION DPRIME *
 *******************
 *
 * computes d-prime based on hit and false alarm rate
 *
 */
int dprime_(hr, far1, d)

real *hr, *far1, *d;

{
		static real y, z1, z2, y2;
		extern /* Subroutine */ int znorma_();

		if (*hr == (float)1) {
	*hr = (float).99;
		}
		if (*hr == (float)0) {
	*hr = (float).01;
		}
		if (*far1 == (float)1) {
	*far1 = (float).99;
		}
		if (*far1 == (float)0) {
	*far1 = (float).01;
		}
		znorma_(hr, &z1, &y);
		znorma_(far1, &z2, &y2);
		*d = z1 - z2;
		return 0;
} /* dprime_ */

/* Subroutine */ int znorma_(p, zz, y)
real *p, *zz, *y;
{
		/* Initialized data */

		static real c1 = (float)2.515517;
		static real c2 = (float).802853;
		static real c3 = (float).010328;
		static real c4 = (float)1.432788;
		static real c5 = (float).189269;
		static real c6 = (float).001308;

		/* System generated locals */
		real r_1;

		/* Builtin functions */
		double log(), sqrt(), r_sign(), exp();

		/* Local variables */
		static real t, ab, cc, dd, pq, sz;

		if ((r_1 = *p - (float).5) < (float)0.) {
	goto L11;
		} else if (r_1 == 0.0) {
	goto L10;
		} else {
	goto L9;
		}
L9:
		pq = (float)1. - *p;
		sz = (float)1.;
		goto L12;
L10:
		*zz = (float)0.;
		goto L99;
L11:
		pq = *p;
		sz = (float)-1.;
L12:
		t = sqrt(log(pq) * (-2.0));
		ab = ((c6 * t + c5) * t + c4) * t + (float)1.;
		cc = (c3 * t + c2) * t + c1;
		dd = cc / ab;
		*zz = t - dd;
		*zz = r_sign(zz, &sz);
L99:
		*y = exp(-(double)(*zz) * (*zz) / (float)2.) * (float).39894228;
		return 0;
}                                                                /* znorma_ */

double r_sign(a, b)
real *a, *b;
{
	if (*b<0.0) return(-1.0*fabs(*a));
	else return(fabs(*a));
}