[refactor] move tests to tests/, add pytest functions and main() entrypoints

- Moved src/tests/ → tests/
- Added test_* functions with assertions to script-style test files
- Added main() to each so IDEs offer it as a separate run target from pytest
- Fixed cupy_test.py: remove spurious x_gpu += x_cpu, fix duplicate xlabel

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

tests/README.md

@@ -0,0 +1,134 @@
+# Tests
+
+Run any test from the `src/tests/` directory:
+
+```bash
+cd src/tests
+python test_<name>.py
+```
+
+---
+
+## Self-contained tests
+
+These generate their own data and run without external files.
+
+| File | RBM type | What it tests |
+|---|---|---|
+| `test_binary_autoencoder.py` | BB-RBM | Auto-encoder on synthetic bars-and-stripes binary images (8×8). Trains, reconstructs, and plots original vs. reconstruction. |
+| `test_model.py` | BB-RBM (4 layers) | Deep auto-encoder on random binary data (1024-dim). Prints per-pattern reconstruction error. |
+| `test_xor.py` | BB-RBM | Auto-encoder on XOR-like 3-bit patterns using the older `Layer` API. |
+| `test_gaussian_autoencoder.py` | GB-RBM | Auto-encoder on synthetic Gaussian blob images (8×8). Trains, reconstructs, and plots original vs. reconstruction. |
+| `test_linear.py` | GB-RBM or GG-RBM | GB/GG-RBM on random continuous data (3000-dim). Toggle `do_gaussian_hidden` at the top of the file. |
+
+---
+
+## Tests requiring external image data
+
+These load images from a fixed path on disk.
+
+| File | RBM type | What it tests |
+|---|---|---|
+| `test_faces_sub_image.py` | GB-RBM | Auto-encoder on 32×32 patches from Caltech WebFaces (`/media/jens/cifs/bilder/MachineVision/Caltech_WebFaces/`). Uses `SubImage` with 50% overlap (stride=16); reconstruction blends overlapping patches by pixel-averaging. Supports RGB and grayscale modes. CLI args: `--grayscale`, `--load_model`, `--do_train`, `--l1_lambda`, `--num_epochs`, `--mini_batch_size`. After training: shows learned weight filters, patch-level reconstructions, and a full overlap-blended image reconstruction. Results summary (grayscale, 256 hidden units, 1000 epochs): L1=0.01 → MSE 0.168, MAE 0.214, structured edge/texture filters; L1=0.05 → MSE 0.195, MAE 0.228, sparse blob detectors. **L1=0.01 is the recommended sweet spot.** Mini-batch size has little effect on final quality (mini_batch=200 vs 1000 give equivalent MSE) but smaller batches amplify the effective L1 penalty per epoch due to more frequent updates. |
+| `test_sub_image.py` | — | Unit test for `SubImage` (patch extraction). No RBM involved. |
+| `test_conv2d.py` | — | Unit test for `conv2d`. No RBM involved. |
+
+---
+
+## Tests requiring external data
+
+These load data from `/home/jens/work/repos/Rbm/` (Armadillo `.dat` files).
+
+| File | RBM type | What it tests |
+|---|---|---|
+| `test_norbs.py` | GB-RBM or GG-RBM | Single-layer auto-encoder on NORB small (96×96 images). Toggle `DO_HIDDEN_GAUSSIAN` at the top of the file. |
+| `test_deep_norbs.py` | GB-RBM + BB-RBM | Two-layer deep auto-encoder: loads pre-trained GB-RBM weights, trains a BB-RBM on top. |
+| `test_learn_norbs_labels.py` | GB-RBM + BB-RBM | Multimodal model combining NORB image features with one-hot label embeddings in a joint hidden layer. |
+| `test_rbm.py` | Any (from `.prj` file) | Full pipeline via `StackFactory`. Takes a project name as CLI argument (`python test_rbm.py <name>`). Uses OpenCV for interactive reconstruction display. |
+
+---
+
+## Tests requiring saved state
+
+These expect a trained model to already exist in `results/`.
+
+| File | Prerequisite |
+|---|---|
+| `test_label.py` | Requires a saved `Label` fitter. Trains a binary label encoder (vocabulary of 100 integers → 5×5 binary codes) and visualises encoded labels. Set `do_train = True` on first run. |
+| `test_learn_encoded_labels.py` | Requires a saved one-hot label fitter. Trains a BB-RBM on one-hot encoded integer labels and plots hidden representation, reconstruction, and original. Set `do_train = True` on first run. |
+
+---
+
+## Implementation notes
+
+### Bug fixes applied (see `rbm/train.py`, `rbm/matrix.py`, `rbm/label.py`)
+
+#### `train.py` — momentum reset every epoch (critical)
+
+`grad_zero()` was called inside the epoch loop, zeroing the momentum accumulator before every gradient update. For the default full-batch case (one mini-batch per epoch), momentum had zero effect — every step was just `lr * dw`. Moved `grad_zero()` outside the loop so momentum accumulates correctly across epochs. **GB-RBM reconstruction MSE dropped from ~0.22 to ~0.014 on the Gaussian blob test.**
+
+#### `train.py` — `cd_gaussian_binary` negative phase (two fixes)
+
+1. **Spurious Gaussian noise** — `h_probs_neg` was computed from `data_neg + N(0,1)` instead of `data_neg`. Noise injected variance into every gradient step.
+2. **Soft hidden states for reconstruction** — `data_neg` was computed from `h_probs_pos` (probabilities) instead of sampled binary states, biasing the fantasy particle toward the mean.
+
+#### `train.py` — `cd_gaussian_gaussian` negative phase (same pattern)
+
+Same two bugs as `cd_gaussian_binary`: noise was added to `data_neg` before computing `h_probs_neg`, and `h_probs_pos` (mean) was used to drive the negative reconstruction instead of a Gaussian sample. Fixed to use sampled `h` for reconstruction and clean means for weight updates.
+
+#### `matrix.py` — `rms_error` wrong denominator
+
+`d_err_squared[1]` (row 1 of the matrix) was used instead of `d_err_squared.shape[1]` (column count). Dead code but now correct.
+
+#### `label.py` — CuPy compatibility and incomplete implementation
+
+- `_dec_binary` used Python's `reversed()` on a CuPy array (unsupported) — replaced with `np.flip()`.
+- `label2vec_onehot` returned a Python list instead of a stacked array — now uses `np.stack()`.
+- `vec2label_onehot` was unimplemented (returned `None`) — now uses `np.argmax()`.
+
+#### `entity.py` — L1 and L2 regularisation gradients wrong
+
+`grad_compute` computed scalar norms and multiplied by `W`, giving gradients that scaled with the total weight magnitude rather than the correct per-element derivatives:
+- **L1**: gradient should be `λ·sign(W)`, not `λ·‖W‖₁·W`
+- **L2**: gradient should be `2λ·W`, not `λ·‖W‖₂²·W`
+
+The combined `(l1+l2)·W` term also incorrectly mixed both penalties into one. Fixed to apply each gradient independently. `weight_decay` (the correct `λW` form of L2) is kept as a separate term.
+
+#### `image.py` — `normalize()` fragile and verbose
+
+The `repeat`/`reshape` approach was replaced with `keepdims=True` broadcasting and a `std < 1e-8` guard to prevent NaN on constant patches. The duplicate normalization logic in `test_faces_sub_image.py:load_image_patches` was removed in favour of calling `normalize()` directly.
+
+#### `train.py` — dead `cd_jens` function removed
+
+`cd_jens` was never called by `train()` or anywhere else. It contained multiple bugs (`prob()` applied to raw data, `entity.forward()` called twice). Removed entirely.
+
+#### `train.py` — `cd_gaussian_binary` redundant visible bias subtraction
+
+`dbv = sum(data_pos - b_v) - sum(data_neg - b_v)` — the `b_v` terms cancel algebraically, making the subtraction redundant and misleading. Simplified to `sum(data_pos) - sum(data_neg)`.
+
+#### `train.py` — `cd_binary_gaussian` Gibbs loop used means instead of samples
+
+For CD-k > 1, the Gibbs loop must sample both visible and hidden units to continue the Markov chain. Previously the loop used the visible mean (`prob(...)` instead of `sample(prob(...))`) and never sampled the Gaussian hidden units. Fixed to respect `do_rao_blackwell`: when False, visible is sampled and Gaussian noise is added to hidden at each step. CD-1 (default) is unaffected since the loop runs 0 times.
+
+#### `train.py` — no data shuffling between epochs
+
+Mini-batches were always sliced in the same order every epoch, introducing systematic gradient bias for mini-batch training. A random permutation is now applied to the full batch at the start of each epoch.
+
+#### `train.py` — `if` chains instead of `elif/else` for CD function selection
+
+All four `if` statements were evaluated even after a match, and an unrecognised entity type would leave `cd_func = None`, crashing later with a cryptic `TypeError`. Replaced with `elif/else` that raises `ValueError` immediately.
+
+#### `entity.py` — L1/L2 regularisation divided by mini-batch size
+
+`state_adjust(grad, k=1/N)` applied the `1/N` scale to the entire gradient, including the regularisation terms. The CD gradient is a batch sum so `1/N` normalisation is correct; L1/L2/weight-decay penalties are per-weight and batch-size independent. Dividing them by `N` made the effective regularisation strength `lambda / N`, so e.g. `l1_lambda=0.5` with `mini_batch_size=1000` was equivalent to an effective L1 of `0.0005` — too small to have any visible effect on filters. Fixed by moving the regularisation update into `state_adjust` where it is applied at full strength before the `1/N` CD update.
+
+---
+
+## RBM type reference
+
+| Type | Visible | Hidden | `EntityParams` |
+|---|---|---|---|
+| BB-RBM | Binary | Binary | `EntityParams()` |
+| GB-RBM | Gaussian | Binary | `EntityParams(do_gaussian_visible=True)` |
+| BG-RBM | Binary | Gaussian | `EntityParams(do_gaussian_hidden=True)` |
+| GG-RBM | Gaussian | Gaussian | `EntityParams(do_gaussian_visible=True, do_gaussian_hidden=True)` |

