/*****************************************************************************
  FILE           : rcc_learn.c
  SHORTNAME      : 
  SNNS VERSION   : 3.2

  PURPOSE        : Functions of RCC
  NOTES          :

  AUTHOR         : Michael Schmalzl
  DATE           : 5.2.92

  CHANGED BY     : Guenter Mamier
  IDENTIFICATION : @(#)rcc_learn.c	1.12 4/7/94
  SCCS VERSION   : 1.12
  LAST CHANGE    : 4/7/94

             Copyright (c) 1990-1994  SNNS Group, IPVR, Univ. Stuttgart, FRG

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

#include <stdio.h>
#include <math.h>
#include <time.h>  
#include <memory.h>
#include <malloc.h>
#include <values.h>

#include "kr_typ.h"      /*  Kernel Types and Constants  */
#include "kr_const.h"    /*  Constant Declarators for SNNS-Kernel  */
#include "kr_def.h"      /*  Default Values  */
#include "kernel.h"      /*  kernel function prototypes  */
#include "kr_mac.h"      /*  Kernel Macros   */
#include "random.h"

#include "kr_ui.h"
#include "cc_type.h"
#include "cc_mac.h"
#include "rcc_learn.ph"
#include "cc_rcc.h"
#include "kr_newpattern.h"


/*****************************************************************************
  FUNCTION : rcc_searchRecurrentLinks

  PURPOSE  : Searches for recurrent links.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
krui_err rcc_searchRecurrentLinks(void)
{
    int CurrentUnit;
    struct Unit *unit_ptr;
    float linkValue;
 
    FOR_ALL_UNITS(unit_ptr){
	if(IS_HIDDEN_UNIT(unit_ptr)){
	    KernelErrorCode = 
		krui_setCurrentUnit(CurrentUnit = GET_UNIT_NO(unit_ptr)); 
	    ERROR_CHECK;
	    if(kr_isConnected(CurrentUnit,&linkValue)) {

		/* delete recurrent link */
		KernelErrorCode = krui_deleteLink(); 
		ERROR_CHECK;

		/* generate recurent link */
		KernelErrorCode = krui_createLink(CurrentUnit,linkValue); 
		ERROR_CHECK;
	    }
	}
    }
    return(KRERR_NO_ERROR);
}



/*****************************************************************************
  FUNCTION : rcc_test

  PURPOSE  : Tests wether to continue learning or not.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static int rcc_test(int StartPattern, int EndPattern, float maxPixelError)
{
    int p,o,pat,sub;
    int start, end;
    Patterns out_pat;
    struct Unit *unitPtr;


    /* compute the necessary sub patterns */

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    if(KernelErrorCode != KRERR_NO_ERROR)
        return (KernelErrorCode);

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    for(p=start; p<=end;p++){
	kr_getSubPatternByNo(&pat,&sub,p);
	rcc_propagateNetForward(pat,sub);

	out_pat = kr_getSubPatData(pat,sub,OUTPUT,NULL);
	FOR_ALL_OUTPUT_UNITS(unitPtr,o){
	    if((fabs(*(out_pat++) - unitPtr->Out.output))>maxPixelError){
		return(CONTINUE_LEARNING);
	    } 
	}
    }
    return(STOP_LEARNING);
}



/*****************************************************************************
  FUNCTION : rcc_propagateNetForward

  PURPOSE  : Propagates a pattern forward through the net.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_propagateNetForward(int PatternNo, int sub_pat_no)
{
    register struct Unit  *inputUnitPtr,*outputUnitPtr,*hiddenUnitPtr;
    register Patterns  in_pat;
    register int dummy;

    in_pat = kr_getSubPatData(PatternNo,sub_pat_no,INPUT,NULL);

    FOR_ALL_INPUT_UNITS(inputUnitPtr,dummy){
	if(inputUnitPtr->out_func == OUT_IDENTITY) {
	    inputUnitPtr->Out.output = inputUnitPtr->act = *in_pat++;
	}else{
	    inputUnitPtr->Out.output = 
		(*inputUnitPtr->out_func) (inputUnitPtr->act = *in_pat++);
	}
    }

    FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,dummy) {
	hiddenUnitPtr->lln = reset[PatternNo];
	hiddenUnitPtr->act = (*hiddenUnitPtr->act_func)(hiddenUnitPtr);
	if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
	    hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	}else{
	    hiddenUnitPtr->Out.output = 
		(*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	}
    }

    FOR_ALL_OUTPUT_UNITS(outputUnitPtr,dummy) {
	outputUnitPtr->lln = reset[PatternNo];
	outputUnitPtr->act = (*outputUnitPtr->act_func) (outputUnitPtr);
	if(outputUnitPtr->out_func == OUT_IDENTITY) {
	    outputUnitPtr->Out.output = outputUnitPtr->act;
	}else{
	    outputUnitPtr->Out.output = 
		(*outputUnitPtr->out_func) (outputUnitPtr->act);
	}
    }
}



/*****************************************************************************
  FUNCTION : rcc_calculateOutputUnitError

  PURPOSE  : Calculates the error of the output units and stores it in the
             array OutputUnitError.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_calculateOutputUnitError(int StartPattern, int EndPattern)
{
    register struct Unit  *inputUnitPtr,*outputUnitPtr,*hiddenUnitPtr;
    register Patterns  in_pat,out_pat;
    register int dummy,o,p;
    int start, end;
    int pat, sub;

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    for(p=start; p<=end;p++){

	kr_getSubPatternByNo(&pat,&sub,p);
	in_pat = kr_getSubPatData(pat,sub,INPUT,NULL);
	out_pat = kr_getSubPatData(pat,sub,OUTPUT,NULL);

	FOR_ALL_INPUT_UNITS(inputUnitPtr,dummy){
	    if(inputUnitPtr->out_func == OUT_IDENTITY) {
		inputUnitPtr->Out.output = inputUnitPtr->act = *in_pat++;
	    }else{
		inputUnitPtr->Out.output = 
		    (*inputUnitPtr->out_func) (inputUnitPtr->act = *in_pat++);
	    }
	}

	FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,dummy) {
	    hiddenUnitPtr->lln = reset[p];
	    hiddenUnitPtr->act = (*hiddenUnitPtr->act_func) (hiddenUnitPtr);
	    if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
		hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	    }else{
		hiddenUnitPtr->Out.output = 
		    (*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	    }
	}
    
	FOR_ALL_OUTPUT_UNITS(outputUnitPtr,o) {
	    outputUnitPtr->lln = reset[p];
	    outputUnitPtr->act = (*outputUnitPtr->act_func)(outputUnitPtr);
	    if(outputUnitPtr->out_func == OUT_IDENTITY) {
		outputUnitPtr->Out.output = outputUnitPtr->act;
	    }else{
		outputUnitPtr->Out.output = 
		    (*outputUnitPtr->out_func) (outputUnitPtr->act);
	    }
	    OutputUnitSumError[o] += 
		(OutputUnitError[p][o] =  
		       (outputUnitPtr->Out.output-(*out_pat++))*
		       ((*outputUnitPtr->act_deriv_func)(outputUnitPtr)+0.1)); 
	} 
    } 
}



/*****************************************************************************
  FUNCTION : rcc_calculateSpecialUnitActivation

  PURPOSE  : Calculates covariance between the special units and the output
             units and stores it in the array CorBetweenSpecialActAndOutError.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_calculateSpecialUnitActivation(int StartPattern, int EndPattern)
{
    register struct Unit  *inputUnitPtr,*specialUnitPtr,*outputUnitPtr,*hiddenUnitPtr;
    register Patterns  in_pat;
    register int dummy,o,s,p;
    int pat,sub;
    int start,end;
  
    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    for(p=start; p<=end;p++){

	kr_getSubPatternByNo(&pat,&sub,p);
	in_pat = kr_getSubPatData(pat,sub,INPUT,NULL);
 
	FOR_ALL_INPUT_UNITS(inputUnitPtr,dummy){
	    if(inputUnitPtr->out_func == OUT_IDENTITY) {
		inputUnitPtr->Out.output = inputUnitPtr->act = *in_pat++;
	    }else{
		inputUnitPtr->Out.output = 
		    (*inputUnitPtr->out_func) (inputUnitPtr->act = *in_pat++);
	    }
	}

	FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,dummy) {
	    hiddenUnitPtr->lln = reset[p];
	    hiddenUnitPtr->act = (*hiddenUnitPtr->act_func) (hiddenUnitPtr);
	    if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
		hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	    }else{
		hiddenUnitPtr->Out.output = 
		    (*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	    }
	}
    
	FOR_ALL_SPECIAL_UNITS(specialUnitPtr,s) {
	    specialUnitPtr->lln = reset[p];
	    specialUnitPtr->act = (*specialUnitPtr->act_func) (specialUnitPtr);
	    if(specialUnitPtr->out_func == OUT_IDENTITY) {
		specialUnitPtr->Out.output = specialUnitPtr->act;
	    }else{
		specialUnitPtr->Out.output = 
		    (*specialUnitPtr->out_func) (specialUnitPtr->act);
	    }
	    SpecialUnitSumAct[s] += 
		SpecialUnitAct[p][s] = 
		    specialUnitPtr->Out.output;
	}
    } 
  
    for(p=start; p<=end;p++){
	FOR_ALL_SPECIAL_UNITS(specialUnitPtr,s) {
	    FOR_ALL_OUTPUT_UNITS(outputUnitPtr,o) {
		CorBetweenSpecialActAndOutError[s][o] += 
		    SpecialUnitAct[p][s] * OutputUnitError[p][o];
	    }
	}
    }    
}



/************* begin rprop routines *********************/


