# Generated by scripts/export_project_setup.py from _includes/linalg/week7_setup.md, _includes/projects/week7_setup.md
# Requires NumPy and Matplotlib. No network access is needed.
# Ferdige forsøksdata og visningshjelpere til SVD-uken.
# Bildedata lastes lokalt; selve matematikken ligger i SVD, trunkering og uskarphetstransformasjonen.
# Du trenger ikke endre det kodede portrettet for å undersøke de matematiske valgene.

import numpy as np
import matplotlib.pyplot as plt
import base64, zlib

# Original: https://wiki.math.ntnu.no/_media/imax3011/2025h/mm.gif
# Gråtoner, skalert til 96 x 96 for korte beregninger i nettleseren.
_portrait_data = "eJxtmndYW1eWwHE2maw32cl+k29SZpJNNpsyk5nUSeKJHfducMMFG2NcKDbNNNMxvfcOBkxHgAEhOqJKoIJ67713CXDbP/feJ4FJds8n6cHnp98597R77jNajVqlUMjlMrlUIpXKEJEjHwogSvenUqXR6TRqJbxPJpWIxRKRWCKBV7eIxSLkJZEAArwffkWlVqmUKq0WXpAvyqRuBUAdQof3gDtUCqXW6lzfeLJmNWhVCjnku8mil3yAd/OlwDoFYpRKDcgqrUajVgOFL/lQvVwuBxegQK3RaA32J8+fPX/x4sXTjXWXzaCRS8ViDxTwpFt8kQTyIQRRABYAjNNAPlgAogB4xQNWQNs1wHqN3uLcePLsOeQ/e7Luspj1erAMGUAB/PYFbHlLKpW73QoXr4UKVIgCuRTBQ74KUa/W6IwWm8Ppevp0w2W1ORxWs16ulEnEKoNOA+KwPQDb+Zt2qgBJq9Nu8WFQEc0qJDwajd4EqC6b1WzUG5UipUoi5FJW6VQyiSWQwyVIEO/8X/6moeCl1cIIAE9ApwO6Eok+tF2j1ut1Oi1fqVapQRboNHIFi82j4ohUoIPJZAhBqGRu3KZfZFAfjIDc7Wdgqlbr1qBERKFA3K9QqrVqidph0SolbIbG7HBqRTyV1ibX6DUSsVwgVgpZbK5YKpdshVQmRVJTLkWsh96AFkO+BipAXIQoRD5Vap3OtuEwyvhClca+5jSwMaNLdIbL5rIaNEalRKkUyEV8iVy6zSUyJKpyxHQPB5iqgXyNBkYTpron69VqvfOJy6xTaQz2tTWHTkRbHGof7Va7nHq5UC6hCCR8tYQllXkUeFwCqXLFyxSEokEEuBgIrDjIBuHQrTmtNqvRYHI4HHajlL44g6mqL52wmMV0IoWxQBNQpXK+UAEDAHyO2I6UpHsFm3wVYj+iAFKhCqhDrdZarVq1VmewOV0Og0pMxw/19+dmFuQz5ey5sZnFWQIPx5WL+VKoAPE+UrRymHmICzb5AOoJMBAYCpBKYBlKk0WjMhikHAXQomVTF+YnW8qSQtPD2laIIyOjU9NEPoEvlYokMiR7POnoLinlS9/AJNduikbj/tCBBRhsBtPamkk41T7LpEkYxKmhof70e36Rfvc60Y/7+qbn8Vw6qF0BaArAO1KRSLbJUyHJrVSo3L/C/rZNA3QWaGEai0FncdlVtLHm4oczo0uLmIEuVGO0/z3vwKKHXY/6p3A4jkApFrEZDK5EIWLRmHyxVObOC9V2QdwDasjDV0MNOoVcD36yWcR4VOH9xPv5zegxVHdP30h6RPrhkOLa9rrxGTyBL1ZIhDwuj89mkOlMcOXzJQoVUvRQhTtRkA4DBSjx1LFWrxLLZVK1zaZhLXTnBZ+6dbewFdXcicIsNJaXX05vamhvmJ5apImB75livUNDp6rsZo2Yz+VyuFLgDvVLQXLSjff4B1aAWqviiY0WLQc31ZaXcO5w4P1CVOvj0WXSFLo2u76ro7drcGiZy2dT+ASGymaQKE1GvU4hZrI4VIYYJshLPnSQhw8/QUtk0vVr6wq2WKVkTHfWlT2IjTt7IqoQhV7CU9mk+dqWzsmJsY72PiKHT54jjOL4cp1KIuDwhRIhi0Sn4slSmCVuuka1xXeLVsMjDK/KNHoRk8eZn+gqzUnJj4i5cSa9AU2gMRgSPrGrb5RHWuwbw5BY3OmpURRFKpfx6YRl4iqHxyBTiSt4klTtSXYNUlAazRYe9EpGQ2b3ktwmYbCoNEJzfnZCTFxcYEzRo1k2j8YRyTiYFbmMR1rlrTKp+KbHPf3LWo2cuTw/2ts/vcphrC4vEeZJcrBVwzrSqt0lhYA9oqgPf9BFFEqZPKVKPJNyMejW3ciY7JJhupiFJ8sVEgboaHIaj0em4rsKGuoGZS6jVkidfNzbOzhBYFDxWOz0CA8po21Zv40v67x01C9jVMmTWYymoUuHfAJDAirLMHialITCCEVKAYkjorAoyzRyS0puTf0YUwzaGhu/hB9Ho6fwOBx2dKRnXoPY70a7A7vpfx2n/4zPiduPpAKxzqIsv37i6PG7uTUN0/M4Hq61b5XLGJ6fRVMYBOJo7I2A27H3m2awlGUug0mfGuztQI1OjY6ie1sf8XVbTESH236gS6WTYe7suxmYzJIIZDpye9rf9lxOedQ7PIAdX1xAoyiCPvTE48Ul0vJsRc1jDKpjZIlKoc4sEInTj1E9gwN9mO4BdH9T7qJRD6H6LbfoPQkKiovZGXG9sHRKwufK5ubyP/N+UIdZJSwuERcnRgm85Ukmj8vhUEgUGpNGxK2wuRwBf74f09OJGZtBo9pRfaNoVE3ysEW3nQ9+0LujoNUoyKP1VR2jMwwhk0mgF36bOEHjSuRCoYRFWqYwiAq9UiYTCUArYNJWyWQajSEWz5W2NnZxePNDNY3d6NnxhvJ741bA1xsMHg2Qv6lKmhdfh2MLBRQWhzBHKQuaVIH7jCajViljUHgaYIhGKZWKBXwWk0al0ugsHh6TV9Iwzlxs669pbCcuDJbX3JtUykGXAUOS3uCGI3yoQquuOR054VxXLa6yiKNTjdU8k9lqsdjMJuhSoxHYpFMrpGKhgMtk0BkUCpNFnGnJqn88VF7Z0tMzLFzsLSrKS6pi6xC6wb0G5OV+qwlZgUnDTPLIKJUwOV7VLwfTlc1ut1ptFqvVYgQa9GByEwkFHDaDRqeQaHQiYbSsqrksqrBhcJLAxXbmZsQduMYyQjB0kfuK+Ade1JLxwuiwrB5s9+T8zGDRlM4Apjew/zrsYBlms9kEFqCUw77MYoAhiLDKWF4YrcyouHo+bwRDZVGHcpOTww4cJpmNW+bDqycOBr1GNF+UFnenqratcWR0qJeiM5qB+UCDw2a1ABUWk0GnlklEAjadRiYS8GQyiTtVmPjg8rXy4dllIrY2NjIxePdhvBWhb64CCDQffGrESzlBN69fPptS0DKOxqwqAdJmQ/g24CC4AJ1aLhELWHTKCn5uFkdawfWXZSTGRGRXDk1OthfEhkZc3e+zYJWpDTBcyBI8C0AySrBUGRpRFLjP2y+zuv7hMMNkMYMB1OmwgzBYoQajXqOQAvtBM1ucnpwlEvGopvLC8NC7Bf1jw40ZEZE3fU/eJFlYUqOHjyjQuflGnZDWExUx0xd6ws/3SkxxFUpsNpltTqfTbgMaLGDMBYmkkvJ5LDoZhx0fH8cuzk0MDnSk3w3J7VyYRxUXpIcF+N7jaRfFJrdjPHxPmkoF9NaY0NHVxohb5/YevXozDac3WO1OYL/dhpgPVqpVSAVcFo2wPD+OwYyMjY9PLgykx8QVN2MXp7q7Cu4GBqTxmSMSEzAXWYJ+q9T0ShoXVxF76z4elRV/7dD+3ceTZ9UW4HygwIEoMBtAlwL7OpdJJeKx4+gh9FBHXXNfZVx+YUNLR8/SxNhIeej1nPn+Xo0JwpGS2axfrVZOYQ5Xl0WfSRnuq0+LuOp9p5ejt4HshAuwQQUWk14DjhcCNpOyvDQ7gX482Ntc1VQcG1JaXT9dnthQ1VVbURyXXJiLBUUPxZNDnkanpC8Pj1REnzieghquj7sR85ijNNpAejpdiIdsZpDXaqVMLADth7SMHRvu7eloKCuIj86pKG+j9gXmVbQUZ2Zk3Y8uFZlM7vAajS/7v15Ln5rsLcwI8Pa94r3fv2ZJIFcbLFZgvwtJUPAVk0Ehk/JJM+jh0YVZTH93R2tNeX5mYkZp3QCWNNeU9yA6PuFBfFC/0W29wZM7nkaqo0wNdE0154T6/fSNfyOBJ5Cq9EaTxe60wdQHqQOaAxeHqcmO8798M6YN3dPaXFtZmp+altODwgyBEWm0PDYnMjzmHsFo3NSwjW/QEedbaNiRzoTAC365wyuPWwcWqUKVwWwzw84DD1GMnvzY4Ov+l+8lx/r5l3fXV5UUPUhJyS99PDM3UltQ14JuTUq5GZw4KDUAB5kQJ3l6NFiKUU/qHWMuj5cEB18LKumvPu17PTypvHNsiQv+VaeHnROTHRFV1NJLYZEXutPOJVQW5uampMcV5pegZvuLg6LSippbE2LComt5RkA3w3gZPbWFxILQRaVOlEZH+Z29VFDqG12eHh96IyileEoJB0elREDpqql6vEDhiMRCJh6VeiY8OzMjJTkp4879/PLcgsIHSRWFdeVROfdicEZIB3yTybDZg4DglpmjbWlhd6+f3JsR9GCFycSPt+aml08IRFKlWsKjYPu6u/tHsAsECoM02VOfciU8IT7pQW51aHBcZlVDSXpcUmZWQ0Ka//lJM+SbPPxNvGmFMduQlxaXcvfC8ehg9BKNTcai6srqhpZoPKlayiGM9zQ1t6GGBlFdKDTqYWVeQlhUfEp2Y2tkcExhd29FUWJsXHJBYsi1QILZbT8QpJG63xbqfF9hSlZ0qJ+vz50MFBo7i11aQtclp7bNM6QaOWu+v7G5d5YpFTMnytNzC3JSY6PvxaUUNvUkx2YUjZFGG9PCI0ITUyKDS1Sg8Vo2+UY332g0c6Z6UoNjszNunzsV9aC8vnNRufbMqZopSukhS3UK1mxnXe+iQCYVS+iTVTEJqfcjomJiUwof9idEpxalFk3PVd64HRB2/96DVasVUWCBLcsTW5BNZtlYbXFR7N2K8pCz90KjMrpIK9OYWbKKWlS9rDApWPOdKBxlMOniiYA7yfkFOfkpMbHxcSkFjW2p4cmZFYl3+mfjfc9eu3evWmkF1W5DdFiQPufOJotpqrRibigjPa3geph/YFxG0Pfvv7X7RuboSgdeY1Wy5sfZ+PtHrv79tT/sePW9X25FJxVXZadlFrfVpyQllpc0NaZM9p71CbgV2ywz25Ft2wK2bRMUM1gL2Ers9JLE6e6MnOqaiFi/vafef+Nvn//5zsxsEWluWW9TMJfIysa8HmxdyLn9f/J67dKjmenuqoKCsvaq5LTCzKre4YdVuDB/P/+Ifi6g26ECK9g1YByQlZitVn1zVFVX4vXG3oSU/LN/8PogvriLAI627EWcxiShr9DkQj5jcRrdXnf7h8M5E5TFodaqosravNDA4pK4xrmFypHoC1cuHkXzQT+3w45oe8m3wh1cvZSe15pyJwiTWzzSf/CTpCEiU67RK/lLK2qdcJXEECukPCZpbmYeOzg8jCUujnVWF1RXFvl538ktii8Z7a4695cff/lnlwr2W+giZNezIKOBBaxFTBT05/aWRJ57VPZwca5hkCngSrSgZ7JxZLWaTyExBPCcKxdwhRK5VMRhEcbay3LKqrK8vS+HNnbE59SlfvT2sc+8u82OTT7YNWxwfELMty/x7Pq+2pK4G3kN7QSeGBQe8JlZI6QuUlUyDpVM54kVag3YI7VqlULMY6xgHuZnFRWEnzns658zWJ+aF/T2nz4+eQdrgzMTVOEePNyLsJjtc4tmu6qxrKWt5WEbWa41wzusOhkHP0uR8OmrRBpbpFBpwXYEn7eIubSFgers9JxMvz3nYvvHJ9oy4s/85x9Phj8Q2+CO53ALEgQbkkp2xSBZSM9PG6CiaxpWpMjeYjerBavzWLqIQyXgiHShQgvHGXDukvAZhOmu4rSklJiL+8JRJAYTk3L72p7DyZGdMmTHRgRGAC4ACoCpaIS+sLgmqrCxal6sBbOPzWpU8gjjsyIpl7I4uUARKfRgJjKAIYXHIs8PN+fcC4uKuH4nt/zRGLYu9saV4MjUuEXhmhOKA7499iOxAGOmTVF+I22Mox8rGGcpjDAqBiVvBb0gNShos6OzOJ5SZwLu1CuFrFXceHdZ4m3/GzeDY5ITskszkwvP+QUH7slrUT8FO7ab73Ru40ORZAWHN3H15MJOglAL9kabUcEn9y3KzWr8QP/4JF6hNdocVoOCT1ue7GvMvHPpwhW/G1HR4QkZFSUhZy5cPPztf394lPfE5dwU+6bxiLjEZfmVj8aZlMqKaYYSbI12i1qIb+yj0pY6Syr7hwdAPZisFr2CT5kfbiuN9z/pc8Y3KCw6OjajvSTk0t6rR957xcsr9sWaC4jbSx4BEQe/uazdD2fw4Mhcltu9ItKCEFt0opHk2xVVOdc/++RmdQ5OKFbrDFo5jzzVUZsRdunkqfN+abdSCx5kVVYk+ATf+K9Xf+flFfI/62trHhWQ7HBsxsPpMixOMfn41vy8KgxdqjWCBihDnXr9C/+A777be9g/QyjmCBVgxGUtDD0sSgw6d+Kcb1Bm4oUk1MjjzvzinJ93vvFvb384+WwN8hEdbvYW3+lcl8zhpmqzirM68FyVzmgxaMgxgfs+33eruG1kmCaXckGVcel4TGdlTkLQ+VO+54sqjv75y1u51TVluYdf+/TtD/aEmtZdEO3mu6kuzxuIZbG0OCPuTmLVBF0K0kWvlgqZ2I4ePFcolUlFIt4qoW8JN95Tl592P+Ty2bPHMwt3/cfO//7oYEykzzteXm98HTD71O0eZBHQfuCotc2Ig5rQ1t29diYwLrdnhSNTg7FRIhaK5GDul4C5gcdlUxfaJuaH2irzMhPDAi6c/vIXnz2/P3dh5wdXbhz9vZfXH8NIG+vA8cjL5fHPVrSdSGUTon1PHbiU2TBDFciVCik4ksrgAxcRncbh0lepZMLSZF9rVVluavTty6c//s/Xdn/y959ePRG///2dXn9u0j9B4G7q2po7qi7ny4x1Wo3V568cikosx5A5YjFfJBh71DeyNFTTkH0nugm3ukqnrUwPdT18WJoVF3zd588H3z7x1ee/ez8y78irXl5xa+trW6BtbtmqZzjqO1DexxqyYnN7cAw+a3U88J+79v7tp2/+8sPXX3z+8Y/ZS0wOAzc50NVak5cSErD/9dN+X+z69Mv9fzr4kdfOUKlrmy+2eWU73+6YOfrdcd/E7M45Kh3X/u2bb/3xq2sXA4oj7379zjtv7Ejn8nirc2NDPQ3FyaHXv9zx10++/u6Tiwe8vLzebQTf/RX/V3i3Arhx0lP3fP5ZVGEHlkxdjvHa8crlbr1IsiFVVAR8/tFxvIDHXV0cG+quK7offPEPXq95nfr+rTOfeHn9ZcpuNTu2Mzct3ioCxHqHRcZM++KfcYWtM1QOD3/pqz8EYxgKmWSuL/LsoeujYBNmk+cnRgaaCuJu+by7w+uVIx99e/4dr/fGZGqrHdrssdtTvwjSU8WILodVSurOvxVU0DjF4EtlOZeC0hv7+IzRR8khvhfLBEKugLkyPTra35QbEXDs01d27PB6/evdv9txDtU1v4bIZm5Cvt29VTo2FwP5dvXjjIy9n2bUjVN4QnH7D3tv1/XUtXfn+J85lkYXcNlsKn4KM9hRmXj74sH3Xn3dy+uVDz/0ev0mqhm17uEjbDffLY7NGLgfB5Da00OOhVf0LdJYbEL4B/+4XfewDx1/7EwFlcWgs+grsyMDHTXZ4QHHfnn7lTff9PL6/ZteHxQ9riC5XirwaNnqzVvNFNmWZ9NjwsMLGtB4EpUyFvbl53uvX9l3KGOKwaJR6TTy3Fh/R0NeUvjFAz++9f07//7Wjp1vvP7lgyqc1eXpCvCC9J+1l8Z7FuBAfjeAofhGYm7LxDKFjp9sTr1zOzRrjMFms+k0Cnlpcrijvig17sahfd9+fHD/vr9/8tX+/3r1x5HNmt1sPlAcv+K78xO2JGVt2K0bqRUoLJXNJBIpDNiWJSI+l00hLc2gO+uKU+4Hnzyw/+vvT8Zlhnv/svczr3eZLjsSV09v9vAdW3h303AiHc/poEVc8gtNbxjCsUQsptQAZhqzWiYWskk4LKa3viT9fujF04d2/2PP3cqixPt3D3331lezdrvb51vwdch3/IbvTl27vi06Mjwus2WaLhQKZPAgb9LIpEIGfhrT01iSFh9y6dzxo/v3nM0qLWhBdYXeDC6dkdu3HAThQH7Nf1l6wH4jp7O2PikytxvPEYJN0WQDu7pEyiPNDLRWZydEBl48f+rkiSOHLufUdU6z2a33Sx4uq4F/Nl2/7sa/bHa/6kbgR4NwoLrqQUJ8MZoqEIplcpVep5aLKdiB6vykiJuXzp/3OXnyyMk9PoWPsCyZjljzsIkCj/qutW32Q/624Dq3NT6XTbXc8aij+e71qkkqW8BlssRKBY+A6apIjQo463P2jI/3qeMnTv9wqvoxWW6waDAr1D6jczt/7Tf2/1aF066XS2bT/xpY3j81tzg3s4KfH++szo0Lu3j0uPcZb58zp0+fu/qNL5bAUVkc9uUlPNridG2D/9Y/2xrey2Xou/7y/rHYwoq66vKaqsKM+IiAq2eOHjvh4+3jc/7CtVuh34cI+BIjaDLS+RWJZ6d1bfGB/B++p1/Adr3x7OmG8MTOf/ns0AX/gCtXLl++4ut98OTRo8dOnTp5+syVoJj4qJ8zDEql2e6w6RfENpg+L2vr/+M7Nk8GcNx98uL5s3Vl+Mevvb7z3a9/+PHnPYeOHj10+Bjge5865RsQkZqTnXKkxKFRWeDTGxLDChe+5vq1/Q7Hr4y3IY8K4UTqfPbi+fNnLkvFkT+++tq/vvHeX7//affeAwcPATl88vSlkKSc8tqa3HO1a1o1mIOtJu6KyeEuX9dLOuA7t20uDvd5Bor9yXPAf7rm7A3Z9+7OV/7lrU/+/vPuXw4cOnhg/4Ejp/3vZZW3dHQ9qvQbcKk18ARnli+rHJvtYRt/q/PD+dyOnDTg7daNZ89eAP88tfQnhR5873dvvvPFtz/v2bf/4MF9e/efuBKTX48aRqNR9aEMm8aAPMvQkRmezN7Eb6xvbGz5x3Mcg3R4Hna+QP4i4flT8cPcrKBd77/zwd/+sWvXrt17fv75lyOX44ubQc7OTg5W5VvNWuQ5iUlLJ1gdvzEeKED4DvecDo+rEG8yO54h/OfP10ktlaXpV//x8YeffvnNjz/u+mnX3kN+MaWtA9hlPH5muGhyzaBHnlUZtDKuHoTXXVcIGrw2Njb9A4+S8JhqMpuMRvsGNB7yDbM9DTVlSZf2ffvNdz/8CMzfe+JKfHn74NTyysoSdqhc6wBs5FGbVqOQWV/6B6CfPNmA9rv3dYebDx+1GW3rz9z8F09F3XmPWhpLEoMungJp733yxPmgBw2o4SkciUxYwLaPbYBzrUlvhE/X9Fq5yIFsu2vA9YjxW3zkEaoVxgmYYln3/D3Ii/+xN353eagf9agsMzkmMjLqXmRcRnXP0NgMbmV1lTiLydWsg5Dp9cizPD046OuBpQh/S+ybYxtMHYA32Vzr60+fI95/8aznQ69bI4OY4Y6GspKyioqKmqYe9PTs/DKBQqMQZwtDn6w57cj/PIAvGvRqhQI+zXe5/fMEyMbT/wW0HKz6"
portrait = np.frombuffer(zlib.decompress(base64.b64decode(_portrait_data)), dtype=np.uint8).reshape(96,96).astype(float)/255

