/* FFNet.c
 *
 * Copyright (C) 1997-2002 David Weenink
 *
 * This program 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 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 djmw 20020712 GPL header
*/

#include "FFNet_Matrix.h"
#include "Matrix_extensions.h"
#include "TableOfReal_extensions.h"
#include "Pattern.h"
#include "Collection.h"
#include "Categories.h"

static int bookkeeping (FFNet me);

#include "oo_DESTROY.h"
#include "FFNet_def.h"
#include "oo_COPY.h"
#include "FFNet_def.h"
#include "oo_EQUAL.h"
#include "FFNet_def.h"
#include "oo_WRITE_ASCII.h"
#include "FFNet_def.h"
#include "oo_WRITE_BINARY.h"
#include "FFNet_def.h"
#include "oo_READ_ASCII.h"
#include "FFNet_def.h"
#include "oo_READ_BINARY.h"
#include "FFNet_def.h"
#include "oo_DESCRIPTION.h"
#include "FFNet_def.h"

/****** non-linearities ****************************************************/

static double sigmoid (I, double x, double *deriv)
{
    double act = NUMsigmoid (x);
	(void) void_me;
    *deriv = act * ( 1.0 - act);
    return act;
}


/* ******************* cost functions ****************************************/

static double minimumSquaredError (I, const float target[])
{
    iam (FFNet); long i, k = my nNodes-my nOutputs+1;
    double cost = 0.0;
    
    for (i=1; i <= my nOutputs; i++, k++)
    {
		double e = my error[k] = target[i] - my activity[k];
		cost += e * e;
    }
    return 0.5 * cost;
}
/* E = - sum (i=1; i=nPatterns; sum (k=1;k=nOutputs; t[k]*ln (o[k]) + (1-t[k])ln (1-o[k])))*/
/* dE/do[k] = (1-t[k])/ (1-o[k]) - t[k]/o[k] */
/* werkt niet bij (grote?) netten */
static double minimumCrossEntropy (I, const float target[])
{
    iam (FFNet); long i, k = my nNodes-my nOutputs+1;
    double cost = 0.0;
    
    for (i=1; i <= my nOutputs; i++, k++)
    {
    	double t1 = 1.0 - target[i], a1 = 1.0 - my activity[k];
   		cost -= target[i] * log (my activity[k]) + t1 * log (a1);
		my error[k] = t1 / a1 - target[i] / my activity[k];
    }
    return cost;
}


/* *********************************************************************** */

static int bookkeeping (FFNet me)
{
    long j, i, n = my nUnitsInLayer[0] + 2, nWeights = 0;
    long firstNodeInPrevious = 1, lastWeightInPrevious = 0;
    
    my nNodes = my nUnitsInLayer[0];
    for (i=1; i <= my nLayers; i++)
    {
		my nNodes += my nUnitsInLayer[i] + 1;
        nWeights += my nUnitsInLayer[i] * (my nUnitsInLayer[i-1] + 1);
    }
    if (my nWeights > 0 && my nWeights != nWeights) return 0;
    my nWeights = nWeights;
    /* FFNet_readAscii already fills my w */ 
	if ( (! my w && ! (my w = NUMdvector (1, my nWeights))) ||
        ! (my activity     = NUMdvector (1, my nNodes)) ||
        ! (my isbias       = NUMlvector (1, my nNodes)) ||
        ! (my nodeFirst    = NUMlvector (1, my nNodes)) ||
        ! (my nodeLast     = NUMlvector (1, my nNodes)) ||
        ! (my wFirst       = NUMlvector (1, my nNodes)) ||
        ! (my wLast        = NUMlvector (1, my nNodes)) ||
        ! (my wSelected    = NUMlvector (1, my nWeights)) ||
        ! (my error        = NUMdvector (1, my nNodes)) ||
        ! (my deriv  	   = NUMdvector (1, my nNodes)) ||
        ! (my dwi          = NUMdvector (1, my nWeights)) ||
        ! (my dw           = NUMdvector (1, my nWeights))) return 0;
	my nInputs = my nUnitsInLayer[0];
	my nOutputs = my nUnitsInLayer[my nLayers];
    my isbias[my nInputs + 1] = my activity[my nInputs + 1] = 1;
    for (j=1; j <= my nLayers; j++,  n++)
    {
        for (i=1; i <= my nUnitsInLayer[j]; i++, n++)
        {
            my isbias[n] = 0;
            my nodeFirst[n] = firstNodeInPrevious;
            my nodeLast[n] = my nodeFirst[n] + my nUnitsInLayer[j-1];
            my wFirst[n] = lastWeightInPrevious + (i-1)* (my nUnitsInLayer[j-1] + 1) + 1;
            my wLast[n] = my wFirst[n] + my nUnitsInLayer[j-1];
        }
        if (j != my nLayers) my isbias[n] = my activity[n] = 1;
        lastWeightInPrevious = my wLast[n-1];
        firstNodeInPrevious += my nUnitsInLayer[j-1] + 1;
     }
     FFNet_selectAllWeights (me);
     return 1;
}