/*****************************************************************************
  FUNCTION : rcc_RPO_trainNet

  PURPOSE  : Minimize the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_RPO_trainNet(int maxNoOfErrorUpdateCycles, float minErrorChange,
                             int outPatience,int StartPattern, int EndPattern, 
			     float epsilonMinus, float epsilonPlus, float dummy, 
			     float **ParameterOutArray, int *NoOfOutParams)
{
    int m,p,counter=0;
    int sub, pat;
    int start, end;
    float oldNetError;
    static float OutParameter[1];

    *NoOfOutParams = 1;
    *ParameterOutArray = OutParameter;


    SumSqError = 0.0;
    cc_initOutputUnits();

    /* compute the necessary sub patterns */

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    if(KernelErrorCode != KRERR_NO_ERROR)
        return;

    /* give oldNetError a meaningful initial value */
    NET_ERROR(OutParameter)=FLOAT_MAX;
    do {
	oldNetError = NET_ERROR(OutParameter);
	for(m=0;m<outPatience;m++) {
	    NET_ERROR(OutParameter) = 0.0;
	    SumSqError = 0.0;
	    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
	    start = kr_AbsPosOfFirstSubPat(StartPattern);
	    end   = kr_AbsPosOfFirstSubPat(EndPattern);
	    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
	    for(p=start; p<=end;p++){
		kr_getSubPatternByNo(&pat,&sub,p);
		rcc_propagateNetForward(pat,sub);
		NET_ERROR(OutParameter) +=rcc_RPO_propagateNetBackward(pat,sub);
	    }
	    rcc_RPO_updateNet(epsilonMinus,epsilonPlus,dummy);
	    if(cc_printOnOff) {
		printf("Epoch: %d NetError: %f \n",++counter,
		       NET_ERROR(OutParameter));
	    } 
	    if((maxNoOfErrorUpdateCycles--) == 0) {
		return;
	    }
	}
    } while(fabs(oldNetError-NET_ERROR(OutParameter)) >= 
	    (minErrorChange * oldNetError));
}



/*****************************************************************************
  FUNCTION : rcc_RPO_propagateNetBackward

  PURPOSE  : Calculate the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static  float rcc_RPO_propagateNetBackward(int PatternNo, int sub_pat_no)
{
    struct Link   *LinkPtr;
    struct Site   *site_ptr;
    struct Unit   *OutputUnitPtr;
    Patterns      out_pat;
    float         error,sum_error,devit;
    int dummy;

    sum_error = 0.0;    
    out_pat = kr_getSubPatData(PatternNo,sub_pat_no,OUTPUT,NULL);

    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,dummy){
	devit =  OutputUnitPtr->Out.output - *(out_pat++);

	sum_error += devit * devit;
	error = devit * ((*OutputUnitPtr->act_deriv_func)(OutputUnitPtr) + 0.1);
	SumSqError += error * error;
 
	BIAS_CURRENT_SLOPE(OutputUnitPtr) += error;

	if (UNIT_HAS_DIRECT_INPUTS(OutputUnitPtr)) {  
	    FOR_ALL_LINKS(OutputUnitPtr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}else{
	    FOR_ALL_SITES_AND_LINKS(OutputUnitPtr,site_ptr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}
    }
    return(sum_error); 
}



/*****************************************************************************
  FUNCTION : rcc_RPO_updateNet

  PURPOSE  : Update the weights of the output units with rprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_RPO_updateNet(float epsilonMinus, float epsilonPlus, 
			      float dummy)
{
    struct Unit *OutputUnitPtr;
    struct Link *LinkPtr;
    float bias_previousSlope,bias_currentSlope,bias_lastWeightChange,
          bias_weightChange=0.0;
    float ln_previousSlope,ln_currentSlope,ln_lastWeightChange,
          ln_weightChange=0.0;
    int o;

    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,o) {
	bias_previousSlope    = BIAS_PREVIOUS_SLOPE(OutputUnitPtr);
	bias_currentSlope     = BIAS_CURRENT_SLOPE(OutputUnitPtr);
	bias_lastWeightChange = 
	    (BIAS_LAST_WEIGHT_CHANGE(OutputUnitPtr) == 0.0) ? 
		(1.0) : (BIAS_LAST_WEIGHT_CHANGE(OutputUnitPtr));
   
	if(bias_currentSlope != 0.0) {
	    if(bias_previousSlope == 0.0) {
		bias_weightChange = 
		    fabs(bias_lastWeightChange) * SIGN(bias_currentSlope);
	    }else if(bias_previousSlope > 0.0) {  
		if(bias_currentSlope > 0.0) {
		    bias_weightChange = epsilonPlus * bias_lastWeightChange;  
		}else if(bias_currentSlope < 0.0) {
		    bias_weightChange = -epsilonMinus * bias_lastWeightChange;
		}
	    }else if(bias_previousSlope < 0.0) {
		if(bias_currentSlope < 0.0) {
		    bias_weightChange = epsilonPlus * bias_lastWeightChange;
		}else if(bias_currentSlope > 0.0) {
		    bias_weightChange = -epsilonMinus * bias_lastWeightChange;
		}
	    }else{
		bias_weightChange = 1.0 * SIGN(bias_currentSlope);
	    }

	    if(fabs(bias_weightChange) < 0.00001) {
		bias_weightChange = 0.00001 * SIGN(bias_weightChange);
	    }
	    if(fabs(bias_weightChange) > 10.0) {
		bias_weightChange = 10.0 * SIGN(bias_weightChange);
	    }
   
	    OutputUnitPtr->bias -= 
		(BIAS_LAST_WEIGHT_CHANGE(OutputUnitPtr) = bias_weightChange); 
	    BIAS_PREVIOUS_SLOPE(OutputUnitPtr) = bias_currentSlope;
	    BIAS_CURRENT_SLOPE(OutputUnitPtr) = 0.0;
	}

	FOR_ALL_LINKS(OutputUnitPtr,LinkPtr) {
   
	    ln_previousSlope    = LN_PREVIOUS_SLOPE(LinkPtr);
	    ln_currentSlope     = LN_CURRENT_SLOPE(LinkPtr);
	    ln_lastWeightChange = 
		(LN_LAST_WEIGHT_CHANGE(LinkPtr) == 0.0) ? 
		    (1.0) : (LN_LAST_WEIGHT_CHANGE(LinkPtr));

	    if(ln_currentSlope != 0.0) {
		if(ln_previousSlope == 0.0) {
		    ln_weightChange = 
			fabs(ln_lastWeightChange) * SIGN(ln_currentSlope);
		}
		else if(ln_previousSlope > 0.0) {  
		    if(ln_currentSlope > 0.0) {
			ln_weightChange = epsilonPlus * ln_lastWeightChange;  
		    }
		    else if(ln_currentSlope < 0.0) {
			ln_weightChange = -epsilonMinus * ln_lastWeightChange;
		    }
		}else if(ln_previousSlope < 0.0) {
		    if(ln_currentSlope < 0.0) {
			ln_weightChange = epsilonPlus * ln_lastWeightChange;
		    }else if(ln_currentSlope > 0.0) {
			ln_weightChange = -epsilonMinus * ln_lastWeightChange;
		    }
		}else{
		    ln_weightChange = 1.0 * SIGN(ln_currentSlope);
		}
    
		if(fabs(ln_weightChange) < 0.00001) {
		    ln_weightChange = 0.00001 * SIGN(ln_weightChange);
		}
		if(fabs(ln_weightChange) > 10) {
		    ln_weightChange = 10 * SIGN(ln_weightChange);
		}

		LinkPtr->weight -= 
		    LN_LAST_WEIGHT_CHANGE(LinkPtr) = ln_weightChange; 
		LN_PREVIOUS_SLOPE(LinkPtr)   = ln_currentSlope;
		LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	    }
	}
    }
}

 
/*****************************************************************************
  FUNCTION : rcc_RPS_trainNet

  PURPOSE  : Maximize the covariance of the special units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_RPS_trainNet(int maxNoOfCovarianceUpdateCycles, 
			     float minCovarianceChange,int specialPatience, 
			     int StartPattern, int EndPattern, 
			     float epsilonMinus, float epsilonPlus, 
			     float dummy, int MaxSpecialUnitNo)
{
    int m,counter=0;
    float oldHighScore,newHighScore=0.0;

    cc_initErrorArrays();
    rcc_calculateOutputUnitError(StartPattern,EndPattern);

    do {
	oldHighScore = newHighScore;
	for(m=0;m<specialPatience;m++) { 
	    counter++;
	    rcc_calculateSpecialUnitActivation(StartPattern,EndPattern);
	    newHighScore = 
		rcc_RPS_propagateNetBackward(StartPattern,EndPattern,counter);
	    rcc_RPS_updateNet(epsilonMinus,epsilonPlus,dummy);
	    cc_initActivationArrays(); 
	    if((maxNoOfCovarianceUpdateCycles--) == 0) {
		return;
	    }
	}
    } while(fabs(newHighScore-oldHighScore) >= 
	    (minCovarianceChange * oldHighScore));
}



/*****************************************************************************
  FUNCTION : rcc_RPS_propagateNetBackward

  PURPOSE  : Calculate the special unit with maximum covariance and return it.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static float rcc_RPS_propagateNetBackward(int StartPattern, int EndPattern, 
					  int counter)
{
    float change=0.0,bestSpecialUnitScore,actPrime,recurrentLinkWeight,dsum;
    int s,o,p,n,h;
    struct Unit *SpecialUnitPtr,*OutputUnitPtr,*hiddenUnitPtr;
    struct Link *LinkPtr,*recurrentLinkPtr;
    int start, end;


    bestSpecialUnitScore = 
	cc_calculateCorrelation(StartPattern,EndPattern,counter); 

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    n = end - start + 1;
    for(p=start; p<=end;p++){

	cc_initInputUnitsWithPattern(p);
	FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,h) { 
	    hiddenUnitPtr->lln = reset[p];
	    hiddenUnitPtr->act = (*hiddenUnitPtr->act_func)(hiddenUnitPtr);
	    if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
		hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	    }else{
		hiddenUnitPtr->Out.output = 
		    (*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	    }
	}
   
	FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,s) {
	    change = 0.0;
	    SpecialUnitPtr->act = SpecialUnitAct[p][s];
	    actPrime = (*SpecialUnitPtr->act_deriv_func)(SpecialUnitPtr);
	    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,o) {
		change -= CorBetweenSpecialActAndOutError[s][o] *
		    ((OutputUnitError[p][o] - OutputUnitSumError[o]/n)/
		     SumSqError);
	    }

	    counter = 0;
	    GET_RECURRENT_LINK(SpecialUnitPtr,recurrentLinkPtr);
	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    recurrentLinkWeight = recurrentLinkPtr->weight; 
	    if(!reset[p]) {	/* calculate slope of the recurrent link */
		dsum = actPrime * (SpecialUnitAct[p-1][s] + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(recurrentLinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum; 
	    }

	    FOR_ALL_NOT_RECURRENT_LINKS(SpecialUnitPtr,LinkPtr) {
		if(reset[p]) {
		    linkArray[s][counter] = 0.0;
		}
		dsum = actPrime * (LinkPtr->to->Out.output + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(LinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum;
	    }

	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    dsum = actPrime * (1 + (recurrentLinkWeight*linkArray[s][counter]));
	    BIAS_CURRENT_SLOPE(SpecialUnitPtr) += change * dsum;
	    linkArray[s][counter] = dsum;
	}
    }
    return(bestSpecialUnitScore);
}



