/**********************************************************************
 * nnet.c
 *
 * (C) 2000 J. Ahrensfeld
 *
 **********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "nnfile.h"
#include "nntypes.h"
#include "nnet.h"

/**********************************************************************/
//#define COOL_FEEDFORWARD
//#define COOL_BACKPROP

/**********************************************************************/
char *unitStr[] =
{
	"Bias", "Input", "Hidden", "Output"
};

/**********************************************************************/
UINT32 NetInit(struct _sNNET *pObj, NETFILE *pNetData, UINT32 netID)
{
	UINT32 i, j, k, layCnt, numUnits;

	pObj->nUnits = 0;
	pObj->ID = netID;
	pObj->pUnit = NULL;
	pObj->pLayer = NULL;
	pObj->nLayers = 0;
	pObj->sse = 0.0;
	pObj->pInputID = 0;
	pObj->pHiddenID = 0;
	pObj->pOutputID = 0;
	pObj->nIn = 0;
	pObj->nOut = 0;
	pObj->nHidden = 0;
	pObj->nnType = pNetData->nntype;
	pObj->initwrange = pNetData->initwrange;

	pObj->nLayers = 2 + pNetData->nHiddenLayers;

	pObj->pLayer = (LAYER*)nnMalloc(pObj->pLayer, pObj->nLayers*sizeof(LAYER));

	layCnt = 0;
	srand(12345);
	rand((UINT16)time(NULL));

	/* Init Units */
	NetAddUnit(pObj, UTYPE_BIAS, 0);			/* Bias Unit has index 0 */

	for(i=0; i < (UINT32)pNetData->LInfo[0].nNeurons; i++)
		NetAddUnit(pObj, UTYPE_INPUT, 0);		/* Input Units */

	for(layCnt++; layCnt < pObj->nLayers-1; layCnt++)
		for(i=0; i < (UINT32)pNetData->LInfo[layCnt].nNeurons; i++)
			NetAddUnit(pObj, UTYPE_HIDDEN, pNetData->LInfo[layCnt].NeuronType);/* Hidden Units */

	for(i=0; i < (UINT32)pNetData->LInfo[layCnt].nNeurons; i++)
		NetAddUnit(pObj, UTYPE_OUTPUT, pNetData->LInfo[layCnt].NeuronType);	/* Output Units */

	j = 0;
	for(i=0; i < pObj->nLayers; i++)
	{
		numUnits = pNetData->LInfo[i].nNeurons;
		if (i==0)
			numUnits++;		// +1 for bias unit

		LayerInit(&pObj->pLayer[i], &pObj->pUnit[j], numUnits, i);
		j += numUnits;
	}

	for(k=0; k < pObj->nUnits; k++)
	{
		pObj->pUnit[k].ppRecurr = (FLOAT64**)realloc(pObj->pUnit[k].ppRecurr, pObj->nUnits*sizeof(FLOAT64*));
		pObj->pUnit[k].ppRecurrOld = (FLOAT64**)realloc(pObj->pUnit[k].ppRecurrOld, pObj->nUnits*sizeof(FLOAT64*));
		for (i=0; i < pObj->nUnits; i++)
		{
			pObj->pUnit[k].ppRecurr[i] = (FLOAT64*)malloc(pObj->nUnits*sizeof(FLOAT64));
			pObj->pUnit[k].ppRecurrOld[i] = (FLOAT64*)malloc(pObj->nUnits*sizeof(FLOAT64));
			memset(pObj->pUnit[k].ppRecurr[i], 0, (pObj->nUnits)*sizeof(FLOAT64));
			memset(pObj->pUnit[k].ppRecurrOld[i], 0, (pObj->nUnits)*sizeof(FLOAT64));
		}
	}

	return EXIT_SUCCESS;
}