static void info (I)
{
    iam (FFNet); long i;
    Melder_info ("I have %ld layer%s.\nA total of %ld units divided as:",
    	my nLayers, my nLayers> 1 ? " (s)": "", FFNet_numberOfUnits (me));
    Melder_info ("\t%ld in layer %ld (output%s)", my nUnitsInLayer[my nLayers],
    	my nLayers, my outputsAreLinear ? ", linear" : "");
    for (i=my nLayers-1; i >= 1; i--)
    	Melder_info ("\t%ld in layer %ld (hidden)", my nUnitsInLayer[i], i);
    Melder_info ("\t%ld in layer 0 (input)", my nUnitsInLayer[0]);
    Melder_info ("This net has %ld weights of which %ld are selected.",
    	my nWeights, FFNet_dimensionOfSearchSpace (me));
    Melder_info ("\n\nnode  isbias nodeFirst nodeLast wFirst wLast");
	for (i=1; i <= my nNodes; i++) Melder_info ("%5ld: %7ld %10ld %9ld %7ld %6ld", i, 
        my isbias[i], my nodeFirst[i], my nodeLast[i], my wFirst[i], my wLast[i]);
} 


class_methods (FFNet, Data)
	class_method_local (FFNet, destroy)
	class_method_local (FFNet, copy)
	class_method_local (FFNet, equal)
	class_method_local (FFNet, writeAscii)
	class_method_local (FFNet, writeBinary)
	class_method_local (FFNet, readAscii)
	class_method_local (FFNet, readBinary)
	class_method_local (FFNet, description)
    class_method (info)
class_methods_end


int FFNet_init (FFNet me, long nLayers, const long nUnitsInLayer[], Any outputCategories,
	int outputsAreLinear, int nonLinearityType)
{
    my nLayers = nLayers;
    if (! (my nUnitsInLayer = NUMlvector_copy (nUnitsInLayer, 0, my nLayers)) ||
    	! bookkeeping (me)) return 0;
	if (outputCategories && ! (my outputCategories = Data_copy (outputCategories))) return 0;
    my outputsAreLinear = outputsAreLinear;
    FFNet_reset (me, 0.1);
    FFNet_setCostFunction (me, FFNet_COST_MSE);
    my nonLinearityType = nonLinearityType;
    FFNet_setNonLinearity (me, nonLinearityType);
    return 1;
}

Any FFNet_create (long nLayers, const long nUnitsInLayer[], Any outputCategories,
	int outputsAreLinear, int nonLinearityType)
{
    FFNet me = new (FFNet);
    if (! me || ! FFNet_init (me, nLayers, nUnitsInLayer, outputCategories,
    	outputsAreLinear, nonLinearityType)) forget (me);
    return me;
}

void FFNet_setNonLinearity (FFNet me, int nonLinearityType)
{
	my nonLinearityType = nonLinearityType;
	my nonLinearity = sigmoid;
	my nlClosure = NULL;
}

void FFNet_setCostFunction (FFNet me, int costType)
{
	my costFunctionType = costType;
	if (costType == 2) my costFunction = minimumCrossEntropy;
	else my costFunction = minimumSquaredError;
	my cfClosure = NULL;
}
       
void FFNet_reset (FFNet me, double wrange)
{
    long i;
    for (i=1; i <= my nWeights; i++)
        if (my wSelected[i]) my w[i] = NUMrandomUniform (-wrange, wrange);
    for (i=1; i <= my nNodes; i++)
        my activity[i] = (my isbias[i] ? 1.0 : 0.0);
    my accumulatedCost = 0.0;
    forget (my minimizer);  
}