/*****************************************************************************
  FUNCTION : rcc_RPS_updateNet

  PURPOSE  : Update the weights of the special units with rprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_RPS_updateNet(float epsilonMinus,float epsilonPlus, float dummy)
{
    struct Unit *specialUnitPtr;
    struct Link *LinkPtr;
    float bias_previousSlope,bias_currentSlope,bias_lastWeightChange,
          bias_weightChange=0.0;
    float ln_previousSlope,ln_currentSlope,ln_lastWeightChange,
          ln_weightChange=0.0;
    int s;

    FOR_ALL_SPECIAL_UNITS(specialUnitPtr,s) {
	bias_previousSlope    = BIAS_PREVIOUS_SLOPE(specialUnitPtr);
	bias_currentSlope     = BIAS_CURRENT_SLOPE(specialUnitPtr);
	bias_lastWeightChange = 
	    (BIAS_LAST_WEIGHT_CHANGE(specialUnitPtr) == 0.0) ? 
		(1.0) : (BIAS_LAST_WEIGHT_CHANGE(specialUnitPtr));
   
	if(bias_currentSlope != 0.0) {
	    if(bias_previousSlope == 0.0) {
		bias_weightChange = 
		    fabs(bias_lastWeightChange) * SIGN(bias_currentSlope);
	    }else if(bias_previousSlope > 0.0) {  
		if(bias_currentSlope > 0.0) {
		    bias_weightChange = epsilonPlus * bias_lastWeightChange;  
		}else if(bias_currentSlope < 0.0) {
		    bias_weightChange = -epsilonMinus * bias_lastWeightChange;
		}
	    }else if(bias_previousSlope < 0.0) {
		if(bias_currentSlope < 0.0) {
		    bias_weightChange = epsilonPlus * bias_lastWeightChange;
		}else if(bias_currentSlope > 0.0) {
		    bias_weightChange = -epsilonMinus * bias_lastWeightChange;
		}
	    }else {
		bias_weightChange = 1.0 * SIGN(bias_currentSlope);
	    }

	    if(fabs(bias_weightChange) < 0.00001) {
		bias_weightChange = 0.00001 * SIGN(bias_weightChange);
	    }
	    if(fabs(bias_weightChange) > 10.0) {
		bias_weightChange = 10.0 * SIGN(bias_weightChange);
	    }
  
	    specialUnitPtr->bias -= 
		(BIAS_LAST_WEIGHT_CHANGE(specialUnitPtr) = bias_weightChange); 
	    BIAS_PREVIOUS_SLOPE(specialUnitPtr) = bias_currentSlope;
	    BIAS_CURRENT_SLOPE(specialUnitPtr) = 0.0;
	}

	FOR_ALL_LINKS(specialUnitPtr,LinkPtr) {
   
	    ln_previousSlope    = LN_PREVIOUS_SLOPE(LinkPtr);
	    ln_currentSlope     = LN_CURRENT_SLOPE(LinkPtr);
	    ln_lastWeightChange = 
		(LN_LAST_WEIGHT_CHANGE(LinkPtr) == 0.0) ? 
		    (1.0) : (LN_LAST_WEIGHT_CHANGE(LinkPtr));

	    if(ln_currentSlope != 0.0) {
		if(ln_previousSlope == 0.0) {
		    ln_weightChange = 
			fabs(ln_lastWeightChange) * SIGN(ln_currentSlope);
		}else if(ln_previousSlope > 0.0) {  
		    if(ln_currentSlope > 0.0) {
			ln_weightChange = epsilonPlus * ln_lastWeightChange;  
		    }else if(ln_currentSlope < 0.0) {
			ln_weightChange = -epsilonMinus * ln_lastWeightChange;
		    }
		}else if(ln_previousSlope < 0.0) {
		    if(ln_currentSlope < 0.0) {
			ln_weightChange = epsilonPlus * ln_lastWeightChange;
		    }else if(ln_currentSlope > 0.0) {
			ln_weightChange = -epsilonMinus * ln_lastWeightChange;
		    }
		}else {
		    ln_weightChange = 1.0 * SIGN(ln_currentSlope);
		}
    
		if(fabs(ln_weightChange) < 0.00001) {
		    ln_weightChange = 0.00001 * SIGN(ln_weightChange);
		}
		if(fabs(ln_weightChange) > 10) {
		    ln_weightChange = 10 * SIGN(ln_weightChange);
		}

		LinkPtr->weight -= 
		    LN_LAST_WEIGHT_CHANGE(LinkPtr) = ln_weightChange; 
		LN_PREVIOUS_SLOPE(LinkPtr)   = ln_currentSlope;
		LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	    }
	}
    }
}
/************* end rprop routines *********************/



