#import <objc/Object.h>
#import "ART_LTM.h"
#import "ART_Func.h"

#define LTM_NODE_RESET 0
#define LTM_NODE_AVAIL 1

#define ART_READ(x) NXRead(stream,&x,sizeof(x))
#define ART_WRITE(x) NXWrite(stream,&x,sizeof(x))

#define ART_WRITE_ASCII(x) NXPrintf(stream,"%7.5f ",x)
#define ART_WRITE_SLASH NXPrintf(stream,"\n")




@implementation ART_LTM
/*
     NXZone *zone;        // Next Heap Zone

     int   M,N;           // field dimensions
     float **BU,**TD,     // learning arrays
           *oldTD,*oldBU,
           d,h;           // learning rate and integation stepsize
     int   *Avail,        // flag array for notation of reset nodes
           estab=1;       // acount of established nodes

*/

+ new
{
     self=[super new];
     return self;
}

- saveinto:(NXStream *)stream as:(int)t
{
     int i,j;
     switch(t) {
       case ART_LTM_BIN:
         ART_WRITE(estab);
         for (i=0; i<M; i++)
           for (j=0; j<estab; j++)
             ART_WRITE(BU[i][j]);
         for (j=0; j<estab; j++)
           for (i=0; i<M; i++)
             ART_WRITE(TD[j][i]);
         break;
       case ART_LTM_ASCII:
         NXPrintf(stream,"TD %d X ",estab);
         NXPrintf(stream,"%d \n",M);
         for (j=0; j<estab; j++) {
           for (i=0; i<M; i++)
             ART_WRITE_ASCII(TD[j][i]);
           ART_WRITE_SLASH;
           }
         NXPrintf(stream,"BU %d X ",estab);
         NXPrintf(stream,"%d \n",M);
         for (j=0; j<estab; j++) {
           for (i=0; i<M; i++)
             ART_WRITE_ASCII(BU[i][j]);
           ART_WRITE_SLASH;
           }

         break;
       }
     return self;
}

- loadfrom:(NXStream *)stream
{
     int i,j;
     ART_READ(estab);
     for (i=0; i<M; i++)
       for (j=0; j<estab; j++)
         ART_READ(BU[i][j]);
     for (j=0; j<estab; j++)
       for (i=0; i<M; i++)
         ART_READ(TD[j][i]);
     return self;
}
// init LTM array and additionals
- init:(int)M_dim of:(int)N_dim
{
     int i,j;
     [super init];
     zone=[self zone];
     M=M_dim;
     N=N_dim;

     // Bottum up Matrix which transforms patterns
     // from F2 to F2  BU[i][j]
     BU=(float **) NXZoneMalloc(zone,M*sizeof(float *));
     for (i=0; i<M; i++)
       BU[i]=(float *) NXZoneMalloc(zone,N*sizeof(float));
     // Top down Matrix which transforms patterns
     // from F2 to F1  TD[j][i]
     TD=(float **) NXZoneMalloc(zone,N*sizeof(float *));
     for (j=0; j<N; j++)
       TD[j]=(float *) NXZoneMalloc(zone,M*sizeof(float));
     oldBU=NXZoneMalloc(zone,M*sizeof(float));
     oldTD=NXZoneMalloc(zone,M*sizeof(float));

     Avail=NXZoneMalloc(zone,N*sizeof(int));
     return self;
}