UINT32 NetWire(struct _sNNET *pObj, UINT32 nnType)
{
	UINT32 layer, i, j, firstUnitInLayer_i, firstUnitInLayer_j;


	switch(nnType)
	{
		case NNTYPE_MLP:

		/* Wire bias unit */
		layer=0;
		for(j=pObj->pLayer[layer].nUnits; j < pObj->nUnits; j++)
			NetAddConnection(pObj,0, j);

		/* Wire input units */
		firstUnitInLayer_i = 1;
		for (i=0; i < pObj->pLayer[layer].nUnits-1; i++)
		{
			firstUnitInLayer_j = pObj->pLayer[layer+1].ppUnit[0]->globalUnitID;
			for(j=0; j < pObj->pLayer[layer+1].nUnits; j++)
			{
				NetAddConnection(pObj,i+firstUnitInLayer_i, j+firstUnitInLayer_j);
			}

		}

		/* Wire hidden units */
		layer++;
		for (layer; layer < pObj->nLayers-1; layer++)
		{
			firstUnitInLayer_i = pObj->pLayer[layer].ppUnit[0]->globalUnitID;
			for (i=0; i < pObj->pLayer[layer].nUnits; i++)
			{
				firstUnitInLayer_j = pObj->pLayer[layer+1].ppUnit[0]->globalUnitID;
				for(j=0; j < pObj->pLayer[layer+1].nUnits; j++)
				{
					NetAddConnection(pObj,i+firstUnitInLayer_i, j+firstUnitInLayer_j);
				}

			}
		}
		break;

		case NNTYPE_RECURRENT:

		/* Wire bias unit */
		/* Wire input units */
		firstUnitInLayer_i = 1;
		for (i=0; i < pObj->pLayer[0].nUnits; i++)
		{
			for(j=pObj->pLayer[0].nUnits; j < pObj->nUnits; j++)
			{
				NetAddConnection(pObj,i, j);
			}

		}
		/* Wire hidden units */
		for (i=pObj->pLayer[0].nUnits; i < pObj->nUnits; i++)
		{
			for(j=pObj->pLayer[0].nUnits; j < pObj->nUnits; j++)
			{
				NetAddConnection(pObj,i, j);
			}

		}
		break;

		default:
		break;
	}
	return EXIT_SUCCESS;
}

UINT32 NetFree(struct _sNNET *pObj)
{
	UINT32 i;

	for (i=0; i < pObj->nUnits; i++)
		UnitFree(&pObj->pUnit[i]);

	for (i=0; i < pObj->nLayers; i++)
		LayerFree(&pObj->pLayer[i]);

	nnFree(pObj->pUnit);
	nnFree(pObj->pLayer);

	return EXIT_SUCCESS;
}

UINT32 NetAddUnit(struct _sNNET *pObj, UINT32 unitType, UINT32 FType)
{
	UINT32 uid;

	uid = pObj->nUnits;
	pObj->pUnit = (UNIT*)realloc(pObj->pUnit, (uid+1)*sizeof(UNIT));

	switch(unitType)
	{
		case UTYPE_INPUT:
			pObj->pInputID = (UINT32*)realloc(pObj->pInputID, (pObj->nIn+1)*sizeof(UINT32));
			pObj->pInputID[pObj->nIn] = uid;
			pObj->nIn++;
			break;

		case UTYPE_HIDDEN:
			pObj->pHiddenID = (UINT32*)realloc(pObj->pHiddenID, (pObj->nHidden+1)*sizeof(UINT32));
			pObj->pHiddenID[pObj->nHidden] = uid;
			pObj->nHidden++;
			break;

		case UTYPE_OUTPUT:
			pObj->pOutputID = (UINT32*)realloc(pObj->pOutputID, (pObj->nOut+1)*sizeof(UINT32));
			pObj->pOutputID[pObj->nOut] = uid;
			pObj->nOut++;
			break;
	}

	UnitInit(&pObj->pUnit[uid], unitType, FType, uid);
	pObj->nUnits++;

	return uid;
}