/************* begin quickprop routines *********************/

/*****************************************************************************
  FUNCTION : rcc_QPO_trainNet

  PURPOSE  : Minimize the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_QPO_trainNet(int maxNoOfErrorUpdateCycles, float minErrorChange,
                             int outPatience, int StartPattern, int EndPattern, 
			     float epsilon, float mu,float decay, 
			     float **ParameterOutArray,int *NoOfOutParams)
{
    int m,p,pat,sub,counter=0;
    float oldNetError;
    static float OutParameter[1];
    int start, end;

    *NoOfOutParams = 1;
    *ParameterOutArray = OutParameter;

    SumSqError = 0.0;
    cc_initOutputUnits();

    /* compute the necessary sub patterns */

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    if(KernelErrorCode != KRERR_NO_ERROR)
        return;

    /* give oldNetError a meaningful initial value */

    NET_ERROR(OutParameter)=FLOAT_MAX;
    do {
	oldNetError = NET_ERROR(OutParameter);
	for(m=0;m<outPatience;m++) {
	    NET_ERROR(OutParameter) = 0.0;
	    SumSqError = 0.0;
	    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
	    start = kr_AbsPosOfFirstSubPat(StartPattern);
	    end   = kr_AbsPosOfFirstSubPat(EndPattern);
	    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
	    for(p=start; p<=end;p++){
		kr_getSubPatternByNo(&pat,&sub,p);
		rcc_propagateNetForward(pat,sub);
		NET_ERROR(OutParameter) +=rcc_QPO_propagateNetBackward(pat,sub);
	    }
	    rcc_QPO_updateNet(epsilon,mu,decay);
	    if(cc_printOnOff) {
		printf("Epoch: %d NetError: %f \n",++counter,
		       NET_ERROR(OutParameter));
	    }
	    if((maxNoOfErrorUpdateCycles--) == 0) {
		return;
	    }
	}
    } while(fabs(oldNetError-NET_ERROR(OutParameter)) >= 
	    (minErrorChange * oldNetError));
}



/*****************************************************************************
  FUNCTION : rcc_QPO_propagateNetBackward

  PURPOSE  : Calculate the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static  float rcc_QPO_propagateNetBackward(int PatternNo, int sub_pat_no)
{
    struct Link   *LinkPtr;
    struct Site   *site_ptr;
    struct Unit   *OutputUnitPtr;
    Patterns      out_pat;
    float         error,sum_error,devit;
    int dummy;

    sum_error = 0.0;    
    out_pat = kr_getSubPatData(PatternNo,sub_pat_no,OUTPUT,NULL);

    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,dummy){
	devit =  OutputUnitPtr->Out.output - *(out_pat++);

	sum_error += devit * devit;
	error = devit * ((*OutputUnitPtr->act_deriv_func)(OutputUnitPtr) + 0.1);
	SumSqError += error * error;
 
	BIAS_CURRENT_SLOPE(OutputUnitPtr) += error;

	if (UNIT_HAS_DIRECT_INPUTS(OutputUnitPtr)) {  
	    FOR_ALL_LINKS(OutputUnitPtr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}else {
	    FOR_ALL_SITES_AND_LINKS(OutputUnitPtr,site_ptr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}
    }
    return(sum_error); 
}



/*****************************************************************************
  FUNCTION : rcc_QPO_updateNet

  PURPOSE  : Update the weights of the output units with quickprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_QPO_updateNet(float epsilon, float mu, float decay)
{
    struct Unit *OutputUnitPtr;
    struct Link *LinkPtr;
    float shrinkFactor=mu/(mu+1);
    float bias_previousSlope,bias_currentSlope,bias_lastWeightChange,
          bias_weightChange;
    float ln_previousSlope,ln_currentSlope,ln_lastWeightChange,
          ln_weightChange;
    int o;

    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,o) {
	bias_previousSlope = BIAS_PREVIOUS_SLOPE(OutputUnitPtr);
	bias_currentSlope = 
	    BIAS_CURRENT_SLOPE(OutputUnitPtr) + decay * OutputUnitPtr->bias;
	bias_lastWeightChange = BIAS_LAST_WEIGHT_CHANGE(OutputUnitPtr);

	bias_weightChange = 0.0;

	if(bias_previousSlope > 0.0) {
	    if(bias_currentSlope > 0.0) {
		bias_weightChange -= epsilon * bias_currentSlope;
	    }
	    if(bias_currentSlope >= (shrinkFactor * bias_previousSlope)) {
		bias_weightChange += mu * bias_lastWeightChange;
	    }else {
		bias_weightChange += (bias_lastWeightChange*bias_currentSlope)/ 
		    (bias_previousSlope - bias_currentSlope); 
	    }
	}else if(bias_previousSlope < 0.0) {
	    if(bias_currentSlope < 0.0) {
		bias_weightChange -= epsilon * bias_currentSlope;
	    }
	    if(bias_currentSlope <= (shrinkFactor * bias_previousSlope)) {
		bias_weightChange += mu * bias_lastWeightChange;
	    }else {
		bias_weightChange += (bias_lastWeightChange*bias_currentSlope)/ 
		    (bias_previousSlope - bias_currentSlope); 
	    }
	}else {
	    bias_weightChange -= epsilon * bias_currentSlope;
	}
	OutputUnitPtr->bias += 
	    (BIAS_LAST_WEIGHT_CHANGE(OutputUnitPtr) = bias_weightChange); 
	BIAS_PREVIOUS_SLOPE(OutputUnitPtr) = bias_currentSlope;
	BIAS_CURRENT_SLOPE(OutputUnitPtr) = 0.0;

	FOR_ALL_LINKS(OutputUnitPtr,LinkPtr) {
   
	    ln_weightChange = 0.0;
	    ln_previousSlope = LN_PREVIOUS_SLOPE(LinkPtr);
	    ln_currentSlope = LN_CURRENT_SLOPE(LinkPtr) + decay*LinkPtr->weight;
	    ln_lastWeightChange = LN_LAST_WEIGHT_CHANGE(LinkPtr);
   
	    if(ln_previousSlope > 0.0) {
		if(ln_currentSlope > 0.0) {
		    ln_weightChange -= epsilon * ln_currentSlope;
		}
		if(ln_currentSlope >= (shrinkFactor * ln_previousSlope)) {
		    ln_weightChange += mu * ln_lastWeightChange;
		}else {
		    ln_weightChange += (ln_lastWeightChange * ln_currentSlope) / 
			(ln_previousSlope - ln_currentSlope); 
		}
	    }else if(ln_previousSlope < 0.0) {
		if(ln_currentSlope < 0.0) {
		    ln_weightChange -= epsilon * ln_currentSlope;
		}
		if(ln_currentSlope <= (shrinkFactor * ln_previousSlope)) {
		    ln_weightChange += mu * ln_lastWeightChange;
		}else {
		    ln_weightChange += (ln_lastWeightChange*ln_currentSlope) / 
			(ln_previousSlope - ln_currentSlope); 
		}
	    }else {
		ln_weightChange -= epsilon * ln_currentSlope;
	    }
   
	    LinkPtr->weight += 
		LN_LAST_WEIGHT_CHANGE(LinkPtr) = 
		    ln_weightChange; 
	    LN_PREVIOUS_SLOPE(LinkPtr)   = ln_currentSlope;
	    LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	}
    }
}



/*****************************************************************************
  FUNCTION : rcc_QPS_trainNet

  PURPOSE  : Maximize the covariance between the special units and the
             output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_QPS_trainNet(int maxNoOfCovarianceUpdateCycles, 
			     float minCovarianceChange, int specialPatience, 
			     int StartPattern, int EndPattern,
			     float epsilon, float mu, float decay, 
			     int MaxSpecialUnitNo)
{
    int m,counter=0;
    float oldHighScore,newHighScore=0.0;

    cc_initErrorArrays();
    rcc_calculateOutputUnitError(StartPattern,EndPattern);

    do {
	oldHighScore = newHighScore;
	for(m=0;m<specialPatience;m++) { 
	    counter++;
	    rcc_calculateSpecialUnitActivation(StartPattern,EndPattern);
	    newHighScore = 
		rcc_QPS_propagateNetBackward(StartPattern,EndPattern,counter);
	    rcc_QPS_updateNet(epsilon,mu,decay);
	    cc_initActivationArrays(); 
	    if((maxNoOfCovarianceUpdateCycles--) == 0) {
		return;
	    }
	}
    } while(fabs(newHighScore-oldHighScore) >= 
	    (minCovarianceChange * oldHighScore));
}



/*****************************************************************************
  FUNCTION : rcc_QPS_propagateNetBackward

  PURPOSE  : Calculate the special unit with the maximum covariance and return 
             it.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static float rcc_QPS_propagateNetBackward(int StartPattern, int EndPattern, 
					  int counter)
{
    float change=0.0,bestSpecialUnitScore,actPrime,recurrentLinkWeight,dsum;
    int s,o,p,n,h;
    int start, end;
    struct Unit *SpecialUnitPtr,*OutputUnitPtr,*hiddenUnitPtr;
    struct Link *LinkPtr,*recurrentLinkPtr;

 
    bestSpecialUnitScore = 
	cc_calculateCorrelation(StartPattern,EndPattern,counter); 

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    n = end - start + 1;
    for(p=start; p<=end;p++){

	cc_initInputUnitsWithPattern(p);
	FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,h) { 
	    hiddenUnitPtr->lln = reset[p];
	    hiddenUnitPtr->act = (*hiddenUnitPtr->act_func)(hiddenUnitPtr);
	    if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
		hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	    }else {
		hiddenUnitPtr->Out.output = 
		    (*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	    }
	}
   
	FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,s) {
	    change = 0.0;
	    SpecialUnitPtr->act = SpecialUnitAct[p][s];
	    actPrime = (*SpecialUnitPtr->act_deriv_func)(SpecialUnitPtr);
	    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,o) {
		change -= CorBetweenSpecialActAndOutError[s][o] *
		    ((OutputUnitError[p][o] - OutputUnitSumError[o]/n)/
		     SumSqError);
	    }

	    counter = 0;
	    GET_RECURRENT_LINK(SpecialUnitPtr,recurrentLinkPtr);
	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    recurrentLinkWeight = recurrentLinkPtr->weight; 
	    if(!reset[p]) {	/* calculate slope of the recurrent link */
		dsum = actPrime * (SpecialUnitAct[p-1][s] + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(recurrentLinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum; 
	    }

	    FOR_ALL_NOT_RECURRENT_LINKS(SpecialUnitPtr,LinkPtr) {
		if(reset[p]) {
		    linkArray[s][counter] = 0.0;
		}
		dsum = actPrime * (LinkPtr->to->Out.output + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(LinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum;
	    }

	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    dsum = actPrime * (1 + (recurrentLinkWeight*linkArray[s][counter]));
	    BIAS_CURRENT_SLOPE(SpecialUnitPtr) += change * dsum;
	    linkArray[s][counter] = dsum;
	}
    }
    return(bestSpecialUnitScore);
}



