/*************************/
/*                       */
/*      ko_work.c        */
/*                       */
/*      Alan Dix         */
/*                       */
/*      February 1993    */
/*                       */
/*************************/

#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "kohonen.h"

#ifdef ANSI
double prox_discrete(int x1,int y1,int x2,int y2,double diam);
double prox_manhatten(int x1,int y1,int x2,int y2,double diam);
double prox_mexican(int x1,int y1,int x2,int y2,double diam);
double (*prox)(int x1,int y1,int x2,int y2,double diam) = prox_manhatten;
static void init_serial(int len);
static int  next_serial(int len);
static void init_rand(int len);
static int  next_rand(int len);
void (*init_index)(int len) = init_serial;
int  (*next_index)(int len) = next_serial;
#else
double prox_discrete();
double prox_manhatten();
double prox_mexican();
double (*prox)() = prox_manhatten;

static void init_serial();
static int  next_serial();
static void init_rand();
static int  next_rand();
void (*init_index)() = init_serial;
int  (*next_index)() = next_serial;
#endif

double prox_discrete(x1,y1,x2,y2,diam)
    int x1, y1, x2, y2;
    double diam;
{
    (void)diam;  /* don't use diam */
    if ( (x1==x2) && (y1==y2) ) return 1.0;
    else return 0.0;
}


double prox_manhatten(x1,y1,x2,y2,diam)
    int x1, y1, x2, y2;
    double diam;
{
    if ( (x1==x2) && (y1==y2) ) return 1.0;
    if ( abs(x1-x2) > diam ) return 0.0;
    if ( abs(y1-y2) > diam ) return 0.0;
    else return 1.0;
}

double prox_mexican(x1,y1,x2,y2,diam)
    int x1, y1, x2, y2;
    double diam;
{
    int dist = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2);
    return 1.0 / ( 1.0 + dist / (diam * diam) );
}

int set_prox(name)
    char *name;
{
    if ( strcmp(name,"discrete") == 0 ) {
        prox = prox_discrete;
        return 1;
        }
    if ( strcmp(name,"manhatten") == 0 ) {
        prox = prox_manhatten;
        return 1;
        }
    if ( strcmp(name,"mexican") == 0 ) {
        prox = prox_mexican;
        return 1;
        }
    return 0;
}

static serial_ct;
        
static void init_serial(len)
    int len;
{
    (void)len;  /* don't use len */
    serial_ct = 0;
}

static int next_serial(len)
    int len;
{
    serial_ct++;
    if ( serial_ct >= len ) serial_ct = 0;
    return serial_ct;
}

/* static int max_rand; */

static void init_rand(len)
    int len;
{
    (void)len;  /* do nout */
}

static int next_rand(len)
    int len;
{
    return rand() % len;
}

int set_index(name)
    char *name;
{
    if ( strcmp(name,"serial") == 0 ) {
        init_index = init_serial;
        next_index = next_serial;
        return 1;
        }
    if ( strcmp(name,"rand") == 0 ) {
        init_index = init_rand;
        next_index = next_rand;
        return 1;
        }
    return 0;
}


double  find_closest_data(vec,close_j,data)
    Ds vec;
    int *close_j;
    data_p data;
{
    int j, i;
    int cj=0;
    double cdist = 1.0e90;
    for ( j=0; j<data->wid; j++ ) {
        double d = 0.0;
        for ( i=0; i<data->wid; i++ ) {
            double d1 = vec[i] - data->val[j][i];
            d += d1 * d1;
            }
        if ( d < cdist ) {
            cdist = d;
            cj = j;
            }
        }
    *close_j = cj;
    return cdist;
}


double  find_closest_code(vec,close_x,close_y,code)
    Ds vec;
    int *close_x;
    int *close_y;
    code_p code;
{
    int x, y, i;
    int cx=0, cy=0;
    double cdist = 1.0e90;
    for ( x=0; x<code->X; x++ ) {
        for ( y=0; y<code->Y; y++ ) {
            double d = 0.0;
            for ( i=0; i<code->wid; i++ ) {
                double d1 = vec[i] - code->val[x][y][i];
                d += d1 * d1;
                }
            if ( d < cdist ) {
                cdist = d;
                cx = x; cy = y;
                }
            }
        }
    *close_x = cx;
    *close_y = cy;
    return cdist;
}


int update(vec,alpha,diam,code)
    Ds vec;
    double alpha;
    double diam;
    code_p code;
{
    int close_x, close_y, x, y, i;
    (void)find_closest_code(vec,&close_x,&close_y,code);
    for ( x=0; x<code->X; x++ ) {
        for ( y=0; y<code->Y; y++ ) {
            double p = prox(close_x,close_y,x,y,diam);
            if ( p < 0.0001 ) continue;
            for ( i=0; i<code->wid; i++ ) {
                code->val[x][y][i] +=
                              alpha * p * ( vec[i] - code->val[x][y][i] );
                }
            }
        }
    return 1;
}