// free all allocated memory
- free
{
      int i,j;
      zone=[self zone];
      for (i=0; i<M; i++)
        NXZoneFree(zone,BU[i]);
      NXZoneFree(zone,BU);

      for (j=0; j<N; j++)
        NXZoneFree(zone,TD[j]);
      NXZoneFree(zone,TD);

      NXZoneFree(zone,oldTD);
      NXZoneFree(zone,oldBU);
      NXZoneFree(zone,Avail);
      
      return [super free];
}
// set the central parameter
-set_d:(float)learnrate h:(float)stepsize
{
     d=( learnrate > 1.0 ) ? 0.9 : ( ( learnrate < 0.0 ) ? 0.9 : learnrate ) ;
     h=stepsize;
     return self;
}
// give away the weights of a single F2 node
- zero_template:(int)template with:(float)BU_scale
{
     int i;
     for (i=0; i<M; i++)
       BU[i][template]=BU_scale/((1.0-d)*sqrt((float)M));
     for (i=0; i<M; i++)
       TD[template][i]=0.0;
     return self;
}
// give away all weights
- set_zero_with:(float)BU_scale
{
     int j;
     for (j=0; j<N; j++)
       [self zero_template:j with:BU_scale];
     estab=1;
     return self;
}
// transform a pattern from F1 to F2, seems like a korrelation
- (int)korrelate_with:(float *)P to:(float *)Y
{
    int i,j;
    for (j=0; j < estab; j++) {
      Y[j]=0.0;
      if ( Avail[j] )
        for (i=0; i<M; i++)
	  Y[j] += BU[i][j]*P[i];
      }
     return estab-1;
}
// retrun a pointer to the TD array of M from a F2 node
- (float *)TD_of:(int)j
{
  return TD[j];
}
// returns a pointer to the BU array of N from a F1 node
- (float *)BU_of:(int)i
{
  return BU[i];
}
// set a vector of dim to the TD[j] array of LTM
- set:(float *)Vec of:(int)dim to_TD_of:(int)j
{
  int i;
  for (i=0;i<dim;i++) TD[j][i]=Vec[i];
  return self;
}
// set a vector of dim to the BU[i] array of LTM
- set:(float* )Vec of:(int)dim to_BU_of:(int)i
{
  int j;
  for (j=0;j<dim;j++) BU[i][j]=Vec[j];
  return self;

}
// declaire a Node as stable as iis norm excceds the thresh
// and if the node is already learned
- (BOOL)is:(int)Node stable:(float)thresh
{
   if (Node>=(estab-1)) return NO;
   return ( (thresh < (L2_norm(TD[Node],M)*(1-d)) ) ? YES:NO );
}

// methods concern the reset state of the F2 nodes
// prepare for a new input pattern
-(int)newPattern
{
  int j;
  for (j=0; j < (estab); j++) Avail[j]=LTM_NODE_AVAIL;
  return estab-1;
}
// reset a certain node in case of the rho mismatch
-(int)resetNode:(int)j
{
  Avail[j]=LTM_NODE_RESET;
  return estab-1;
}
- (BOOL)NodeAvail:(int)i;
{
      return ( Avail[i]==LTM_NODE_AVAIL ) ? YES:NO;
}
// returns the amount of etablished F2 nodes
-(int)estabNodes
{
  return estab-1;
}
// set the estab var to a new value
- setestabNodes:(int)newestab
{
     estab=newestab+1;
     return self;
}

// perform the learning facility of the LTM array
- (float)learnPattern:(float *)Learn_Vec toNode:(int)F2_on
{
     int i;
     float BUdiff=0.0,TDdiff=0.0;
     float tmp=0.0;
     float delta1, delta2, delta3, delta4;
     float zt1,zt2;

     if (F2_on >= (estab-1) ) estab++;
     for (i=0; i<M; i++) {
       oldBU[i]=BU[i][F2_on];
       oldTD[i]=TD[F2_on][i];

       zt1=d; /* Obacht */
       zt2=Learn_Vec[i];

       // Runge Kutta calculation of the learn DGL
       tmp=BU[i][F2_on];
       delta1 = h * (zt1*(zt2-(tmp)));
       delta2 = h *  zt1*(zt2-(tmp+ delta1/2.0));
       delta3 = h *  zt1*(zt2-(tmp+ delta2/2.0));
       delta4 = h *  zt1*(zt2-(tmp+ delta3));
       BU[i][F2_on]=(float)(tmp+(delta1+2.0*delta2+2.0*delta3+delta4)/6.0);

       tmp=TD[F2_on][i];       
       delta1 = h * (zt1*(zt2-(tmp)));
       delta2 = h *  zt1*(zt2-(tmp+ delta1/2.0));
       delta3 = h *  zt1*(zt2-(tmp+ delta2/2.0));
       delta4 = h *  zt1*(zt2-(tmp+ delta3));
       TD[F2_on][i]=(float)(tmp+(delta1+2.0*delta2+2.0*delta3+delta4)/6.0);

 //      BU[i][F2_on]=RungeKutta(i,BU[i][F2_on],h);
 //      TD[F2_on][i]=RungeKutta(i,TD[F2_on][i],h);

       BUdiff += fabs(oldBU[i]-BU[i][F2_on]);
       }
     TDdiff=vec_diff(TD[F2_on],oldTD,M);
     return (BUdiff+TDdiff);
}
/*
// functions of the LTM learning facility
float dTD(int i, float z)
{
     return(zt1*(zt2-z));
}

float dBU(int i, float z)
{
     return(zt1*(zt2-z));
}

float RungeKutta( int i,float x,float h)
{
     float delta1, delta2, delta3, delta4;
     delta1 = h * dTD(i,x);
     delta2 = h * dTD(i,x + delta1/2.0);
     delta3 = h * dTD(i,x + delta2/2.0);
     delta4 = h * dTD(i,x + delta3);
     return(x+(delta1+2.0*delta2+2.0*delta3+delta4)/6.0);
}
*/
@end