void NetPrint(struct _sNNET *pObj)
{
	UINT32 layer, i, j, id, nAxonTo, firstUnitInLayer_i;

	printf("Connection Table\n");

	/* Wire hidden units */
	for (layer=0; layer < pObj->nLayers; layer++)
	{
		firstUnitInLayer_i = pObj->pLayer[layer].ppUnit[0]->globalUnitID;
		printf("*** LAYER %d ***\n",layer);
		for (i=0; i < pObj->pLayer[layer].nUnits; i++)
		{
			id = pObj->pUnit[i+firstUnitInLayer_i].globalUnitID;
			nAxonTo = pObj->pUnit[id].nAxonTo;
			if (!nAxonTo)
				continue;
			printf("Axon of %s unit %d.%d (ID #%d) connected to \n",unitStr[pObj->pUnit[id].UType],pObj->pUnit[id].layerID,pObj->pUnit[id].unitID, id);

			for (j=0; j < nAxonTo; j++)
			{
				id = pObj->pUnit[i+firstUnitInLayer_i].pAxonTo[j].unitID;
				printf("%s unit %d.%d (ID #%d)\n",unitStr[pObj->pUnit[id].UType],pObj->pUnit[id].layerID,pObj->pUnit[id].unitID,id);
			}

		}
	}
	printf("\n");

}

UINT32 LayerInit(struct _sLAYER *pObj, struct _sUNIT *pUnit, UINT32 nUnits, UINT32 LayID)
{
	UINT32 i;
	pObj->ID = LayID;
	pObj->nUnits = nUnits;

	pObj->ppUnit = (UNIT**) malloc(pObj->nUnits*sizeof(UNIT*));

	for(i=0; i < pObj->nUnits; i++)
	{
		pUnit[i].layerID = LayID;
		pUnit[i].unitID = i;
		pObj->ppUnit[i] = &pUnit[i];
	}

	return EXIT_SUCCESS;
}

UINT32 LayerFree(struct _sLAYER *pObj)
{
	free(pObj->ppUnit);

	return EXIT_SUCCESS;
}


UINT32 UnitInit(struct _sUNIT *pObj, UINT32 type, UINT32 FType, UINT32 ID)
{
	pObj->globalUnitID = ID;
	pObj->layerID = 0;
	pObj->unitID = 0;
	pObj->UType = type;

	pObj->nAxonTo = 0;
	pObj->nWeights = 0;

	pObj->pAxonFrom = (UINT32*)malloc(sizeof(UINT32));
	pObj->pAxonTo = (WTINFO*)malloc(sizeof(WTINFO));
	pObj->pWeight = (FLOAT64*)malloc(sizeof(FLOAT64));
	pObj->pSlope = (FLOAT64*)malloc(sizeof(FLOAT64));
	pObj->pdWeight = (FLOAT64*)malloc(sizeof(FLOAT64));
	pObj->ppRecurr = (FLOAT64**)malloc(sizeof(FLOAT64*));
	*pObj->ppRecurr = (FLOAT64*)malloc(sizeof(FLOAT64));
	pObj->ppRecurrOld = (FLOAT64**)malloc(sizeof(FLOAT64*));
	*pObj->ppRecurrOld = (FLOAT64*)malloc(sizeof(FLOAT64));

	*pObj->pAxonFrom = UNIT_ID_INPUT;
	pObj->pAxonTo->unitID = UNIT_ID_OUTPUT;
	*pObj->pWeight = 1.0;
	*pObj->pSlope = 0.0;
	*pObj->pdWeight = 0.0;
	**pObj->ppRecurr = 0.12345;
	**pObj->ppRecurrOld = 0.12345;

	pObj->net   = 0.0;
	pObj->error = 0.0;
	pObj->delta = 0.0;
	pObj->axon  = 0.0;
	pObj->lastAxon = 0.0;

	pObj->F = 0;
	pObj->Fd = 0;
	pObj->FType = 0;

	if (type==UTYPE_BIAS)
	{
		pObj->axon  = 1.0;
		pObj->lastAxon  = pObj->axon;
		return EXIT_SUCCESS;
	}

	if (type == UTYPE_INPUT)
		return EXIT_SUCCESS;

	UnitSetFunc(pObj, FType);

	return EXIT_SUCCESS;
}

