1. #!/usr/bin/env python
  2. from __future__ import division
  3. __author__ = 'Horea Christian'
  5. from matplotlib.ticker import MultipleLocator, FuncFormatter, FormatStrFormatter
  6. import mdp
  7. import matplotlib.pyplot
  8. import numpy as np
  9. from loadsums import pca_loadsum
  10. from pylab import figure, xlabel, ylabel, show, title, legend
  12. from load_vdaq_conds import data_type, bytes_per_pixel
  13. from get_data import data_name, lenheader, nframesperstim, framewidth, frameheight
  14. nframesperstim = np.array(nframesperstim)
  16. stim_start = 20
  17. stim_length = 5
  18. outfile = data_name + "-pca(80, 8).npy"
  20. timepts = np.array(nframesperstim)
  21. img_range_zero = np.arange(0, nframesperstim)
  22. img_range_one = np.arange(nframesperstim, nframesperstim * 2)
  24. fimg_pca_one = pca_loadsum(data_name, img_range_one, True, data_type, framewidth, frameheight, lenheader, bytes_per_pixel)
  25. fimg_pca_zero = pca_loadsum(data_name, img_range_zero, True, data_type, framewidth, frameheight, lenheader, bytes_per_pixel)
  26. fimg_pca = fimg_pca_one / fimg_pca_zero
  27. #fimg_pca_zero_norm = (fimg_pca_zero - np.mean(fimg_pca_zero)) / np.std(fimg_pca_zero)
  28. #result = matplotlib.mlab.PCA(fimg_pca).project(fimg_pca)
  29. #result = mdp.fastica(fimg_pca, white_comp=6)
  32. inp = fimg_pca.T # where fimg_pca is 120x(504^2)
  33. pca = mdp.nodes.PCANode(output_dim=0.97, svd=True, reduce=True) # keep only 97% of the input variance
  34. pca.train(inp)
  35. pca.stop_training()
  36. result = pca(inp)
  37. #result = np.load(outfile)
  38. result_img = pca.execute(inp,0)
  39. print np.shape(result)
  40. print np.shape(result[:,0])
  41. print result[:,0]
  42. print np.sqrt(len(result[:,0]))
  43. print np.shape(result[0])
  44. print np.sqrt(len(result[0]))
  45. print np.shape(result_img)
  46. print result_img
  48. outfile = data_name + "-pca" + str(np.shape(result)) + ".npy"
  49. np.save(outfile, result)
  52. fig = figure(facecolor='#eeeeee')
  54. x0=fig.add_subplot(111)
  55. x0.plot(result[:,0], 'k', alpha = 0.12, label="Component 1")
  56. x0.plot(result[:,1], 'r', alpha = 0.7, label="Component 2")
  57. x0.plot(result[:,2], 'g', alpha = 0.7, label="Component 3")
  58. #x0.plot(result[:,3], 'b', alpha = 0.7, label="Component 4")
  59. #x0.plot(result[:,4], 'c', alpha = 0.7, label="Component 5")
  60. #x0.plot(result[:,0] - result[:,1], 'g', alpha = 0.7, label="Component 2/Component 1")
  61. #x0.hlines(y=min(result), xmin=0, xmax = stim_start, colors = '#96D0FF', label="Baseline")
  62. #x0.hlines(y=min(result), xmin=stim_start, xmax = stim_start+stim_length, colors = '#0F23FF', label="Stimulus")
  63. #x0.axhline(y=0, color = 'k', alpha = 0.12)
  64. ylabel('R/R$_{baseline}$', fontsize='12')
  65. xlabel('Time [s]', fontsize='12')
  66. x0.xaxis.set_major_formatter(FuncFormatter(lambda x,i: "%.0f" % (x/10)))
  67. x0.xaxis.set_minor_locator(MultipleLocator(5))
  68. x0.yaxis.set_major_formatter(FormatStrFormatter('%.5f'))
  69. legend(bbox_to_anchor=(0., 1.02), loc=3, ncol=5, mode="expand", borderaxespad=0.)
  70. x0.set_xlim(0, timepts-1)
  72. #C = result[:,0]
  73. #height = np.sqrt(len(C))
  74. #a = np.reshape(C,(height,height))
  75. #x0 = fig.add_subplot(3,3,1)
  76. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  77. #cbar = fig.colorbar(cx0)
  78. #cbar.set_label('Reflectance change')
  79. #cbar.ax.tick_params(direction='out')
  80. #
  81. #C = result[:,1]
  82. #a = np.reshape(C,(height,height))
  83. #x0 = fig.add_subplot(3,3,2)
  84. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  85. #cbar = fig.colorbar(cx0)
  86. #cbar.set_label('Reflectance change')
  87. #cbar.ax.tick_params(direction='out')
  88. #
  89. #C = result[:,2]
  90. #a = np.reshape(C,(height,height))
  91. #x0 = fig.add_subplot(3,3,3)
  92. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  93. #cbar = fig.colorbar(cx0)
  94. #cbar.set_label('Reflectance change')
  95. #cbar.ax.tick_params(direction='out')
  97. #C = result[:,3]
  98. #a = np.reshape(C,(height,height))
  99. #x0 = fig.add_subplot(3,3,4)
  100. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  101. #cbar = fig.colorbar(cx0)
  102. #cbar.set_label('Reflectance change')
  103. #cbar.ax.tick_params(direction='out')
  104. #
  105. #C = result[:,4]
  106. #a = np.reshape(C,(height,height))
  107. #x0 = fig.add_subplot(3,3,5)
  108. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  109. #cbar = fig.colorbar(cx0)
  110. #cbar.set_label('Reflectance change')
  111. #cbar.ax.tick_params(direction='out')
  112. #
  113. #C = result[:,5]
  114. #a = np.reshape(C,(height,height))
  115. #x0 = fig.add_subplot(3,3,6)
  116. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  117. #cbar = fig.colorbar(cx0)
  118. #cbar.set_label('Reflectance change')
  119. #cbar.ax.tick_params(direction='out')
  120. #
  121. #C = result[:,6]
  122. #a = np.reshape(C,(height,height))
  123. #x0 = fig.add_subplot(3,3,7)
  124. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  125. #cbar = fig.colorbar(cx0)
  126. #cbar.set_label('Reflectance change')
  127. #cbar.ax.tick_params(direction='out')
  128. #
  129. #C = result[:,7]
  130. #a = np.reshape(C,(height,height))
  131. #x0 = fig.add_subplot(3,3,8)
  132. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  133. #cbar = fig.colorbar(cx0)
  134. #cbar.set_label('Reflectance change')
  135. #cbar.ax.tick_params(direction='out')
  136. #
  137. #C = result[:,8]
  138. #a = np.reshape(C,(height,height))
  139. #x0 = fig.add_subplot(3,3,9)
  140. #cx0 = x0.imshow(a, interpolation='bilinear', origin='upper', extent=None)
  141. #cbar = fig.colorbar(cx0)
  142. #cbar.set_label('Reflectance change')
  143. #cbar.ax.tick_params(direction='out')
  145. show()