Necessary Resources

Letter Image Files

EasyBMP Resources

To Compile, Use These Commands:

g++ -o -g NeuralNet.cpp

g++ -o neuralnet EasyBMP.o NeuralNet.o

neuralnet

With Large L, Segmentation Faults Can Be Avoided Using:

limit stacksize unlimited

view sourceprint?
001#include <stdio.h>    //
002#include <stdlib.h>   //
003#include <ctime>      // needed for compatible  sleep function
004#include <sstream>    // needed for incorporation of variables into file names
005#include "EasyBMP.h"  // needed for BMP input/output
006#include <iostream>   //
007#include <string>     //
008#include <vector>     //
009#include <math.h>     // needed for floor()
010#include <fstream>
011//#include <windows.h>  // needed for Sleep()
012 
013using namespace std;
014using std::cout;
015using std::endl;
016 
017ofstream hammergy("hammergy.dat", ios::app);
018ofstream flips("flips.dat");
019 
020// Changeable parameters
021float       Tstart = 0.0;                   // Start temperature for annealing (put 0 for no annealing)
022int const   L = 40;                         // Size of neural net
023int const   MC_Steps = 1000;                // Number of MC-steps
024float const RN_seed = 0.4f;                 // RNG seed
025float       spindmg = 0.2;                  // Initial spin damage
026int         teachOrtho  = 0;                // Number of pseudoorthogonal images teached
027bool        restrictJbool = false;          // If true restricts J-elements to 1 and -1
028float       bdamage = 0.0f;                 // Amount of random braindamage 0 for no braindamage;
029char*       initname = "A_40";              // Name of initial image
030 
031float       J_lifetime = 0.0;
032BMP         current_BMP;
033int         Step = 0;      
034char        image_directory[] = "Letter Images/";
035char        output_directory[] = "Run Output/";
036 
037 
038void MonteCarloStep(int [L][L], float [L][L][L][L], float T, int[L][L]);
039float CalculateEnergy(int [L][L], float [L][L][L][L]);
040float CalculateEnergyChange_SN (int [L][L], float [L][L][L][L], int, int);
041void InitializeSpins(int [L][L]);
042void InitializeSpins(int [L][L], char []);
043void InitializeJ(int, float [L][L][L][L]);
044void TeachJ(float [L][L][L][L], float, int*);
045void TeachJ(float [L][L][L][L], float, char []);
046void InputImage(char*, int[L][L]);
047void PrintSpins(int [L][L]);
048void PrintSpins_Linux(int [L][L]);
049int* loadImage(char []);
050void outputBMP(int [L][L], char []);
051void restrictJ(float [L][L][L][L]);
052void brainDamageRec(int, int, int, int, float [L][L][L][L],int [L][L]);
053void brainDamageRan(float, float [L][L][L][L]);
054void createOrthogonal(int, float [L][L][L][L]);
055void flipS(float spindmg, int S[L][L]);
056float Hamming(int [L][L], int [L][L]);
057//void sleep(unsigned int);
058 
059int main(void) {
060  float T; 
061  // Seed Random Numbers
062  srand((int)(RAND_MAX * RN_seed));
063 
064  // Create Spin State and J Matrices
065  //    x  y
066  int S[L][L] = {0};
067  //       x  y
068  int refS[L][L] = {0};
069  //      ix iy jx jy
070  float J[L][L][L][L] = {0.0};
071 
072  // Initialize Spin State and J Matrices
073  InitializeSpins(refS, "black");  // Initializes Reference Spin State
074  InitializeSpins(S, initname);              // Initializes Spins to an "A" in italics in a different, serif font
075  //InitializeSpins(S);                      // Initializes Spins Randomly
076  //InitializeJ(2, J);                       // Initializes J-elements randomly
077 
078  // "Teach" J new letters, to increase the weight of the letter, increase parameter 2
079  TeachJ(J, 1.0f, "A");
080  TeachJ(J, 1.0f, "B");
081  TeachJ(J, 1.0f, "C");
082  TeachJ(J, 1.0f, "D");
083  TeachJ(J, 1.0f, "E");
084  TeachJ(J, 1.0f, "F");
085  //TeachJ(J, 1.0f, "G");
086  //TeachJ(J, 1.0f, "H");  
087  //TeachJ(J, 1.0f, "I");
088  //TeachJ(J, 1.0f, "J");  
089     
090  createOrthogonal(teachOrtho,J);
091 
092  // randomize initial spin state
093  flipS(spindmg, S);
094 
095  // Restrict J
096  if(restrictJbool)
097    restrictJ(J);
098 
099  brainDamageRan(bdamage,J);
100   
101  // MC Run
102  // Print Initial Spin State
103  PrintSpins_Linux(S);
104 
105  // Perform MC Steps and Print Spin State
106  for (Step = 0; Step < MC_Steps; Step++) {
107      PrintSpins_Linux(S);                      // <-- In windows, use PrintSpins(), in linux use PrintSpins_Linux()
108      T = (Tstart/MC_Steps)*(MC_Steps - Step);         
109      MonteCarloStep(S, J, T, refS);
110      outputBMP(S, "Trial1");
111  }
112  return 0;
113}
114 
115// Chooses LxL random neurons to switch
116// If the energy decreases, the neuron is switched, if not it stays
117// <param S> Spin state of neurons
118// <param J> J matrix for neuron interactions
119void MonteCarloStep(int S[L][L], float J[L][L][L][L], float T, int refS[L][L]) {
120  int x, y;
121  float dE;
122  int counter;
123  counter = 0;
124  for (int i = 0; i < L*L; i++) {
125      // Choose Random Neuron to Switch
126      x = (int)(L * (rand() / (RAND_MAX + 1.0)));
127      y = (int)(L * (rand() / (RAND_MAX + 1.0)));
128 
129      // Calculate Energy Change
130      dE = CalculateEnergyChange_SN(S, J, x, y);
131 
132      if (dE < 0.0) {
133          S[x][y] = -S[x][y];
134          counter ++;
135      } else {
136          if (exp(-dE/T) > (rand() / (RAND_MAX + 1.0))) {
137              S[x][y] = -S[x][y];
138              counter ++;
139          }
140      }
141  }
142   
143  hammergy << Hamming(S, refS) << ' ' << CalculateEnergy(S, J) << endl;
144  flips << counter << "\n" << endl;
145   
146}
147 
148// Calculates the total instantaneous energy
149// <param S> Spin states of neurons
150// <param J> J matrix for neuron interactions
151// <return>  Returns the energy, E
152float CalculateEnergy(int S[L][L], float J[L][L][L][L]) {
153  int i, j;
154  float E = 0;
155 
156  for (i = 0; i < L*L; i++) {
157    for (j = i + 1; j < L*L; j++) {
158      if (!(i%L == j%L && (int)floor(1.0*i/L) == (int)floor(1.0*j/L))) {
159        E -= J[i%L][(int)floor(1.0*i/L)][j%L][(int)floor(1.0*j/L)] * S[i%L][(int)floor(1.0*i/L)] * S[j%L][(int)floor(1.0*j/L)];
160      }
161    }
162  }
163 
164  return E;
165}
166 
167// Calculates the Change in Energy for a Single Neuron Switch
168// <param S> Spin states of neurons
169// <param J> J matrix for neuron interactions
170// <param x> x coordinate for neuron to switch
171// <param y> y coordinate for neuron to switch
172// <return>  Returns the change in energy, dE
173float CalculateEnergyChange_SN (int S[L][L], float J[L][L][L][L], int x, int y) {
174  int j;
175  float dE = 0;
176 
177  for (j = 0; j < L*L; j++) {
178    if (!(x == j%L && y == (int)floor(1.0*j/L))) {
179      dE -= (float)J[x][y][j%L][(int)floor(1.0*j/L)] * (float)-S[x][y] * (float)S[j%L][(int)floor(1.0*j/L)];
180    }
181  }
182 
183  return dE;
184}
185 
186// Initializes the Spins
187// <param S> Spin states of neurons
188void InitializeSpins(int S[L][L]) {
189  for (int x = 0; x < L; x++) {
190    for (int y = 0; y < L; y++) {
191      S[x][y] = (int)floor(rand() / (RAND_MAX + 1.0) + 0.5);
192      if (S[x][y] == 0) S[x][y] = -1;
193    }
194  }
195}
196 
197// Initializes the Spins
198// <param S> Spin states of neurons
199// <param name> File name in the image directory to load
200void InitializeSpins(int S[L][L], char name[]) {
201  std::string filename;
202  std::stringstream filename_out;
203 
204  filename_out << image_directory << name << ".bmp";
205  filename = filename_out.str();
206 
207  int* image = loadImage((char*)filename.c_str());
208 
209  for (int i = 0; i < L*L; i++) {
210    S[i%L][(int)floor(1.0*i/L)] = image[i];
211  }
212}
213 
214// Initializes the J Matrix
215// <param selection> Chooses how to initialize J
216//   case 1 - All values 0
217//   case 2 - Random Initialization
218// <param J> J matrix for neuron interactions
219void InitializeJ(int selection, float J[L][L][L][L]) {
220  int ix, iy, jx, jy;
221 
222  switch (selection) {
223    // Blank
224  case 1:
225    break;
226    // Random Initialization between -5 and 5
227  case 2:
228    for (ix = 0; ix < L; ix++) {
229      for (iy = 0; iy < L; iy++) {
230        for (jx = 0; jx < L; jx++) {
231          for (jy = 0; jy < L; jy++) {
232            J[ix][iy][jx][jy] = (float)(rand() / (RAND_MAX + 1.0) * 10.0 - 5.0);
233          }
234        }
235      }
236    }
237    break;
238  }
239}
240 
241 
242// "Teaches" new J values based on an input image
243// <param exposure_time> the weight of the image on the J matrix
244// <param image> the image to affect the J matrix
245void TeachJ(float J[L][L][L][L], float exposure_time, int* image) {
246  for (int ix = 0; ix < L; ix++) {
247    for (int iy = 0; iy < L; iy++) {
248      for (int jx = 0; jx < L; jx++) {
249        for (int jy = 0; jy < L; jy++) {
250          J[ix][iy][jx][jy] = J_lifetime / (J_lifetime + exposure_time) * J[ix][iy][jx][jy] +
251              exposure_time / (J_lifetime + exposure_time) * image[ix + iy * L] * image[jx + jy * L];
252        }
253      }
254    }
255  }
256 
257  J_lifetime += exposure_time;
258}
259 
260// "Teaches" new J values based on a bmp image
261// <param exposure_time> the weight of the image on the J matrix
262// <param name> The name of the file to be opened in the image_directory (e.g. "A"
263//              with L=10 would open the file at "Letter Images/A_10.bmp"
264void TeachJ(float J[L][L][L][L], float exposure_time, char name[]) {
265  std::string filename;
266  std::stringstream filename_out;
267 
268  filename_out << image_directory << name << "_" << L << ".bmp";
269  filename = filename_out.str();
270 
271  int* image = loadImage((char*)filename.c_str());
272 
273  for (int ix = 0; ix < L; ix++) {
274    for (int iy = 0; iy < L; iy++) {
275      for (int jx = 0; jx < L; jx++) {
276        for (int jy = 0; jy < L; jy++) {
277          J[ix][iy][jx][jy] = (J_lifetime / (J_lifetime + exposure_time)) * J[ix][iy][jx][jy] +
278              (exposure_time / (J_lifetime + exposure_time)) * image[ix + iy * L] * image[jx + jy * L];
279        }
280      }
281    }
282  }
283 
284  J_lifetime += exposure_time;
285}
286 
287// Prints the spin states to the terminal
288// <param S> Spin states of neurons
289void PrintSpins(int S[L][L]) {
290  system("cls");
291  for (int x = 0; x < L; x++) {
292    for (int y = 0; y < L; y++) {
293      if (S[x][y] == -1) {
294        cout << " ";
295      } else {
296        cout << "#";
297      }
298    }
299    cout << endl;
300  }
301}
302 
303// Prints the current spin matrix
304void PrintSpins_Linux(int S[L][L]) {
305  system("clear");
306  for (int x = 0; x < L; x++) {
307    for (int y = 0; y < L; y++) {
308      if (S[x][y] == -1) {
309        cout << " ";
310      } else {
311        cout << "#";
312      }
313    }
314    cout << endl;
315  }
316}
317 
318// Loads an image from the given filename and returns a binary 2D array
319// <param filename> Relative path of image
320int* loadImage(char filename[]) {
321  BMP image;
322  image.ReadFromFile(filename);
323     
324  int* S_image = new int[L*L];
325 
326  for(int i = 0; i < L*L; i++) {
327    if(image((int)floor(1.0*i/L),i%L)->Red + image((int)floor(1.0*i/L),i%L)->Blue + image((int)floor(1.0*i/L),i%L)->Green < 382)
328      S_image[i] = 1;
329    else
330      S_image[i] = -1;
331  }
332 
333  return S_image;
334}
335 
336// Prints the current spin state to a .bmp image file in the output directory
337// <param S> Spin states of neurons
338// <param name> File name of current run (e.g. "test" in monte carlo step 1
339//              will produce the file "Run Output/test_1.bmp"
340void outputBMP(int S[L][L], char name[]) {
341  std::string filename;
342  std::stringstream filename_out;
343 
344  filename_out << output_directory << "TEMP" << "_" << Step << ".bmp";
345  filename = filename_out.str();
346  current_BMP.SetSize(L,L);
347 
348  for(int x = 0;  x < L; x++) {
349    for(int y = 0; y < L; y++) {
350      if(S[x][y] == -1) {
351        current_BMP(y,x)->Red = 255;
352        current_BMP(y,x)->Green = 255;
353        current_BMP(y,x)->Blue = 255;
354      } else {
355        current_BMP(y,x)->Red = 0;
356        current_BMP(y,x)->Green = 0;
357        current_BMP(y,x)->Blue = 0;
358      }
359    }
360  }
361 
362  current_BMP.WriteToFile((char*)filename.c_str());
363}
364 
365// Restricts elements in J to +1 or -1
366void restrictJ(float J[L][L][L][L]){
367  for (int ix = 0; ix < L; ix++) {
368    for (int iy = 0; iy < L; iy++) {
369      for (int jx = 0; jx < L; jx++) {
370        for (int jy = 0; jy < L; jy++) {
371          if(J[ix][iy][jx][jy]>=0)
372            J[ix][iy][jx][jy]=1;
373          else
374            J[ix][iy][jx][jy]=-1;
375        }
376      }
377    }
378  }
379}
380 
381// Simulates killing neurons in rectangle starting at (sx,sy) and extending dx and dy
382void brainDamageRec(int sx, int sy, int dx, int dy, float J[L][L][L][L], int S[L][L]){
383  for (int x = sx; x < sx+dx; x++) {
384    for (int y = sy; y < sy+dy; y++) {
385      S[y][x] = 1;
386      for (int ix = 0; ix < L; ix++) {
387        for (int iy = 0; iy < L; iy++) {
388          J[y][x][iy][ix] = 0;
389          J[iy][ix][y][x] = 0;
390        }
391      }
392    }
393  }
394}
395 
396// Sets elements of J to 0 with a probability of rate
397void brainDamageRan(float rate, float J[L][L][L][L]){
398  if(rate != 0){
399    for (int ix = 0; ix < L; ix++) {
400      for (int iy = 0; iy < L; iy++) {
401        for (int jx = 0; jx < L; jx++) {
402          for (int jy = 0; jy < L; jy++) {
403            if((rand() / (RAND_MAX + 1.0))<rate)
404              J[ix][iy][jx][jy]=0;
405          }
406        }
407      }
408    }
409  }
410}
411 
412// Creates and teaches pseudoorthogonal images
413void createOrthogonal(int n,float J[L][L][L][L]){
414  int * S1;
415  S1 = new int [L*L];
416  for(int i=0; i<n; i++){
417    for (int x = 0; x < L*L; x++) {
418      if((rand() / (RAND_MAX + 1.0)) < 0.5)
419        S1[x] = 1;
420      else
421        S1[x] = -1;
422    }
423 
424    TeachJ(J,1.0f,S1);
425  }
426}
427 
428// Sets initial spindamage
429void flipS(float spindmg, int S[L][L]) {
430  if(spindmg != 0){
431    for (int x = 0; x < L; x++) {
432      for (int y = 0; y < L; y++) {
433        if((rand() / (RAND_MAX + 1.0))<spindmg)
434      S[x][y] = int(((rand() % 2)-0.5)*2);
435      }
436    }
437  }
438}
439 
440// Calculates Hamming Distance
441float Hamming(int S[L][L], int refS[L][L]) {
442  float hamming = 0.0f;
443 
444  for (int x = 0; x < L; x++) {
445    for (int y = 0; y < L; y++) {
446      hamming += (1.0 / L*L) * pow(refS[x][y] - S[x][y], 2.0);
447    }
448  }
449  return hamming;
450}