UINT32 UnitSetFunc(struct _sUNIT *pObj, UINT32 FuncType)
{
	switch(FuncType)
	{
		case Lin:
		pObj->F = Linear;
		pObj->Fd = dLinear;
		pObj->FType = Lin;
		break;

		case Tanh:
		pObj->F = Tanh2;
		pObj->Fd = dTanh;
		pObj->FType = Tanh;
		break;

		case ASig:
		pObj->F = ASigmoid;
		pObj->Fd = dSigmoid;
		pObj->FType = ASig;
		break;

		case SSig:
		pObj->F = SSigmoid;
		pObj->Fd = dSSigmoid;
		pObj->FType = SSig;
		break;

		case Sig:
		pObj->F = Sigmoid;
		pObj->Fd = dSigmoid;
		pObj->FType = Sig;
		break;

		case FSig:
		pObj->F = FSigmoid;
		pObj->Fd = dFSigmoid;
		pObj->FType = FSig;
		break;

		default:
		pObj->F = Sigmoid;
		pObj->Fd = dSigmoid;
		pObj->FType = Sig;
		break;
	}
	return EXIT_SUCCESS;
}

UINT32 NetAddConnection(struct _sNNET *pObj, UINT32 ID_i, UINT32 ID_j)
{
	UINT32 *pFrom = NULL, i, k;
	FLOAT64 *pWeight = NULL, *pdWeight = NULL;
	WTINFO *pTo = NULL;

	UNIT *pUnit_i, *pUnit_j;

	if (ID_j==0 || ID_j >= pObj->nUnits)
		return EXIT_FAILURE;

	if (ID_i >= pObj->nUnits)
		return EXIT_FAILURE;

	pUnit_i = &pObj->pUnit[ID_i];
	pUnit_j = &pObj->pUnit[ID_j];

	for(i=0; i < pUnit_i->nAxonTo; i++)
	{
		if (pUnit_i->pAxonTo[i].unitID == ID_j)
		{
			printf("Axon of unit #%d already connected to unit #%d\n",ID_i,ID_j);
			return EXIT_FAILURE;
		}
	}

	/* Unit j */
	pUnit_j->pAxonFrom = realloc(pUnit_j->pAxonFrom, (pUnit_j->nWeights + 1)*sizeof(UINT32));
	pUnit_j->pWeight = realloc(pUnit_j->pWeight,(pUnit_j->nWeights + 1)*sizeof(FLOAT64));
	pUnit_j->pSlope = realloc(pUnit_j->pSlope,(pUnit_j->nWeights + 1)*sizeof(FLOAT64));
	pUnit_j->pdWeight = realloc(pUnit_j->pdWeight,(pUnit_j->nWeights + 1)*sizeof(FLOAT64));

	pUnit_j->pAxonFrom[pUnit_j->nWeights] = pUnit_i->globalUnitID;
	pUnit_j->pWeight[pUnit_j->nWeights] = 2.0*(0.5 - (FLOAT64)rand()/(FLOAT64)RAND_MAX)*pObj->initwrange;
	pUnit_j->pdWeight[pUnit_j->nWeights] = 0.0;
	pUnit_j->pSlope[pUnit_j->nWeights] = 0.0;

	/* Unit i */
	pUnit_i->pAxonTo = realloc(pUnit_i->pAxonTo, (pUnit_i->nAxonTo + 1)*sizeof(WTINFO));
	pUnit_i->pAxonTo[pUnit_i->nAxonTo].unitID = pUnit_j->globalUnitID;
	pUnit_i->pAxonTo[pUnit_i->nAxonTo].weightID = pUnit_j->nWeights;
	pUnit_i->nAxonTo++;
	pUnit_j->nWeights++;

	return EXIT_SUCCESS;
}
UINT32 NetRTRLFeedForward(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pOut)
{
	UINT32 k, i, id;
	UNIT * pUnit_i, *pUnit_k;

	// Input units (one pattern per training interval)
	for (k=0; k < pObj->nIn; k++)
	{
		id = pObj->pInputID[k];
		pObj->pUnit[id].axon = pIn[k];
	}
	// ************* FEED FORWARD *********************
	// Calc net(t+1)
	// Hidden Units (update every time step)
	for (k=0; k < pObj->nHidden; k++)
	{
		id = pObj->pHiddenID[k];
		pUnit_k = &pObj->pUnit[id];
		pUnit_k->net =0.0;

		for (i =0; i < pUnit_k->nWeights; i++)
		{
			pUnit_i = &pObj->pUnit[pUnit_k->pAxonFrom[i]];
			pUnit_k->net += pUnit_k->pWeight[i]*pUnit_i->axon;
		}
	}
	// Output units (update every time step)
	for (k=0; k < pObj->nOut; k++)
	{
		id = pObj->pOutputID[k];
		pUnit_k = &pObj->pUnit[id];
		pUnit_k->net =0.0;

		for (i =0; i < pUnit_k->nWeights; i++)
		{
			pUnit_i = &pObj->pUnit[pUnit_k->pAxonFrom[i]];
			pUnit_k->net += pUnit_k->pWeight[i]*pUnit_i->axon;
		}

	}

	/// Axon(t+1)
	for (k=0; k < pObj->nUnits; k++)
	{
		pObj->pUnit[k].lastAxon = pObj->pUnit[k].axon;
		if (pObj->pUnit[k].F)
		{
			pObj->pUnit[k].axon = pObj->pUnit[k].F(pObj->pUnit[k].net);
		}
	}

	for (k=0; k < pObj->nOut; k++)
	{
		pUnit_k = &pObj->pUnit[pObj->pOutputID[k]];
		pOut[k] = pUnit_k->axon;
	}

	return EXIT_SUCCESS;
}
UINT32 NetRTRL(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pOut, FLOAT64 *pTarget, FLOAT64 alpha, FLOAT64 eta, UINT32 nPattern, UINT32 trainInterval)
{

	UINT32 j, k, i, t, id, p, n, uid_i, uid_j;
	UNIT * pUnit_i, *pUnit_j, *pUnit_k, *pUnit_n;
	FLOAT64 sum, err;

	pObj->sse = 0;

	for (p=0; p < nPattern; p++)
	{
		// Input units (one pattern per training interval)
		NetRTRLFeedForward(pObj, &pIn[p*pObj->nIn], &pOut[p*pObj->nOut]);
		
		// delta(t)
		for (j=0; j < pObj->nUnits; j++)
		{

			uid_j = j;
			pUnit_j = &pObj->pUnit[uid_j];

			for (i=0; i < pObj->nUnits; i++)
			{

				uid_i = i;
				pUnit_i = &pObj->pUnit[uid_i];

				for (k=0; k < pObj->nUnits; k++)
				{
					pUnit_k = &pObj->pUnit[k];

					sum = 0.0;
					for (n=0; n < pUnit_k->nWeights; n++)
					{
						pUnit_n = &pObj->pUnit[pUnit_k->pAxonFrom[n]];
						sum += pUnit_k->pWeight[n]*pUnit_n->ppRecurrOld[uid_i][uid_j];

					}
					if (pUnit_k->Fd)
					{
						if (uid_i == k)
							pUnit_k->ppRecurr[uid_i][uid_j] = pUnit_k->Fd(pUnit_k->axon)*(sum+pUnit_j->lastAxon);
						else
							pUnit_k->ppRecurr[uid_i][uid_j] = pUnit_k->Fd(pUnit_k->axon)*sum;

					}
				}
			}
		}

		for (k=0; k < pObj->nUnits; k++)
			for (i=0; i < pObj->nUnits; i++)
				for (j=0; j < pObj->nUnits; j++)
					pObj->pUnit[k].ppRecurrOld[i][j] = pObj->pUnit[k].ppRecurr[i][j];

		// dW(t)
		for (k=0; k < pObj->nOut; k++)
		{
			pUnit_k = &pObj->pUnit[pObj->pOutputID[k]];
			err = pTarget[p*pObj->nOut+k] - pUnit_k->axon;
			pUnit_k->error = err;
//				pUnit_k->error = log10(fabs((1 + err)/(1 - err)));
//			err = log10(fabs((1 + err)/(1 - err)));
			pObj->sse += err*err;

			for (i=0; i < pObj->nUnits; i++)
			{
				uid_i = i;
				pUnit_i = &pObj->pUnit[i];
				for (j=0; j < pUnit_i->nWeights; j++)
				{
					uid_j = pUnit_i->pAxonFrom[j];
					pUnit_j = &pObj->pUnit[uid_j];
 					pUnit_i->pdWeight[j] += alpha*err*pUnit_k->ppRecurr[uid_i][uid_j] + (1.0-alpha)*pUnit_i->pdWeight[j];

				}
			}
		}
		// After training interval adjust Hidden and Output Units
		for (i=0; i < pObj->nUnits; i++)
		{
			pUnit_i = &pObj->pUnit[i];
			for (j=0; j < pUnit_i->nWeights; j++)
			{
				pUnit_i->pWeight[j] += eta*pUnit_i->pdWeight[j];
				pUnit_i->pdWeight[j] = 0;
			}
		}
	}
	pObj->sse *= 0.5;
	return EXIT_SUCCESS;

}