/*****************************************************************************
  FUNCTION : rcc_QPS_updateNet

  PURPOSE  : Update the weights of the special units with quickprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_QPS_updateNet(float epsilon, float mu, float decay)
{
    struct Unit *SpecialUnitPtr;
    struct Link *LinkPtr;
    float shrinkFactor=mu/(mu+1);
    float bias_previousSlope,bias_currentSlope,bias_lastWeightChange,
          bias_weightChange;
    float ln_previousSlope,ln_currentSlope,ln_lastWeightChange,ln_weightChange;
    int o;

    FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,o) {
	bias_previousSlope = BIAS_PREVIOUS_SLOPE(SpecialUnitPtr);
	bias_currentSlope = BIAS_CURRENT_SLOPE(SpecialUnitPtr) + 
	                    decay * SpecialUnitPtr->bias;
	bias_lastWeightChange = BIAS_LAST_WEIGHT_CHANGE(SpecialUnitPtr);

	bias_weightChange = 0.0;

	if(bias_previousSlope > 0.0) {
	    if(bias_currentSlope > 0.0) {
		bias_weightChange -= epsilon * bias_currentSlope;
	    }
	    if(bias_currentSlope >= (shrinkFactor * bias_previousSlope)) {
		bias_weightChange += mu * bias_lastWeightChange;
	    }else {
		bias_weightChange += (bias_lastWeightChange*bias_currentSlope)/ 
		    (bias_previousSlope - bias_currentSlope); 
	    }
	}else if(bias_previousSlope < 0.0) {
	    if(bias_currentSlope < 0.0) {
		bias_weightChange -= epsilon * bias_currentSlope;
	    }
	    if(bias_currentSlope <= (shrinkFactor * bias_previousSlope)) {
		bias_weightChange += mu * bias_lastWeightChange;
	    }else {
		bias_weightChange += (bias_lastWeightChange*bias_currentSlope)/ 
		    (bias_previousSlope - bias_currentSlope); 
	    }
	}else {
	    bias_weightChange -= epsilon * bias_currentSlope;
	}
	SpecialUnitPtr->bias += 
	    (BIAS_LAST_WEIGHT_CHANGE(SpecialUnitPtr) = bias_weightChange); 
	BIAS_PREVIOUS_SLOPE(SpecialUnitPtr) = bias_currentSlope;
	BIAS_CURRENT_SLOPE(SpecialUnitPtr) = 0.0;

	FOR_ALL_LINKS(SpecialUnitPtr,LinkPtr) {
   
	    ln_weightChange = 0.0;
	    ln_previousSlope = LN_PREVIOUS_SLOPE(LinkPtr);
	    ln_currentSlope = LN_CURRENT_SLOPE(LinkPtr) + decay*LinkPtr->weight;
	    ln_lastWeightChange = LN_LAST_WEIGHT_CHANGE(LinkPtr);
   
	    if(ln_previousSlope > 0.0) {
		if(ln_currentSlope > 0.0) {
		    ln_weightChange -= epsilon * ln_currentSlope;
		}
		if(ln_currentSlope >= (shrinkFactor * ln_previousSlope)) {
		    ln_weightChange += mu * ln_lastWeightChange;
		}else {
		    ln_weightChange += (ln_lastWeightChange * ln_currentSlope) / 
			(ln_previousSlope - ln_currentSlope); 
		}
	    }else if(ln_previousSlope < 0.0) {
		if(ln_currentSlope < 0.0) {
		    ln_weightChange -= epsilon * ln_currentSlope;
		}
		if(ln_currentSlope <= (shrinkFactor * ln_previousSlope)) {
		    ln_weightChange += mu * ln_lastWeightChange;
		}else {
		    ln_weightChange += (ln_lastWeightChange * ln_currentSlope) / 
			(ln_previousSlope - ln_currentSlope); 
		}
	    }else {
		ln_weightChange -= epsilon * ln_currentSlope;
	    }
   
	    LinkPtr->weight += LN_LAST_WEIGHT_CHANGE(LinkPtr) = ln_weightChange; 
	    LN_PREVIOUS_SLOPE(LinkPtr)   = ln_currentSlope;
	    LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	}
    } 
}


/******************* end quickprop routines *************************/


/************* begin Backprop routines *********************/