/***** OPERATION: ***********************************************************/
/* step 1 */
void FFNet_propagate (FFNet me, const float input[], float output[])
{
    double act; long i, j, k = 1;
    long nNodes = my outputsAreLinear ? my nNodes - my nOutputs : my nNodes;

	/* clamp input pattern on the network */
	for (i=1; i <= my nUnitsInLayer[0]; i++) my activity[i] = input[i];
    	
    /* on hidden units use activation function */
    for (i=my nUnitsInLayer[0]+2; i <= nNodes; i++)
    {
        if (my isbias[i]) continue;
        for (act=0.0, j=my nodeFirst[i]; j <= my nodeLast[i]; j++, k++)
            act += my w[k] * my activity[j];
        my activity[i] = my nonLinearity (me, act, & my deriv[i]);
    }
    /* on output units use another activation function */
    if (my outputsAreLinear)
    {
    	for (i=nNodes+1; i <= my nNodes; i++)  
    	{
        	if (my isbias[i]) continue;
        	for (act=0.0, j=my nodeFirst[i]; j <= my nodeLast[i]; j++, k++)
            	act += my w[k] * my activity[j];
			my activity[i] = act;
			my deriv[i] = 1.0;
    	}
    }
    if (output) for (k=my nNodes-my nOutputs+1, i=1; i <= my nOutputs; i++, k++)
    	output[i] = my activity[k];  
}


double FFNet_computeError (FFNet me, const float target[])
{
	long i, j, k; double cost;
	
	/* compute error at output layer */
	cost = my costFunction (me, target);
    for (i=1; i <= my nNodes - my nOutputs; i++) my error[i] = 0.0;
    /* backpropagation of errors from output to first hidden layer*/
    for (i=my nNodes; i > my nInputs+1; i--)
    {
        if (my isbias[i]) continue;
        my error[i] *= my deriv[i];
        if (my nodeFirst[i] > my nInputs + 1)
        {
            for (k=my wFirst[i], j=my nodeFirst[i]; j <= my nodeLast[i]-1; j++, k++)
	        	my error[j] += my error[i] * my w[k];
        }
    }
    return cost;
}

void FFNet_computeDerivative (FFNet me)
{
    long i, j, k = 1;
    for (i=my nInputs+2; i <= my nNodes; i++)
    {
        if ( ! my isbias[i]) for (j=my nodeFirst[i]; j <= my nodeLast[i]; j++, k++)
            my dwi[k] = - my error[i] * my activity[j];
    }
}

/******* end operation ******************************************************/

long FFNet_getWinningUnit (FFNet me, int labeling)
{
    long i, pos = 1, k = my nNodes - my nOutputs;
	if (labeling == 2) /* stochastic */
	{
		double sum = 0, random;
		for (i=1; i <= my nOutputs; i++) sum += my activity[k+i];
		random = NUMrandomUniform (0, sum);
		for (pos=my nOutputs; pos >= 2; pos--) if (random > (sum -= my activity[k+pos])) break;
	}
	else /* winner-takes-all */
	{
		double max = my activity[k+1];
		for (i=2; i <= my nOutputs; i++) if (my activity[k+i] > max)
		{ max = my activity[k+i]; pos = i; }
	}
	return pos;
}

void FFNet_propagateToLayer (FFNet me, const float input[], float activity[], long layer)
{
	long i, k = 0;
	Melder_assert (activity);
	FFNet_propagate (me, input, NULL);
	for (i=0; i < layer; i++) k += my nUnitsInLayer[i] + 1;
	for (i=1; i <= my nUnitsInLayer[layer]; i++) activity[i] = my activity[k+i];
}

void FFNet_selectAllWeights (FFNet me)
{
	long i;
	for (i=1; i <= my nWeights; i++) my wSelected[i] = 1;
	my dimension = my nWeights;
}

long FFNet_dimensionOfSearchSpace (FFNet me)
{
	long i, n = 0;
	for (i=1; i <= my nWeights; i++) if (my wSelected[i]) n++;
	return n;
}

void FFNet_selectBiasesInLayer (FFNet me, long layer)
{
	long i, node = my nUnitsInLayer[0]+1;
	if (layer < 1 || layer > my nLayers) return;
	for (i=1; i <= my nWeights; i++) my wSelected[i] = 0;
	for (i=1; i < layer; i++) node += my nUnitsInLayer[i] + 1;
	for (i=node+1; i <= node+my nUnitsInLayer[layer]; i++)
		my wSelected[my wLast[i]] = 1;
	my dimension = my nUnitsInLayer[layer];
}

void FFNet_weightConnectsUnits (FFNet me, long index, long *fromUnit, long *toUnit, 
    long *layer)
{
    long i = 1, np, nw = my nUnitsInLayer[1] * (my nInputs + 1);

    Melder_assert (index > 0 && index <= my nWeights);
    while (index > nw)
    { i++; nw += ( np = my nUnitsInLayer[i] * (my nUnitsInLayer[i-1] + 1)); }
    if (i > 1) index -= nw - np;
    *fromUnit = index % (my nUnitsInLayer[i-1] + 1);
    *toUnit = (index - 1) / (my nUnitsInLayer[i-1] + 1) + 1;
    *layer = i;
}

