/*  koh.c  */

/*  Return codes:
 *  0 - ok
 *  1 - error
 *  2 - meerdere entries op 1 veld
 */

#define __NO_MATH_INLINES

#ifdef __WIN32__
#define my_PATH_SEP '\\'
#else
#define my_PATH_SEP '/'
#endif

#ifdef __MSDOS__
#ifndef __COMPACT__
#error Memory model COMPACT required
#endif /* __COMPACT__ */
#include <alloc.h>
#include <dir.h>
#endif /* __MSDOS__ */

#include <ctype.h>
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <values.h>

#ifndef MAXFLOAT
#include <float.h>
#define MAXFLOAT FLT_MAX
#endif

#define BUFSIZE 2048

#define MIN(A, B) ((A) < (B) ? (A) : (B))

typedef struct {
    int
        x1,
        y1,
        x2,
        y2;
     float
        f;
} LINK;

LINK
    *link;
int
    link_n,
    link_max;

typedef struct
{
    char
	*name;
    float
	*v;
    int
	x,
        y,
        used,
        cl1,
        cl2,
        cluster;
} VEC;

typedef struct
{
    float
	f;
    int
	x,
	y;
    char
	c;
} LINE;

VEC
    *vec = NULL;

LINE
    *line;

float
    ***map,
    **inc,
    **vdiff;

int
    *cluster,
    epoch,
    max_epoch = 100,
    vec_max = 0,
    vec_n = 0,
    vec_size,
    map_max = 0,
    map_size = 0,
    verbose = 1,
    input_line;

char
    *programname,
    buffer [BUFSIZE],
    buf1 [BUFSIZE],
    buf2 [BUFSIZE],
    buf3 [BUFSIZE],
    *infile = NULL,
    *baseout = NULL,
    *no_mem_buffer,
    Out_of_memory [] = "Out of memory";

FILE
    *fp_log;

void
    get_programname (char const *argv0),
    syntax (void),
    inc_mapsize (void),
    process_args (int argc, char *argv []),
    setup (void),
    read_files (void),
    update_nn (int i),
    errit (char const *format, ...),
    *m_malloc (size_t size),
    *m_realloc (void *block, size_t size),
    hussel (void),
    Heapsort (
	void *base,
	size_t nelem,
	size_t width,
	int (*fcmp)(void *, void *)
    );
char
    *get_arg (int *argc, char *argv [], int *index),
    *m_strdup (const char *s),
    *to_ps (char const *string);
int
    save_files (void),
    vecpos_cmp (void *, void *),
    line_cmp (void *, void *),
    getline (FILE *fp, int required, char const *filename);
float
    square_vecdiff (float *v1, float *v2);
FILE
    *r_fopen (char const *filename),
    *w_fopen (char const *filename);
VEC
    *store_vec (char *name);

