/*
 * File: formant6.c
 *
 * (c) Peter Kleiweg
 *     Wed Feb  7 16:59:48 2007
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2,
 * or (at your option) any later version.
 *
 */

#define my_VERSION "0.01"

#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 <dir.h>
#endif  /* __MSDOS__  */
#include <ctype.h>
#include <errno.h>
#include <float.h>
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#define BUFSIZE 4095

typedef enum { FALSE = 0, TRUE } BOOL_;

typedef struct {
    double
        f [6];
} F6_;

typedef struct {
    int
        n_strings,
	max_strings;
    F6_
        *f6;
} ITEM_;

typedef struct {
    char
        *s;
    int
        i;
} LABEL_;

typedef struct {
    double
        f;
    int
        n;
} DIFF_;

ITEM_
    *items = NULL;

LABEL_
    *labels = NULL;

DIFF_
    **diff;

int
    rowlen = 0,
    *row1 = NULL,
    *row2 = NULL,
    lineno,
    n_locations = 0,
    arg_c;

FILE
    *fp_in = NULL,
    *fp_out = NULL;

char
    qq [] = "|/-\\",
    *label_file = NULL,
    *output_file = "formant6.out",
    buffer [BUFSIZE + 1],
    **arg_v,
    *programname,
    *no_mem_buffer,
    out_of_memory [] = "Out of memory";

char const
    *filename;

BOOL_
    verbose = TRUE;

void
    open_read (char const *s),
    read_labelfile (void),
    get_programname (char const *argv0),
    process_args (void),
    errit (char const *format, ...),
    syntax (void),
    *s_malloc (size_t size),
    *s_realloc (void *block, size_t size);
char
    *getline (BOOL_ required),
    *get_arg (void),
    *s_strdup (char const *s);
char const
    *unquote (char *s);
int
    kgv (int l1, int l2),
    cmpi (void const *, void const *),
    cmps (void const *, void const *),
    searchcmp (void const *key, void const *p),
    getlocation (void);
double
    diffun (F6_ *, F6_ *);