void FFNet_nodeToUnitInLayer (FFNet me, long node, long *unit, long *layer)
{
    long i = 0, nn = my nUnitsInLayer[0] + 1;
    Melder_assert (node > 0 && node <= my nNodes);
    while (node > nn) { nn += my nUnitsInLayer[++i] + 1; }
    if (i > 0) node -= nn - (my nUnitsInLayer[i] + 1);
    *unit = node % (my nUnitsInLayer[i] + 1);
    *layer = i;
}

long FFNet_numberOfWeights (FFNet me) { return my nWeights; }

long FFNet_numberOfLayers (FFNet me) {  return my nLayers; }

long FFNet_numberOfUnits (FFNet me) { return my nNodes - my nLayers; }

long FFNet_numberOfHiddenLayers (FFNet me) { return my nLayers - 1; }

long FFNet_numberOfUnitsInLayer (FFNet me, int layer)
{
    if (layer > my nLayers || layer < 0) return 0;
    return my nUnitsInLayer[layer];  
}

double FFNet_getMinimum (FFNet me)
{
	return my minimizer ? Minimizer_getMinimum (my minimizer) : NUMundefined;
}

void FFNet_drawTopology (FFNet me, Graphics g)
{
    long i, maxNumOfUnits = my nUnitsInLayer[0];
    int dxIsFixed = 1;
    double r, dx, dy = 1.0 / (my nLayers + 1);

    for (i=1; i <= my nLayers; i++)
	if (my nUnitsInLayer[i] > maxNumOfUnits) maxNumOfUnits = my nUnitsInLayer[i];
    dx = 1.0 / maxNumOfUnits;
    r = dx / 10;
    Graphics_setInner (g);
    Graphics_setWindow (g, 0.0, 1.0, 0.0, 1.0); 
    for (i=0; i <= my nLayers; i++)
    {
		long j, k;
		double dx2 = dx, x2WC, y2WC = dy / 2 + i * dy;
		double x2 = (maxNumOfUnits - my nUnitsInLayer[i] + 1) * dx2/2;
		/* draw the units */
		if (! dxIsFixed) { dx2 = 1.0 / my nUnitsInLayer[i]; x2 = dx2 / 2; }
		if (i == 0)
		{
	    	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_TOP);
	    	for (x2WC = x2, j=1; j <= my nInputs; j++)
	    	{
				Graphics_arrow (g, x2WC, y2WC - r - dy/4, x2WC, y2WC - r);
				x2WC += dx2;
	    	}
		}  
		Graphics_setColour (g, Graphics_RED);
		for (x2WC = x2, j=1; j <= my nUnitsInLayer[i]; j++)
		{
	   	 	Graphics_circle (g, x2WC, y2WC, r); x2WC += dx2;
		}
		Graphics_setColour (g, Graphics_BLACK);
		if (i > 0)
		{
	    	double dx1 = dx, x1 = (maxNumOfUnits - my nUnitsInLayer[i-1] + 1) * dx1 / 2;
	   		double y1WC = y2WC - dy;
	    	if (! dxIsFixed) { dx1 = 1.0 / my nUnitsInLayer[i-1]; x1 = dx1 / 2; }
	    	x2WC = x2;	
	    	for (j=1; j <= my nUnitsInLayer[i]; j++)
	    	{
				double x1WC = x1;
				for (k=1; k <= my nUnitsInLayer[i-1]; k++)
				{
		    		double xd = x2WC - x1WC;
		    		double cosa = xd / sqrt (xd*xd + dy * dy);
		    		double sina = dy / sqrt (xd*xd + dy * dy);
		    		Graphics_line (g, x1WC + r * cosa, y1WC + r * sina,
					     x2WC - r * cosa, y2WC - r * sina);
		    		x1WC += dx1;
				}
				x2WC += dx2;
	    	}
		}
		if (i == my nLayers)
		{
	    	x2WC = x2;
	    	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_BOTTOM);
	    	for (j=1; j <= my nOutputs; j++)
	    	{
				Graphics_arrow (g, x2WC, y2WC + r, x2WC, y2WC + r + dy / 4);
				if (my outputCategories) Categories_drawItem (my outputCategories, g, j, x2WC, y2WC + dy/4); 
				x2WC += dx2;
	    	}
		}
    }
    Graphics_unsetInner (g);
}