int main (int argc, char *argv [])
{
    int
	i,
	j,
	x,
	y,
	xx,
	yy,
	incr,
	size,
        startsize,
        retval;
    float
	min,
	diff,
        alfa;

    no_mem_buffer = (char *)malloc (1024);

    get_programname (argv [0]);

    process_args (argc, argv);

    if (! infile)
        errit ("Naam input-file ontbreekt");
    if (! baseout)
	baseout = "koh_out";

    read_files ();

    setup ();

    sprintf (buffer, "%s.log", baseout);
    fp_log = w_fopen (buffer);

    sprintf (
	buffer,
	"Mapsize: %ix%i\n"
	"Vectorlength: %i\n"
	"Number of vectors: %i\n"
	"Number of epochs: %i\n",
	map_max,
	map_max,
	vec_size,
	vec_n,
	max_epoch
    );
    fprintf (fp_log, buffer);
    fflush (fp_log);
    if (verbose) {
	printf (buffer);
	fflush (stdout);
    }

    startsize = map_size;
    for (epoch = 1; epoch <= max_epoch; epoch++) {
        alfa = 1.0 - 0.95 * ((float) epoch) / ((float) max_epoch);
/*      alfa = 1.0 / (float) epoch;     */
	size = (int)(
	    2.0 +
	    ((float)(epoch - 1)) *
            ((float)(map_max - startsize)) /
            ((float)(max_epoch / 2))
	);
	if (size > map_size && map_size < map_max) {
	    inc_mapsize ();
	    incr = 1;
	} else
	    incr = 0;
	if ((epoch == 1)
	 || (max_epoch < 10)
	 || (epoch % (max_epoch / 10) == 0)
	 || (epoch == max_epoch)
	 || (incr)) {
	    fprintf (
		fp_log,
		"epoch: %i\tmap: %ix%i\talfa: %f\n",
		epoch,
		map_size,
		map_size,
		alfa
	    );
	    fflush (fp_log);
	}
	if (verbose) {
	    printf (
		" epoch: %i  map: %ix%i  alfa: %f\r",
		epoch,
		map_size,
		map_size,
		alfa
	    );
	    fflush (stdout);
	}
	hussel ();
	for (i = 0; i < vec_n; i++) {
	    min = MAXFLOAT;
	    for (x = 1; x <= map_size; x++)
		for (y = 1; y <= map_size; y++) {
                    diff = square_vecdiff (map [x][y], vec [i].v);
		    if (diff < min) {
			min = diff;
			xx = x;
			yy = y;
		    }
		}
	    for (x = (xx <= 1) ? 1 : (xx - 1);
		   x <= xx + 1 && x <= map_size;
		     x++)
		for (y = (yy <= 1) ? 1 : (yy - 1);
		       y <= yy + 1 && y <= map_size;
			 y++)
		    for (j = 0; j < vec_size; j++)
			map [x][y][j] +=
			    alfa
			  * (vec [i].v [j] - map [x][y][j])
			  * ((x == xx) ? 1.0 : 0.6)
			  * ((y == yy) ? 1.0 : 0.6);
	}
    }

    if (verbose)
	printf ("\a\n");

    retval = save_files ();

    fclose (fp_log);

    return retval;
}

void inc_mapsize ()
{
    int
	i,
	x,
	y,
	x1,
	x2,
	y1,
	y2;
    float
	xp,
	xd,
	yp,
	yd;

    map_size++;

    for (i = 0; i < vec_size; i++) {

	for (x = 1; x < map_size; x++)
	    for (y = 1; y < map_size; y++)
		inc [x][y] = map [x][y][i];

	for (x = 1; x <= map_size; x++) {
	    xp = 1.0 +
		 ((float) (x - 1)) /
		 ((float) map_size) *
		 ((float) (map_size - 1));
	    x1 = floor (xp);
	    x2 = ceil (xp);
	    if (x2 == map_size)
		x2 = x1;
	    xd = xp - x1;
	    for (y = 1; y <= map_size; y++) {
		yp = 1.0 +
		     ((float) (y - 1)) /
		     ((float) map_size) *
		     ((float) (map_size - 1));
		y1 = floor (yp);
		y2 = ceil (yp);
		if (y2 == map_size)
		    y2 = y1;
		yd = yp - y1;
		map [x][y][i] =
		    inc [x1][y1] * (1.0 - xd) * (1.0 - yd) +
		    inc [x2][y1] * xd * (1.0 - yd) +
		    inc [x1][y2] * (1.0 - xd) * yd +
		    inc [x2][y2] * xd * yd;
	    }
	}
    }
}


void read_files ()
{
    FILE
	*fp;
    VEC
	*v;
    int
	i;

    fp = r_fopen (infile);
    getline (fp, 1, infile);
    if (sscanf (buffer, "%i", &vec_size) != 1)
	errit (
	    "file \"%s\", line %i\nVectorsize expected",
	    infile,
	    input_line
	);
    while (getline (fp, 0, NULL)) {
	if (sscanf (buffer, " %s", buf1) != 1)
	    errit ("file \"%s\", line %i", infile, input_line);
	v = store_vec (buf1);
	v->v = (float *) m_malloc (vec_size * sizeof (float));
	for (i = 0; i < vec_size; i++) {
	    getline (fp, 1, infile);
	    if (sscanf (buffer, "%f", &(v->v[i])) != 1)
		errit (
		    "file \"%s\", line %i\n"
		    "Missing value for vector \'%s\'",
		    infile,
		    input_line,
		    v->name
		);
	}
    }
    fclose (fp);
}