tests/cupy_test.py

@@ -0,0 +1,107 @@
+import time
+
+import numpy as np
+import cupy as cp
+
+from matplotlib import pyplot as plt
+
+def warmup_cupy():
+	a = cp.random.random((10,10))
+	b = cp.random.random((10,10))
+	_ = cp.dot(a, b)
+
+def calc_cpu(a: np.array, b: np.array):
+	return np.dot(a, b)
+
+def calc_gpu(a: np.array, b: np.array):
+	return cp.dot(a, b)
+
+def perf_test_2():
+	### Numpy and CPU
+	print("Creation")
+	s = time.time()
+	x_cpu = np.ones((1000, 1000, 1000))
+	e = time.time()
+	time_cpu = e - s
+	s = time.time()
+	x_gpu = cp.ones((1000, 1000, 1000))
+	cp.cuda.Stream.null.synchronize()
+	e = time.time()
+	time_gpu = e - s
+	print(f"Speedup: {time_cpu/time_gpu}")
+
+	### Numpy and CPU
+	print("Multiply with constant")
+	s = time.time()
+	x_cpu *= 5
+	e = time.time()
+	time_cpu = e - s
+	s = time.time()
+	x_gpu *= 5
+	cp.cuda.Stream.null.synchronize()
+	e = time.time()
+	time_gpu = e - s
+	print(f"Speedup: {time_cpu/time_gpu}")
+
+	### Numpy and CPU
+	print("Multiply with constant, Square, Square")
+	s = time.time()
+	x_cpu *= 5
+	x_cpu *= x_cpu
+	x_cpu += x_cpu
+	e = time.time()
+	time_cpu = e - s
+	s = time.time()
+	x_gpu *= 5
+	x_gpu *= x_gpu
+	x_gpu += x_gpu
+	cp.cuda.Stream.null.synchronize()
+	e = time.time()
+	time_gpu = e - s
+	print(f"Speedup: {time_cpu/time_gpu}")
+
+def perf_test_1():
+	n_values = []
+	np_times = []
+	cp_times = []
+	ratio_times = []
+
+	for N in range(100, 5100, 100):
+		a_np = np.random.rand(1024, N)
+		b_np = np.random.rand(N, 1024)
+
+		a_cp = cp.asarray(a_np)
+		b_cp = cp.asarray(b_np)
+
+		start_time = time.time()
+		c_np = calc_cpu(a_np, b_np)
+		end_time = time.time()
+		numpy_time = end_time - start_time
+
+		start_time = time.time()
+		c_cp = calc_gpu(a_cp, b_cp)
+		c_np = cp.asnumpy(c_cp)
+		end_time = time.time()
+		cupy_time = end_time - start_time
+
+		n_values.append(N)
+		np_times.append(numpy_time)
+		cp_times.append(cupy_time)
+		ratio_times.append(numpy_time/cupy_time)
+
+	plt.plot(n_values, np_times, label='numpy time')
+	plt.plot(n_values, cp_times, label='cupy time')
+	plt.plot(n_values, ratio_times, label='numpy/cupy time')
+	plt.xlabel("Matrix size [N]")
+	plt.ylabel("Time [s]")
+	plt.title("Performance numpy vs cupy")
+	plt.grid()
+	plt.legend()
+	plt.show()
+
+def main():
+	warmup_cupy()
+	perf_test_2()
+
+if __name__ == "__main__":
+	main()

tests/test_binary_autoencoder.py

@@ -0,0 +1,86 @@
+import random
+import matplotlib.pyplot as plt
+
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np
+
+N = 8
+N_VIS = N * N
+N_HID = 32
+N_CASES = 400
+
+
+class BinaryAEModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity(
+			(N_VIS, N_HID),
+			EntityParams(),
+			TrainingParams(learning_rate=0.1, momentum=0.9, num_epochs=5000, do_rao_blackwell=True)
+		)
+
+	def forward(self, x: Mat) -> Mat:
+		return self.unit1.forward(x)
+
+	def reconstruct(self, x: Mat) -> Mat:
+		return self.unit1.reconstruct(x)
+
+
+def make_bars_and_stripes(n_cases: int, n: int = N) -> Mat:
+	data = np.zeros((n_cases, n * n), dtype=np.float64)
+	for i in range(n_cases):
+		img = np.zeros((n * n,), dtype=np.float64)
+		if random.random() > 0.5:
+			row = random.randint(0, n - 1)
+			img[row * n:(row + 1) * n] = 1.0
+		else:
+			col = random.randint(0, n - 1)
+			img[col::n] = 1.0
+		data[i] = img
+	return data
+
+
+def test_binary_autoencoder():
+	model = BinaryAEModel("binary_autoencoder_pytest", "results")
+	model.unit1.training_params.num_epochs = 500
+	model.init(0.01)
+	train_batch = make_bars_and_stripes(N_CASES)
+	model.train(train_batch)
+	mae = float(np.mean(np.abs(train_batch - np.array([model.reconstruct(model.forward(x)) for x in train_batch]))))
+	assert mae < 0.4, f"Reconstruction MAE {mae:.4f} too high"
+
+
+def main():
+	model = BinaryAEModel("binary_autoencoder", "results")
+	model.init(0.01)
+
+	train_batch = make_bars_and_stripes(N_CASES)
+	model.train(train_batch)
+	model.save()
+
+	test_batch = make_bars_and_stripes(10)
+	n_show = len(test_batch)
+
+	fig, axes = plt.subplots(2, n_show, figsize=(n_show * 1.5, 3))
+	axes[0, 0].set_ylabel("Original")
+	axes[1, 0].set_ylabel("Recon")
+
+	total_error = 0.0
+	for i, inp in enumerate(test_batch):
+		h = model.forward(inp)
+		recon = model.reconstruct(h)
+		total_error += float(np.mean(np.abs(inp - recon)))
+
+		axes[0, i].imshow(np.asnumpy(np.reshape(inp, (N, N))), cmap='gray', vmin=0, vmax=1)
+		axes[0, i].axis('off')
+		axes[1, i].imshow(np.asnumpy(np.reshape(recon, (N, N))), cmap='gray', vmin=0, vmax=1)
+		axes[1, i].axis('off')
+
+	print(f"Mean reconstruction error: {total_error / n_show:.4f}")
+	plt.tight_layout()
+	plt.show()
+
+
+if __name__ == "__main__":
+	main()