/*****************************************************************************
  FUNCTION : rcc_BPO_trainNet

  PURPOSE  : Minimize the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_BPO_trainNet(int maxNoOfErrorUpdateCycles, float minErrorChange,
                             int outPatience, int StartPattern, int EndPattern, 
			     float eta, float mu, float fse, 
			     float **ParameterOutArray, int *NoOfOutParams)
{
    int m,p,counter=0;
    float oldNetError;
    static float OutParameter[1];
    float dummy=0.0;
    int start,end;
    int pat,sub;

    *NoOfOutParams = 1;
    *ParameterOutArray = OutParameter;

    SumSqError = 0.0;
    cc_initOutputUnits();

    /* compute the necessary sub patterns */
    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    if(KernelErrorCode != KRERR_NO_ERROR)
        return;

    /* give oldNetError a meaningful initial value */
    NET_ERROR(OutParameter)=FLOAT_MAX;
    do {
	oldNetError = NET_ERROR(OutParameter);
	for(m=0;m<outPatience;m++) {
	    NET_ERROR(OutParameter) = 0.0;
	    SumSqError = 0.0;
	    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
	    start = kr_AbsPosOfFirstSubPat(StartPattern);
	    end   = kr_AbsPosOfFirstSubPat(EndPattern);
	    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
	    for(p=start; p<=end;p++){
		kr_getSubPatternByNo(&pat,&sub,p);
		rcc_propagateNetForward(pat,sub);
		NET_ERROR(OutParameter) += 
		    rcc_BPO_propagateNetBackward(pat,sub,fse);
	    }
	    rcc_BPO_updateNet(eta,mu,dummy);
	    if(cc_printOnOff){
		printf("Epoch: %d NetError: %f \n",++counter,
		       NET_ERROR(OutParameter));
	    }
	    if((maxNoOfErrorUpdateCycles--) == 0) {
		return;
	    }
	}
    } while(fabs(oldNetError-NET_ERROR(OutParameter)) >= 
	    (minErrorChange * oldNetError));
}



/*****************************************************************************
  FUNCTION : rcc_BPO_propagateNetBackward

  PURPOSE  : Calculate the error of the output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static  float rcc_BPO_propagateNetBackward(int PatternNo, int sub_pat_no,
					   float fse)
{
    struct Link   *LinkPtr;
    struct Site   *site_ptr;
    struct Unit   *OutputUnitPtr;
    Patterns      out_pat;
    float         error,sum_error,devit;
    int dummy=0.0;

    sum_error = 0.0;    
    out_pat = kr_getSubPatData(PatternNo,sub_pat_no,OUTPUT,NULL);

    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,dummy){
	devit =  OutputUnitPtr->Out.output - *(out_pat++);

	sum_error += devit * devit;
	error = devit * ((*OutputUnitPtr->act_deriv_func)(OutputUnitPtr) + fse);
	SumSqError += error * error;

	BIAS_CURRENT_SLOPE(OutputUnitPtr) += error;

	if (UNIT_HAS_DIRECT_INPUTS(OutputUnitPtr)) {  
	    FOR_ALL_LINKS(OutputUnitPtr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}else {
	    FOR_ALL_SITES_AND_LINKS(OutputUnitPtr,site_ptr,LinkPtr) {
		LN_CURRENT_SLOPE(LinkPtr) += error * LinkPtr->to->Out.output;
	    }
	}
    }
    return(sum_error); 
}



/*****************************************************************************
  FUNCTION : rcc_BPO_updateNet

  PURPOSE  : Update the weights of the output units with backprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_BPO_updateNet(float eta, float mu, float dummy)
{
    struct Unit *outputUnitPtr;
    struct Link *LinkPtr;
    float bias_weightChange;
    float ln_weightChange;
    int s;

    FOR_ALL_OUTPUT_UNITS(outputUnitPtr,s) {
	bias_weightChange = BIAS_CURRENT_SLOPE(outputUnitPtr) * eta + 
	                    BIAS_PREVIOUS_SLOPE(outputUnitPtr) * mu;
	outputUnitPtr->bias -= bias_weightChange; 
	BIAS_PREVIOUS_SLOPE(outputUnitPtr) = BIAS_CURRENT_SLOPE(outputUnitPtr);
	BIAS_CURRENT_SLOPE(outputUnitPtr) = 0.0;

	FOR_ALL_LINKS(outputUnitPtr,LinkPtr) {
	    ln_weightChange = LN_CURRENT_SLOPE(LinkPtr) * eta + 
		              LN_PREVIOUS_SLOPE(outputUnitPtr) * mu;
	    LinkPtr->weight -= ln_weightChange; 
	    LN_PREVIOUS_SLOPE(outputUnitPtr) = LN_CURRENT_SLOPE(outputUnitPtr);
	    LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	}
    }
}



/*****************************************************************************
  FUNCTION : rcc_BPS_trainNet

  PURPOSE  : Maximize the covariance between the special units and the 
             output units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_BPS_trainNet(int maxNoOfCovarianceUpdateCycles, 
			     float minCovarianceChange, int specialPatience, 
			     int StartPattern, int EndPattern, float eta, 
			     float mu, float dummy1, int MaxSpecialUnitNo)
{
    int m,counter=0;
    float oldHighScore,newHighScore=0.0;

    cc_initErrorArrays();
    rcc_calculateOutputUnitError(StartPattern,EndPattern);

    do {
	oldHighScore = newHighScore;
	for(m=0;m<specialPatience;m++) { 
	    counter++;
	    rcc_calculateSpecialUnitActivation(StartPattern,EndPattern);
	    newHighScore = rcc_BPS_propagateNetBackward(StartPattern,EndPattern,
							counter);
	    rcc_BPS_updateNet(eta,mu,dummy1);
	    cc_initActivationArrays(); 
	    if((maxNoOfCovarianceUpdateCycles--) == 0) {
		return;
	    }    
	}
    } while(fabs(newHighScore-oldHighScore) >= 
	    (minCovarianceChange * oldHighScore));
}



/*****************************************************************************
  FUNCTION : rcc_BPS_propagateNetBackward

  PURPOSE  : Calculate the special unit with the maximum covariance and return 
             it.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static float rcc_BPS_propagateNetBackward(int StartPattern, int EndPattern, 
					  int counter)
{
    float change=0.0,bestSpecialUnitScore,actPrime,recurrentLinkWeight,dsum;
    int s,o,p,n,h;
    struct Unit *SpecialUnitPtr,*OutputUnitPtr,*hiddenUnitPtr;
    struct Link *LinkPtr,*recurrentLinkPtr;
    int start,end;

    bestSpecialUnitScore = 
	cc_calculateCorrelation(StartPattern,EndPattern,counter); 

    KernelErrorCode = kr_initSubPatternOrder(StartPattern,EndPattern);
    start = kr_AbsPosOfFirstSubPat(StartPattern);
    end   = kr_AbsPosOfFirstSubPat(EndPattern);
    end  += kr_NoOfSubPatPairs(EndPattern) - 1;
    n = end - start +1;
    for(p=start; p<=end;p++){

	cc_initInputUnitsWithPattern(p);
	FOR_ALL_HIDDEN_UNITS(hiddenUnitPtr,h) { 
	    hiddenUnitPtr->lln = reset[p];
	    hiddenUnitPtr->act = (*hiddenUnitPtr->act_func)(hiddenUnitPtr);
	    if(hiddenUnitPtr->out_func == OUT_IDENTITY) {
		hiddenUnitPtr->Out.output = hiddenUnitPtr->act;
	    }else {
		hiddenUnitPtr->Out.output = 
		    (*hiddenUnitPtr->out_func) (hiddenUnitPtr->act);
	    }
	}
   
	FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,s) {
	    change = 0.0;
	    SpecialUnitPtr->act = SpecialUnitAct[p][s];
	    actPrime = (*SpecialUnitPtr->act_deriv_func)(SpecialUnitPtr);
	    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,o) {
		change -= CorBetweenSpecialActAndOutError[s][o] *
		   ((OutputUnitError[p][o]-OutputUnitSumError[o]/n)/SumSqError);
	    }

	    counter = 0;
	    GET_RECURRENT_LINK(SpecialUnitPtr,recurrentLinkPtr);
	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    recurrentLinkWeight = recurrentLinkPtr->weight; 
	    if(!reset[p]) {	/* calculate slope of the recurrent link */
		dsum = actPrime * (SpecialUnitAct[p-1][s] + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(recurrentLinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum; 
	    }

	    FOR_ALL_NOT_RECURRENT_LINKS(SpecialUnitPtr,LinkPtr) {
		if(reset[p]) {
		    linkArray[s][counter] = 0.0;
		}
		dsum = actPrime * (LinkPtr->to->Out.output + 
				   (recurrentLinkWeight*linkArray[s][counter]));
		LN_CURRENT_SLOPE(LinkPtr) += change * dsum;
		linkArray[s][counter++] = dsum;
	    }

	    if(reset[p]) {
		linkArray[s][counter] = 0.0;
	    }
	    dsum = actPrime * (1 + (recurrentLinkWeight*linkArray[s][counter]));
	    BIAS_CURRENT_SLOPE(SpecialUnitPtr) += change * dsum;
	    linkArray[s][counter] = dsum;
	}
    }
    return(bestSpecialUnitScore);
}