int save_files ()
{
    int
	i,
	j,
        k,
        p1,
        p2,
	x,
	y,
	xx,
	yy,
	l,
	retval,
        worst;
    float
        d,
	f,
	min,
	maxmin,
	diff,
	maxdiff,
	som,
        vecmin,
        vecmax;
    FILE
	*fp,
	*fp_ps;

    retval = 0;

    sprintf (buf3, "%s.map", baseout);
    fp = w_fopen (buf3);
    fprintf (
	fp,
	"# Mapsize\n"
	"\t%i\n"
	"# Vectorlength\n"
	"\t%i\n"
	"# Mapvectors\n",
	map_size,
        vec_size
    );
    for (x = 1; x <= map_size; x++)
        for (y = 1; y <= map_size; y++) {
            fprintf (fp, "# X: %i  Y: %i\n", x, y);
            for (i = 0; i < vec_size; i++)
                fprintf (fp, "\t%f\n", map [x][y][i]);
        }
    fclose (fp);

    sprintf (buf3, "%s.dif", baseout);
    fp = w_fopen (buf3);
    fprintf (fp, "# Mapsize\n%i\n# Differences\n", map_size);
    l = 0;
    for (y = 1; y <= map_size; y++)
        for (x = 1; x < map_size; x++) {
            f = sqrt (square_vecdiff (map[x][y], map[x+1][y])) / sqrt (vec_size);
	    fprintf (
		fp,
		"%f\t# %i,%i - %i,%i\n",
                f,
		x,
		y,
		x+1,
		y
	    );
            line [l].f = f;
            line [l].x = x;
            line [l].y = y;
            line [l++].c = 'r';
        }
    for (x = 1; x <= map_size; x++)
        for (y = 1; y < map_size; y++) {
            f = sqrt (square_vecdiff (map[x][y], map[x][y+1])) / sqrt (vec_size);
	    fprintf (
		fp,
		"%f\t# %i,%i - %i,%i\n",
                f,
		x,
		y,
		x,
		y+1
	    );
            line [l].f = f;
            line [l].x = x;
            line [l].y = y;
            line [l++].c = 'u';
        }
    fclose (fp);

    Heapsort (line, l, sizeof (LINE), line_cmp);
    maxdiff = line [l - 1].f;

    sprintf (buf3, "%s.ps", baseout);
    fp_ps = w_fopen (buf3);
    fprintf (
        fp_ps,
        "%c!PS-Adobe-3.0 EPSF-3.0\n"
        "%c%cBoundingBox: 99 99 %i %i\n"
        "%c%cCreator: %s, (c) P. Kleiweg 1995 - 1999\n"
        "%c%cTitle: %s\n"
        "%c%cEndComments\n",
        '%', '%', '%',
        101 + 30 * map_size,
        101 + 30 * map_size,
        '%', '%', programname,
        '%', '%', infile,
        '%', '%'
    );
    fprintf (
	fp_ps,
	"\n"
        "40 dict begin\n"
        "\n"
	"/X0 100 def\n"
	"/Y0 100 def\n"
	"/DX 30 def\n"
	"/DY 30 def\n"
	"\n"
	"/Times-Roman findfont 10 scalefont setfont\n"
	"\n"
        "/SHOW_GRAY true def\n"
        "\n"
        "/SHOW_LINKS true def\n"
        "/ANGLE -20 def\n"
	"\n"
    );
    fprintf (
        fp_ps,
	"/center {\n"
	"  dup stringwidth pop\n"
	"  2 div neg 0 rmoveto\n"
	"} bind def\n"
	"\n"
     );
    fprintf (
        fp_ps,
	"/box {\n"
	"  /mapsize exch def\n"
	"  0.2 setlinewidth\n"
	"  0 setgray\n"
	"  X0 Y0 moveto\n"
	"  DX mapsize mul 0 rlineto\n"
	"  0 DY mapsize mul rlineto\n"
	"  DX mapsize mul neg 0 rlineto\n"
	"  closepath\n"
        "  stroke\n"
        "} bind def\n"
        "\n"
    );
    fprintf (
        fp_ps,
        "/col {\n"
        "  SHOW_GRAY {\n"
        "    %f div\n"
        "    dup 1 gt {\n"
        "      pop 1\n"
        "    } if\n"
        "    1 exch sub\n"
        "    setgray\n"
        "    2 setlinewidth\n"
        "  } {\n"
        "    pop\n"
        "    0 setgray\n"
        "    0.2 setlinewidth\n"
        "  } ifelse\n"
        "} bind def\n"
        "\n",
        maxdiff * 0.7
    );
    fprintf (
        fp_ps,
        "/r {\n"
        "  col\n"
        "  /y exch def\n"
        "  /x exch def\n"
        "  X0 x DX mul add\n"
        "  Y0 y 1 sub DY mul add moveto\n"
        "  0 DY rlineto\n"
        "  stroke\n"
        "} bind def\n"
        "\n"
    );
    fprintf (
        fp_ps,
        "/u {\n"
        "  col\n"
        "  /y exch def\n"
        "  /x exch def\n"
        "  X0 x 1 sub DX mul add\n"
	"  Y0 y DY mul add moveto\n"
        "  DX 0 rlineto\n"
        "  stroke\n"
	"} bind def\n"
        "\n"
    );
    fprintf (
        fp_ps,
        "/linkmax 0 def\n"
        "/l {\n"
        "  /dd exch def\n"
        "  /y2 exch def\n"
        "  /x2 exch def\n"
        "  /y1 exch def\n"
        "  /x1 exch def\n"
        "  SHOW_LINKS {\n"
        "    [ dd dup mul linkmax dup mul div 5 mul\n"
        "      dup 6 exch sub\n"
        "      exch ] 0 setdash\n"
        "    0.2 setlinewidth\n"
        "    0 setgray\n"
        "      x1 x2 gt\n"
        "      x1 x2 eq y1 y2 gt and\n"
        "    or {\n"
        "      /x1 x2 /x2 x1 def def\n"
        "      /y1 y2 /y2 y1 def def\n"
        "    } if\n"
        "    x1 0.5 sub DX mul X0 add\n"
        "    y1 0.5 sub DY mul Y0 add moveto\n"
        "      x1 x2 sub 1 eq\n"
        "      x2 x1 sub 1 eq\n"
        "      y1 y2 sub 1 eq\n"
        "      y2 y1 sub 1 eq\n"
        "    or or or {\n"
        "      x2 0.5 sub DX mul X0 add\n"
        "      y2 0.5 sub DY mul Y0 add lineto\n"
        "      stroke\n"
        "    } {\n"
        "      /xx0 x1 0.5 sub DX mul X0 add def\n"
        "      /yy0 y1 0.5 sub DY mul Y0 add def\n"
        "      /xx3 x2 0.5 sub DX mul X0 add def\n"
        "      /yy3 y2 0.5 sub DY mul Y0 add def\n"
        "      /alfa yy3 yy0 sub xx3 xx0 sub atan def\n"
        "      /len\n"
        "        xx0 xx3 sub dup mul yy0 yy3 sub dup mul add sqrt\n"
        "        1 2 ANGLE cos mul add\n"
        "        div\n"
        "      def\n"
        "      /xx1 xx0 len alfa ANGLE add cos mul add def\n"
        "      /yy1 yy0 len alfa ANGLE add sin mul add def\n"
        "      /xx2 xx3 len alfa 180 add ANGLE sub cos mul add def\n"
        "      /yy2 yy3 len alfa 180 add ANGLE sub sin mul add def\n"
        "      xx1 yy1 xx2 yy2 xx3 yy3 curveto stroke\n"
        "    } ifelse\n"
        "    [ ] 0 setdash\n"
        "  } if\n"
        "} bind def\n"
	"\n"
    );
    fprintf (
        fp_ps,
        "/n {\n"
        "  /s exch def\n"
        "  /y exch def\n"
        "  /x exch def\n"
        "  0 setgray\n"
        "  X0 x DX mul add Y0 y DY mul add moveto\n"
        "  DX 2 div neg DY 2 div neg SHIFT sub rmoveto\n"
        "  s center show\n"
        "} bind def\n"
	"\n"
    );
    fprintf (
        fp_ps,
        "/SHIFT\n"
        "gsave\n"
        "  newpath\n"
        "  0 0 moveto\n"
        "  (-) false charpath\n"
        "  pathbbox\n"
        "grestore\n"
        "exch pop add exch pop\n"
        "2 div\n"
        "def\n"
        "\n"
    );

    maxmin = som = 0.0;
    for (i = 0; i < vec_n; i++) {
        min = MAXFLOAT;
	for (x = 1; x <= map_size; x++)
	    for (y = 1; y <= map_size; y++) {
                diff = sqrt (square_vecdiff (map [x][y], vec [i].v));
                if (diff < min) {
                    min = diff;
		    xx = x;
		    yy = y;
		}
	    }
        vec [i].x = xx;
        vec [i].y = yy;
        som += min;
        if (min > maxmin) {
            maxmin = min;
            worst = i;
        }
    }
    fprintf (
        fp_log,
        "Badness: %f\n"
        "Worst:   %f (%s)\n",
        som,
        maxmin,
        vec [worst].name
    );
    if (verbose)
	printf (
            "Badness: %f\n"
            "Worst:   %f (%s)\n",
            som,
            maxmin,
            vec [worst].name
        );
    Heapsort (vec, vec_n, sizeof (VEC), vecpos_cmp);

    for (i = 0; i < l; i++)
        fprintf (
            fp_ps,
            "%i %i %f %c\n",
            line [i].x,
            line [i].y,
            line [i].f,
            line [i].c
        );

    if (verbose) {
        printf ("Creating minimal spanning tree");
        fflush (stdout);
    }
    sprintf (buf3, "%s.spn", baseout);
    fp = w_fopen (buf3);
    
    for (i = 0; i < vec_n; i++) {
        vdiff [i][i] = 0.0;
        for (j = 0; j < i; j++)
            vdiff [j][i] = vdiff [i][j] = sqrt (square_vecdiff (vec [i].v, vec [j].v));
    }

    link_n = 0;
    link_max = vec_n;
    link = (LINK *) m_malloc (link_max * sizeof (LINK));
    
    if (vec_max < 2 * vec_n - 1) {
        vec_max = 2 * vec_n - 1;
        vec = (VEC *) m_realloc (vec, vec_max * sizeof (VEC));
    }
    for (i = 0; i < vec_n; i++) {
         vec [i].used = 0;
         vec [i].cl1 = vec [i].cl2 = -1;
         vec [i].cluster = i;
    }

    for (i = vec_n; i < 2 * vec_n - 1; i++) {
        vec [i].used = 0;
        d = MAXFLOAT;
        for (j = 0; j < i; j++)
            if (! vec [j].used)
  	        for (k = 0; k < j; k++)
	            if ((! vec [k].used) && (vdiff [j][k] < d)) {
                        p1 = j;
                        p2 = k;
                        d = vdiff [j][k];
	            }
        vec [p1].used = vec [p2].used = 1;
        vec [i].cl1 = p1;
        vec [i].cl2 = p2;
        vec [i].cluster = i;
        if (i == vec_n)
            vecmin = d;

        for (j = 0; j < vec_n; j++)
            if (vec [j].cluster == p1)
                for (k = 0; k < vec_n; k++)
		    if (vec [k].cluster == p2 &&
                        vdiff [j][k] == d &&
                        (vec [j].x != vec [k].x || vec [j].y != vec [k].y)
                    ) {
                        fprintf (
                            fp,
                            "%i %i %i %i %f\n",
                            vec [j].x,
                            vec [j].y,
                            vec [k].x,
                            vec [k].y,
                            d
                        );
			if (link_n == link_max) {
                            link_max += vec_n;
                            link = (LINK *) m_realloc (link, link_max * sizeof (LINK));
			}
                        link [link_n].x1 = vec [j].x;
                        link [link_n].y1 = vec [j].y;
                        link [link_n].x2 = vec [k].x;
                        link [link_n].y2 = vec [k].y;
                        link [link_n++].f = d;
		    }

        update_nn (i);

        if (verbose) {
            printf (".");
            fflush (stdout);
        }
    }
    vecmax = d;

    fprintf (
        fp_ps,
        "/linkmin %f def\n"
        "/linkmax %f def\n",
        vecmin,
        vecmax
    );

    for (i = 0; i < link_n; i++)
        fprintf (
            fp_ps,
            "%i %i %i %i %f l\n",
            link [i].x1,
            link [i].y1,
            link [i].x2,
            link [i].y2,
            link [i].f
        );

    if (verbose)
        printf ("\n");

    sprintf (buf3, "%s.top", baseout);
    fp = w_fopen (buf3);
    fprintf (fp, "# Mapsize\n%i\n# Positions\n", map_size);
    for (i = 0; i < vec_n; i++) {
        fprintf (fp, "%8i%8i\t%s\n", vec [i].x, vec [i].y, vec [i].name);
        fprintf (fp_ps, "%i %i (%s) n\n", vec [i].x, vec [i].y, to_ps (vec [i].name));
        if (i > 0 && vec [i].x == vec [i - 1].x && vec [i].y == vec [i - 1].y) {
            retval = 2;
            fprintf (fp_log, "%s == %s\n", vec [i - 1].name, vec [i].name);
            if (verbose)
                printf ("%s == %s\n", vec [i - 1].name, vec [i].name);
        }
    }
    fclose (fp);

    fprintf (fp_ps, "\n%i box\nend\nshowpage\n%c%cEOF\n", map_size, '%', '%');
    fclose (fp_ps);

    return retval;
}

