/***************************************************************************************************** 
 	Code: The PINN Class that creates, trains and evaluates the PINN memory.
	Author: Dmitry O. Gorodnichy, 1993-2003
	Copyright: The author grants the right to use the code, 
		provided the appropriate references are acknowledged.

	For the description of all parameters in the code, and for the detailed description of the 
	properties and performance of the PINN memory see http://www.cv.iit.nrc.ca/~dmitry/pinn 
	and the following papers.
	
[1] D.O. Gorodnichy. "A Way to Improve Error Correction Capability of Hopfield Associative Memory in the 
Case Of Saturation", Proc. of HELNET 94-95 Intern. Workshop on Neural Networks, Vol. I/II, pp.198-212, VU University Press, Amsterdam.

[2] D.O. Gorodnichy and A.M. Reznik. "Increasing Attraction of Pseudo-Inverse Autoassociative Networks", 
Neural Processing Letters, volume 5, issue 2, pp. 123-127, 1997. 

[3] D.O. Gorodnichy and A.M. Reznik. "Static and Dynamic Attractors of Autoassociative Neural Networks", 
Proc. of 9th Intern. Conf. on Image Analysis and Processing (ICIAP'97), Florence, Vol. II, pp. 238-245, 1997. 

[3] Dmitry O. Gorodnichy. "The Optimal Value of Self-Connection", CD-ROM Proc. of Intern. Joint Conf. on 
Neural Networks (IJCNN'99), Washington DC, July 1999 - "Best Presentation" award.
	
Details of the code:
	
- Memorization of the patterns is done using the recursive Greville formula (obtained in [1]), 
	which is very fast and suitable for on-fly memory update. 
- Recognition of the patterns is done using the flood-fill neuro-processing technique (described in [2])
 	which yields very fast recognition suitable for on-line recognition.
 	
Instructions: 
	Under Visial C++, the code compiles and runs as it is. 
	- It creates M random binary samples of size N and stores them in a file.
	- It memorizes the patterns in the network, 
		after which is computes the network's remaining capacity
	- Then it retrieves from memory the same patterns corrupted by nNoise percentage of noise,
		and prints the result of recognition: remaining noise (nHemming) and projection residue (E)
	
	The size of the network (N), the sample size (M), and the percentage of noise (nNoise) can be changed 
	in main(). 	The examples of suitable values are: 
	N=1000, M=300, nNoise=10% - for regular network
	N=1000, M=300, nNoise=30% - for desaturated network
	
***********************************************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// Utility functions
int compareSamples(signed char *Y, signed char *V, int N);
bool writeRandomBinarySamples(int nNetSize=100, int nSamples=50, char *sFileName="pinn-samples.txt", double dPercentageOfOnes=0.2, int nRndSeed=1);
void getSample(signed char **V, int N, char *string);
void addNoise(signed char *Y, signed char **V, int N, int nNoise);


class CPINN  
{
public:
	CPINN();
	virtual ~CPINN();
	void createMemory(int N=100);
	void deleteMemory();
	bool putInMemory(signed char *V);
	void analyzeMemory(); // examine how filled the memory is 
	void setDesaturationLevel(double dDesaturationLevel);
	void desaturateMemory();
	void setMemoryToDormantState(); // in which all neurons are in unexcited (-1) state
	int  retrieveFromMemory(signed char *Y0, signed char **Yattractor);

	bool bCreated;
	int m_N; 
	int m_M; // number of stored prototypes in network's memory 
	float 
		**m_C,	// matrix C of weightes (aka synaptic weights) of the network
		*m_S,	// S = CY, postsynaptic potentials (PSP) of the network
		*m_S0,
		m_E;	// E = || CY - Y ||
	signed char 
		*m_Y;	// state (aka potentials, or neurons, or neuron states) of the network
	double m_D;
	int *m_aExcited; // buffer which stored indices of excited neurons
	int m_nExcitedNow, m_nExcitedLast;  // number of excited neurons at current and previous state of the network
	bool m_bCycleOccured; // redundant, as m_nExcitedNow>0 => m_bCycleOccured
	int m_nIter;
	double m_dAveCii, m_dAveCij2, m_dAveAbsCij2;

};

CPINN::CPINN()
{
	bCreated = false;
}
CPINN::~CPINN()
{
	deleteMemory();
}

void CPINN::deleteMemory()
{
	if (bCreated)
	{
		delete m_Y;
		delete m_S;
		delete m_S0;
		delete m_aExcited;
		for (int i=0; i<m_N; i++)
			delete m_C[i];
		delete  m_C;
		bCreated = false;
	}
}


void CPINN::createMemory(int N)
{
	m_M = 0; // No pattern stored
	m_N = N;
    m_C = new float*[N];
	for (int i=0; i<N; i++)
	{
		m_C[i] = new float[N];
		for (int j=0; j<N; j++)
			m_C[i][j]=0;
	}
	m_Y = new signed char[N];
	m_S = new float[N];
	m_S0 = new float[N];
	m_aExcited = new int[N];
	m_nExcitedNow=0; 
	m_nExcitedLast=0;
	bCreated = true;
	m_D = 1.0; //0.15; // The most optimal desaturation value for most cases
}


bool CPINN::putInMemory(signed char *V)
{
	int i,j;
	for(m_E=m_N,i=0; i<m_N; i++)
	{
		for(m_S[i]=0,j=0; j<m_N; j++)
			m_S[i] += m_C[j][i]*V[j];
		m_E -= m_S[i]*V[i];
	}
	if (m_E < 0.001)
		return false; // Either the memory is full OR the same(similar) pattern has been already stored

    for(i=0;i<m_N;i++)
		for(j=0;j<m_N;j++)
			m_C[i][j]+=(V[i]-m_S[i])*(V[j]-m_S[j])/m_E; // Notice that m_C is symetric.

	m_M ++; // One pattern has been added to memory
	return true;
}

void CPINN::analyzeMemory() // examine how filled the memory is 
{
	double X,Y,Z; int i, j;

// Compute the actual values
	for(X=0,Y=0,Z=0, i=0;i<m_N;i++)
    {
		for(j=0;j<m_N;j++)
		{
			Y+=m_C[i][j]*m_C[i][j];  Z+=fabs(m_C[i][j]);
		}
		X+=m_C[i][i];
    }
  
  /* Average values of weights Cii and Cij */
  X=X/m_N; 
  Y=Y/m_N/m_N; 
  Z=(Z/m_N/m_N)*(Z/m_N/m_N); 
  printf("\nC:  <Cii> = %.5f,\t   <Cij^> = %.5f  vs  <|Cij|>^=%.5f\n", X,Y,Z);
  m_dAveCii = X; 
  m_dAveCij2 = Y; 
  m_dAveAbsCij2 = Z;

 // Compute approximations of Average values of weights Cii and Cij by formula from Ref [2] 
  X=(double)m_M/m_N; 
  Y=(double)m_M*(m_N-m_M)/m_N/m_N/m_N;  
  printf("vs   M/N  = %.5f,\tM(N-M)/N3 = %.5f \n", X,Y);

}     