UINT32 NetTrain(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pOut, FLOAT64 *pTarget, FLOAT64 eta, FLOAT64 alpha, UINT32 len)
{
	UINT32 j, k, i, p, n;
	UNIT *pUnit_j, *pUnit_k;
	FLOAT64 slope;

	n=0;
	pObj->sse =0.0;
	for (p=0; p < len; p++)
	{
		NetFeedForward(pObj, &pIn[p*(pObj->nIn)], &pOut[n]);
		NetBackPropErr(pObj, &pIn[p*(pObj->nIn)], &pTarget[n]);

		/* Adjust Hidden and Output Units */
		for (k=pObj->nIn+1; k < pObj->nUnits; k++)
		{
			pUnit_k = &pObj->pUnit[k];
			slope = 0.0;
			for (i =0; i < pUnit_k->nWeights; i++)
			{
				j = pUnit_k->pAxonFrom[i];
				pUnit_j = &pObj->pUnit[j];
				pUnit_k->pdWeight[i] = alpha*pUnit_j->axon*pUnit_k->delta + (1.0-alpha)*pUnit_k->pdWeight[i];
				pUnit_k->pWeight[i] += eta*pUnit_k->pdWeight[i];
			}
		}
		n++;
	}
	pObj->sse *= 0.5;
	return EXIT_SUCCESS;
}