void update_nn (int i)
{
    int
        j,
        cl1,
        cl2;

    cl1 = vec [i].cl1;
    cl2 = vec [i].cl2;
    vdiff [i][i] = 0.0;
    for (j = 0; j < i; j++)
        vdiff [i][j] = vdiff [j][i] = MIN (vdiff [cl1][j], vdiff [cl2][j]);

    for (j = 0; j < i; j++)
        if (vec [j].cluster == vec [i].cl1 ||
            vec [j].cluster == vec [i].cl2
        )
            vec [j].cluster = i;
}

void syntax ()
{
    fprintf (
	stderr,
        "\nKohonen Map Generator\n"
        "With minimal spanning tree and visible cluster boundaries\n"
        "(c) P. Kleiweg 1995 - 1999\n\n"
        "Usage: %s -i string [-o string]\n"
        "\t[-m int] [-s int] [-e int] [-q]\n\n"
        "-i: input file\n"
        "-o: name without extension of output files\n"
        "-m: map size\n"
        "-s: initial map size\n"
        "-e: number of epochs\n"
        "-q: quiet\n\n",
	programname
    );
    exit (1);
}

void setup ()
{
    int
	i,
	j,
	k;
    float
	f,
	max,
	min;

    srand ((unsigned int) time (NULL));

    max = -MAXFLOAT;
    min = MAXFLOAT;
    for (i = 0; i < vec_n; i++)
	for (j = 0; j < vec_size; j++) {
	    if (vec [i].v[j] < min)
		min = vec [i].v[j];
	    if (vec [i].v[j] > max)
		max = vec [i].v[j];
	}

    if (map_max < 1)
	map_max = sqrt (vec_n * 2) + 1;
    if (map_max < 2)
	map_max = 2;
    if (map_size < 2 || map_size > map_max)
        map_size = (vec_size == 1) ? map_max : 2;

    map = (float ***) m_malloc ((map_max + 1) * sizeof (float **));
    for (i = 1; i <= map_max; i++) {
	map [i] = (float **) m_malloc ((map_max + 1) * sizeof (float *));
	for (j = 1; j <= map_max; j++) {
	    map [i][j] = (float *) m_malloc (vec_size * sizeof (float));
	    for (k = 0; k < vec_size; k++) {
		f = ((float)rand ()) / ((float)RAND_MAX);    /*   0 ... 1 */
		f = f * 2.0 - 1.0;                           /*  -1 ... 1 */
		f = f * (max - min) / 10.0;
		f += (min + max) / 2.0;    /* 1/5 halverwege min en max   */
		map [i][j][k] = f;
	    }
	}
    }

    line = (LINE *) m_malloc (
	2 * map_max * (map_max - 1) * sizeof (LINE)
    );

    vdiff = (float **) m_malloc ((2 * vec_n - 1) * sizeof (float *));
    for (i = 0; i < 2 * vec_n - 1; i++)
        vdiff [i] = (float *) m_malloc ((2 * vec_n - 1) * sizeof (float));

    inc = (float **) m_malloc ((map_max + 1) * sizeof (float *));
    for (i = 1; i <= map_max; i++)
	inc [i] = (float *) m_malloc ((map_max + 1) * sizeof (float));

    cluster = (int *) m_malloc (vec_n * sizeof (int));
}

