|
10 | 10 | "cell_type": "code", |
11 | 11 | "source": "%matplotlib inline\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport time\nfrom smlr import Surrogate, StrengthDataset, StrengthSample\nfrom smlr.metrics import normalized_l2\n", |
12 | 12 | "metadata": {}, |
13 | | - "id": "5342efbb" |
| 13 | + "id": "5342efbb", |
| 14 | + "outputs": [], |
| 15 | + "execution_count": null |
14 | 16 | }, |
15 | 17 | { |
16 | 18 | "cell_type": "markdown", |
|
22 | 24 | "cell_type": "code", |
23 | 25 | "source": "def strength_function(params, energy):\n \"\"\"Multi-peak strength function with parameter-dependent features.\n \n Parameters control:\n - params[0]: Controls peak 1 position and peak 2 amplitude\n - params[1]: Controls peak 2 position and global width\n \"\"\"\n p1, p2 = params\n \n # Peak 1: Low energy\n e1 = 5 + 4 * p1\n s1 = 0.8\n \n # Peak 2: High energy, amplitude depends on p1\n e2 = 15 + 5 * p2\n s2 = 0.6 + 0.4 * p1\n \n # Peak 3: Small satellite\n e3 = 25 + 2 * p2\n s3 = 0.2\n \n # Width depends on p2\n gamma = 1.5 + 0.8 * p2\n \n def lorentzian(e, e0, g, s):\n return s * g / np.pi / ((e - e0)**2 + g**2)\n \n return (lorentzian(energy, e1, gamma, s1) +\n lorentzian(energy, e2, gamma, s2) +\n lorentzian(energy, e3, gamma, s3))\n\n# Energy grid\nenergy = np.linspace(0, 40, 250)\n\n# Training data: 6x6 grid in [0,1]^2\ntrain_p1 = np.linspace(0, 1, 6)\ntrain_p2 = np.linspace(0, 1, 6)\n\ntrain_samples = []\nfor p1 in train_p1:\n for p2 in train_p2:\n params = np.array([p1, p2])\n spectrum = strength_function(params, energy)\n train_samples.append(StrengthSample(params, energy, spectrum))\n\ntrain_dataset = StrengthDataset(train_samples)\nprint(f\"Training samples: {len(train_samples)}\")\n\n# Test data: 10 random interpolation points\nrng = np.random.default_rng(42)\ntest_interp_params = rng.uniform(0.1, 0.9, size=(10, 2))\n\n# Extrapolation test: points outside training range\ntest_extrap_params = np.array([\n [1.1, 0.5], # p1 extrapolated\n [0.5, 1.1], # p2 extrapolated\n [1.1, 1.1], # Both extrapolated\n [-0.1, 0.5], # Negative p1\n [0.5, -0.1], # Negative p2\n])\n\nprint(f\"Interpolation test points: {len(test_interp_params)}\")\nprint(f\"Extrapolation test points: {len(test_extrap_params)}\")\n", |
24 | 26 | "metadata": {}, |
25 | | - "id": "20c73ab3" |
| 27 | + "id": "20c73ab3", |
| 28 | + "outputs": [], |
| 29 | + "execution_count": null |
26 | 30 | }, |
27 | 31 | { |
28 | 32 | "cell_type": "markdown", |
|
34 | 38 | "cell_type": "code", |
35 | 39 | "source": "def evaluate_interpolation(model, test_params, energy):\n \"\"\"Compute errors on test points.\"\"\"\n errors = []\n for params in test_params:\n result = model.predict(params, energy)\n truth = strength_function(params, energy)\n errors.append(normalized_l2(result.spectrum, truth, energy))\n return np.array(errors)\n\n# Configure models to compare\nconfigs = {\n 'Regression (n=3)': {'backend': 'regression', 'n_components': 3, 'width_mode': 'global', 'random_state': 42},\n 'Regression (n=5)': {'backend': 'regression', 'n_components': 5, 'width_mode': 'global', 'random_state': 42},\n 'Regression (n=7)': {'backend': 'regression', 'n_components': 7, 'width_mode': 'global', 'random_state': 42},\n 'PMM (poles=5)': {'backend': 'pmm', 'n_poles': 5, 'retain': 0.8},\n 'PMM (poles=7)': {'backend': 'pmm', 'n_poles': 7, 'retain': 0.8},\n 'PMM (poles=10)': {'backend': 'pmm', 'n_poles': 10, 'retain': 0.8},\n}\n\nresults = {}\nfor name, cfg in configs.items():\n backend = cfg.pop('backend')\n \n # Time training\n t_start = time.time()\n model = Surrogate(backend, **cfg)\n model.fit(train_dataset)\n train_time = time.time() - t_start\n \n # Evaluate\n interp_errors = evaluate_interpolation(model, test_interp_params, energy)\n \n results[name] = {\n 'model': model,\n 'train_time': train_time,\n 'interp_mean': np.mean(interp_errors),\n 'interp_max': np.max(interp_errors),\n 'interp_errors': interp_errors,\n }\n \n print(f\"{name:20s}: mean={np.mean(interp_errors):.4f}, max={np.max(interp_errors):.4f}, time={train_time:.3f}s\")\n", |
36 | 40 | "metadata": {}, |
37 | | - "id": "542fcb5b" |
| 41 | + "id": "542fcb5b", |
| 42 | + "outputs": [], |
| 43 | + "execution_count": null |
38 | 44 | }, |
39 | 45 | { |
40 | 46 | "cell_type": "code", |
41 | 47 | "source": "# Visualize interpolation comparison\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\n\n# Left: Box plot of errors\nax = axes[0]\nreg_errors = [results[k]['interp_errors'] for k in results if 'Regression' in k]\npmm_errors = [results[k]['interp_errors'] for k in results if 'PMM' in k]\nreg_labels = [k for k in results if 'Regression' in k]\npmm_labels = [k for k in results if 'PMM' in k]\n\npositions_reg = np.arange(len(reg_errors)) * 2\npositions_pmm = np.arange(len(pmm_errors)) * 2 + 0.8\n\nbp1 = ax.boxplot(reg_errors, positions=positions_reg, widths=0.6, patch_artist=True)\nbp2 = ax.boxplot(pmm_errors, positions=positions_pmm, widths=0.6, patch_artist=True)\n\nfor patch in bp1['boxes']:\n patch.set_facecolor('lightblue')\nfor patch in bp2['boxes']:\n patch.set_facecolor('lightcoral')\n\nax.set_xticks((positions_reg + positions_pmm) / 2)\nax.set_xticklabels(['Low', 'Med', 'High'])\nax.set_xlabel('Model Complexity')\nax.set_ylabel('Normalized L2 Error')\nax.set_title('Interpolation Error Distribution')\nax.legend([bp1['boxes'][0], bp2['boxes'][0]], ['Regression', 'PMM'])\n\n# Right: Training time comparison\nax = axes[1]\nnames = list(results.keys())\ntimes = [results[k]['train_time'] for k in names]\ncolors = ['steelblue' if 'Regression' in k else 'coral' for k in names]\nax.barh(names, times, color=colors)\nax.set_xlabel('Training Time (s)')\nax.set_title('Training Time Comparison')\n\nplt.tight_layout()\nfig\n", |
42 | 48 | "metadata": {}, |
43 | | - "id": "f9a59c06" |
| 49 | + "id": "f9a59c06", |
| 50 | + "outputs": [], |
| 51 | + "execution_count": null |
44 | 52 | }, |
45 | 53 | { |
46 | 54 | "cell_type": "markdown", |
|
52 | 60 | "cell_type": "code", |
53 | 61 | "source": "# Pick best of each type based on interpolation\nbest_reg = 'Regression (n=5)'\nbest_pmm = 'PMM (poles=7)'\n\n# Evaluate extrapolation\nprint(\"Extrapolation Errors:\")\nprint(\"-\" * 50)\nfor name in [best_reg, best_pmm]:\n model = results[name]['model']\n extrap_errors = []\n for params in test_extrap_params:\n try:\n result = model.predict(params, energy)\n truth = strength_function(params, energy)\n err = normalized_l2(result.spectrum, truth, energy)\n except:\n err = np.nan\n extrap_errors.append(err)\n \n results[name]['extrap_errors'] = np.array(extrap_errors)\n print(f\"{name}: {extrap_errors}\")\n", |
54 | 62 | "metadata": {}, |
55 | | - "id": "e406eabb" |
| 63 | + "id": "e406eabb", |
| 64 | + "outputs": [], |
| 65 | + "execution_count": null |
56 | 66 | }, |
57 | 67 | { |
58 | 68 | "cell_type": "code", |
59 | 69 | "source": "# Visualize extrapolation at one point\nextrap_point = np.array([1.1, 0.5])\ntruth = strength_function(extrap_point, energy)\n\nfig, axes = plt.subplots(1, 2, figsize=(12, 4))\n\nfor ax, name in zip(axes, [best_reg, best_pmm]):\n model = results[name]['model']\n try:\n result = model.predict(extrap_point, energy)\n pred = result.spectrum\n err = normalized_l2(pred, truth, energy)\n \n ax.plot(energy, truth, 'k-', lw=2, label='Ground truth')\n ax.plot(energy, pred, 'r--', lw=2, label=f'{name}')\n ax.set_title(f'{name}\\nExtrap. error = {err:.4f}')\n except Exception as e:\n ax.text(0.5, 0.5, f'Error: {e}', transform=ax.transAxes)\n ax.set_title(f'{name}\\nFailed')\n \n ax.set_xlabel('Energy')\n ax.set_ylabel('Strength')\n ax.legend()\n ax.set_xlim(0, 40)\n\nplt.tight_layout()\nfig\n", |
60 | 70 | "metadata": {}, |
61 | | - "id": "92d09437" |
| 71 | + "id": "92d09437", |
| 72 | + "outputs": [], |
| 73 | + "execution_count": null |
62 | 74 | }, |
63 | 75 | { |
64 | 76 | "cell_type": "markdown", |
|
70 | 82 | "cell_type": "code", |
71 | 83 | "source": "def compute_sum_rule(spectrum, energy):\n \"\"\"Compute integrated spectral weight.\"\"\"\n return np.trapezoid(spectrum, energy)\n\n# Compute sum rules for training data\ntrain_sum_rules = [compute_sum_rule(s.strength, s.energy) for s in train_dataset.samples]\nmean_sum_rule = np.mean(train_sum_rules)\nprint(f\"Training data mean sum rule: {mean_sum_rule:.4f}\")\n\n# Check sum rule preservation on test points\nfig, ax = plt.subplots(figsize=(8, 5))\n\nfor name in [best_reg, best_pmm]:\n model = results[name]['model']\n predicted_srs = []\n true_srs = []\n \n for params in test_interp_params:\n result = model.predict(params, energy)\n truth = strength_function(params, energy)\n \n predicted_srs.append(compute_sum_rule(result.spectrum, energy))\n true_srs.append(compute_sum_rule(truth, energy))\n \n relative_error = (np.array(predicted_srs) - np.array(true_srs)) / np.array(true_srs) * 100\n \n results[name]['sum_rule_error'] = relative_error\n ax.scatter(range(len(relative_error)), relative_error, label=name, s=100, alpha=0.7)\n\nax.axhline(0, color='k', linestyle='--', alpha=0.5)\nax.set_xlabel('Test point index')\nax.set_ylabel('Sum rule error (%)')\nax.set_title('Sum Rule Preservation')\nax.legend()\nax.set_ylim(-10, 10)\nfig\n", |
72 | 84 | "metadata": {}, |
73 | | - "id": "39da1378" |
| 85 | + "id": "39da1378", |
| 86 | + "outputs": [], |
| 87 | + "execution_count": null |
74 | 88 | }, |
75 | 89 | { |
76 | 90 | "cell_type": "markdown", |
|
82 | 96 | "cell_type": "code", |
83 | 97 | "source": "# Estimate parameter counts\ndef estimate_params(backend, config, n_train, n_dim=2):\n \"\"\"Rough estimate of effective parameters.\"\"\"\n if backend == 'regression':\n nc = config.get('n_components', 5)\n # Poles + strengths + widths, each regressed from params\n return nc * 3 * (n_dim + 1) # Linear regression coefficients\n else: # pmm\n np_ = config.get('n_poles', 10)\n # Matrix elements scale as n_poles^2 * n_dim\n return np_**2 + np_ * n_dim\n\n# Summary plot\nfig, ax = plt.subplots(figsize=(8, 5))\n\nmarkers = {'Regression': 'o', 'PMM': 's'}\ncolors_backend = {'Regression': 'steelblue', 'PMM': 'coral'}\n\nfor name, res in results.items():\n backend = 'Regression' if 'Regression' in name else 'PMM'\n config = configs[name] if name in configs else {}\n n_params = estimate_params(backend.lower(), config, len(train_samples))\n \n ax.scatter(n_params, res['interp_mean'], \n marker=markers[backend], c=colors_backend[backend],\n s=150, alpha=0.8, label=name if 'n=3' in name or 'poles=5' in name else '')\n ax.annotate(name.split('(')[1].rstrip(')'), \n (n_params, res['interp_mean']),\n textcoords='offset points', xytext=(5, 5), fontsize=8)\n\nax.set_xlabel('Estimated Number of Parameters')\nax.set_ylabel('Mean Interpolation Error')\nax.set_title('Parameter Efficiency')\nax.legend(['Regression', 'PMM'], loc='upper right')\nax.set_yscale('log')\nfig\n", |
84 | 98 | "metadata": {}, |
85 | | - "id": "7e369c08" |
| 99 | + "id": "7e369c08", |
| 100 | + "outputs": [], |
| 101 | + "execution_count": null |
86 | 102 | }, |
87 | 103 | { |
88 | 104 | "cell_type": "markdown", |
|
94 | 110 | "cell_type": "code", |
95 | 111 | "source": "# Generate summary table\nimport pandas as pd\n\nsummary_data = []\nfor name in [best_reg, best_pmm]:\n res = results[name]\n backend = 'Regression' if 'Regression' in name else 'PMM'\n \n summary_data.append({\n 'Backend': backend,\n 'Configuration': name.split('(')[1].rstrip(')'),\n 'Mean Interp. Error': f\"{res['interp_mean']:.4f}\",\n 'Max Interp. Error': f\"{res['interp_max']:.4f}\",\n 'Training Time (s)': f\"{res['train_time']:.3f}\",\n 'Sum Rule RMSE (%)': f\"{np.sqrt(np.mean(res['sum_rule_error']**2)):.2f}\",\n })\n\nprint(\"=\" * 70)\nprint(\"BACKEND COMPARISON SUMMARY\")\nprint(\"=\" * 70)\nfor d in summary_data:\n print(f\"\\n{d['Backend']} ({d['Configuration']}):\")\n for k, v in d.items():\n if k not in ['Backend', 'Configuration']:\n print(f\" {k}: {v}\")\n", |
96 | 112 | "metadata": {}, |
97 | | - "id": "4d38320a" |
| 113 | + "id": "4d38320a", |
| 114 | + "outputs": [], |
| 115 | + "execution_count": null |
98 | 116 | } |
99 | 117 | ], |
100 | 118 | "metadata": { |
|
0 commit comments