tests/test_conv2d.py

@@ -0,0 +1,12 @@
+from image.sub_image import conv2d
+from rbm.matrix import np
+
+
+if __name__ == "__main__":
+	A = np.random.rand(3, 3)
+	K = np.random.rand(3, 3)
+
+	y = conv2d(A, K, (1,1))
+	print(y)
+
+	print("Test: [passed]")

tests/test_deep_norbs.py

@@ -0,0 +1,61 @@
+import os
+import matplotlib.pyplot as plt
+
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np, read_armadillo
+
+class TestModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity((96*96, 333), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False), TrainingParams(learning_rate=0.001, momentum=0.9, num_epochs=1000), enable_training=False)
+		self.unit2 = Entity((333, 256), EntityParams(), TrainingParams(learning_rate=0.1, momentum=0.9, num_epochs=10000, do_rao_blackwell=True))
+
+	def forward(self, x: Mat):
+		x = self.unit1.forward(x)
+		x = self.unit2.forward(x)
+		return x
+
+	def backward(self, x: Mat):
+		x = self.unit2.reconstruct(x)
+		x = self.unit1.reconstruct(x)
+		return x
+
+if __name__ == "__main__":
+	work_dir = "results"
+	prj_name = "deep_norb_small_16h_v2"
+	prj_root = "/home/jens/work/repos/Rbm"
+
+	# Create model
+	model = TestModel(prj_name, "results")
+
+	# Init state
+	model.init(0.1)
+
+	# load state
+	model.load()
+
+	# Load train data
+	train_batch = read_armadillo(os.path.join(prj_root, f"norb_small_16h_v2.training.dat"))
+
+	# Load test data
+	test_batch = read_armadillo(os.path.join(prj_root, f"norb_small_16h_v2.test.dat"))
+
+	# Train
+	model.train(train_batch)
+
+	# save state
+	model.save()
+
+	fig, axes = plt.subplots(1, len(test_batch), figsize=(12, 3))
+	for index, inp in enumerate(test_batch):
+		out_normalized = model.backward(model.forward(inp))
+		img = 2*(out_normalized + 0.5)
+		img = np.reshape(img, (96, 96))
+		axes[index].imshow(np.asnumpy(img))
+		axes[index].axis('off')
+
+	plt.show()
+
+
+

tests/test_faces_sub_image.py