UINT32 NetTrain2(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pOut, FLOAT64 *pError, FLOAT64 eta, FLOAT64 alpha, UINT32 len)
{
	UINT32 j, k, i, p, n;
	UNIT *pUnit_j, *pUnit_k;

	n=0;
	for (p=0; p < len; p++)
	{
		NetBackPropErr2(pObj, &pIn[p*(pObj->nIn)], &pError[n]);

		/* Adjust Hidden and Output Units */
		for (k=pObj->nIn+1; k < pObj->nUnits; k++)
		{
			pUnit_k = &pObj->pUnit[k];
			for (i =0; i < pUnit_k->nWeights; i++)
			{
				j = pUnit_k->pAxonFrom[i];
				pUnit_j = &pObj->pUnit[j];
				pUnit_k->pdWeight[i] = eta*pUnit_j->axon*pUnit_k->delta;
				pUnit_k->pWeight[i] += alpha*pUnit_k->pdWeight[i];
			}
		}
		n++;
	}
	return EXIT_SUCCESS;
}

UINT32 NetFeedForward(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pOut)
{
	UINT32 j, k, m, i;
	UNIT *pUnit_j, *pUnit_k, **ppU;

	m=0;
	/* Input Layer */
	for (k=1; k < pObj->nIn+1; k++)
	{
		pUnit_k = &pObj->pUnit[k];
		pUnit_k->axon = pIn[m++];
	}
#ifndef COOL_FEEDFORWARD

	/* Input, Hidden and Output Layers */
	for (k=pObj->pLayer[0].nUnits; k < pObj->nUnits; k++)
	{
		pUnit_k = &pObj->pUnit[k];
		pUnit_k->net =0.0;
		for (i =0; i < pUnit_k->nWeights; i++)
		{
			j = pUnit_k->pAxonFrom[i];
			pUnit_j = &pObj->pUnit[j];
			pUnit_k->net += pUnit_k->pWeight[i]*pUnit_j->axon;
		}
		pUnit_k->axon = pUnit_k->F(pUnit_k->net);
	}
#else
#pragma message ("nnet.c: Using cool feed forward computing.")
	/********************************************************/
	/* TEST: don't care for layer organization */
	/********************************************************/
	for (i=0; i < pObj->nUnits; i++)
		pObj->pUnit[i].net = 0.0;

	for (j=0; j < pObj->nUnits; j++)
	{
		pUnit_j = &pObj->pUnit[j];
		if (pUnit_j->F)
			pUnit_j->axon = pUnit_j->F(pUnit_j->net);

		for (i=0; i < pUnit_j->nAxonTo; i++)
		{
			k = pUnit_j->pAxonTo[i].unitID;
			wid = pUnit_j->pAxonTo[i].weightID;
			pUnit_k = &pObj->pUnit[k];
			pUnit_k->net += pUnit_k->pWeight[wid]*pUnit_j->axon;
		}
	}
#endif
	/********************************************************/
	ppU = pObj->pLayer[pObj->nLayers-1].ppUnit;
	for (j = 0; j < pObj->nOut; j++)
	{
		pOut[j] = (*ppU)->axon;
		ppU++;
	}
	return EXIT_SUCCESS;
}