# Like bildedimensjoner gjør parameterbudsjettet sammenlignbart på tvers av mønstre.
def challenge_images():
    n = 96
    text = np.zeros((n,n))
    # Fire enkle bokstavformer: H, I, L, I; ingen skrifttypeavhengighet.
    for col in (8, 30, 52, 74):
        text[24:72, col:col+5] = 1
    text[24:72, 20:25] = 1
    text[46:51, 8:25] = 1
    text[67:72, 52:69] = 1
    return {'portrett': portrait.copy(), 'tekst': text,
            'diagonal': np.eye(n),
            'tilfeldig': np.random.default_rng(7).random((n,n))}

# Samme gråtoneskala i alle bilder; visning endrer ikke matrisene brukt til feilberegning.
def show_images(images):
    fig, axes = plt.subplots(1, len(images), figsize=(3*len(images),3), squeeze=False)
    for ax, (name, value) in zip(axes[0], images.items()):
        ax.imshow(value, cmap='gray', vmin=0, vmax=1)
        ax.set_title(name); ax.axis('off')
    plt.tight_layout(); plt.show()

# Visningshjelper: behold k SVD-ledd.
def rank_image(A, k):
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    return (U[:, :k]*s[:k]) @ Vt[:k, :]

# Lag kjent signal, uskarpe data og én støyrealisasjon for et kontrollert inversjonsforsøk.
def blur_problem(noise_level=.005, seed=17):
    n = 80
    t = np.linspace(0, 1, n)
    truth = ((t>.2)&(t<.4)).astype(float) + .6*np.exp(-((t-.72)/.06)**2)
    # Hver rad beskriver hvordan én resultatkoordinat beregnes fra signalverdier i nærliggende posisjoner.
    H = np.exp(-.5*((t[:,None]-t[None,:])/.035)**2)
    # Radsum 1 gjør at et konstant signal forblir konstant etter glatting.
    H /= H.sum(axis=1, keepdims=True)
    # Dette er data uten støy; støyen legges til etter at signalet er glattet.
    clean = H @ truth
    observed = clean + noise_level*np.random.default_rng(seed).standard_normal(n)
    return t, H, truth, clean, observed

