From ca1d4bffe97647ec9016a1016df6ca904513b44f Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Fri, 15 Aug 2025 12:12:21 +0200 Subject: [PATCH] Add MATLAB comparison framework for LFM kernel validation --- .../2025-08-15_matlab-comparison-framework.md | 76 +++++ scripts/compare_with_matlab.py | 318 ++++++++++++++++++ 2 files changed, 394 insertions(+) create mode 100644 backlog/features/2025-08-15_matlab-comparison-framework.md create mode 100644 scripts/compare_with_matlab.py diff --git a/backlog/features/2025-08-15_matlab-comparison-framework.md b/backlog/features/2025-08-15_matlab-comparison-framework.md new file mode 100644 index 00000000..f9db5675 --- /dev/null +++ b/backlog/features/2025-08-15_matlab-comparison-framework.md @@ -0,0 +1,76 @@ +--- +id: "matlab-comparison-framework" +title: "Create MATLAB comparison framework for LFM kernel validation" +status: "In Progress" +priority: "High" +created: "2025-08-15" +last_updated: "2025-08-15" +owner: "Neil Lawrence" +github_issue: "" +dependencies: "lfm-kernel-code-review" +tags: +- lfm +- kernel +- validation +- matlab +- comparison +--- + +# Create MATLAB comparison framework for LFM kernel validation + +## Description +Create a comprehensive comparison framework to validate our GPy LFM kernel implementation against the MATLAB reference implementation in GPmat. + +## Background +- We have analyzed the complete MATLAB implementation (SIM, DISIM kernels) +- Need to validate our GPy implementation against the reference +- Comparison framework will ensure mathematical correctness and numerical accuracy +- Will help catch implementation errors and validate parameter handling + +## Implementation Tasks +- [x] Create MATLAB comparison script (`scripts/compare_with_matlab.py`) +- [ ] Test MATLAB script with existing GPmat installation +- [ ] Create standard test cases for SIM and DISIM kernels +- [ ] Implement GPy computation integration in comparison script +- [ ] Add parameter validation and constraint testing +- [ ] Create visualization tools for comparison results +- [ ] Add cross-kernel computation validation +- [ ] Document comparison methodology and tolerance standards + +## Test Cases to Implement +- [ ] Basic SIM kernel with standard parameters +- [ ] SIM kernel with fast decay and no delay +- [ ] Basic DISIM kernel with hierarchical parameters +- [ ] Edge cases (zero delay, extreme decay values) +- [ ] Multi-output scenarios +- [ ] Cross-kernel computations (SIM × RBF, DISIM × SIM) + +## Acceptance Criteria +- [ ] MATLAB comparison script runs successfully +- [ ] Standard test cases produce consistent results +- [ ] Comparison framework can detect implementation errors +- [ ] Tolerance standards defined and documented +- [ ] Results visualization and reporting implemented +- [ ] Framework integrated into development workflow + +## Implementation Notes +- Use scipy.io for loading MATLAB .mat files +- Support both MATLAB and Octave as reference implementations +- Implement robust error handling for missing dependencies +- Create standardized test data sets for reproducible comparisons +- Consider numerical precision differences between platforms + +## Related +- Backlog: lfm-kernel-code-review +- Backlog: implement-lfm-kernel-core +- MATLAB Implementation: ~/lawrennd/GPmat/matlab/ + +## Progress Updates + +### 2025-08-15 +Created initial MATLAB comparison framework: +- Implemented `MATLABComparison` class with automatic MATLAB/Octave detection +- Created test case generation for SIM and DISIM kernels +- Added result comparison with tolerance checking +- Framework ready for integration with GPy implementation +- Script can generate MATLAB code dynamically for different test cases diff --git a/scripts/compare_with_matlab.py b/scripts/compare_with_matlab.py new file mode 100644 index 00000000..d0a1b697 --- /dev/null +++ b/scripts/compare_with_matlab.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python3 +""" +Comparison script for LFM kernel implementation against MATLAB reference. + +This script generates test cases and compares results between GPy and MATLAB +implementations to validate our LFM kernel implementation. +""" + +import numpy as np +import subprocess +import tempfile +import os +import json +from pathlib import Path + + +class MATLABComparison: + """Compare GPy LFM kernel results with MATLAB reference.""" + + def __init__(self, matlab_path=None): + """Initialize with path to MATLAB/Octave executable.""" + self.matlab_path = matlab_path or self._find_matlab() + self.temp_dir = tempfile.mkdtemp() + + def _find_matlab(self): + """Try to find MATLAB or Octave executable.""" + # Try common MATLAB paths + matlab_paths = [ + '/Applications/MATLAB_R2023b.app/bin/matlab', # macOS + '/usr/local/bin/matlab', # Linux + 'matlab', # In PATH + ] + + # Try Octave as fallback + octave_paths = [ + '/usr/local/bin/octave', # macOS + '/usr/bin/octave', # Linux + 'octave', # In PATH + ] + + for path in matlab_paths + octave_paths: + try: + result = subprocess.run([path, '--version'], + capture_output=True, text=True, timeout=5) + if result.returncode == 0: + print(f"Found MATLAB/Octave: {path}") + return path + except (subprocess.TimeoutExpired, FileNotFoundError): + continue + + raise RuntimeError("Could not find MATLAB or Octave executable") + + def create_matlab_script(self, test_case): + """Create MATLAB script for kernel computation.""" + script = f""" +% MATLAB script for LFM kernel comparison +% Generated test case: {test_case['name']} + +% Add GPmat to path (assuming it's in ~/lawrennd/GPmat) +addpath('~/lawrennd/GPmat/matlab'); + +% Test parameters +{self._matlab_params(test_case)} + +% Create kernel +{self._matlab_kernel_creation(test_case)} + +% Compute kernel matrix +{self._matlab_kernel_computation(test_case)} + +% Save results +results = struct(); +results.K = K; +results.params = {self._matlab_params_struct(test_case)}; +results.test_case = '{test_case['name']}'; + +save('{os.path.join(self.temp_dir, 'matlab_results.mat')}', 'results'); +fprintf('MATLAB computation completed\\n'); +""" + return script + + def _matlab_params(self, test_case): + """Generate MATLAB parameter definitions.""" + params = test_case.get('params', {}) + lines = [] + for key, value in params.items(): + if isinstance(value, (int, float)): + lines.append(f"{key} = {value};") + elif isinstance(value, list): + lines.append(f"{key} = [{', '.join(map(str, value))}];") + return '\n'.join(lines) + + def _matlab_kernel_creation(self, test_case): + """Generate MATLAB kernel creation code.""" + kernel_type = test_case.get('kernel_type', 'sim') + if kernel_type == 'sim': + return """ +% Create SIM kernel +kern = kernCreate(t, 'sim'); +kern.decay = decay; +kern.delay = delay; +kern.variance = variance; +kern.inverseWidth = inverseWidth; +kern = kernParamInit(kern); +""" + elif kernel_type == 'disim': + return """ +% Create DISIM kernel +kern = kernCreate(t, 'disim'); +kern.decay = decay; +kern.di_decay = di_decay; +kern.variance = variance; +kern.di_variance = di_variance; +kern.inverseWidth = inverseWidth; +kern.rbf_variance = rbf_variance; +kern = kernParamInit(kern); +""" + else: + raise ValueError(f"Unknown kernel type: {kernel_type}") + + def _matlab_kernel_computation(self, test_case): + """Generate MATLAB kernel computation code.""" + return """ +% Compute kernel matrix +K = kernCompute(kern, t); +""" + + def _matlab_params_struct(self, test_case): + """Generate MATLAB parameter structure.""" + params = test_case.get('params', {}) + param_str = ', '.join([f"'{k}', {v}" for k, v in params.items()]) + return f"struct({param_str})" + + def run_matlab_computation(self, test_case): + """Run MATLAB computation for given test case.""" + script = self.create_matlab_script(test_case) + + # Write script to temporary file + script_path = os.path.join(self.temp_dir, 'compute_kernel.m') + with open(script_path, 'w') as f: + f.write(script) + + # Run MATLAB + try: + if 'octave' in self.matlab_path.lower(): + cmd = [self.matlab_path, '--no-gui', '--silent', script_path] + else: + cmd = [self.matlab_path, '-batch', f"run('{script_path}')"] + + result = subprocess.run(cmd, capture_output=True, text=True, timeout=30) + + if result.returncode != 0: + print(f"MATLAB error: {result.stderr}") + return None + + # Load results + import scipy.io + results_path = os.path.join(self.temp_dir, 'matlab_results.mat') + if os.path.exists(results_path): + matlab_data = scipy.io.loadmat(results_path) + return matlab_data['results'][0, 0] + else: + print("MATLAB results file not found") + return None + + except subprocess.TimeoutExpired: + print("MATLAB computation timed out") + return None + except Exception as e: + print(f"Error running MATLAB: {e}") + return None + + def create_test_cases(self): + """Create standard test cases for comparison.""" + test_cases = [ + { + 'name': 'sim_basic', + 'kernel_type': 'sim', + 'params': { + 'decay': 1.0, + 'delay': 0.1, + 'variance': 1.0, + 'inverseWidth': 1.0 + }, + 't': np.linspace(0, 5, 20).reshape(-1, 1) + }, + { + 'name': 'sim_fast_decay', + 'kernel_type': 'sim', + 'params': { + 'decay': 2.0, + 'delay': 0.0, + 'variance': 1.0, + 'inverseWidth': 0.5 + }, + 't': np.linspace(0, 3, 15).reshape(-1, 1) + }, + { + 'name': 'disim_basic', + 'kernel_type': 'disim', + 'params': { + 'decay': 1.0, + 'di_decay': 0.5, + 'variance': 1.0, + 'di_variance': 1.0, + 'inverseWidth': 1.0, + 'rbf_variance': 1.0 + }, + 't': np.linspace(0, 5, 20).reshape(-1, 1) + } + ] + return test_cases + + def compare_results(self, matlab_results, gpy_results, test_case): + """Compare MATLAB and GPy results.""" + if matlab_results is None or gpy_results is None: + return {'status': 'error', 'message': 'Missing results'} + + # Extract kernel matrices + matlab_K = matlab_results['K'] + gpy_K = gpy_results['K'] + + # Basic shape check + if matlab_K.shape != gpy_K.shape: + return { + 'status': 'error', + 'message': f'Shape mismatch: MATLAB {matlab_K.shape} vs GPy {gpy_K.shape}' + } + + # Compute differences + abs_diff = np.abs(matlab_K - gpy_K) + rel_diff = abs_diff / (np.abs(matlab_K) + 1e-10) + + comparison = { + 'status': 'success', + 'test_case': test_case['name'], + 'shapes_match': True, + 'max_abs_diff': float(np.max(abs_diff)), + 'mean_abs_diff': float(np.mean(abs_diff)), + 'max_rel_diff': float(np.max(rel_diff)), + 'mean_rel_diff': float(np.mean(rel_diff)), + 'matlab_shape': matlab_K.shape, + 'gpy_shape': gpy_K.shape + } + + # Check if results are close enough + tolerance = 1e-6 + comparison['within_tolerance'] = comparison['max_abs_diff'] < tolerance + + return comparison + + def cleanup(self): + """Clean up temporary files.""" + import shutil + shutil.rmtree(self.temp_dir, ignore_errors=True) + + +def main(): + """Main comparison function.""" + print("LFM Kernel MATLAB Comparison") + print("=" * 40) + + # Initialize comparison framework + try: + comparator = MATLABComparison() + except RuntimeError as e: + print(f"Error: {e}") + print("Please ensure MATLAB or Octave is installed and in PATH") + return + + # Create test cases + test_cases = comparator.create_test_cases() + + results = [] + + for test_case in test_cases: + print(f"\nTesting: {test_case['name']}") + print("-" * 20) + + # Run MATLAB computation + print("Running MATLAB computation...") + matlab_results = comparator.run_matlab_computation(test_case) + + if matlab_results is None: + print("MATLAB computation failed") + continue + + # TODO: Run GPy computation (when implemented) + print("GPy computation not yet implemented") + gpy_results = None + + # Compare results + if gpy_results is not None: + comparison = comparator.compare_results(matlab_results, gpy_results, test_case) + results.append(comparison) + + if comparison['status'] == 'success': + print(f"✓ Shapes match: {comparison['shapes_match']}") + print(f"✓ Max abs diff: {comparison['max_abs_diff']:.2e}") + print(f"✓ Within tolerance: {comparison['within_tolerance']}") + else: + print(f"✗ Error: {comparison['message']}") + else: + print("Skipping comparison (GPy not implemented yet)") + + # Save results + results_file = 'matlab_comparison_results.json' + with open(results_file, 'w') as f: + json.dump(results, f, indent=2) + + print(f"\nResults saved to: {results_file}") + + # Cleanup + comparator.cleanup() + + +if __name__ == "__main__": + main()