int train_code(rlen,alpha0,diam0,K,data,code)
    long rlen;
    double alpha0;
    double diam0;
    double K;
    data_p data;
    code_p code;
{
    long to_do = rlen;
    int i;
    long t = 0;
    double diam  = diam0;
    double alpha = alpha0;
    double lambda = exp(-K);
    if ( data->wid != code->wid ) return -1;
    if ( rlen <= 0 ) return 0;
    for ( i=0; i<data->len && to_do>=0; i++, to_do-- ) {
        update( data->val[i], alpha, diam, code );
        }
    init_index(data->len);
    for ( t=0; t<to_do; t++, i++ ) {
        int index = next_index(data->len);
        if ( i >= data->len ) {
            i = 0;
            diam  = diam * lambda;
            alpha = alpha * lambda;
            }
        update( data->val[index], alpha, diam, code );
        }
    return 1;
}

int find_label(lab, labs)
    char *lab;
    label_p labs;
{
    int i;
    if ( lab == NULL ) return 0;
    for ( i=0; i<labs->nos; i++ ) {
        if ( strcmp(lab,labs->labels[i]) == 0 ) return i;
        }
    return -1;
}


label_p gen_label_list(data)
    data_p data;
{
    label_p labs = (label_p)malloc(sizeof(label_t));
    int  j;
    labs->nos = 1;
    labs->labels = (Css)malloc(labs->nos*sizeof(Cs));
    labs->labels[0] = strdup("");
    for ( j=0; j<data->len; j++ ) {
        if ( find_label(data->labels[j], labs) >= 0 ) continue;
        labs->nos++;
        labs->labels = (Css)realloc(labs->labels,labs->nos*sizeof(Cs));
        labs->labels[labs->nos-1] = strdup(data->labels[j]);
        }
    return labs;
}

code_p label_code_stats(data,labs,code)
    data_p  data;
    label_p labs;
    code_p  code;
{
    int x, y,j,i;
    code_p stats = alloc_code(labs->nos, code->X, code->Y);

    if ( data->wid != code->wid ) return NULL;
    
    for ( x=0; x<stats->X; x++ ) {
        for ( y=0; y<stats->Y; y++ ) {
            for ( i=0; i<stats->wid; i++ ) {
                stats->val[x][y][i] = 0.0;
                }
            }
        }
    for ( j=0; j<data->len; j++ ) {
        int lab = find_label(data->labels[j], labs);
        if ( lab < 0 ) return NULL;
        (void)find_closest_code(data->val[j],&x,&y,code);
        stats->val[x][y][lab] += 1.0;
        }
    return stats;
}

int label_code_prop(stats,labs,code,hi,lo)
    code_p  stats;
    label_p labs;
    code_p  code;
    double  hi;
    double  lo;
{
    int x, y,i;

    for ( x=0; x<stats->X; x++ ) {
        for ( y=0; y<stats->Y; y++ ) {
            Ds vec = stats->val[x][y];
            double sum = 0.0;
            double fst = 0.0, snd = 0.0;
            int lab = 0;

            for ( i=0; i<stats->wid; i++ ) {
                sum += vec[i];
                }
            for ( i=0; i<stats->wid; i++ ) {
                vec[i] = vec[i] / sum;
                if ( vec[i] >= fst ) {
                    snd = fst;
                    fst = vec[i];
                    lab = i;
                    }
                else if ( vec[i] >= snd ) {
                    snd = vec[i];
                    }
                } /* for i */
            if ( fst >= hi && snd < lo ) {
                stats->labels[x][y] = strdup(labs->labels[lab]);
                code->labels[x][y]  = strdup(labs->labels[lab]);
                }
            else {
                stats->labels[x][y] = strdup(labs->labels[0]);
                code->labels[x][y]  = strdup(labs->labels[0]);
                }
            } /* for Y */
        } /* for X */
    return 1;
}

int label_code(stats,labs,code,hi,lo)
    code_p  stats;
    label_p labs;
    code_p  code;
    double  hi;
    double  lo;
{
    int x, y,i;

    for ( x=0; x<stats->X; x++ ) {
        for ( y=0; y<stats->Y; y++ ) {
            Ds vec = stats->val[x][y];
            double fst = 0.0, snd = 0.0, sum = 0.0;
            int lab = 0;

            for ( i=1; i<stats->wid; i++ ) {
                sum += vec[i];
                if ( vec[i] >= fst ) {
                    snd = fst;
                    fst = vec[i];
                    lab = i;
                    }
                else if ( vec[i] >= snd ) {
                    snd = vec[i];
                    }
                } /* for i */
            if ( fst >= hi*sum && snd < lo*sum ) {
                stats->labels[x][y] = strdup(labs->labels[lab]);
                code->labels[x][y]  = strdup(labs->labels[lab]);
                }
            else {
                stats->labels[x][y] = strdup(labs->labels[0]);
                code->labels[x][y]  = strdup(labs->labels[0]);
                }
            } /* for Y */
        } /* for X */
    return 1;
}