# Ferdige forsøksverktøy. Rangregelen bruker bare data og oppgitt støynorm.
# Feil mot fasiten beregnes etterpå og påvirker ikke rangvalget.

def _choose_svd_rank(residuals, delta, tau):
    if delta <= 0 or tau < 1:
        raise ValueError('Bruk positivt støynivå og tau >= 1.')
    admissible = np.flatnonzero(residuals <= tau*delta)
    return int(admissible[0]) if admissible.size else None


def image_svd_trial(noise_level=.12, seed=23, tau=1.05, image_name='portrett'):
    if noise_level <= 0:
        raise ValueError('Bruk positivt støynivå.')
    C = challenge_images()[image_name]
    E = noise_level*np.random.default_rng(seed).standard_normal(C.shape)
    Y = C+E
    U,s,Vt = np.linalg.svd(Y,full_matrices=False)
    delta = np.linalg.norm(E,'fro')
    # Ingen klipping av data eller rekonstruksjon før feilberegningen.
    residuals = np.sqrt(np.r_[np.cumsum(s[::-1]**2)[::-1],0.])
    k = _choose_svd_rank(residuals,delta,tau)
    # Ortonormale koordinater gir de to leddene i prosjektets feilidentitet.
    signal_coordinates, noise_coordinates = U.T@C, U.T@E
    signal_energy = np.sum(signal_coordinates**2,axis=1)
    noise_energy = np.sum(noise_coordinates**2,axis=1)
    lost = np.r_[np.cumsum(signal_energy[::-1])[::-1],0.]
    kept = np.r_[0.,np.cumsum(noise_energy)]
    errors = np.sqrt(lost+kept)
    Yk = (U[:,:k]*s[:k])@Vt[:k,:]
    actual_error = np.linalg.norm(Yk-C,'fro')
    print(f'Frø {seed}; delta = {delta:.6g}; tau = {tau}; valgt rang = {k}')
    print('Alle normer i tabellen er absolutte Frobeniusnormer.')
    print('                  Rang    Dataavvik    Feil mot rent bilde')
    print(f'Valgt             {k:4d}    {np.linalg.norm(Y-Yk):9.4g}    {actual_error:9.4g}')
    print(f'Ubehandlet Y      {len(s):4d}    {0.:9.4g}    {delta:9.4g}')
    print('Feil² direkte:',actual_error**2,
          '= tapt signal² + beholdt støy²:',lost[k],'+',kept[k])
    best = int(np.argmin(errors))
    print(f'Beste rang mot fasit (bare etterpå): {best}; feil = {errors[best]:.6g}')
    ranks = np.arange(len(s)+1)
    fig,ax = plt.subplots(figsize=(7,4))
    ax.plot(ranks,residuals,label='Avvik mot Y')
    ax.plot(ranks,errors,label='Feil mot C')
    ax.plot(ranks,np.sqrt(lost),'--',label='Tapt signal')
    ax.plot(ranks,np.sqrt(kept),'--',label='Beholdt støy')
    ax.axhline(tau*delta,color='gray',linestyle=':',label='Tillatt dataavvik')
    ax.axvline(k,color='black',linestyle=':')
    ax.set(xlabel='Antall komponenter k',ylabel='Absolutt Frobeniusnorm')
    ax.legend(); fig.tight_layout(); plt.show()
    # To bilder i bredden, samme gråtoneskala; figurene endrer ikke dataene.
    fig,axes = plt.subplots(1,2,figsize=(7,3.5))
    for ax,M,title in zip(axes,[Y,Yk],['Data med støy',f'Regelens valg: k={k}']):
        ax.imshow(M,cmap='gray',vmin=0,vmax=1); ax.set_title(title); ax.axis('off')
    fig.tight_layout(); plt.show()
    return dict(k=k,delta=delta,residuals=residuals,errors=errors,
                lost_squared=lost,noise_squared=kept,singular_values=s,
                C=C,E=E,Y=Y,reconstruction=Yk)