UINT32 NetBackPropErr(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pTarget)
{
	UINT32 j, k, i, target, numLayers, layer;
	UNIT **ppUnit_j, **ppUnit_k, *pUnit_j, *pUnit_k;
	FLOAT64 err;


	numLayers = pObj->nLayers;

	/* Zero Errors of all units */
	for (j = pObj->pLayer[0].nUnits; j < pObj->nUnits; j++)
	{
		pObj->pUnit[j].error = 0;
	}

	/* Calc Error Output Layer */
	target=0;
	ppUnit_j = pObj->pLayer[numLayers-1].ppUnit;

	for (j=0; j < pObj->pLayer[numLayers-1].nUnits; j++)
	{
		pUnit_j = *(ppUnit_j++);
		err = pTarget[target++] - pUnit_j->axon;
		pUnit_j->error = log10(fabs((1 + err)/(1 - err)));
		pUnit_j->delta = pUnit_j->error*pUnit_j->Fd(pUnit_j->axon);
		pObj->sse += err*err;
	}

#ifndef COOL_BACKPROP

	/* Calc Error Hidden Layer */
	for (layer = numLayers-1; layer > 1; layer--)
	{
		ppUnit_k = pObj->pLayer[layer].ppUnit;
		for (k =0; k < pObj->pLayer[layer].nUnits; k++)
		{
			pUnit_k = *(ppUnit_k++);
			for (i=1; i < pUnit_k->nWeights; i++)
			{
				j = pUnit_k->pAxonFrom[i];
				pUnit_j = &pObj->pUnit[j];
				if (pUnit_j->UType == UTYPE_HIDDEN)
				{
					pUnit_k->pSlope[j] += pUnit_k->delta*pUnit_k->axon;
					pUnit_j->error += pUnit_k->pWeight[i]*pUnit_k->delta;
					pUnit_j->delta = pUnit_j->error*pUnit_j->Fd(pUnit_j->axon);
				}
			}
		}
	}

#else
#pragma message ("nnet.c: Using cool error back propagation.")
#error ("nnet.c: Implement cool error back propagation.")
#endif

return EXIT_SUCCESS;
}