float square_vecdiff (float *v1, float *v2)
{
    int
	i;
    float
	err,
        sum;

    sum = 0.0;
    for (i = 0; i < vec_size; i++) {
	err = v1 [i] - v2 [i];
        sum += err * err;
    }
    return sum;
}

void process_args (int argc, char *argv [])
{
    int
	index;

    if (argc == 1)
        syntax ();

    for (index = 1; index < argc; index++) {
	if (argv [index][0] == '-') {
	    switch (argv [index][1]) {
		case 'q':
		    verbose = 0;
		    break;
		case 'e':
		    if (sscanf (
			    get_arg (&argc, argv, &index),
			    "%i",
			    &max_epoch) != 1)
			errit ("Wrong number of epochs");
		    break;
		case 'm':
		    if (sscanf (
			    get_arg (&argc, argv, &index),
			    "%i",
			    &map_max) != 1)
			errit ("Wrong mapsize");
		    break;
		case 's':
		    if (sscanf (
			    get_arg (&argc, argv, &index),
			    "%i",
			    &map_size) != 1)
			errit ("Wrong initial mapsize");
		    break;
		case 'i':
		    infile = m_strdup (get_arg (&argc, argv , &index));
		    break;
		case 'o':
		    baseout = m_strdup (get_arg (&argc, argv , &index));
#ifdef __MSDOS__
		    fnsplit (baseout, NULL, NULL, buf1, buf2);
		    if (! buf1 [0])
			errit ("Missing file name: \"%s\"", baseout);
		    if (buf2 [0])
			errit ("No file extension allowed: \"%s\"", baseout);
#endif
		    break;
		default:
                    errit ("Invalid option '%s'", argv [index]);
	    }
	} else
            errit ("Invalid parameter \"%s\"", argv [index]);
    }
}