def signal_svd_trial(noise_level=.005, seed=17, tau=1.05):
    if noise_level <= 0:
        raise ValueError('Bruk positivt støynivå.')
    t,H,truth,clean,b = blur_problem(noise_level=noise_level,seed=seed)
    eta = b-clean
    U,s,Vt = np.linalg.svd(H,full_matrices=False)
    delta = np.linalg.norm(eta)
    beta,alpha,e = U.T@b,Vt@truth,U.T@eta
    cutoff = max(H.shape)*np.finfo(float).eps*s[0]
    cap = int(np.sum(s>cutoff))
    ranks = np.arange(cap+1)
    # Resterende datakoordinater teller med også når rangsøket stopper ved cap.
    residual_formula = np.sqrt(np.r_[np.cumsum(beta[::-1]**2)[::-1],0.])[:cap+1]
    k = _choose_svd_rank(residual_formula,delta,tau)
    print(f'Frø {seed}; delta = {delta:.6g}; tau = {tau}; største tillatte rang = {cap}')
    if k is None:
        print('Ingen numerisk tillatt rang oppfyller residualkravet.')
        return dict(k=None,delta=delta,cap=cap)
    # Hver kolonne i solutions er en rekonstruksjon; første kolonne er x_0=0.
    increments = Vt[:cap,:].T*(beta[:cap]/s[:cap])
    solutions = np.column_stack([np.zeros(H.shape[1]),np.cumsum(increments,axis=1)])
    residuals = np.linalg.norm(b[:,None]-H@solutions,axis=0)
    errors = np.linalg.norm(solutions-truth[:,None],axis=0)
    amplified = np.r_[0.,np.cumsum((e[:cap]/s[:cap])**2)]
    lost = np.r_[np.cumsum(alpha[::-1]**2)[::-1],0.][:cap+1]
    print('Alle normer i tabellen er absolutte euklidske normer.')
    print('                  Rang    Residual    Rekonstruksjonsfeil')
    for label,j in [('Regelens valg',k),('Fast referanse',5)]:
        if j <= cap:
            print(f'{label:17s} {j:4d}    {residuals[j]:9.4g}    {errors[j]:9.4g}')
    print('Feil² direkte:',errors[k]**2,
          '= forsterket støy² + tapt signal²:',amplified[k],'+',lost[k])
    print('Residual² direkte / koordinatformel:',residuals[k]**2,residual_formula[k]**2)
    best = int(np.argmin(errors))
    print(f'Beste tillatte rang mot fasit (bare etterpå): {best}; feil = {errors[best]:.6g}')
    fig,axes = plt.subplots(2,1,figsize=(7,7))
    axes[0].semilogy(ranks,np.maximum(residuals,1e-16),label='Residual mot b')
    axes[0].semilogy(ranks,np.maximum(errors,1e-16),label='Feil mot kjent signal')
    axes[0].axhline(tau*delta,color='gray',linestyle=':',label='Tillatt residual')
    axes[0].axvline(k,color='black',linestyle=':')
    axes[0].set(xlabel='Antall beholdte retninger k',ylabel='Absolutt norm (logaritmisk)')
    axes[1].plot(t,truth,label='Kjent signal')
    axes[1].plot(t,solutions[:,k],label=f'Regelens valg: k={k}')
    if cap >= 5:
        axes[1].plot(t,solutions[:,5],'--',label='Fast referanse: k=5')
    axes[1].set(xlabel='Posisjon',ylabel='Signalverdi')
    for ax in axes: ax.legend(); ax.grid(alpha=.25)
    fig.tight_layout(); plt.show()
    return dict(k=k,delta=delta,cap=cap,residuals=residuals,errors=errors,
                residual_formula=residual_formula,lost_squared=lost,
                noise_squared=amplified,singular_values=s,
                H=H,b=b,truth=truth,eta=eta,reconstruction=solutions[:,k])