void CPINN::setDesaturationLevel(double dDesaturationLevel)
{
	m_D = dDesaturationLevel;
}
void CPINN::desaturateMemory()
{
	for(int i=0;i<m_N;i++)
		m_C[i][i]=m_D*m_C[i][i];
}

// By default, all neurons are not excited, i.e. are "-1"
void CPINN::setMemoryToDormantState() // in which all neurons are in unexcited (-1) state
{
	int i,j;
	for(i=0;i<m_N;i++)
		for(j=0, m_S0[i]=0;j<m_N;j++)
		  m_S0[i]-=m_C[j][i];
}
 

// Recognition stage. Retrieval from memory. Testing/examining the network.
//
// Using the flood-fill neuro-processing technique described in Ref [3], 
// which is based on the fact that in brain only a small number of neurons 
// is actually excited (i.e. are "+1"), and therefore should be processed.

int CPINN::retrieveFromMemory(signed char *Y0, signed char **Yattractor)
{	
	int i;
	bool bAttractorAchieved;
	int nOscilatingNeurons;
	// Initialize state of the network to Y0, by setting the buffer of '+1' neurons
	for(m_nExcitedNow=0, i=0; i<m_N; i++)
	{
		m_S[i]=m_S0[i];
		m_Y[i] = Y0[i];
		if (Y0[i]>0) 
		    m_aExcited[m_nExcitedNow++]=+i+1; 
	}		

	// Iterate until the network converges to an attractor, either static or dynamic(cycle).
	for(m_nIter=0, bAttractorAchieved=false; !bAttractorAchieved ;m_nIter++)
    {
		for(int e=0; e<m_nExcitedNow; e++) // Trace synapces of all excited neurons and update PSP of the connected neurons
			if (m_aExcited[e]<0)
				for(int j=0;j<m_N;j++)
					m_S[j] -= 2*m_C[-m_aExcited[e]-1][j]; // -1, because it starts from '+/-1' not '0'
			else
				for(int j=0;j<m_N;j++)
					m_S[j] += 2*m_C[m_aExcited[e]-1][j];
      
		m_nExcitedLast=m_nExcitedNow;
		m_nExcitedNow=0;
		nOscilatingNeurons=0;
		for(int i=0; i<m_N; i++) // Update the state of the network
		{
			if ( (m_S[i]<0) && (m_Y[i]>0) )
			{ 
				m_Y[i]=-1; 
				if (m_aExcited[m_nExcitedNow]==i+1)  // cycle check: see whether this neuron was excited in previous iteration
						nOscilatingNeurons++;	// if yes, update nOscilatingNeurons
					m_aExcited[m_nExcitedNow++]=-i-1; // update buffer
				}
			else if  ( (m_S[i]>=0) && (m_Y[i]<0) ) // another option: (m_S[i]>=0)
			{
				m_Y[i]=+1; 
				if (m_aExcited[m_nExcitedNow]==-i-1)  // for check of a cycle 
					nOscilatingNeurons++;
				m_aExcited[m_nExcitedNow++]=+i+1;
			}
			else
				i=i;
		}
		if (m_nExcitedNow==0) 
		{
			m_bCycleOccured = false; // it is a static attractor
			bAttractorAchieved = true;
		}
		if (m_nExcitedNow == m_nExcitedLast && nOscilatingNeurons==m_nExcitedNow) 
		{
			m_bCycleOccured = true;// it is a cycle  
			bAttractorAchieved = true;
		}	
	} // end of iteration loop
	
// Check how close the retrieved pattern is to those which were stored: m_E = ||CY-Y||^2

	double E;
	for(m_E=m_N, i=0, E=0; i<m_N; i++)
	{
		(*Yattractor)[i] = m_Y[i];
		m_E -= m_S[i]*m_Y[i]; //  equivalent to E += (m_S[i]-m_V[i])*(m_S[i]-m_V[i]);
	}

	return 1;
}

  
void main(void)
{
	CPINN pinn;
	int M	=	200, 
		N	=	1000, 
		nNoise= 300; // (int)N/10;
		
	FILE *f, *fout;
	char sFileName[] = "pinn-samples.txt";
	char sPattern[1002]; // make sure it is larger than N
//	sPattern = new char[N+1];
	signed char *V, *Y;
	V = new signed char[N]; 
	Y = new signed char[N]; 

	writeRandomBinarySamples(N, M, sFileName, 0.3, 0);
	pinn.createMemory(N);


// Memorize patterns

	printf("\nLearning patterns (N=%i, M=%i): \n", N,M);
	f=fopen(sFileName, "r");
	fgets(sPattern, 100, f);  /* first two comment lines */
	fgets(sPattern, 100, f);  /* first two comment lines */
	for (int i=0;i<M;i++)
	{
		if ( ! fgets(sPattern, N+1, f) )
			break;
		fgetc(f);
		getSample(&V, N, sPattern);
		pinn.putInMemory(V);
		printf("%i ",i);
//		pinn.getRemainedCapacity();
	}
	fclose (f);

#if 1

// Analyze memory and desaturate it if needed
	pinn.analyzeMemory();

	printf("\t<Cii>^2=%.4f, <Cij^2>=%.4f,  D=%.2f\n", 
		pinn.m_dAveCii*pinn.m_dAveCii, pinn.m_dAveCij2, pinn.m_D);  

/***
  Memory recognizes well when it is not saturated. 
  Whether the memory is saturated or not can be told by looking at ratio <Cii> / <Cij>, 
  which increases as the memory learns new patterns.
  It is advisable to desaturate the memory when  <Cii> / <Cij>  approaches   1/N.
***/
//  Uncomment these two lines to see the improvement due to the desaturation !
//
//	pinn.setDesaturationLevel(0.15);
//	pinn.desaturateMemory();

	pinn.setMemoryToDormantState();

// Recognize(Retrieve) patterns from 

	printf("\nRetrieval: \n\tN=%i, M=%i, Noise=%i\n", N,M,nNoise);  
	printf("# \tHemming\t|CV-V|\t\t#Iter \t#Oscilating\n");  

	fout=fopen("pinn-res.txt", "w");
	fprintf(fout, "Retrieval from Noise = Noise=%i\n", nNoise);  
	fprintf(fout, "# \tHemming\t|CV-V|\t\t#Iter \t#Oscilations\n");  

	f=fopen(sFileName, "r");
	fgets(sPattern, 100, f);   /* first comment line */
	fgets(sPattern, 100, f);   /* first comment line */
	

	for (i=0;i<M;i++)
	{
		if ( ! fgets(sPattern, N+1, f) )
			break;
		fgetc(f);
		getSample(&Y, N, sPattern);		
	    addNoise(Y, &V, N, nNoise); 
//		V = pinn.retrieveFromMemory(Y);
		pinn.retrieveFromMemory(V,&V);
		int nHemming = compareSamples(Y,V,N);
		printf("%i \t%i \t%+.6f \t%i \t%i\n", i, nHemming, pinn.m_E, pinn.m_nIter, pinn.m_nExcitedNow);  
//		printf("%i ",i);

		fprintf(fout, "%i \t%i \t%+.6f \t%i \t%i\n", i, nHemming, pinn.m_E, pinn.m_nIter, pinn.m_nExcitedNow);  
	}
	fclose (f);
	fclose (fout);
			
#endif

	delete Y;		
	delete V;
//	delete sPattern;
	M=M;

//	pinn.deleteMemory();
	return;
} // main()