/*****************************************************************************
  FUNCTION : rcc_BPS_updateNet

  PURPOSE  : Update the weights of the special units with backprop.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static void rcc_BPS_updateNet(float eta, float mu, float dummy)
{
    struct Unit *specialUnitPtr;
    struct Link *LinkPtr;
    float bias_weightChange;
    float ln_weightChange;
    int s;

    FOR_ALL_SPECIAL_UNITS(specialUnitPtr,s) {
	bias_weightChange = BIAS_CURRENT_SLOPE(specialUnitPtr) * eta + 
	                    BIAS_PREVIOUS_SLOPE(specialUnitPtr) * mu;
	specialUnitPtr->bias -= bias_weightChange; 
	BIAS_PREVIOUS_SLOPE(specialUnitPtr) =BIAS_CURRENT_SLOPE(specialUnitPtr);
	BIAS_CURRENT_SLOPE(specialUnitPtr) = 0.0;

	FOR_ALL_LINKS(specialUnitPtr,LinkPtr) {
	    ln_weightChange = LN_CURRENT_SLOPE(LinkPtr) * eta + 
		              LN_PREVIOUS_SLOPE(specialUnitPtr) * mu;
	    LinkPtr->weight -= ln_weightChange; 
	    LN_PREVIOUS_SLOPE(specialUnitPtr) =LN_CURRENT_SLOPE(specialUnitPtr);
	    LN_CURRENT_SLOPE(LinkPtr)    = 0.0;
	}
    }
}


/******************* end packprop routines *************************/



/*****************************************************************************
  FUNCTION : rcc_generateSpecialUnits

  PURPOSE  : Generates the special units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static krui_err rcc_generateSpecialUnits(int MaxSpecialUnitNo, int OutputXMax, 
					 int type)
{
    int i,selector;
    struct Unit *UnitPtr;
    int CurrentUnit;

    for(i=0;i<MaxSpecialUnitNo;i++) {
	if(type==RANDOM){
	    selector = i % (NO_OF_ACT_FUNCS - 1);
	}else {
	    selector = type;
	}
	KernelErrorCode = 
	    kr_unitSetTType(CurrentUnit=kr_makeDefaultUnit(),SPECIAL); 
	ERROR_CHECK;

	KernelErrorCode = 
	    krui_setUnitActFunc(CurrentUnit,rcc_actFuncArray[selector]);
	ERROR_CHECK;  

	UnitPtr = kr_getUnitPtr(CurrentUnit); 
	ERROR_CHECK;

	SET_UNIT_XPOS(UnitPtr,OutputXMax+3);
	SET_UNIT_YPOS(UnitPtr,2+i);
	KernelErrorCode = krui_setCurrentUnit(CurrentUnit); 
	ERROR_CHECK;

	/* links between special units and input units */
	FOR_ALL_UNITS(UnitPtr){
	    if((IS_INPUT_UNIT(UnitPtr) || IS_HIDDEN_UNIT(UnitPtr)) && 
	       UNIT_IN_USE(UnitPtr)) {
		KernelErrorCode = 
		    krui_createLink(GET_UNIT_NO(UnitPtr),
				    cc_generateRandomNo(RCC_MAX_VALUE)); 
		ERROR_CHECK;
	    }
	}
	KernelErrorCode = krui_createLink(CurrentUnit,
					  cc_generateRandomNo(RCC_MAX_VALUE));
	ERROR_CHECK;
    }  
    return(KRERR_NO_ERROR);
}



/*****************************************************************************
  FUNCTION : rcc_generateHiddenUnit 

  PURPOSE  : Generates the hidden units.
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
static krui_err rcc_generateHiddenUnit(int maxYPosOfHiddenUnit, 
				       int *xPosOfLastInsertedHiddenUnit,
                                       int *yPosOfLastInsertedHiddenUnit)
{
    int CurrentUnit,dummy,NewHiddenUnit,s;
    struct Unit *HiddenUnitPtr,*SpecialUnitPtr,*OutputUnitPtr;
    FlintType linkValue=0.0;


    FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,s){
	KernelErrorCode = 
	    krui_setCurrentUnit(CurrentUnit = GET_UNIT_NO(SpecialUnitPtr)); 
	ERROR_CHECK;

	if(krui_isConnected(CurrentUnit)) {
	    if(SpecialUnitPtr == bestSpecialUnitPtr) {
		linkValue = ((struct Link *)(SpecialUnitPtr->sites))->weight;
	    }
	    KernelErrorCode = krui_deleteLink(); /* delete recurent link */
	    ERROR_CHECK;
	}else {
	    KernelErrorCode = KRERR_CC_ERROR5;
	    ERROR_CHECK;
	} 
    }

    CurrentUnit = 
	KernelErrorCode = 
	    kr_copyUnit(ONLY_INPUTS,GET_UNIT_NO(bestSpecialUnitPtr));
    if(KernelErrorCode < 0) {
	ERROR_CHECK;
    }
    KernelErrorCode = KRERR_NO_ERROR;

    KernelErrorCode = kr_unitSetTType(CurrentUnit,HIDDEN); 
    ERROR_CHECK;

    HiddenUnitPtr = kr_getUnitPtr(CurrentUnit); 
    ERROR_CHECK;

    KernelErrorCode = 
	cc_setHiddenUnit(HiddenUnitPtr,maxYPosOfHiddenUnit,
			 xPosOfLastInsertedHiddenUnit,
			 yPosOfLastInsertedHiddenUnit);
    ERROR_CHECK;

    KernelErrorCode = krui_setCurrentUnit(CurrentUnit); 
    ERROR_CHECK;

    /* generate recurent link */
    KernelErrorCode = krui_createLink(CurrentUnit,linkValue); 
    ERROR_CHECK;

    NewHiddenUnit = CurrentUnit;
    /* generate Links between output unit and new hidden unit */
    FOR_ALL_OUTPUT_UNITS(OutputUnitPtr,dummy){
	CurrentUnit = GET_UNIT_NO(OutputUnitPtr);
	KernelErrorCode = krui_setCurrentUnit(CurrentUnit); 
	ERROR_CHECK;

	KernelErrorCode = krui_createLink(NewHiddenUnit,0.0); 
	ERROR_CHECK;
    }
 
    FOR_ALL_SPECIAL_UNITS(SpecialUnitPtr,s){
	KernelErrorCode = 
	    krui_setCurrentUnit(CurrentUnit = GET_UNIT_NO(SpecialUnitPtr)); 
	ERROR_CHECK;
	
	KernelErrorCode = 
	    krui_createLink(NewHiddenUnit,cc_generateRandomNo(RCC_MAX_VALUE));
	ERROR_CHECK;

	/* generate recurent link */
	KernelErrorCode = 
	    krui_createLink(CurrentUnit,cc_generateRandomNo(RCC_MAX_VALUE)); 
	ERROR_CHECK;
    }

    KernelErrorCode = kr_topoSort(TOPOLOGICAL_RCC);
    ERROR_CHECK;    

    if(CC_TEST) {
	cc_printUnitArray();
	cc_printTopoPtrArray();
    }

    KernelErrorCode = cc_setPointers();
    ERROR_CHECK;

    NetModified = FALSE;
    return(KRERR_NO_ERROR);
}