UINT32 NetBackPropErr2(struct _sNNET *pObj, FLOAT64 *pIn, FLOAT64 *pError)
{
	UINT32 j, k, i, error, numLayers, layer;
	UNIT **ppUnit_j, **ppUnit_k, *pUnit_j, *pUnit_k;
	FLOAT64 err;


	numLayers = pObj->nLayers;

	/* Zero Errors of all units */
	for (j = pObj->pLayer[0].nUnits; j < pObj->nUnits; j++)
	{
		pObj->pUnit[j].error = 0;
	}

	/* Calc Error Output Layer */
	pObj->sse =0.0;
	error=0;
	ppUnit_j = pObj->pLayer[numLayers-1].ppUnit;

	for (j=0; j < pObj->pLayer[numLayers-1].nUnits; j++)
	{
		pUnit_j = *(ppUnit_j++);
		err = pError[error++];
		pUnit_j->error = log10(fabs((1 + err)/(1 - err)));
		pUnit_j->delta = pUnit_j->error*pUnit_j->Fd(pUnit_j->axon);
		pObj->sse += err*err;
	}

	/* Calc Error Hidden Layer */
	for (layer = numLayers-1; layer > 1; layer--)
	{
		ppUnit_k = pObj->pLayer[layer].ppUnit;
		for (k =0; k < pObj->pLayer[layer].nUnits; k++)
		{
			pUnit_k = *(ppUnit_k++);
			for (i=1; i < pUnit_k->nWeights; i++)
			{
				j = pUnit_k->pAxonFrom[i];
				pUnit_j = &pObj->pUnit[j];
				pUnit_j->error += pUnit_k->pWeight[i]*pUnit_k->delta;
				pUnit_j->delta = pUnit_j->error*pUnit_j->Fd(pUnit_j->axon);
			}
		}
	}
	return EXIT_SUCCESS;
}

UINT32 NetWeightUpd(struct _sNNET *pObj, FLOAT64 *pIn, UINT32 len)
{
	return EXIT_SUCCESS;
}

UINT32 UnitFree(struct _sUNIT *pObj)
{
	return EXIT_SUCCESS;
}

void* nnMalloc(void *pBuffer, UINT32 size)
{

	if (pBuffer == NULL)
		pBuffer = malloc(size);
	else
		pBuffer = NULL;

	return pBuffer;
}

void* nnFree(void *pBuffer)
{

	if (pBuffer != NULL)
		free(pBuffer);

	pBuffer = NULL;

	return pBuffer;
}

/*************************************************************************/
#define K_SIGM	1.0
#define K_TANH	1.0

FLOAT64 Sigmoid(FLOAT64 x)
{
	return 1.0 / (1.0 + exp(-K_SIGM*x));
}

FLOAT64 dSigmoid(FLOAT64 x)
{
	return x*(1.0-x);
}

FLOAT64 Linear(FLOAT64 x)
{
  return x;
}

FLOAT64 dLinear(FLOAT64 x)
{
  return 1;
}

FLOAT64 Tanh2(FLOAT64 x)
{
	FLOAT64 e1 = exp(2*K_TANH*x);
	e1 = (e1-1.0)/(e1+1.0);
	return e1;
}

FLOAT64 dTanh2(FLOAT64 x)
{
	FLOAT64 e1 = exp(2*K_TANH*x);
	FLOAT64 e11 = e1+1.0;

	return 4*K_TANH*e1 / (e11*e11);
}

FLOAT64 ASigmoid(FLOAT64 x)
{
	return 1.0 / (1.0 + exp(-K_SIGM*x));
}

FLOAT64 SSigmoid(FLOAT64 x)
{
	return - 0.5 + 1.0 / (1.0 + exp(-x));
}

FLOAT64 dSSigmoid(FLOAT64 x)
{
	FLOAT64 e1 = exp(K_SIGM*x);
	FLOAT64 e11 = 1.0 + e1;
	return K_SIGM*e1 / (e11*e11);
}

FLOAT64 FSigmoid(FLOAT64 x)
{
	FLOAT64 kx = K_SIGM*x;
	return kx / (1.0 + fabs(kx));
}

FLOAT64 dFSigmoid(FLOAT64 x)
{
	return fabs(x*(1.0-x));
}

FLOAT64 dTanh(FLOAT64 x)
{
	return 1.0-x*x;
}

FLOAT64 NullFunc(FLOAT64 x)
{
	return 0;
}
/**************************************************************************/