/*** Utility functions ***
********************************************************************************************/

int compareSamples(signed char *Y, signed char *V, int N)
{
	for(int nHemming=0,i=0;i<N;i++) 
		if (Y[i]!=V[i]) 
			nHemming++;
	return nHemming;
}

bool writeRandomBinarySamples(int N, int M, char *sFileName, double dPercentageOfOnes, int nRndSeed)
{
    FILE *f;
	double dRnd;
    f=fopen(sFileName, "w");
	if (f==NULL) return false; 

	char * sSample = new char[N+1];

    fprintf(f, "Prototype set: N=%i, M=%i, L=%.2f. Seed=%i\n\n",N, M, dPercentageOfOnes, nRndSeed);
    srand( (unsigned)nRndSeed ); 
    for (int m=0;m<M;m++)
    {
		for (int i=0;i<N;i++)
		{
			dRnd = (double)rand()/RAND_MAX;
			if (dRnd < dPercentageOfOnes)
				sSample[i] = '1';
			else 
				sSample[i] = '0';
		}
		sSample[i]='\0'; 
      	fprintf(f, "%s\n", sSample);
      }
    fclose(f); 

	delete sSample;
	return true;
}

void getSample(signed char **V, int N, char *string)
{
    for(int i=0;i<N;i++)
		(*V)[i]=('1'-string[i])*2 - 1;
}

void addNoise(signed char *Y, signed char **V, int N, int nNoise)
{
	bool *bY = new bool [N];
	for(int i=0;i<N;i++)  
		bY[i]=false;
	for(i=0;i<nNoise; )  
    {  
		int kk=rand() % N;
		if (bY[kk] == false)
		{
			bY[kk]=true;  i++;
		}
    }
	for(i=0;i<N;i++)  
		if (bY[i] == true) 
			(*V)[i] = - Y[i];
		else
			(*V)[i] = Y[i];

	delete bY;
}

/*********************************************************************************************
*** Utility functions ***/