char *get_arg (int *argc, char *argv [], int *index)
{
    int
	i;
    i = 1;
    while (argv [*index][++i]) {
	if (! isspace ((unsigned char) argv [*index][i]))
	    return argv [*index] + i;
    }
    if (*index == *argc - 1)
	errit ("Missing argument for \"%s\"", argv [*index]);
    return argv [++*index];
}

void *m_realloc (void *block, size_t size)
{
    void
	*p;

    p = realloc (block, size);
    if (! p) {
        free (no_mem_buffer);
	errit (Out_of_memory);
    }
    return p;
}

void *m_malloc (size_t size)
{
    void
	*p;

    p = malloc (size);
    if (! p) {
        free (no_mem_buffer);
	errit (Out_of_memory);
    }
    return p;
}

char *m_strdup (const char *s)
{
    char
	*p = strdup (s);
    if (! p) {
        free (no_mem_buffer);
	errit (Out_of_memory);
    }
    return p;
}

void errit (char const *format, ...)
{
    va_list
	list;

    fprintf (stderr, "\nError %s: ", programname);

    va_start (list, format);
    vfprintf (stderr, format, list);

    fprintf (stderr, "\n\n");

    exit (1);
}

int getline (FILE *fp, int required, char const *filename)
/* Lees een regel in
 * Plaats in buffer
 * Negeer lege regels en regels die beginnen met #
 */
{
    int
	i;

    for (;;) {
	if (fgets (buffer, BUFSIZE, fp) == NULL) {
	    if (required)
		errit ("Unexpected end of file in \"%s\"", filename);
	    else
		return 0;
	}
	input_line++;
	i = strlen (buffer);
	while (i && isspace ((unsigned char) buffer [i - 1]))
	    buffer [--i] = '\0';
	i = 0;
	while (buffer [i] && isspace ((unsigned char) buffer [i]))
	    i++;
	if (buffer [i] == '#')
	    continue;
	if (buffer [i]) {
	    memmove (buffer, buffer + i, strlen (buffer) + 1);
	    return 1;
	}
    }
}