@@ -0,0 +1,235 @@
+import os
+import numpy
+import cv2
+import matplotlib.pyplot as plt
+
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from image.sub_image import SubImageExtract, normalize
+from rbm.matrix import Mat, np, convert
+from rbm.status import CheckpointStatus
+
+DATA_DIR  = '/media/jens/cifs/bilder/MachineVision/Caltech_WebFaces/'
+PATCH     = 8
+STRIDE    = 4    # 50% overlap; set equal to PATCH for non-overlapping
+GRAYSCALE = False  # reassigned in __main__ when --grayscale is set
+N_CH      = 1 if GRAYSCALE else 3
+N_VIS     = N_CH * PATCH * PATCH   # 1024 grayscale / 3072 colour
+N_HID     = 64
+N_IMAGES  = 50
+
+
+class TestModel(Model):
+	def __init__(self, name: str, work_dir: str = '.', l1_lambda: float = 0.0,
+	             num_epochs: int = 1000, mini_batch_size: int = 1000):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity(
+			(N_VIS, N_HID),
+			EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False),
+			TrainingParams(learning_rate=0.001, momentum=0.9, num_epochs=num_epochs,
+			               mini_batch_size=mini_batch_size, l1_lambda=l1_lambda)
+		)
+
+	def forward(self, x: Mat) -> Mat:
+		return self.unit1.forward(x)
+
+	def reconstruct(self, x: Mat) -> Mat:
+		return self.unit1.reconstruct(x)
+
+
+def pad_to_stride(img: numpy.ndarray) -> numpy.ndarray:
+	"""Pad so (H - PATCH) and (W - PATCH) are divisible by STRIDE (required by SubImageExtract.check)."""
+	h, w = img.shape[:2]
+	ph = (-h) % STRIDE if h >= PATCH else PATCH - h
+	pw = (-w) % STRIDE if w >= PATCH else PATCH - w
+	return numpy.pad(img, ((0, ph), (0, pw), (0, 0)), mode='constant')
+
+
+def extract_patches_numpy(img_chw: numpy.ndarray) -> numpy.ndarray:
+	"""Extract overlapping patches (stride=STRIDE) from (C,H,W) array; returns (N, C*P*P)."""
+	C, H, W = img_chw.shape
+	rows = [img_chw[:, x:x+PATCH, y:y+PATCH].flatten()
+	        for x in range(0, H - PATCH + 1, STRIDE)
+	        for y in range(0, W - PATCH + 1, STRIDE)]
+	return numpy.stack(rows, axis=0)
+
+
+def _img_to_chw(img: numpy.ndarray) -> numpy.ndarray:
+	"""Convert cv2 BGR image to (C, H, W) float64 array, respecting GRAYSCALE flag."""
+	if GRAYSCALE:
+		gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY).astype(numpy.float64) / 255.0
+		return gray[numpy.newaxis]                                    # (1, H, W)
+	img_f = img[:, :, ::-1].astype(numpy.float64) / 255.0           # BGR→RGB
+	return numpy.ascontiguousarray(numpy.transpose(img_f, (2, 0, 1)))  # (3, H, W)
+
+
+def load_patches(data_dir: str, n_images: int = N_IMAGES, start: int = 0):
+	"""Load images, extract overlapping patches, return normalised cupy array."""
+	all_patches = []
+	files = sorted(f for f in os.listdir(data_dir) if f.endswith('.jpg'))[start:start + n_images]
+	for fname in files:
+		img = cv2.imread(os.path.join(data_dir, fname))
+		if img is None or img.shape[0] < PATCH or img.shape[1] < PATCH:
+			continue
+		img = pad_to_stride(img)
+		img_chw = _img_to_chw(img)
+		all_patches.append(extract_patches_numpy(img_chw))
+
+	patches_np = numpy.concatenate(all_patches, axis=0)
+	std = numpy.std(patches_np, axis=1)
+	patches_np = patches_np[std > 0.01]
+
+	# Normalize on CPU; return numpy — train() converts mini-batches to GPU on the fly
+	mean = patches_np.mean(axis=1, keepdims=True)
+	std2 = patches_np.std(axis=1, keepdims=True)
+	return (patches_np - mean) / numpy.where(std2 < 1e-8, 1.0, std2)
+
+
+def load_image_patches(path: str):
+	"""Extract normalised patches from a single image via SubImageExtract; return (patches, nx_steps, ny_steps)."""
+	sub = SubImageExtract(PATCH, PATCH, STRIDE, STRIDE)
+	img = cv2.imread(path)
+	img = pad_to_stride(img)
+	h, w = img.shape[:2]
+	nx_steps = (h - PATCH) // STRIDE + 1
+	ny_steps = (w - PATCH) // STRIDE + 1
+
+	img_nchw = np.array(_img_to_chw(img)[numpy.newaxis])
+	patches_nchw = sub(img_nchw)
+	patches = patches_nchw.reshape(-1, N_VIS)
+	return normalize(patches), nx_steps, ny_steps
+
+
+def assemble_overlapping(patch_list, nx_steps: int, ny_steps: int) -> numpy.ndarray:
+	"""Reconstruct image from overlapping patches by averaging contributions per pixel."""
+	H = (nx_steps - 1) * STRIDE + PATCH
+	W = (ny_steps - 1) * STRIDE + PATCH
+	accum = numpy.zeros((N_CH, H, W), dtype=numpy.float64)
+	count = numpy.zeros((1,    H, W), dtype=numpy.float64)
+
+	for k, p in enumerate(patch_list):
+		x = (k // ny_steps) * STRIDE
+		y = (k %  ny_steps) * STRIDE
+		arr = convert(p).reshape(N_CH, PATCH, PATCH)
+		accum[:, x:x+PATCH, y:y+PATCH] += arr
+		count[0, x:x+PATCH, y:y+PATCH] += 1.0
+
+	img = (accum / numpy.maximum(count, 1)).transpose(1, 2, 0).squeeze()  # (H,W,C) or (H,W)
+	lo, hi = img.min(), img.max()
+	return numpy.clip((img - lo) / (hi - lo + 1e-8), 0, 1)
+
+
+def show_filters(entity: Entity, rows: int = 8, cols: int = 16):
+	W = convert(entity.state.w_hv)
+	fig, axes = plt.subplots(rows, cols, figsize=(cols * 1.2, rows * 1.2))
+	fig.suptitle(f'Learned GB-RBM filters  ({N_HID} hidden units, {PATCH}×{PATCH} RGB patches)', fontsize=10)
+	for j in range(rows * cols):
+		ax = axes[j // cols][j % cols]
+		ax.axis('off')
+		if j >= N_HID:
+			continue
+		filt = W[:, j].reshape(N_CH, PATCH, PATCH).transpose(1, 2, 0).squeeze()
+		lo, hi = filt.min(), filt.max()
+		ax.imshow(numpy.clip((filt - lo) / (hi - lo + 1e-8), 0, 1),
+		          cmap='gray' if GRAYSCALE else None)
+	plt.tight_layout()
+
+
+def show_reconstructions(model: TestModel, patches, n_show: int = 10):
+	fig, axes = plt.subplots(2, n_show, figsize=(n_show * 1.5, 3.5))
+	axes[0, 0].set_ylabel('Original')
+	axes[1, 0].set_ylabel('Recon')
+	fig.suptitle('Patch reconstructions (normalised display)', fontsize=10)
+	err_total = 0.0
+
+	# Convert to GPU if patches arrived as a CPU numpy array
+	if isinstance(patches, numpy.ndarray):
+		patches = np.array(patches)
+
+	cmap = 'gray' if GRAYSCALE else None
+
+	def to_img(p):
+		a = convert(p).reshape(N_CH, PATCH, PATCH).transpose(1, 2, 0).squeeze()
+		return numpy.clip((a - a.min()) / (a.max() - a.min() + 1e-8), 0, 1)
+
+	for i in range(n_show):
+		inp   = patches[i]
+		recon = model.reconstruct(model.forward(inp))
+		err_total += float(np.mean(np.abs(inp - recon)))
+		axes[0, i].imshow(to_img(inp),   cmap=cmap); axes[0, i].axis('off')
+		axes[1, i].imshow(to_img(recon), cmap=cmap); axes[1, i].axis('off')
+	print(f'Mean patch reconstruction MAE: {err_total / n_show:.4f}')
+	plt.tight_layout()
+
+
+def show_image_reconstruction(model: TestModel, img_path: str):
+	patches, nx_steps, ny_steps = load_image_patches(img_path)
+	recons = [model.reconstruct(model.forward(patches[i])) for i in range(patches.shape[0])]
+
+	orig_grid  = assemble_overlapping([patches[i] for i in range(len(recons))], nx_steps, ny_steps)
+	recon_grid = assemble_overlapping(recons, nx_steps, ny_steps)
+
+	fig, axes = plt.subplots(1, 2, figsize=(12, 6))
+	cmap = 'gray' if GRAYSCALE else None
+	axes[0].imshow(orig_grid,  cmap=cmap); axes[0].set_title('Original (normalised)');        axes[0].axis('off')
+	axes[1].imshow(recon_grid, cmap=cmap); axes[1].set_title('Reconstructed (overlap blended)'); axes[1].axis('off')
+	fig.suptitle(os.path.basename(img_path), fontsize=9)
+	plt.tight_layout()
+
+
+if __name__ == '__main__':
+	from argparse import ArgumentParser
+	ap = ArgumentParser()
+	ap.add_argument('--load_model', type=lambda s: s.lower() != 'false', default=True,
+	                help='Load saved model weights (default: true)')
+	ap.add_argument('--do_train',   type=lambda s: s.lower() != 'false', default=False,
+	                help='Train the model (default: false)')
+	ap.add_argument('--grayscale',  action='store_true', default=False,
+	                help='Use single-channel grayscale patches (default: false)')
+	ap.add_argument('--l1_lambda',  type=float, default=0.0,
+	                help='L1 regularisation strength (default: 0.0)')
+	ap.add_argument('--num_epochs',      type=int,   default=1000,
+	                help='Number of training epochs (default: 1000)')
+	ap.add_argument('--mini_batch_size', type=int,   default=1000,
+	                help='Mini-batch size (default: 1000)')
+	ap.add_argument('--n_images',        type=int,   default=N_IMAGES,
+	                help=f'Number of training images (default: {N_IMAGES})')
+	args = ap.parse_args()
+
+	if args.grayscale:
+		GRAYSCALE = True
+		N_CH      = 1
+		N_VIS     = PATCH * PATCH   # 1024
+
+	prj_name = 'faces_sub_image_gray' if args.grayscale else 'faces_sub_image'
+	work_dir = 'results'
+
+	model = TestModel(prj_name, work_dir, l1_lambda=args.l1_lambda,
+	                  num_epochs=args.num_epochs, mini_batch_size=args.mini_batch_size)
+	model.init(0.001)
+
+	if args.load_model:
+		model.load()
+
+	if args.do_train:
+		print(f'Loading {args.n_images} training images (stride={STRIDE})...')
+		train_patches = load_patches(DATA_DIR, args.n_images)
+		print(f'Training on {train_patches.shape[0]} patches ({N_VIS}-dim each)')
+		model.train(train_patches, status=CheckpointStatus(save_fn=model.save))
+		model.save()
+
+	# Show learned weight filters
+	show_filters(model.unit1)
+
+	# Show patch-level reconstructions on held-out images
+	print('Loading test patches...')
+	test_patches = load_patches(DATA_DIR, n_images=10, start=N_IMAGES)
+	if test_patches.shape[0] < 10:
+		test_patches = load_patches(DATA_DIR, n_images=10)
+	show_reconstructions(model, test_patches)
+
+	# Show full image reconstruction with overlap blending
+	test_img = os.path.join(DATA_DIR, sorted(os.listdir(DATA_DIR))[60])
+	show_image_reconstruction(model, test_img)
+
+	plt.show()

tests/test_gaussian_autoencoder.py

@@ -0,0 +1,87 @@
+import random
+import matplotlib.pyplot as plt
+
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from image.sub_image import normalize
+from rbm.matrix import Mat, np
+
+N = 8
+N_VIS = N * N
+N_HID = 32
+N_CASES = 400
+
+
+class GaussianAEModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity(
+			(N_VIS, N_HID),
+			EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False),
+			TrainingParams(learning_rate=0.001, momentum=0.9, num_epochs=3000)
+		)
+
+	def forward(self, x: Mat) -> Mat:
+		return self.unit1.forward(x)
+
+	def reconstruct(self, x: Mat) -> Mat:
+		return self.unit1.reconstruct(x)
+
+
+def make_gaussian_blobs(n_cases: int, n: int = N) -> Mat:
+	xs = np.arange(n, dtype=np.float64)
+	ys = np.arange(n, dtype=np.float64)
+	xx, yy = np.meshgrid(xs, ys)
+	data = np.zeros((n_cases, n * n), dtype=np.float64)
+	for i in range(n_cases):
+		cx = random.uniform(1.0, n - 2.0)
+		cy = random.uniform(1.0, n - 2.0)
+		img = np.exp(-((xx - cx) ** 2 + (yy - cy) ** 2) / 4.0)
+		data[i] = img.flatten()
+	return data
+
+
+def test_gaussian_autoencoder():
+	model = GaussianAEModel("gaussian_autoencoder_pytest", "results")
+	model.unit1.training_params.num_epochs = 1500
+	model.init(0.01)
+	train_batch = normalize(make_gaussian_blobs(N_CASES))
+	model.train(train_batch)
+	for x in train_batch[:10]:
+		recon = model.reconstruct(model.forward(x))
+		assert recon.size == x.size
+
+
+def main():
+	model = GaussianAEModel("gaussian_autoencoder", "results")
+	model.init(0.01)
+
+	train_batch = normalize(make_gaussian_blobs(N_CASES))
+	model.train(train_batch)
+	model.save()
+
+	test_batch = normalize(make_gaussian_blobs(10))
+	n_show = len(test_batch)
+
+	fig, axes = plt.subplots(2, n_show, figsize=(n_show * 1.5, 3))
+	axes[0, 0].set_ylabel("Original")
+	axes[1, 0].set_ylabel("Recon")
+
+	total_error = 0.0
+	for i, inp in enumerate(test_batch):
+		h = model.forward(inp)
+		recon = model.reconstruct(h)
+		total_error += float(np.mean(np.abs(inp - recon)))
+
+		axes[0, i].imshow(np.asnumpy(np.reshape(inp, (N, N))), cmap='viridis')
+		axes[0, i].axis('off')
+		axes[1, i].imshow(np.asnumpy(np.reshape(recon, (N, N))), cmap='viridis')
+		axes[1, i].axis('off')
+
+	print(f"Mean reconstruction error: {total_error / n_show:.4f}")
+	plt.tight_layout()
+	plt.show()
+
+
+if __name__ == "__main__":
+	main()

tests/test_label.py

@@ -0,0 +1,47 @@
+import matplotlib.pyplot as plt
+
+from rbm.matrix import np
+from label.label import Label
+
+if __name__ == "__main__":
+	voc_size = 100
+	do_train = True
+	work_dir = "../../results"
+	prj_root = "/home/jens/work/repos/Rbm"
+
+	# Create Labeler
+	label_w = 5
+	label_h = 5
+	enc = Label(voc_size, label_w*label_h, Label.EncodingType.Binary, work_dir)
+
+	# Load
+	enc.fitter.load()
+
+	# Train
+	if do_train:
+		train_labels = list(range(voc_size))
+		enc.fit(train_labels)
+		enc.fitter.save()
+
+	# Encode test labels
+	test_labels = np.array([1, 23, 99, 37, 55, 7, 31, 10, 19, 70])
+	encoded, binary_labels = enc.encode(test_labels)
+
+	# Decode
+	decoded_soft = enc.decode(encoded)
+	decode_hard = (decoded_soft > 0.90).astype(int)
+
+	print(f"Soft Error: {np.sum((decoded_soft - binary_labels) ** 2)}")
+	print(f"Hard Error: {np.sum((decode_hard - binary_labels) ** 2)}")
+
+	fig, axes = plt.subplots(1, len(encoded), figsize=(12, 3))
+	for index, inp in enumerate(encoded):
+		img = np.reshape(inp, (label_w, label_h))
+		axes[index].imshow(img)
+		axes[index].axis('off')
+		axes[index].set_title(f'{test_labels[index]}')
+
+	plt.show()
+
+	print("Test: [passed]")
+

tests/test_learn_encoded_labels.py

@@ -0,0 +1,82 @@
+import matplotlib.pyplot as plt
+
+
+from model.model import Model
+from label.label import Label
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np
+
+class LabelLearner(Model):
+	def __init__(self, name: str, dim, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity(dim, EntityParams(do_gaussian_hidden=False), TrainingParams(learning_rate=0.01, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+
+	def forward(self, x: Mat):
+		x = self.unit1.forward(x)
+		return x
+
+	def backward(self, x: Mat):
+		x = self.unit1.reconstruct(x)
+		return x
+
+if __name__ == "__main__":
+	voc_size = 100
+	do_train = False
+	work_dir = "../../results"
+	prj_root = "/home/jens/work/repos/Rbm"
+
+	# Create Labeler
+	label_w = 5
+	label_h = 5
+	enc = Label(voc_size, label_w*label_h, Label.EncodingType.OneHot, work_dir)
+
+	# Load
+	enc.fitter.load()
+
+	# Train
+	if do_train:
+		train_labels = list(range(voc_size))
+		enc.fit(train_labels)
+		enc.fitter.save()
+
+	# Encode test labels
+	test_labels = np.array([1, 23, 99, 37, 55, 7, 31, 10, 19, 70])
+	encoded, _ = enc.encode(test_labels)
+
+	# Create model
+	model = LabelLearner("Label-Learner", (label_w*label_h, 16), work_dir)
+
+	# Init state
+	model.init(0.01)
+
+	# load state
+	model.load()
+
+	# Train
+	model.train(encoded)
+
+	# save state
+	model.save()
+
+	# Plot
+	cmap = 'Grays'
+	fig, axes = plt.subplots(3, len(encoded), figsize=(12, 3))
+	for index, inp in enumerate(encoded):
+		hidden = model.forward(inp)
+		recon = model.backward(hidden)
+		img_hidden = np.reshape(hidden, (4, 4))
+		img_recon = np.reshape(recon, (label_w, label_h))
+		axes[0][index].imshow(img_hidden.get(), cmap=cmap)
+		axes[0][index].axis('off')
+		axes[0][index].set_title(f'{test_labels[index]}')
+		axes[1][index].imshow(img_recon.get(), cmap=cmap)
+		axes[1][index].axis('off')
+		axes[1][index].set_title(f'{test_labels[index]}')
+		axes[2][index].imshow(np.reshape(inp.get(), (label_w, label_h)), cmap=cmap)
+		axes[2][index].axis('off')
+		axes[2][index].set_title(f'{test_labels[index]}')
+
+	plt.show()
+
+	print("Test: [passed]")
+

tests/test_learn_norbs_labels.py

@@ -0,0 +1,92 @@
+import os
+import matplotlib.pyplot as plt
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np, read_armadillo
+from rbm.train import train, Status
+from label.label import Label
+
+class TestModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+#		self.unit_image = Entity((96*96, 16), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False, num_gibbs_samples=3), TrainingParams(learning_rate=0.00001, momentum=0.9, do_rao_blackwell=False, num_epochs=1000))
+		self.unit_image = Entity((96 * 96, 16), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False))
+		self.unit_combine = Entity((32, 64), EntityParams(), TrainingParams(learning_rate=0.01, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+
+	def forward2(self, images: Mat, labels: Mat):
+		h = self.unit_image.forward(images)
+		h = np.ndarray.flatten(h)
+		v_combine = np.concat((h, labels), axis=0)
+		h_res = self.unit_combine.forward(v_combine)
+		return h_res
+
+	def reconstruct(self, h: Mat):
+		v_combine = self.unit_combine.reconstruct(h)
+		v_image = self.unit_image.reconstruct(np.transpose(np.transpose(v_combine)[0:self.unit_image.shape[1]]))
+		v_label = np.transpose(np.transpose(v_combine)[self.unit_image.shape[1]:self.unit_image.shape[1]+self.unit_image.shape[1]])
+		return np.ndarray.flatten(v_image), np.ndarray.flatten(v_label)
+
+	def train2(self, images: Mat, labels: Mat):
+		train(self.unit_image, images, Status())
+		h11 = self.unit_image.forward(images)
+		x2 = np.concat((h11, labels), axis=1)
+		train(self.unit_combine, x2, Status())
+
+if __name__ == "__main__":
+	work_dir = "results"
+	prj_name = "norb_small_16h_v2"
+	prj_root = "/home/jens/work/repos/Rbm"
+
+	# Create model
+	model = TestModel("norb_small_16h_v2", "results")
+
+	# Init state
+	model.init(0.01)
+
+	# load state
+	model.load()
+
+	# Load train data
+	train_batch = read_armadillo(os.path.join(prj_root, f"{prj_name}.training.dat"))
+
+	# Load test data
+	test_batch = read_armadillo(os.path.join(prj_root, f"{prj_name}.test.dat"))
+
+	# Train
+	# Create Labeler
+	enc = Label(20, 16, Label.EncodingType.OneHot, work_dir)
+
+	# Load
+	enc.fitter.load()
+
+	# Encode test labels
+	train_labels = np.array([1, 2, 3, 4, 5])
+	enc.fit(train_labels)
+
+	encoded, _ = enc.encode(train_labels)
+
+	model.train2(train_batch, encoded)
+
+	# save state
+	model.save()
+
+	cmap = 'Grays'
+	fig, axes = plt.subplots(3, len(train_batch), figsize=(12, 3))
+	for index, image in enumerate(train_batch):
+		label_zero = np.zeros(shape=16)
+		image_zero = np.zeros(shape=96*96)
+		image_reconst, labels_reconst = model.reconstruct(model.forward2(image, label_zero))
+#		image_reconst, labels_reconst = model.reconstruct(model.forward2(image_zero, encoded[index, :]))
+		img = 1-2*(image_reconst + 0.5)
+		img = np.reshape(img, (96, 96))
+		axes[0][index].imshow(img, cmap=cmap)
+		axes[0][index].axis('off')
+		axes[1][index].imshow(np.reshape(labels_reconst, (4,4)), cmap=cmap)
+		axes[1][index].axis('off')
+		axes[2][index].imshow(np.reshape(encoded[index, :], (4,4)), cmap=cmap)
+		axes[2][index].axis('off')
+
+	plt.show()
+
+
+

tests/test_linear.py

@@ -0,0 +1,67 @@
+from rbm.matrix import Mat, np
+from rbm.entity import Entity, EntityParams, TrainingParams
+from model.model import Model
+from image.sub_image import normalize
+
+WORK_DIR = "../../results"
+USE_OPTIMIZER = True
+N_VIS = 3000
+N_CASES = 1000
+class LinearModel(Model):
+	def __init__(self, name: str, work_dir: str = '.', do_gaussian_hidden=False):
+		super().__init__(name, work_dir)
+
+		if do_gaussian_hidden:
+			# Hidden gaussian
+			self.unit1 = Entity((N_VIS, 64), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=True),
+								TrainingParams(learning_rate=0.01, momentum=0.9, num_epochs=1000))
+		else:
+			# Hidden binary
+			self.unit1 = Entity((N_VIS, 1000), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False),
+								TrainingParams(learning_rate=0.005, momentum=0.9, num_epochs=1000, mini_batch_size=1000))
+
+	def forward(self, x: Mat):
+		x = self.unit1.forward(x)
+		return x
+
+	def reconstruct(self, x: Mat):
+		x = self.unit1.reconstruct(x)
+		return x
+
+class _SmallLinearModel(Model):
+	def __init__(self):
+		super().__init__("linear_pytest", "results")
+		self.unit1 = Entity((64, 32), EntityParams(do_gaussian_visible=True),
+		                    TrainingParams(learning_rate=0.005, momentum=0.9, num_epochs=100))
+
+	def forward(self, x: Mat):
+		return self.unit1.forward(x)
+
+	def reconstruct(self, x: Mat):
+		return self.unit1.reconstruct(x)
+
+
+def test_linear():
+	model = _SmallLinearModel()
+	model.init(0.1)
+	batch = normalize(np.random.randn(50, 64, dtype=np.float64))
+	model.train(batch)
+	for pattern in batch:
+		recon = model.reconstruct(model.forward(pattern))
+		assert recon.size == pattern.size
+
+
+def main():
+	model = LinearModel("linear", "results")
+	model.init(0.1)
+	model.load()
+	training_batch = np.random.randn(N_CASES, N_VIS, dtype=np.float64)
+	training_batch_norm = normalize(training_batch)
+	model.train(training_batch_norm)
+	model.save()
+	for pattern in training_batch:
+		model.reconstruct(model.forward(pattern))
+
+
+if __name__ == "__main__":
+	main()

tests/test_model.py

@@ -0,0 +1,56 @@
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np
+
+class DeepModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+		self.unit1 = Entity((1024, 333), EntityParams(), TrainingParams(learning_rate=0.1, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+		self.unit2 = Entity((333, 64), EntityParams(), TrainingParams(learning_rate=0.1, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+		self.unit3 = Entity((64, 128), EntityParams(), TrainingParams(learning_rate=0.1, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+		self.unit4 = Entity((128, 128), EntityParams(), TrainingParams(learning_rate=0.1, momentum=0.9, do_rao_blackwell=True, num_epochs=1000))
+
+	def forward(self, x: Mat):
+		x = self.unit1.forward(x)
+		x = self.unit2.forward(x)
+		x = self.unit3.forward(x)
+		x = self.unit4.forward(x)
+		return x
+
+	def backward(self, x: Mat):
+		x = self.unit4.reconstruct(x)
+		x = self.unit3.reconstruct(x)
+		x = self.unit2.reconstruct(x)
+		x = self.unit1.reconstruct(x)
+		return x
+
+def test_deep_model():
+	model = DeepModel("TestModel_pytest", "results")
+	model.unit1.training_params.num_epochs = 100
+	model.unit2.training_params.num_epochs = 100
+	model.unit3.training_params.num_epochs = 100
+	model.unit4.training_params.num_epochs = 100
+	model.init(0.1)
+	batch = (np.random.rand(32, 1024) > 0.5).astype(np.float64)
+	model.train(batch)
+	errors = [float(np.mean((model.backward(model.forward(x)) - x) ** 2)) for x in batch]
+	assert sum(errors) / len(errors) < 0.5, "Deep model reconstruction MSE too high"
+
+
+def main():
+	model = DeepModel("TestModel", "results")
+	model.init(0.1)
+	model.load()
+	batch = (np.random.rand(64, 1024) > 0.5).astype(np.float64)
+	model.train(batch)
+	model.save()
+	for index, inp in enumerate(batch):
+		out = model.backward(model.forward(inp))[0]
+		print(f"- Pattern {index} -------------------------")
+		print(f"Input  : {inp}")
+		print(f"Output : {(out > 0.9).astype(np.float64)}")
+		print(f"Error  : {np.mean((out - inp)**2):0.3f}")
+
+
+if __name__ == "__main__":
+	main()

tests/test_norbs.py

@@ -0,0 +1,71 @@
+import os
+import matplotlib.pyplot as plt
+
+from model.model import Model
+from rbm.entity import Entity, EntityParams, TrainingParams
+from rbm.matrix import Mat, np, read_armadillo
+
+DO_HIDDEN_GAUSSIAN = True
+
+class TestModel(Model):
+	def __init__(self, name: str, work_dir: str = '.'):
+		super().__init__(name, work_dir)
+
+		if DO_HIDDEN_GAUSSIAN:
+			# Hidden gaussian
+			self.unit1 = Entity((96*96, 16), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=True),
+								TrainingParams(learning_rate=0.00001, momentum=0.9, num_epochs=1000))
+		else:
+			# Hidden binary
+			self.unit1 = Entity((96 * 96, 333), EntityParams(do_gaussian_visible=True, do_gaussian_hidden=False),
+								TrainingParams(learning_rate=0.001, momentum=0.9, num_epochs=1000, l2_lambda=0.001))
+
+	def forward(self, x: Mat):
+		x = self.unit1.forward(x)
+		return x
+
+	def reconstruct(self, x: Mat):
+		x = self.unit1.reconstruct(x)
+		return x
+
+if __name__ == "__main__":
+	work_dir = "results"
+	prj_name = "norb_small_16h_v2"
+	prj_root = "/home/jens/work/repos/Rbm"
+
+	# Create model
+	model = TestModel(prj_name, "results")
+
+	# Init state
+	if DO_HIDDEN_GAUSSIAN:
+		model.init(0.01)
+	else:
+		model.init(0.1)
+
+	# load state
+#	model.load()
+
+	# Load train data
+	train_batch = read_armadillo(os.path.join(prj_root, f"{prj_name}.training.dat"))
+
+	# Load test data
+	test_batch = read_armadillo(os.path.join(prj_root, f"{prj_name}.test.dat"))
+
+	# Train
+	model.train(train_batch)
+
+	# save state
+	model.save()
+
+	fig, axes = plt.subplots(1, len(test_batch), figsize=(12, 3))
+	for index, inp in enumerate(test_batch):
+		out_normalized = model.reconstruct(model.forward(inp))
+		img = 2*(out_normalized + 0.5)
+		img = np.reshape(img, (96, 96))
+		axes[index].imshow(np.asnumpy(img))
+		axes[index].axis('off')
+
+	plt.show()
+
+
+

tests/test_rbm.py

@@ -0,0 +1,87 @@
+import os.path
+import cv2 as cv
+from argparse import ArgumentParser
+from stack.factory import StackFactory
+from rbm.status import Status
+from stack.deep import StackDeep
+from rbm.matrix import Mat, np, convert, read_armadillo
+from rbm.entity import Entity
+
+def cv_show(name: str, vec: Mat, shape):
+	img = cv.Mat(convert(np.resize(vec, shape)))
+	img_n = cv.normalize(src=img, dst=None, alpha=255, beta=0, norm_type=cv.NORM_MINMAX, dtype=cv.CV_8U)
+	cv.imshow(f"{name}", img_n)
+
+class MyStatus(Status):
+	def __init__(self, _stack: StackDeep, _batch: Mat):
+		Status.__init__(self, update_interval=10)
+		self.stack = _stack
+		self.batch = _batch
+		self.index = 0
+
+	def on_report(self, entity: Entity, status: dict) -> bool:
+		do_continue = True
+		# print status values
+		Status.print_status(entity, status)
+		# Shape of training vector
+		shape = self.stack.from_index(0).shape[0:2] + (1,)
+
+		# User input
+		max_index = self.batch.shape[0] - 1
+		key = cv.waitKeyEx(1)
+		if key == ord('q'):
+			do_continue = False
+		if key == ord('-'):
+			self.index = max(0, self.index-1)
+		if key == ord('+'):
+			self.index = min(max_index, self.index+1)
+
+		# Show reconstruction
+		img = self.stack.pass_down_up(self.batch[self.index,:])
+		cv_show("img", img, shape)
+
+		return do_continue
+
+
+def main(prj_name: str = "test"):
+	work_dir = "../../results"
+	prj_root = "/home/jens/work/repos/Rbm"
+	prj_path = os.path.join(prj_root, f"{prj_name}.prj")
+	# Create stack from project file
+	stack = StackFactory.from_file(prj_path, work_dir=work_dir)
+
+	# Init state
+	stack.state_init(0.01)
+
+	# Load state
+	stack.state_load()
+
+	# Load train data
+	training_data = read_armadillo(os.path.join(prj_root, f"{prj_name}.training.dat"))
+	try:
+		test_data = read_armadillo(os.path.join(prj_root, f"{prj_name}.test.dat"))
+	except FileNotFoundError:
+		test_data = training_data
+
+	# Prepare status listener
+	my_status = MyStatus(stack, test_data)
+
+	# Train
+	stack.train(training_data, status=my_status)
+
+	# Save state
+	stack.state_save()
+
+if __name__ == "__main__":
+	ap = ArgumentParser()
+	ap.add_argument("name", type=str,
+		default='default',
+		help="Name of project")
+
+	args = ap.parse_args()
+	var_args = vars(args)
+
+	main(var_args["name"])
+	cv.destroyAllWindows()
+
+	print("Test: [passed]")

tests/test_rnn.py

@@ -0,0 +1,99 @@
+"""Tests for StackRnn — shared-weights and unrolled (own-weights) modes."""
+import numpy as _np_cpu
+from stack.rnn import StackRnn
+from rbm.matrix import np, convert
+from rbm.entity import EntityParams, TrainingParams
+from rbm.status import Status
+
+SENSORY_SIZE = 8
+H_SIZE       = 4
+T            = 16    # sequence length
+NUM_SEQ      = 10    # sequences in the batch
+WORK_DIR     = "../../results"
+
+_PARAMS = TrainingParams(learning_rate=0.05, momentum=0.5, num_epochs=5,
+                         do_rao_blackwell=True)
+
+
+def _make_sequences() -> _np_cpu.ndarray:
+    rng  = _np_cpu.random.RandomState(0)
+    base = (rng.rand(NUM_SEQ, SENSORY_SIZE) > 0.5).astype(_np_cpu.float64)
+    return _np_cpu.stack([base] * T, axis=1)   # (NUM_SEQ, T, SENSORY_SIZE)
+
+
+def test_rnn_shared():
+    """Shared-weights mode: one entity reused at every time step."""
+    seqs = _make_sequences()
+    rnn  = StackRnn("test_rnn_shared", WORK_DIR)
+    rnn.append(StackRnn.make_layer("layer0", SENSORY_SIZE, H_SIZE,
+                                   EntityParams(), _PARAMS))
+    rnn.state_init(0.01)
+
+    assert rnn.is_shared
+    assert rnn.sensory_size() == SENSORY_SIZE
+    assert rnn.h_size()       == H_SIZE
+
+    rnn.train(seqs)
+
+    rnn.reset(batch_size=1)
+    for t in range(T):
+        h = rnn.step(np.array(seqs[0, t][None, :]))
+        assert h.shape == (1, H_SIZE), f"bad shape at t={t}"
+
+    recon = rnn.reconstruct(h)
+    assert recon.shape == (1, SENSORY_SIZE)
+
+    print(f"Original : {convert(np.array(seqs[0, -1][None, :]))}")
+    print(f"Recon    : {convert(recon)}")
+    print("test_rnn_shared: [passed]")
+
+
+def test_rnn_unrolled():
+    """Unrolled mode: T entities, one per time step, each with own weights."""
+    seqs = _make_sequences()
+    rnn  = StackRnn("test_rnn_unrolled", WORK_DIR)
+    for layer in StackRnn.make_unrolled(T, SENSORY_SIZE, H_SIZE,
+                                        EntityParams(), _PARAMS):
+        rnn.append(layer)
+    rnn.state_init(0.01)
+
+    assert not rnn.is_shared
+    assert rnn.num_layers() == T
+
+    rnn.train(seqs)
+
+    # Inference — _t advances through all T entities
+    rnn.reset(batch_size=1)
+    for t in range(T):
+        h = rnn.step(np.array(seqs[0, t][None, :]))
+        assert h.shape == (1, H_SIZE), f"bad shape at t={t}"
+        assert rnn._t == t + 1
+
+    recon = rnn.reconstruct(h)
+    assert recon.shape == (1, SENSORY_SIZE)
+
+    # next_entity wraps around modularly
+    rnn._t = T + 3
+    assert rnn.next_entity() is rnn.from_index(3).entity
+
+    print("test_rnn_unrolled: [passed]")
+
+
+def test_rnn_save_load():
+    seqs = _make_sequences()
+    rnn  = StackRnn("test_rnn_sl", WORK_DIR)
+    for layer in StackRnn.make_unrolled(T, SENSORY_SIZE, H_SIZE,
+                                        EntityParams(), _PARAMS):
+        rnn.append(layer)
+    rnn.state_init(0.01)
+    rnn.train(seqs)
+    rnn.state_save()
+    rnn.state_load()
+    print("test_rnn_save_load: [passed]")
+
+
+if __name__ == "__main__":
+    test_rnn_shared()
+    test_rnn_unrolled()
+    test_rnn_save_load()
+    print("All RNN tests passed.")

tests/test_rnn_helper.py

@@ -0,0 +1,91 @@
+import numpy as np
+import pytest
+from stack.rnn_helper import shift_right, shift_left, shift_up, shift_down, ch2idx, idx2ch, VOCAB
+
+M = np.array([
+    [1, 2, 3],
+    [4, 5, 6],
+], dtype=float)
+
+
+def test_shift_right():
+    result = shift_right(M)
+    expected = np.array([
+        [0, 1, 2],
+        [0, 4, 5],
+    ], dtype=float)
+    assert np.array_equal(result, expected)
+    assert result.shape == M.shape
+
+
+def test_shift_left():
+    result = shift_left(M)
+    expected = np.array([
+        [2, 3, 0],
+        [5, 6, 0],
+    ], dtype=float)
+    assert np.array_equal(result, expected)
+    assert result.shape == M.shape
+
+
+def test_shift_up():
+    result = shift_up(M)
+    expected = np.array([
+        [4, 5, 6],
+        [0, 0, 0],
+    ], dtype=float)
+    assert np.array_equal(result, expected)
+    assert result.shape == M.shape
+
+
+def test_shift_down():
+    result = shift_down(M)
+    expected = np.array([
+        [0, 0, 0],
+        [1, 2, 3],
+    ], dtype=float)
+    assert np.array_equal(result, expected)
+    assert result.shape == M.shape
+
+
+def test_ch2idx_known():
+    assert ch2idx(' ') == 0
+    assert ch2idx('A') == 4
+    assert ch2idx('Z') == 29
+    assert ch2idx('0') == 30
+    assert ch2idx('9') == 39
+
+
+def test_idx2ch_known():
+    assert idx2ch(0) == ' '
+    assert idx2ch(4) == 'A'
+    assert idx2ch(29) == 'Z'
+    assert idx2ch(30) == '0'
+    assert idx2ch(39) == '9'
+
+
+def test_ch2idx_invalid():
+    with pytest.raises(KeyError):
+        ch2idx('a')
+
+
+def test_idx2ch_invalid():
+    with pytest.raises(IndexError):
+        idx2ch(len(VOCAB))
+
+
+def test_ch2idx_idx2ch_roundtrip():
+    for i, ch in enumerate(VOCAB):
+        assert ch2idx(ch) == i
+        assert idx2ch(i) == ch
+
+
+if __name__ == "__main__":
+    test_shift_right()
+    test_shift_left()
+    test_shift_up()
+    test_shift_down()
+    test_ch2idx_known()
+    test_idx2ch_known()
+    test_ch2idx_idx2ch_roundtrip()
+    print("All tests passed.")

tests/test_sub_image.py

@@ -0,0 +1,26 @@
+from image.sub_image import SubImageExtract
+from rbm.matrix import np
+
+
+def test_sub_image_rgb():
+	n, c, w, h = 2, 3, 8, 8
+	img = np.reshape(np.linspace(1, n*c*w*h, num=n*c*w*h, dtype=np.float64), (n, c, w, h))
+	sub = SubImageExtract(4, 4, 1, 1)
+	result = sub(img, 3)
+	assert result.shape[1] == 3
+	assert result.shape[2] == 4 and result.shape[3] == 4
+
+
+def test_sub_image_gray():
+	n, w, h = 2, 8, 8
+	img = np.reshape(np.linspace(1, n*w*h, num=n*w*h, dtype=np.float64), (n, w, h))
+	sub = SubImageExtract(4, 4, 1, 1)
+	result = sub(img, 1)
+	assert result.shape[1] == 1
+	assert result.shape[2] == 4 and result.shape[3] == 4
+
+
+if __name__ == "__main__":
+	test_sub_image_rgb()
+	test_sub_image_gray()
+	print("Test: [passed]")

tests/test_xor.py

@@ -0,0 +1,46 @@
+import os.path
+
+from rbm.layer import Layer
+from rbm.status import Status
+from rbm.train import train
+from rbm.matrix import Mat, np
+from rbm.entity import EntityParams, TrainingParams
+
+WORK_DIR = "../../results"
+USE_OPTIMIZER = True
+
+def xor():
+	# Create params
+	entity_params = EntityParams()
+	training_params = TrainingParams()
+	entity_params.do_rao_blackwell = True
+	entity_params.num_gibbs_samples = 3
+
+	# Create layer
+	layer = Layer("Layer_0", (3, 1, 0, 16), entity_params, training_params)
+
+	# Init weights
+	layer.init(0.01)
+
+	# Load weights (if exists)
+	layer.load(os.path.join(WORK_DIR, "xor_layer0_state.npz"))
+
+	# Prepare training data
+	training_batch = Mat([[0,1,1], [0,0,0], [1,1,0], [1,0,1]], dtype=np.float64)
+
+	# Train layer
+	train(layer.entity, training_batch, Status())
+
+	# Save weights
+	layer.save(os.path.join(WORK_DIR, "xor_layer0_state.npz"))
+
+	# Test with test data
+	test_batch = Mat([[0,0,0], [0,1,0], [1,0,0], [1,1,0]], dtype=np.float64)
+	for pattern in test_batch:
+		h = layer.entity.forward(pattern)
+		v = layer.entity.reconstruct(h)
+		print(f"P{pattern} : {v}")
+
+if __name__ == "__main__":
+	xor()
+	print("Test: [passed]")