void FFNet_drawActivation (FFNet me, Graphics g)
{
    long i, node = 1, maxNumOfUnits = my nUnitsInLayer[0];
    int dxIsFixed = 1;
    double r, dx, dy = 1.0 / ( my nLayers + 1);

    Graphics_setInner (g);
    Graphics_setWindow (g, 0.0, 1.0, 0.0, 1.0); 
    for (i=1; i <= my nLayers; i++)
		if (my nUnitsInLayer[i] > maxNumOfUnits) maxNumOfUnits = my nUnitsInLayer[i];
    dx = 1.0 / maxNumOfUnits;
    r = dx / 5;
    for (i=0; i <= my nLayers; i++, node++)
    {
		long j;
		double dx2 = dx, x2WC, y2WC = dy / 2 + i * dy;
		double x2 = (maxNumOfUnits - my nUnitsInLayer[i] + 1) * dx2/2;
		if (! dxIsFixed) { dx2 = 1.0 / my nUnitsInLayer[i]; x2 = dx2 / 2; }
		x2WC = x2;
		for (j=1; j <= my nUnitsInLayer[i]; j++, node++)
		{
	    	double activity = my activity[node];
	    	double radius = r * (fabs (activity) < 0.05 ? 0.05 : r * fabs (activity));
	    	Graphics_setColour (g, activity < 0 ? Graphics_BLACK : Graphics_RED);
	    	if (activity < 0) Graphics_fillCircle (g, x2WC, y2WC, radius);
	    	Graphics_circle (g, x2WC, y2WC, radius);
	    	x2WC += dx2;
		}
    }
    Graphics_setColour (g, Graphics_BLACK);
    Graphics_unsetInner (g);
}

void FFNet_drawWeightsToLayer (FFNet me, Graphics g, int layer, int scaling, int garnish)
{
    Matrix weights;
    
    if (layer < 1 || layer > my nLayers || 
    	! (weights = FFNet_weightsToMatrix (me, layer, 0))) return;
    Matrix_scale (weights, scaling);
    Matrix_drawAsSquares (weights, g, 0, 0, 0, 0, 0);
    if (garnish)
    {
    	double x1WC, x2WC, y1WC, y2WC; char text[30];
		Graphics_inqWindow (g, & x1WC, & x2WC, & y1WC, & y2WC);
		sprintf (text, "Units in layer %ld ->", layer);
		Graphics_textBottom (g, 0, text);
		if (layer == 1) strcpy (text, "Input units ->");
		else sprintf (text, "Units in layer %ld ->", layer-1); 
		Graphics_textLeft (g, 0, text);
		/* how do I find out the current settings ??? */
		Graphics_setTextAlignment (g, Graphics_RIGHT, Graphics_HALF);
		Graphics_setInner (g);
		Graphics_text (g, 0.5, weights->ny, "bias");
		Graphics_unsetInner (g);
    }
    forget (weights);
}

void FFNet_drawCostHistory (FFNet me, Graphics g, long iFrom, long iTo, 
    double costMin, double costMax, int garnish)
{
    if (my minimizer) Minimizer_drawHistory (my minimizer, g, iFrom, iTo,
		costMin, costMax, garnish);
}

Collection FFNet_createIrisExample (long nHidden)
{
	TableOfReal iris = NULL; Collection c = NULL;
	FFNet me = NULL; Pattern thee = NULL; Categories him = NULL, uniq = NULL;
	long i, j, nLayers = 0, nUnitsInLayer[5];
	nUnitsInLayer[0] = 4; /*inputs */
	if (nHidden > 0) nUnitsInLayer[++nLayers] = nHidden;
	nUnitsInLayer[++nLayers] = 3; /* categories */
	if (! (c = Collection_create (classData, 3)) ||
		! (uniq = Categories_sequentialNumbers (3)) ||
		! (me = FFNet_create (nLayers, nUnitsInLayer, uniq, 0, FFNet_NONLIN_SIGMOID)) ||
		! Collection_addItem (c, me) ||
		! (iris = TableOfReal_createIrisDataset ()) ||
		! TableOfReal_split (iris, 0, 0, 0, 0, &thee, &him) ||
		! Collection_addItem (c, thee) ||
		! Collection_addItem (c, him)) goto end;
	/*
		Scale data to interval [0-1]
	*/
	for (i=1; i <= 150; i++)
	{
		for (j=1; j <= 4; j++) thy z[i][j] /= 10.0;
	}
	Thing_setName (me, "iris");
	Thing_setName (thee, "iris");
	Thing_setName (him, "iris"); 
end:
	forget (uniq); forget (iris);
	if (! Melder_hasError()) return c;
	forget (c);
	return NULL;
}