int main (int argc, char *argv [])
{
    int
	argn,
	loc,
	i,
	i1,
	i2,
	j,
	l,
	l1,
	l2,
	n,
	q;
    double
	newsum,
	sum;
    F6_
	f6,
	*f6pi,
	*f6pj;
    BOOL_
	busy;


    no_mem_buffer = (char *) malloc (1024);

    get_programname (argv [0]);

    arg_c = argc;
    arg_v = argv;
    process_args ();

    if (arg_c == 1)
	syntax ();

    if (access (output_file, F_OK) == 0)
        errit ("File exists: \"%s\"", output_file);

    fp_out = fopen (output_file, "w");
    if (! fp_out)
	errit ("Creating file \"%s\": %s", output_file, strerror (errno));

    if (! label_file)
	errit ("Missing label file");

    if (n_locations < 2)
        errit ("Illegal number of locations: %i", n_locations);

    diff = (DIFF_ **) s_malloc ((n_locations + 1) * sizeof (DIFF_ *));
    for (i = 2; i <= n_locations; i++) {
	diff [i] = (DIFF_ *) s_malloc (i * sizeof (DIFF_));
	for (j = 1; j < i; j++) {
	    diff [i][j].f = 0.0;
	    diff [i][j].n = 0;
	}
    }

    read_labelfile ();

    items = (ITEM_ *) s_malloc ((n_locations + 1) * sizeof (ITEM_));
    for (i = 1; i <= n_locations; i++) {
	items [i].max_strings = 0;
	items [i].f6 = NULL;
    }

    q = 0;

    for (argn = 1; argn < arg_c; argn++) {

	open_read (arg_v [argn]);

	if (verbose) {
	    fprintf (stderr, "  %s\r", filename);
	    fflush (stderr);
	}

	for (i = 1; i <= n_locations; i++)
	    items [i].n_strings = 0;


	loc = -1;
	while (getline (FALSE)) {
	    if (buffer [0] == ':')
		loc = getlocation ();
	    else {
		if (loc < 0)
		    errit ("Vlues before label, in file %s, line %li", arg_v [argn], lineno);
		if (sscanf (buffer, "%lf %lf %lf %lf %lf %lf",
			    &(f6.f[0]), &(f6.f[1]), &(f6.f[2]), &(f6.f[3]), &(f6.f[4]), &(f6.f[5])) != 6)
		    errit ("Missing value(s), in file %s, line %li", arg_v [argn], lineno);
		if (items [loc].n_strings == items [loc].max_strings) {
		    items [loc].max_strings += 16;
		    items [loc].f6 = (F6_ *) s_realloc (items [loc].f6, items [loc].max_strings * sizeof (F6_));
		}
		items [loc].f6 [items [loc].n_strings++] = f6;
	    }
	}
	fclose (fp_in);


	for (i = 2; i <= n_locations; i++) {
	    if ((l1 = items [i].n_strings) == 0)
		continue;
	    f6pi = items [i].f6;
	    for (j = 1; j < i; j++) {

		if (verbose) {
		    if (++q == 4)
			q = 0;
		    fprintf (stderr, "%c\r", qq [q]);
		    fflush (stderr);
		}

		if ((l2 = items [j].n_strings) == 0)
		    continue;
		f6pj = items [j].f6;

		l = kgv (l1, l2);
		if (l > rowlen) {
		    rowlen = l;
		    row1 = (int *) s_realloc (row1, rowlen * sizeof (int));
		    row2 = (int *) s_realloc (row2, rowlen * sizeof (int));
		}
		for (n = 0; n < l; n++) {
		    row1 [n] = n % l1;
		    row2 [n] = n % l2;
		}

		sum = 0.0;
		for (n = 0; n < l; n++)
		    sum += diffun (f6pi + row1 [n], f6pj + row2 [n]);

		busy = TRUE;
		while (busy) {

		    busy = FALSE;
		    for (i1 = 0; i1 < l; i1++)
			for (i2 = i1 + 1; i2 < l; i2++) {
			    if (diffun (f6pi + row1 [i1], f6pj + row2 [i1]) + diffun (f6pi + row1 [i2], f6pj + row2 [i2]) >
				diffun (f6pi + row1 [i1], f6pj + row2 [i2]) + diffun (f6pi + row1 [i2], f6pj + row2 [i1]))
				{
				    n = row2 [i1];
				    row2 [i1] = row2 [i2];
				    row2 [i2] = n;
				    busy = TRUE;
				}
			}
		    if (busy) {
			newsum = 0.0;
			for (n = 0; n < l; n++)
			    newsum += diffun (f6pi + row1 [n], f6pj + row2 [n]);
			if (sum > newsum)
			    sum = newsum;
			else
			    busy = FALSE;
		    }
		}

		diff [i][j].f += sum / l;
		diff [i][j].n++;
	    }
	}

	if (verbose) {
	    fprintf (stderr, " \n");
	    fflush (stderr);
	}
    }


    if (label_file)
        qsort (labels, n_locations, sizeof (LABEL_), cmpi);

    fprintf (fp_out, "%i\n", n_locations);
    if (label_file) {
	for (i = 0; i < n_locations; i++)
	    fprintf (fp_out, "%s\n", labels [i].s);
    } else {
	for (i = 1; i <= n_locations; i++)
	    fprintf (fp_out, "%i\n", i);
    }

    for (i = 2; i <= n_locations; i++)
        for (j = 1; j < i; j++)
	    if (diff [i][j].n)
		fprintf (fp_out, "%g\n", diff [i][j].f / diff [i][j].n);
	    else {
		fprintf (fp_out, "NA\n");
		if (verbose)
		    fprintf (stderr, "No value for: %i \\ %i\n", j, i);
	    }

    fclose (fp_out);

    return 0;
}

double diffun (F6_ *fi, F6_ *fj)
{
    int
	i;
    double
	d,
	sum;

    sum = 0.0;
    for (i = 0; i < 6; i++) {
	d = fi->f [i] - fj->f [i];
	sum += d * d;
    }
    return sqrt (sum);
}


/*
 * kgv = kleinst gemeemschappelijk veelvoud
 */
int kgv (int l1, int l2)
{
    int
        ll1,
        ll2;

    ll1 = l1;
    ll2 = l2;
    for (;;) {
        if (ll1 == ll2)
            return ll1;
        while (ll1 < ll2)
            ll1 += l1;
        while (ll2 < ll1)
            ll2 += l2;
    }
}

void read_labelfile ()
{
    int
        i,
        n;

    open_read (label_file);
    labels = (LABEL_ *) s_malloc (n_locations * sizeof (LABEL_));
    for (i = 0; i < n_locations; i++) {
        labels [i].s = NULL;
        labels [i].i = i + 1;
    }
    while (getline (FALSE)) {
        if (sscanf (buffer, "%i%n", &i, &n) < 1)
            errit ("Missing number in \"%s\", line %lu", filename, lineno);
        if (i < 1 || i > n_locations)
            errit ("Number out of range in \"%s\", line %lu", filename, lineno);
        if (labels [i - 1].s)
            errit ("Duplicate label for nr %i in \"%s\", line %lu", i, filename, lineno);
        while (buffer [n] && isspace ((unsigned char) buffer [n]))
            n++;
        if (! buffer [n])
            errit ("Missing label in \"%s\", line %lu", filename, lineno);
        labels [i - 1].s = s_strdup (unquote (buffer + n));
    }
    fclose (fp_in);
    for (i = 0; i < n_locations; i++)
        if (! labels [i].s)
            errit ("Missing number and label for location nr %i in file %s", i + 1, filename);
    qsort (labels, n_locations, sizeof (LABEL_), cmps);
    for (i = 1; i < n_locations; i++)
        if (! strcmp (labels [i - 1].s, labels [i].s))
            errit ("Duplicate label \"%s\" in file %s", labels [i].s, filename);
}