VEC *store_vec (char *name)
{
    if (vec_n == vec_max) {
	vec_max += 256;
	vec = (VEC *) m_realloc (vec, vec_max * sizeof (VEC));
    }

    vec [vec_n].name = m_strdup (name);

    return &(vec [vec_n++]);
}

FILE *r_fopen (char const *filename)
{
    FILE
	*fp;

    fp = fopen (filename, "r");
    if (! fp)
	errit ("Bestand \"%s\" niet gevonden", filename);
    input_line = 0;
    return fp;
}

FILE *w_fopen (char const *filename)
{
    FILE
	*fp;

    fp = fopen (filename, "w");
    if (! fp)
	errit ("Bestand \"%s\" kan niet geopend worden", filename);
    return fp;
}

void get_programname (char const *argv0)
{
#ifdef __MSDOS__
    char
        name [MAXFILE];
    fnsplit (argv0, NULL, NULL, name, NULL);
    programname = strdup (name);
#else
    char
        *p;
    p = strrchr (argv0, my_PATH_SEP);
    if (p)
        programname = strdup (p + 1);
    else
        programname = strdup (argv0);
#endif    
}

void Heapsort
    (
	void *base,
	size_t nelem,
	size_t width,
	int (*fcmp)(void *, void *)
    )
{
    unsigned
	L, R, S, J;
    char
	*A,
	*heapswap;

    if (nelem < 2)
	return;

    A = (char *) base;
    heapswap = (char *) m_malloc (width);

    A -= width;                 /* table: A[1]..        */
    R = nelem;                  /*              ..A[R]  */
    L = (R >> 1) + 1;
    while (R > 1) {

	if (L > 1)
	    /* generate reverse partial ordered tree */
	    --L;
	else {
	    /* delete max from reverse partial ordered list
	       and append in front of ordered list    */
	    memcpy (heapswap, A + width, width);
	    memcpy (A + width, A + R * width, width);
	    memcpy (A + R-- * width, heapswap, width);
	}

	/* A[L+1]..A[R] is reverse partial ordered tree
	   insert A[L] to make A[L]..A[R] reverse partial ordered tree  */
	memcpy (heapswap, A + L * width, width);
	J = L;
	for (;;) {
	    S = J;              /* mother         */
	    J <<= 1;            /* left daughter  */
	    if (J > R)          /* no daughters   */
		break;
	    if ( (J < R) && (fcmp (A + J * width, A + (J + 1) * width) < 0) )
		J++;            /*  right daughter  */
	    if (fcmp (heapswap, A + J * width) >= 0)
		break;
	    memcpy (A + S * width, A + J * width, width);
	}
	memcpy (A + S * width, heapswap, width);
    }

    free (heapswap);
}