/*****************************************************************************
  FUNCTION : LEARN_RecCasCor 

  PURPOSE  : The main routine of RCC
  NOTES    :

  UPDATE   : 5.2.93
******************************************************************************/
krui_err LEARN_RecCasCor(int StartPattern, int EndPattern, 
			 float *ParameterInArray, int NoOfInParams, 
			 float **ParameterOutArray, int *NoOfOutParams)
{
    static int   OldMaxSpecialUnitNo=0,MaxSpecialUnitNo;
    static int   maxNoOfErrorUpdateCycles,maxNoOfCovarianceUpdateCycles,
                 outPatience,specialPatience;
    static int   maxYPosOfHiddenUnit,xPosOfLastInsertedHiddenUnit,
                 yPosOfLastInsertedHiddenUnit;
    static int   specialFuncType,oldSpecialFuncType,outputXMax,learnFunc;
    static float maxPixelError,minErrorChange,minCovarianceChange,
                 param1,param2,param3,param4,param5,param6;
  
    srand48((long)time((long *)0)); 

    if(strcmp(krui_getUpdateFunc(),"RCC_Order")){
	return(KRERR_CC_ERROR10);
    }

    /* Enter this path only  if "ALL" in the remote pannel was pressed */
    if(cc_allButtonIsPressed == 1) { 

	param1 = ParameterInArray[0];
	param2 = ParameterInArray[1];
	param3 = ParameterInArray[2];

	param4 = ParameterInArray[3];
	param5 = ParameterInArray[4];
	param6 = 0.001;
	/*  param6 = ParameterInArray[5]; */
  
	/* cc_data.GLOBAL.pixelError */
	maxPixelError = ParameterInArray[6];

	/* cc_data.GLOBAL.learningFunc */
	learnFunc     = (int)ParameterInArray[7]; 

	/* cc_data.GLOBAL.onOff */ 
	cc_printOnOff = (int)ParameterInArray[8]; 

	/* cc_data.CAND.covarianceChange */
	minCovarianceChange = ParameterInArray[9];

	/* cc_data.CAND.candidatePatience */ 
	specialPatience = ParameterInArray[10]; 

	/* cc_data.CAND.maxNoOfUpdateCycles */
	maxNoOfCovarianceUpdateCycles = ParameterInArray[11]; 

	/* cc_data.CAND.maxNoOfCandUnits */
	MaxSpecialUnitNo = ParameterInArray[12]; 

	/* cc_data.CAND.actFunc */
	specialFuncType = (int)ParameterInArray[13]; 
        
	/* cc_data.OUT.errorChange */
	minErrorChange = ParameterInArray[14];

	/* cc_data.OUT.outputPatience */ 
	outPatience = ParameterInArray[15]; 

	/* cc_data.OUT.maxNoOfUpdateCycles */
	maxNoOfErrorUpdateCycles = ParameterInArray[16]; 

	cc_end = 0;
	cc_cascade = 1;

	cc_compareActFunctions(specialFuncType,RCC);
	rcc_manageResetArray(StartPattern,EndPattern,1);
	rcc_manageLinkArray(MaxSpecialUnitNo,1);

	switch(learnFunc) {
	  case BACKPROP:
	    cc_trainOutputUnits  = rcc_BPO_trainNet;
	    cc_trainSpecialUnits = rcc_BPS_trainNet; 
	    break;
	  case QUICKPROP:
	    cc_trainOutputUnits  = rcc_QPO_trainNet;
	    cc_trainSpecialUnits = rcc_QPS_trainNet;
	    break;
	  case RPROP: 
	    cc_trainOutputUnits  = rcc_RPO_trainNet;
	    cc_trainSpecialUnits = rcc_RPS_trainNet; 
	    break;
	  default: CC_ERROR(KRERR_CC_ERROR3);
	}

	maxYPosOfHiddenUnit =
	    xPosOfLastInsertedHiddenUnit =
		yPosOfLastInsertedHiddenUnit =
		    outputXMax=0;
	KernelErrorCode = 
	    cc_calculateNetParameters(&maxYPosOfHiddenUnit,
				      &xPosOfLastInsertedHiddenUnit,
				      &yPosOfLastInsertedHiddenUnit,
				      &outputXMax);
	ERROR_CHECK;
    }

    if(cc_printOnOff){
	if(strcmp(krui_getUpdateFunc(),"RCC_Order")){
	    return(KRERR_CC_ERROR10);
	}
	if(strcmp(krui_getInitialisationFunc(),"RCC_Weights")){
	    return(KRERR_CC_ERROR11);
	}    
    }


  
    if(cc_end){
	return(KRERR_NO_ERROR);
    }

    if(NetModified || (TopoSortID!=TOPOLOGICAL_RCC) || 
       (LearnFuncHasChanged) || (OldMaxSpecialUnitNo!=MaxSpecialUnitNo) || 
       (cc_update)) {

	OldMaxSpecialUnitNo = MaxSpecialUnitNo;
	oldSpecialFuncType  = specialFuncType;

	KernelErrorCode = rcc_searchRecurrentLinks();
	ERROR_CHECK;

	KernelErrorCode = cc_deleteAllSpecialUnits();
	ERROR_CHECK; 

	KernelErrorCode = rcc_generateSpecialUnits(MaxSpecialUnitNo,outputXMax,
						   specialFuncType);
	ERROR_CHECK;

	KernelErrorCode = kr_topoSort(TOPOLOGICAL_RCC);
	ERROR_CHECK;    

	if(CC_TEST) {
	    cc_printUnitArray();
	    cc_printTopoPtrArray();
	}

	KernelErrorCode = cc_setPointers();
	ERROR_CHECK;

	cc_update = 0;
	LearnFuncHasChanged = 0;
	NetModified = 0;
    }

    if(cc_allButtonIsPressed == 1) {
	KernelErrorCode = cc_freeStorage(StartPattern,EndPattern,0);
	ERROR_CHECK;

	KernelErrorCode = 
	    cc_allocateStorage(StartPattern,EndPattern,MaxSpecialUnitNo);
	ERROR_CHECK;
    }

    if(oldSpecialFuncType != specialFuncType){
	oldSpecialFuncType = specialFuncType;
	cc_changeActFuncOfSpecialUnits(specialFuncType,RCC);
    }

    KernelErrorCode = cc_initSpecialUnitLinks();
    ERROR_CHECK;

    /* For safety, this error should only appear, if someone has changed 
       the program in a wrong way!!!
       */
    if(cc_storageFree){ 
	CC_ERROR(KRERR_CC_ERROR2);    
    }
    if(cc_allButtonIsPressed == 1) {
	cc_allButtonIsPressed = 0;
	if((outPatience != 0) && (maxNoOfErrorUpdateCycles != 0)) {
	    (*cc_trainOutputUnits)(maxNoOfErrorUpdateCycles,minErrorChange,
				   outPatience,StartPattern,EndPattern,param1,
				   param2,param3,ParameterOutArray,
				   NoOfOutParams);
	}
    }

    if(rcc_test(StartPattern,EndPattern,maxPixelError) ==  CONTINUE_LEARNING) {
	if((specialPatience != 0) && (maxNoOfCovarianceUpdateCycles != 0)) {
	    (*cc_trainSpecialUnits)(maxNoOfCovarianceUpdateCycles,
				    minCovarianceChange,specialPatience,
				    StartPattern,EndPattern,param4,
				    param5,param6,MaxSpecialUnitNo);
	}
    }else {			/* Stop learning */
	cc_end = 1;
	return(KRERR_NO_ERROR);
    }

    if((specialPatience != 0) && (maxNoOfCovarianceUpdateCycles != 0)) {
	KernelErrorCode = 
	    rcc_generateHiddenUnit(maxYPosOfHiddenUnit,
				   &xPosOfLastInsertedHiddenUnit,
				   &yPosOfLastInsertedHiddenUnit);
	ERROR_CHECK;
    }

    if((outPatience != 0) && (maxNoOfErrorUpdateCycles != 0)) {
	(*cc_trainOutputUnits)(maxNoOfErrorUpdateCycles,minErrorChange,
			       outPatience,StartPattern,EndPattern,param1,
			       param2,param3,ParameterOutArray,NoOfOutParams); 
    }
    return(KRERR_NO_ERROR);
}