int cmpi (void const *p1, void const *p2)
{
    return ((LABEL_ *)p1)->i - ((LABEL_ *)p2)->i;
}

int cmps (void const *p1, void const *p2)
{
    return strcmp (((LABEL_ *)p1)->s, ((LABEL_ *)p2)->s);
}

int searchcmp (void const *key, void const *p)
{
    return strcmp ((char *) key, ((LABEL_ *)p)->s);
}

/*
 * remove quotes from string, overwrite with result
 */
char const *unquote (char *s)
{
    int
        i,
        j;

    if (s [0] != '"')
        return s;
    j = 0;
    for (i = 1; s [i]; i++) {
        if (s [i] == '"')
            break;
        if (s [i] == '\\') {
            if (s [i + 1])
                s [j++] = s [++i];
        } else
            s [j++] = s [i];
    }
    s [j] = '\0';
    return s;
}

int getlocation ()
{
    char const
        *s;
    LABEL_
        *p;
    int
        i;

    i = 1;
    while (buffer [i] && isspace ((unsigned char) buffer [i]))
        i++;

    s = unquote (buffer + i);

    p = (LABEL_ *) bsearch (s, labels, n_locations, sizeof (LABEL_), searchcmp);
    if (! p)
	errit ("Undefined label \"%s\" in %s, line %li", s, filename, lineno);

    return p->i;
}

/*
 * get line from `fp_in' into `buffer'
 * remove leading and trailing white space
 * skip empty lines and lines starting with #
 */
char *getline (BOOL_ required)
{
    int
        i;

    for (;;) {
        lineno++;
        if (fgets (buffer, BUFSIZE, fp_in) == NULL) {
            if (required)
                errit ("Reading file \"%s\", line %li: %s",
                       filename, lineno, errno ? strerror (errno) : "End of file");
            else
                return NULL;
        }
        i = strlen (buffer);
        while (i > 0 && isspace ((unsigned char) buffer [i - 1]))
            buffer [--i] = '\0';
        i = 0;
        while (buffer [i] && isspace ((unsigned char) buffer [i]))
            i++;
        if (i)
            memmove (buffer, buffer + i, strlen (buffer + i) + 1);
        if (buffer [0] == '\0' || buffer [0] == '#')
            continue;
        return buffer;
    }
}

void open_read (char const *s)
{
    filename = s;
    lineno = 0;
    fp_in = fopen (s, "r");
    if (! fp_in)
        errit ("Opening file \"%s\": %s", s, strerror (errno));
}

void process_args ()
{
    while (arg_c > 1 && arg_v [1][0] == '-') {
        switch (arg_v [1][1]) {
            case 'l':
                label_file = get_arg ();
                break;
            case 'n':
                n_locations = atoi (get_arg ());
                break;
            case 'o':
                output_file = get_arg ();
                break;
            case 'q':
                verbose = FALSE;
                break;
            default:
                errit ("Illegal option '%s'", arg_v [1]);
        }
	arg_c--;
	arg_v++;
    }
}

char *get_arg ()
{
    if (arg_v [1][2])
        return arg_v [1] + 2;

    if (arg_c == 2)
        errit ("Missing argument for '%s'", arg_v [1]);

    arg_v++;
    arg_c--;
    return arg_v [1];
}

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);
}

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 *s_malloc (size_t size)
{
    void
	*p;

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

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

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

char *s_strdup (char const *s)
{
    char
	*s1;

    if (s) {
	s1 = (char *) s_malloc (strlen (s) + 1);
	strcpy (s1, s);
    } else {
	s1 = (char *) s_malloc (1);
	s1 [0] = '\0';
    }
    return s1;
}

void syntax ()
{
    fprintf (
	stderr,
	"\n"
	"Version " my_VERSION "\n"
	"\n"
	"Usage: %s -n number -l filename [-o filename] datafile(s)\n"
	"\n"
	"\t-l : file with location labels\n"
	"\t-n : number of locations\n"
	"\t-o : output file\n"
	"\n",
	programname
    );
    exit (1);
}