int vecpos_cmp (void *p1, void *p2)
{
    VEC
        *v1,
        *v2;

    v1 = (VEC *) p1;
    v2 = (VEC *) p2;

    if (v1->x != v2->x)
        return v1->x - v2->x;

    if (v1->y != v2->y)
        return v1->y - v2->y;

    return strcmp (v1->name, v2->name);
}

int line_cmp (void *p1, void *p2)
{
    float
        f;

    f = ((LINE *)p1)->f - ((LINE *)p2)->f;
    if (f < 0)
        return -1;
    if (f > 0)
        return 1;
    return 0;
}

void hussel ()
{
    char
	*swap;
    int
	i,
	j;

    swap = (char *) m_malloc (sizeof (VEC));
    for (i = vec_n - 1; i > 1; i--) {
	j = rand () % i;
	memcpy (swap, &(vec [i]), sizeof (VEC));
	memcpy (&(vec [i]), &(vec [j]), sizeof (VEC));
	memcpy (&(vec [j]), swap, sizeof (VEC));
    }
    free (swap);
}

char *to_ps (char const *string)
{
    int
        i, j;

    i = j = -1;
    while (string [++i])
        switch (string [i]) {
            case '(':
            case ')':
            case '\\':
                buffer [++j] = '\\';
                /*  no break  */
    	default:
                buffer [++j] = string [i];
	}
    buffer [++j] = '\0';

    return buffer;
}