patchreconst.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348
  1. import numpy as np
  2. import tensorflow as tf
  3. from sklearn.decomposition import SparseCoder
  4. import timeit
  5. PATCH_SIZE = 32
  6. #Size of random sample from entire input space
  7. BATCH_SIZE = 5000
  8. def responses(samples, scBasis, sparsity, nonNeg, icaFilters):
  9. samples = samples.reshape((samples.shape[0], PATCH_SIZE, PATCH_SIZE, 1))
  10. patchCount = samples.shape[0]
  11. v1Simple = np.zeros((patchCount, 11, 11, 3, 12, 2))
  12. v1Complex = np.zeros((patchCount, 4356))
  13. anglesFile = np.zeros((patchCount, 11, 11, 3, 12, 1))
  14. loopCount = int(np.ceil(patchCount / BATCH_SIZE))
  15. with tf.device('/gpu:0'):
  16. for i in range(loopCount):
  17. bStart = i * BATCH_SIZE
  18. bEnd = bStart + BATCH_SIZE
  19. batch = samples[bStart:bEnd, :, :, :]
  20. filtersTensor = tf.convert_to_tensor(np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1)))
  21. v1Responses = tf.nn.conv2d(batch, filtersTensor, [1, 3, 3, 1], 'SAME')
  22. v1Responses = tf.reshape(v1Responses, [v1Responses.shape[0], 11, 11, 3, 12, 2])
  23. pair0, pair1 = tf.split(v1Responses, num_or_size_splits = 2, axis = -1)
  24. angles = tf.atan2(pair1, pair0)
  25. v1CResponses = tf.sqrt(tf.reduce_sum(tf.square(v1Responses), axis = -1))
  26. v1CResponses = tf.reshape(v1CResponses, [v1CResponses.shape[0], -1])
  27. v1Simple[bStart:bEnd, :] = v1Responses[:].numpy()
  28. v1Complex[bStart:bEnd, :] = v1CResponses[:].numpy()
  29. anglesFile[bStart:bEnd, :, :, :, :, :] = angles[:].numpy()
  30. # val = np.min(v1Complex) - 1.0
  31. # v1Complex = v1Complex.reshape((patchCount, 11, 11, 3, 12))
  32. # v1Complex[:, 5, 5, :, :] = val
  33. # v1Complex[:, 6, 5, :, :] = val
  34. # v1Complex[:, 5, 6, :, :] = val
  35. # v1Complex[:, 6, 6, :, :] = val
  36. # v1Complex = v1Complex.reshape((patchCount, -1))
  37. v1ComplexMean = np.mean(v1Complex, axis = 1, keepdims = True)
  38. v1Complex -= v1ComplexMean
  39. forwardPCA = np.load('forwardPCA.npy')
  40. pcaTransformed = np.dot(v1Complex, forwardPCA.T)
  41. sc = SparseCoder(scBasis, transform_algorithm = 'lasso_cd', transform_alpha = sparsity, positive_code = nonNeg)
  42. scResp = sc.transform(pcaTransformed)
  43. icaResp = np.dot(pcaTransformed, icaFilters.T)
  44. icaResp[np.where(icaResp < 0.0)] = 0.0
  45. return scResp, icaResp, v1Simple, v1ComplexMean, anglesFile
  46. def responsesV1(v1, scBasis, sparsity, nonNeg, icaFilters):
  47. patchCount = v1.shape[0]
  48. v1Complex = np.zeros((patchCount, 4356))
  49. anglesFile = np.zeros((patchCount, 11, 11, 3, 12, 1))
  50. loopCount = int(np.ceil(patchCount / BATCH_SIZE))
  51. with tf.device('/gpu:0'):
  52. for i in range(loopCount):
  53. bStart = i * BATCH_SIZE
  54. bEnd = bStart + BATCH_SIZE
  55. v1Responses = v1[bStart:bEnd, :]
  56. filtersTensor = tf.convert_to_tensor(np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1)))
  57. pair0, pair1 = tf.split(v1Responses, num_or_size_splits = 2, axis = -1)
  58. angles = tf.atan2(pair1, pair0)
  59. v1CResponses = tf.sqrt(tf.reduce_sum(tf.square(v1Responses), axis = -1))
  60. v1CResponses = tf.reshape(v1CResponses, [v1CResponses.shape[0], -1])
  61. v1Complex[bStart:bEnd, :] = v1CResponses[:].numpy()
  62. anglesFile[bStart:bEnd, :, :, :, :, :] = angles[:].numpy()
  63. v1ComplexMean = np.mean(v1Complex, axis = 1, keepdims = True)
  64. v1Complex -= v1ComplexMean
  65. forwardPCA = np.load('forwardPCA.npy')
  66. pcaTransformed = np.dot(v1Complex, forwardPCA.T)
  67. sc = SparseCoder(scBasis, transform_algorithm = 'lasso_cd', transform_alpha = sparsity, positive_code = nonNeg)
  68. scResp = sc.transform(pcaTransformed)
  69. icaResp = np.dot(pcaTransformed, icaFilters.T)
  70. icaResp[np.where(icaResp < 0.0)] = 0.0
  71. return scResp, icaResp, v1ComplexMean, anglesFile
  72. def responsesPCAV1C(samples):
  73. samples = samples.reshape((samples.shape[0], PATCH_SIZE, PATCH_SIZE, 1))
  74. patchCount = samples.shape[0]
  75. v1Simple = np.zeros((patchCount, 11, 11, 3, 12, 2))
  76. v1Complex = np.zeros((patchCount, 4356))
  77. anglesFile = np.zeros((patchCount, 11, 11, 3, 12, 1))
  78. loopCount = int(np.ceil(patchCount / BATCH_SIZE))
  79. with tf.device('/gpu:0'):
  80. for i in range(loopCount):
  81. bStart = i * BATCH_SIZE
  82. bEnd = bStart + BATCH_SIZE
  83. batch = samples[bStart:bEnd, :, :, :]
  84. filtersTensor = tf.convert_to_tensor(np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1)))
  85. v1Responses = tf.nn.conv2d(batch, filtersTensor, [1, 3, 3, 1], 'SAME')
  86. v1Responses = tf.reshape(v1Responses, [v1Responses.shape[0], 11, 11, 3, 12, 2])
  87. pair0, pair1 = tf.split(v1Responses, num_or_size_splits = 2, axis = -1)
  88. angles = tf.atan2(pair1, pair0)
  89. v1CResponses = tf.sqrt(tf.reduce_sum(tf.square(v1Responses), axis = -1))
  90. v1CResponses = tf.reshape(v1CResponses, [v1CResponses.shape[0], -1])
  91. v1Simple[bStart:bEnd, :] = v1Responses[:].numpy()
  92. v1Complex[bStart:bEnd, :] = v1CResponses[:].numpy()
  93. anglesFile[bStart:bEnd, :, :, :, :, :] = angles[:].numpy()
  94. # val = np.min(v1Complex) - 1.0
  95. # v1Complex = v1Complex.reshape((patchCount, 11, 11, 3, 12))
  96. # v1Complex[:, 5, 5, :, :] = val
  97. # v1Complex[:, 6, 5, :, :] = val
  98. # v1Complex[:, 5, 6, :, :] = val
  99. # v1Complex[:, 6, 6, :, :] = val
  100. # v1Complex = v1Complex.reshape((patchCount, -1))
  101. v1CCopy = v1Complex.copy()
  102. v1ComplexMean = np.mean(v1Complex, axis = 1, keepdims = True)
  103. v1Complex -= v1ComplexMean
  104. forwardPCA = np.load('forwardPCA.npy')
  105. pcaTransformed = np.dot(v1Complex, forwardPCA.T)
  106. return pcaTransformed, v1CCopy
  107. def responsesCropV1(samples, scBasis, sparsity, nonNeg, icaFilters):
  108. samples = samples.reshape((samples.shape[0], PATCH_SIZE, PATCH_SIZE, 1))
  109. patchCount = samples.shape[0]
  110. v1Simple = np.zeros((patchCount, 11, 11, 3, 12, 2))
  111. v1SimpleCropped = np.zeros((patchCount, 11, 11, 3, 12, 2))
  112. v1Complex = np.zeros((patchCount, 4356))
  113. anglesFile = np.zeros((patchCount, 11, 11, 3, 12, 1))
  114. loopCount = int(np.ceil(patchCount / BATCH_SIZE))
  115. with tf.device('/gpu:0'):
  116. for i in range(loopCount):
  117. bStart = i * BATCH_SIZE
  118. bEnd = bStart + BATCH_SIZE
  119. batch = samples[bStart:bEnd, :, :, :]
  120. filtersTensor = tf.convert_to_tensor(np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1)))
  121. v1Responses = tf.nn.conv2d(batch, filtersTensor, [1, 3, 3, 1], 'SAME')
  122. v1Responses = tf.reshape(v1Responses, [v1Responses.shape[0], 11, 11, 3, 12, 2])
  123. v1Simple[bStart:bEnd, :] = v1Responses[:].numpy()
  124. v1Np = v1Responses.numpy()
  125. v1Np[:, 3:, :, :, :, :] = 0.0
  126. v1Responses = tf.constant(v1Np)
  127. pair0, pair1 = tf.split(v1Responses, num_or_size_splits = 2, axis = -1)
  128. angles = tf.atan2(pair1, pair0)
  129. v1CResponses = tf.sqrt(tf.reduce_sum(tf.square(v1Responses), axis = -1))
  130. v1CResponses = tf.reshape(v1CResponses, [v1CResponses.shape[0], -1])
  131. v1SimpleCropped[bStart:bEnd, :] = v1Responses[:].numpy()
  132. v1Complex[bStart:bEnd, :] = v1CResponses[:].numpy()
  133. anglesFile[bStart:bEnd, :, :, :, :, :] = angles[:].numpy()
  134. v1ComplexMean = np.mean(v1Complex, axis = 1, keepdims = True)
  135. v1Complex -= v1ComplexMean
  136. forwardPCA = np.load('forwardPCA.npy')
  137. pcaTransformed = np.dot(v1Complex, forwardPCA.T)
  138. sc = SparseCoder(scBasis, transform_algorithm = 'lasso_cd', transform_alpha = sparsity, positive_code = nonNeg)
  139. scResp = sc.transform(pcaTransformed)
  140. icaResp = np.dot(pcaTransformed, icaFilters.T)
  141. icaResp[np.where(icaResp < 0.0)] = 0.0
  142. return scResp, icaResp, v1Simple, v1SimpleCropped, v1ComplexMean, anglesFile
  143. def reconstruct(basis, codes, v1cMean, angles):
  144. filters = np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1))
  145. invPCA = np.load('inversePCA.npy')
  146. super = np.dot(codes, basis)
  147. super = np.dot(super, invPCA)
  148. super += v1cMean[:super.shape[0]]
  149. super = super.reshape((super.shape[0], 11, 11, 3, 12, 1))
  150. quadPair0 = super * np.cos(angles)
  151. quadPair1 = super * np.sin(angles)
  152. v1 = np.concatenate((quadPair0, quadPair1), axis = -1)
  153. with tf.device('/gpu:0'):
  154. v1Tensor = tf.convert_to_tensor(v1, dtype = tf.float32)
  155. v1Tensor = tf.reshape(v1Tensor, [v1Tensor.shape[0], 11, 11, -1])
  156. gaborsTensor = tf.convert_to_tensor(filters, dtype = tf.float32)
  157. gaborsTensor = tf.reshape(gaborsTensor, [12, 12, 1, -1])
  158. reconst = tf.nn.conv2d_transpose(v1Tensor, gaborsTensor, [v1Tensor.shape[0], 32, 32, 1], [1, 3, 3, 1], padding = 'SAME')
  159. return reconst.numpy()[:, :, :, 0]
  160. def reconstructV1C(v1C, v1cMean, angles):
  161. filters = np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1))
  162. v1C += v1cMean
  163. v1C = v1C.reshape((v1C.shape[0], 11, 11, 3, 12, 1))
  164. quadPair0 = v1C * np.cos(angles)
  165. quadPair1 = v1C * np.sin(angles)
  166. v1 = np.concatenate((quadPair0, quadPair1), axis = -1)
  167. with tf.device('/gpu:0'):
  168. v1Tensor = tf.convert_to_tensor(v1, dtype = tf.float32)
  169. v1Tensor = tf.reshape(v1Tensor, [v1Tensor.shape[0], 11, 11, -1])
  170. gaborsTensor = tf.convert_to_tensor(filters, dtype = tf.float32)
  171. gaborsTensor = tf.reshape(gaborsTensor, [12, 12, 1, -1])
  172. reconst = tf.nn.conv2d_transpose(v1Tensor, gaborsTensor, [v1Tensor.shape[0], 32, 32, 1], [1, 3, 3, 1], padding = 'SAME')
  173. return reconst.numpy()[:, :, :, 0]
  174. def reconstructPCA(pcaTransformed, v1cMean, angles):
  175. filters = np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1))
  176. invPCA = np.load('inversePCA.npy')
  177. pcaTransformed = np.dot(pcaTransformed, invPCA)
  178. pcaTransformed += v1cMean[:pcaTransformed.shape[0]]
  179. pcaTransformed = pcaTransformed.reshape((pcaTransformed.shape[0], 11, 11, 3, 12, 1))
  180. quadPair0 = pcaTransformed * np.cos(angles)
  181. quadPair1 = pcaTransformed * np.sin(angles)
  182. v1 = np.concatenate((quadPair0, quadPair1), axis = -1)
  183. with tf.device('/gpu:0'):
  184. v1Tensor = tf.convert_to_tensor(v1, dtype = tf.float32)
  185. v1Tensor = tf.reshape(v1Tensor, [v1Tensor.shape[0], 11, 11, -1])
  186. gaborsTensor = tf.convert_to_tensor(filters, dtype = tf.float32)
  187. gaborsTensor = tf.reshape(gaborsTensor, [12, 12, 1, -1])
  188. reconst = tf.nn.conv2d_transpose(v1Tensor, gaborsTensor, [v1Tensor.shape[0], 32, 32, 1], [1, 3, 3, 1], padding = 'SAME')
  189. return reconst.numpy()[:, :, :, 0]
  190. def reconstructV1(v1):
  191. filters = np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1))
  192. with tf.device('/gpu:0'):
  193. v1Tensor = tf.convert_to_tensor(v1, dtype = tf.float32)
  194. v1Tensor = tf.reshape(v1Tensor, [v1Tensor.shape[0], 11, 11, -1])
  195. gaborsTensor = tf.convert_to_tensor(filters, dtype = tf.float32)
  196. reconst = tf.nn.conv2d_transpose(v1Tensor, gaborsTensor, [v1Tensor.shape[0], 32, 32, 1], [1, 3, 3, 1], padding = 'SAME')
  197. return reconst.numpy()[:, :, :, 0]
  198. def reconstructForV1(basis, codes, v1cMean, angles):
  199. filters = np.load('gabors.npy').astype(np.float32).reshape((12, 12, 1, -1))
  200. invPCA = np.load('inversePCA.npy')
  201. super = np.dot(codes, basis)
  202. super = np.dot(super, invPCA)
  203. super += v1cMean[:super.shape[0]]
  204. super = super.reshape((super.shape[0], 11, 11, 3, 12, 1))
  205. quadPair0 = super * np.cos(angles)
  206. quadPair1 = super * np.sin(angles)
  207. v1 = np.concatenate((quadPair0, quadPair1), axis = -1)
  208. return v1