diff --git a/devel/bin/mccode-git-diff-code b/devel/bin/mccode-git-diff-code index c2258ea3bf..1b92ca7dd3 100755 --- a/devel/bin/mccode-git-diff-code +++ b/devel/bin/mccode-git-diff-code @@ -6,14 +6,16 @@ Drop-in replacement for `git diff --name-only BASE [TARGET]` that leaves out A .comp/.instr file is listed only if its "code" differs between BASE and TARGET, where code means: - everything from the DEFINE COMPONENT / DEFINE INSTRUMENT line onward, and - - the %Example: lines of the header (mctest uses them as test definitions), - so a header edit that touches %Example: lines still triggers a test run. + - the %Example: and %Scan: lines of the header (mctest uses them as test + definitions), so a header edit that touches them still triggers a test run. + A %Scan: line includes its {}-enclosed target values, which may span + several header lines; re-wrapping those values is not a change. Other header edits (mcdoc text, %P parameter docs, typos, units) are ignored, and so are whitespace-only edits: indentation, blank lines, re-wrapped lines, spacing around operators. Whitespace between two words (int a vs inta) and the line break after a #-preprocessor line or a line with // are still significant. -For each listed .comp/.instr file, the reason (new, code, %Example) is written +For each listed .comp/.instr file, the reason (new, code, %Example/%Scan) is written to stderr; stdout carries only the file names. Added or renamed files are always listed; deleted files never are (nothing to @@ -27,7 +29,8 @@ import subprocess import sys DEFINE = re.compile(r'^[ \t]*DEFINE[ \t]+(COMPONENT|INSTRUMENT)\b', re.M) -EXAMPLE = re.compile(r'%Example:.*') +# a %Scan: line runs on to the } closing its target values, possibly over several lines +TESTLINE = re.compile(r'%(?:Example:[^\n]*|Scan:[^\n{]*(?:\{[^}]*\}[^\n]*)?)') def git(*args): @@ -65,11 +68,14 @@ def squash(text): def split(text): - ''' (%Example lines, code body) of a .comp/.instr file, whitespace-insensitive. ''' + ''' (%Example/%Scan lines, code body) of a .comp/.instr file, whitespace-insensitive. ''' m = DEFINE.search(text) if not m: return '', squash(text) # unrecognised layout: compare everything - return squash('\n'.join(EXAMPLE.findall(text[:m.start()]))), squash(text[m.start():]) + # the header's * line prefixes inside a multi-line {} value list are layout, not content + tests = [re.sub(r'\{[^}]*\}', lambda v: v.group(0).replace('*', ' '), t) + for t in TESTLINE.findall(text[:m.start()])] + return squash('\n'.join(tests)), squash(text[m.start():]) def main(argv): @@ -85,7 +91,7 @@ def main(argv): why = 'new' else: (ex0, body0), (ex1, body1) = split(old), split(new) - why = ', '.join(w for w, a, b in (('code', body0, body1), ('%Example', ex0, ex1)) if a != b) + why = ', '.join(w for w, a, b in (('code', body0, body1), ('%Example/%Scan', ex0, ex1)) if a != b) if not why: continue sys.stderr.write('mccode-git-diff-code: %s (%s)\n' % (path, why)) diff --git a/mcstas-comps/examples/LLB/LLB_6T2/LLB_6T2.instr b/mcstas-comps/examples/LLB/LLB_6T2/LLB_6T2.instr index 79fbe3c557..a4292ba4fa 100644 --- a/mcstas-comps/examples/LLB/LLB_6T2/LLB_6T2.instr +++ b/mcstas-comps/examples/LLB/LLB_6T2/LLB_6T2.instr @@ -202,7 +202,7 @@ COMPONENT Pos_2 = Arm( AT (0,0,b2mon) RELATIVE Bar_out ROTATED (0,monoc_ang,0) ABSOLUTE -COMPONENT Monoc=Monochromator_curved( +SPLIT 10 COMPONENT Monoc=Monochromator_curved( zwidth=0.015, yheight=0.1, gap=0.0005, NH=7, NV=1, mosaich=25.0, mosaicv=25.0, reflect="HOPG.rfl", transmit="HOPG.trm", diff --git a/mcstas-comps/examples/Risoe/TAS1_C1/TAS1_C1.instr b/mcstas-comps/examples/Risoe/TAS1_C1/TAS1_C1.instr index 807fe531c6..1d290d4335 100644 --- a/mcstas-comps/examples/Risoe/TAS1_C1/TAS1_C1.instr +++ b/mcstas-comps/examples/Risoe/TAS1_C1/TAS1_C1.instr @@ -26,6 +26,9 @@ * detector is at the sample position. * * %Example: PHM=-37.077 Detector: sng_I=0.000429426 +* %Scan: mcrun TAS1_C1.instr PHM=-38.077,-36.077 -N11 -n1e6 Detector: sng_I={ +* 1.7037e-05,5.77737e-05,0.000131759,0.000247282,0.00036353,0.000415802,0.000377557,0.000267179, +* 0.000149367,5.89201e-05,1.85848e-05 } * * %Parameters * PHM: [deg] Monochromator arm angle, aka A1 diff --git a/mcstas-comps/examples/Risoe/TAS1_C1_Tilt/TAS1_C1_Tilt.instr b/mcstas-comps/examples/Risoe/TAS1_C1_Tilt/TAS1_C1_Tilt.instr index a3a9b6e5cf..383263ca1d 100644 --- a/mcstas-comps/examples/Risoe/TAS1_C1_Tilt/TAS1_C1_Tilt.instr +++ b/mcstas-comps/examples/Risoe/TAS1_C1_Tilt/TAS1_C1_Tilt.instr @@ -26,12 +26,16 @@ * detector is at the sample position and there is no analyzer. * * %Example: TTM=-74 Detector: sng_I=0.000405056 +* %Scan: mcrun TAS1_C1_Tilt.instr OMC1=-44.5,55.5 -N21 -n1e6 Detector: sng_I={ +* 7.82478e-05,0.000118605,0.000165766,0.000222136,0.000277334,0.000324033,0.00036983,0.000407361, +* 0.000420118,0.000412682,0.000402763,0.000373695,0.000323706,0.000280374,0.000224993,0.000165908, +* 0.000121172,8.11465e-05,4.92901e-05,2.88144e-05,1.46245e-05 } * * %Parameters * PHM: [deg] Monochromator arm angle, aka A1 * TTM: [deg] Monochromator take-off angle, aka A2 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * * %Link * The McStas User manual diff --git a/mcstas-comps/examples/Risoe/TAS1_Diff_Powder/TAS1_Diff_Powder.instr b/mcstas-comps/examples/Risoe/TAS1_Diff_Powder/TAS1_Diff_Powder.instr index 79907965c5..d59528a67f 100644 --- a/mcstas-comps/examples/Risoe/TAS1_Diff_Powder/TAS1_Diff_Powder.instr +++ b/mcstas-comps/examples/Risoe/TAS1_Diff_Powder/TAS1_Diff_Powder.instr @@ -27,6 +27,9 @@ * sample is a powder and there is no analyzer. * * %Example: PHM=-37.077 Detector: sng_I=4.50892e-07 +* %Scan: mcrun TAS1_Diff_Powder.instr TT=32.92,34.32 -N15 -n1e6 Detector: sng_I={ +* 3.44621e-09,1.42523e-08,3.66193e-08,9.83779e-08,1.88888e-07,3.09967e-07,4.32565e-07,5.50027e-07, +* 5.59591e-07,5.02795e-07,3.84587e-07,2.46907e-07,1.35535e-07,6.38448e-08,2.13895e-08 } * * %Parameters * PHM: [deg] Monochromator arm angle, aka A1 @@ -34,7 +37,7 @@ * TT: [deg] Take-off angle at the sample position, aka A4 * TTA: [deg] Take-off angle at the analyzer position, aka A6 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * C2: [min] Collimator 2 aperture (sample-ana arm) * C3: [min] Collimator 3 aperture (ana-detector arm) * diff --git a/mcstas-comps/examples/Risoe/TAS1_Diff_Slit/TAS1_Diff_Slit.instr b/mcstas-comps/examples/Risoe/TAS1_Diff_Slit/TAS1_Diff_Slit.instr index a73f2c1adf..782b51d707 100644 --- a/mcstas-comps/examples/Risoe/TAS1_Diff_Slit/TAS1_Diff_Slit.instr +++ b/mcstas-comps/examples/Risoe/TAS1_Diff_Slit/TAS1_Diff_Slit.instr @@ -27,13 +27,16 @@ * sample is a slit and there is no analyzer. * * %Example: C1=30 TT=0 Detector: sng_I=2e-5 +* %Scan: mcrun TAS1_Diff_Slit.instr C1=30 TT=-0.6,0.6 -N13 -n1e6 Detector: sng_I={ +* 6.85677e-07,2.11237e-06,4.51499e-06,8.20823e-06,1.39162e-05,1.73145e-05,2.01992e-05,2.05377e-05, +* 1.71222e-05,1.22857e-05,8.94974e-06,5.09308e-06,2.51101e-06 } * * %Parameters * PHM: [deg] Monochromator arm angle, aka A1 * TTM: [deg] Monochromator take-off angle, aka A2 * TT: [deg] Take-off angle at the sample position, aka A4 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * C2: [min] Collimator 2 aperture (sample-ana arm) * C3: [min] Collimator 3 aperture (ana-detector arm) * diff --git a/mcstas-comps/examples/Risoe/TAS1_Diff_Vana/TAS1_Diff_Vana.instr b/mcstas-comps/examples/Risoe/TAS1_Diff_Vana/TAS1_Diff_Vana.instr index ec11525698..10de482760 100644 --- a/mcstas-comps/examples/Risoe/TAS1_Diff_Vana/TAS1_Diff_Vana.instr +++ b/mcstas-comps/examples/Risoe/TAS1_Diff_Vana/TAS1_Diff_Vana.instr @@ -27,6 +27,10 @@ * sample is a vanadium cylinder and there is no analyzer. * * %Example: PHM=-37.077 Detector: sng_I=1.95173e-09 +* %Scan: mcrun TAS1_Diff_Vana.instr TT=31.52,35.52 -N21 -n1e6 Detector: sng_I={ +* 1.8001e-09,1.87932e-09,1.86258e-09,1.89395e-09,1.86904e-09,1.86482e-09,1.85878e-09,1.87728e-09, +* 1.81479e-09,1.84234e-09,1.91997e-09,1.89089e-09,1.82982e-09,1.9059e-09,1.88333e-09,1.775e-09, +* 1.80894e-09,1.99095e-09,1.89532e-09,1.94438e-09,1.94377e-09 } * * %Parameters * PHM: [deg] Monochromator arm angle, aka A1 @@ -34,7 +38,7 @@ * TT: [deg] Take-off angle at the sample position, aka A4 * TTA: [deg] Take-off angle at the analyzer position, aka A6 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * C2: [min] Collimator 2 aperture (sample-ana arm) * C3: [min] Collimator 3 aperture (ana-detector arm) * diff --git a/mcstas-comps/examples/Risoe/TAS1_Powder/TAS1_Powder.instr b/mcstas-comps/examples/Risoe/TAS1_Powder/TAS1_Powder.instr index 22cd1cc7cb..fe207cdca8 100644 --- a/mcstas-comps/examples/Risoe/TAS1_Powder/TAS1_Powder.instr +++ b/mcstas-comps/examples/Risoe/TAS1_Powder/TAS1_Powder.instr @@ -25,6 +25,12 @@ * The sample is a powder and the analyzer is a single plate. * * %Example: PHM=-37.077 Detector: sng_I=2.26723e-07 +* With TTA=0 the detector looks straight through the analyzer, so an OMA scan shows +* a transmission dip where the analyzer is at its Bragg angle: +* %Scan: mcrun TAS1_Powder.instr OMA=-19.45,-15.45 -N21 -n1e6 Detector: sng_I={ +* 4.37937e-07,4.53735e-07,4.60594e-07,4.49352e-07,4.43857e-07,4.36457e-07,4.22756e-07,4.15322e-07, +* 3.49642e-07,2.82077e-07,2.31752e-07,2.48424e-07,3.1068e-07,3.96655e-07,4.27444e-07,4.42215e-07, +* 4.49148e-07,4.60529e-07,4.36828e-07,4.39458e-07,4.52569e-07 } * * %Parameters * PHM: [deg] Monochromator rotation angle, aka A1 @@ -33,7 +39,7 @@ * OMA: [deg] Analyzer rotation angle, aka A5 * TTA: [deg] Take-off angle at the analyzer position, aka A6 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * C2: [min] Collimator 2 aperture (sample-ana arm) * C3: [min] Collimator 3 aperture (ana-detector arm) * diff --git a/mcstas-comps/examples/Risoe/TAS1_Vana/TAS1_Vana.instr b/mcstas-comps/examples/Risoe/TAS1_Vana/TAS1_Vana.instr index a63a3ae25e..53e8e737a1 100644 --- a/mcstas-comps/examples/Risoe/TAS1_Vana/TAS1_Vana.instr +++ b/mcstas-comps/examples/Risoe/TAS1_Vana/TAS1_Vana.instr @@ -25,6 +25,12 @@ * The sample is a vanadium and the analyzer is a single plate. * * %Example: PHM=-37.077 Detector: sng_I=1.11099e-09 +* With TTA=0 the detector looks straight through the analyzer, so an OMA scan shows +* a transmission dip where the analyzer is at its Bragg angle: +* %Scan: mcrun TAS1_Vana.instr OMA=-19.45,-15.45 -N21 -n1e6 Detector: sng_I={ +* 1.81292e-09,1.90855e-09,1.8668e-09,1.86678e-09,1.85029e-09,1.834e-09,1.76611e-09,1.60322e-09, +* 1.33626e-09,1.0796e-09,1.09548e-09,1.2646e-09,1.46433e-09,1.72966e-09,1.83299e-09,1.77587e-09, +* 1.82133e-09,1.94325e-09,1.88317e-09,1.89597e-09,1.94066e-09 } * * %Parameters * PHM: [deg] Monochromator rotation angle, aka A1 @@ -33,7 +39,7 @@ * OMA: [deg] Analyzer rotation angle, aka A5 * TTA: [deg] Take-off angle at the analyzer position, aka A6 * C1: [min] Collimator 1 aperture (mono-sample arm) -* OMC1: [deg] Tilt angle of the Collimator 1 +* OMC1: [arcmin] Tilt angle of the Collimator 1 * C2: [min] Collimator 2 aperture (sample-ana arm) * C3: [min] Collimator 3 aperture (ana-detector arm) * diff --git a/mcstas-comps/examples/TUDelft/SESANS_Delft/SESANS_Delft.instr b/mcstas-comps/examples/TUDelft/SESANS_Delft/SESANS_Delft.instr index a78758bd02..11b97b4657 100644 --- a/mcstas-comps/examples/TUDelft/SESANS_Delft/SESANS_Delft.instr +++ b/mcstas-comps/examples/TUDelft/SESANS_Delft/SESANS_Delft.instr @@ -24,7 +24,11 @@ * * %Example: -y Detector: det_I=2.2537e-07 * -* Scan Example: SESANS_Delft -n 1000000 -N 31 L0=2.165 DL=0.02 By=0,0.0468 +* %Scan: mcrun SESANS_Delft.instr L0=2.165 DL=0.02 By=0,0.0468 -N31 -n1e6 Detector: det_I={ +* 2.35556e-07,2.35296e-07,2.30432e-07,2.25294e-07,2.20061e-07,2.13707e-07,2.11276e-07,2.05654e-07, +* 2.03486e-07,2.00872e-07,1.99138e-07,2.00791e-07,1.99118e-07,1.98852e-07,1.98976e-07,1.99064e-07, +* 1.98964e-07,1.99602e-07,1.98235e-07,2.00064e-07,1.99617e-07,1.99683e-07,1.98427e-07,1.98624e-07, +* 1.97647e-07,1.98342e-07,1.98835e-07,1.99153e-07,1.99708e-07,1.99172e-07,1.98348e-07 } * * %P * : [] diff --git a/mcstas-comps/examples/Templates/Tomography/README.md b/mcstas-comps/examples/Templates/Tomography/README.md index 98a2427c1d..1d1688ce36 100644 --- a/mcstas-comps/examples/Templates/Tomography/README.md +++ b/mcstas-comps/examples/Templates/Tomography/README.md @@ -13,19 +13,22 @@ ```text Instrument to study tomographic imaging by means of the feature of OFF shape samples. +The sample (geometry, an OFF file, default socket.off) is rotated by omega around the +vertical axis, and the transmitted beam is recorded on a 2D detector (monitor). -Example: mcrun Tomography.instr offfile=bunny.off -n1e4 -N18 omega=0,340 -d TomoScan -(Note that to achieve proper statistics for tomographic reconstruction, MUCH higher ncounts +A tomography is a scan of omega over a full rotation, e.g. +mcrun Tomography.instr omega=0,355 -N72 -n1e7 -d TomoScan +(to achieve proper statistics for tomographic reconstruction, MUCH higher ncounts are needed) -Use the provided Matlab tomo_recon.m function (requires imaging toolbox, PGPLOT output data -and a Unix/Mac) in the tools/matlab folder to reconstruct a 3D volume of the object. Use e.g. -isosurface to do thresholding for extraction of the object surface. +Use the provided tomo_recon.py (numpy + matplotlib) in this folder to reconstruct a 3D +volume of the object from the scan directory: python tomo_recon.py TomoScan [--save] ``` ## Examples -- **Test: omega=0 Detector: monitor_I=2.23492e-09** +- **Test: omega=0 Detector: monitor_I=9.37708e-10** +- **Scan: mcrun Tomography.instr omega=0,355 -N72 -n1e6 Detector: monitor_I={72 values}** ## Input parameters @@ -41,9 +44,9 @@ Parameters in **boldface** are required; the others are optional. | div_h | deg | Source horisontal divergence (angular width) | 1e-4 | | source_w | m | Source width | 0.4 | | source_h | m | Source height | 0.2 | -| det_w | m | Detector width | 0.4 | -| det_h | m | Detector height | 0.2 | -| opts | string | Monitor_nD options string | "x bins=80 y bins=40" | +| det_w | m | Detector width | 0.25 | +| det_h | m | Detector height | 0.15 | +| opts | string | Monitor_nD options string | "x bins=128 y bins=64" | ## Links diff --git a/mcstas-comps/examples/Templates/Tomography/Tomography.instr b/mcstas-comps/examples/Templates/Tomography/Tomography.instr index 4d00b1c1b6..f43d44c0a9 100644 --- a/mcstas-comps/examples/Templates/Tomography/Tomography.instr +++ b/mcstas-comps/examples/Templates/Tomography/Tomography.instr @@ -13,16 +13,28 @@ * * %Description * Instrument to study tomographic imaging by means of the feature of OFF shape samples. +* The sample (geometry, an OFF file, default socket.off) is rotated by omega around the +* vertical axis, and the transmitted beam is recorded on a 2D detector (monitor). * -* Example: mcrun Tomography.instr offfile=bunny.off -n1e4 -N18 omega=0,340 -d TomoScan -* (Note that to achieve proper statistics for tomographic reconstruction, MUCH higher ncounts +* A tomography is a scan of omega over a full rotation, e.g. +* mcrun Tomography.instr omega=0,355 -N72 -n1e7 -d TomoScan +* (to achieve proper statistics for tomographic reconstruction, MUCH higher ncounts * are needed) * -* Use the provided Matlab tomo_recon.m function (requires imaging toolbox, PGPLOT output data -* and a Unix/Mac) in the tools/matlab folder to reconstruct a 3D volume of the object. Use e.g. -* isosurface to do thresholding for extraction of the object surface. +* Use the provided tomo_recon.py (numpy + matplotlib) in this folder to reconstruct a 3D +* volume of the object from the scan directory: python tomo_recon.py TomoScan [--save] * -* %Example: omega=0 Detector: monitor_I=2.23492e-09 +* %Example: omega=0 Detector: monitor_I=9.37708e-10 +* %Scan: mcrun Tomography.instr omega=0,355 -N72 -n1e6 Detector: monitor_I={ +* 9.42128e-10,9.25881e-10,9.14476e-10,9.04288e-10,8.94769e-10,8.8652e-10,8.84254e-10,8.86881e-10, +* 8.87397e-10,8.86793e-10,8.80893e-10,8.80094e-10,8.82899e-10,8.84653e-10,8.8854e-10,8.9252e-10, +* 9.01629e-10,9.10865e-10,9.22842e-10,9.09788e-10,8.99241e-10,8.92795e-10,8.87909e-10,8.83615e-10, +* 8.80372e-10,8.81549e-10,8.8362e-10,8.86581e-10,8.88443e-10,8.89103e-10,8.85436e-10,8.86795e-10, +* 8.92286e-10,9.02972e-10,9.15912e-10,9.25596e-10,9.39192e-10,9.28645e-10,9.15208e-10,9.03769e-10, +* 8.93097e-10,8.87378e-10,8.87768e-10,8.88077e-10,8.88221e-10,8.88865e-10,8.84663e-10,8.81642e-10, +* 8.83464e-10,8.85594e-10,8.8875e-10,8.90938e-10,8.98263e-10,9.09941e-10,9.22377e-10,9.11957e-10, +* 8.98938e-10,8.90529e-10,8.85399e-10,8.83245e-10,8.80727e-10,8.78931e-10,8.82134e-10,8.84583e-10, +* 8.84924e-10,8.87263e-10,8.88867e-10,8.87931e-10,8.92295e-10,9.00247e-10,9.15281e-10,9.26476e-10 } * * %Parameters * geometry: [string] Name of the OFF file describing the sample shape @@ -42,7 +54,7 @@ * * %End *******************************************************************************/ -DEFINE INSTRUMENT Tomography(string geometry="socket.off", omega=0, sigma_abs=100, frac_scatt=0, div_v=1e-4, div_h=1e-4, source_w=0.4, source_h=0.2, det_w=0.4, det_h=0.2, string opts="x bins=80 y bins=40") +DEFINE INSTRUMENT Tomography(string geometry="socket.off", omega=0, sigma_abs=100, frac_scatt=0, div_v=1e-4, div_h=1e-4, source_w=0.4, source_h=0.2, det_w=0.25, det_h=0.15, string opts="x bins=128 y bins=64") DEPENDENCY " -DUSE_OFF " TRACE diff --git a/mcstas-comps/examples/Templates/Tomography/tomo_recon.py b/mcstas-comps/examples/Templates/Tomography/tomo_recon.py new file mode 100644 index 0000000000..4cd7e06d93 --- /dev/null +++ b/mcstas-comps/examples/Templates/Tomography/tomo_recon.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 +"""Reconstruct a volume from a Tomography.instr omega scan, e.g. + + mcrun Tomography.instr omega=0,350 -N36 -n1e7 -d TomoScan + python tomo_recon.py TomoScan + +Filtered back-projection (ramp*Hann filter, like Matlab iradon's 'hann'), +one sinogram per detector row. Shows 3 central slices; --save writes vol.npy. +Python port of tools/matlab/tomo_recon.m (PW, 20080620). +""" +import glob, os, sys +import numpy as np +import matplotlib.pyplot as plt + + +def header(path, key): + for line in open(path): + if line.startswith('# ' + key): + return line[len(key) + 2:].strip() + + +def iradon(sino, theta): + """sino: (nbins, nproj), theta in degrees -> (nbins, nbins) slice.""" + n = sino.shape[0] + pad = 1 << int(np.ceil(np.log2(2 * n))) + f = np.fft.fftfreq(pad) + h = 2 * np.abs(f) * (1 + np.cos(2 * np.pi * f)) / 2 + filt = np.fft.ifft(np.fft.fft(sino, pad, axis=0) * h[:, None], axis=0).real[:n] + x = np.arange(n) - (n - 1) / 2 + X, Y = np.meshgrid(x, -x) + img = np.zeros((n, n)) + for p, t in zip(filt.T, np.deg2rad(theta)): + img += np.interp(X * np.cos(t) + Y * np.sin(t), x, p, left=0, right=0) + return img * np.pi / (2 * len(theta)) + + +def main(datadir, save=False): + dat = os.path.join(datadir, 'mccode.dat') + nproj = int(header(dat, 'Numpoints:')) + lo, hi = map(float, header(dat, 'xlimits:').split()) + theta = np.linspace(lo, hi, nproj) + + mons = [] + for j in range(nproj): + f = glob.glob(os.path.join(datadir, str(j), '*.x_y'))[0] + I = np.loadtxt(f) + mons.append(I[:len(I) // 3]) # I block only, rows = y slices + geometry = header(f, 'Param: geometry=') + mons = np.array(mons) # (nproj, nslice, nbins) + pos = mons[mons > 0] + mons = np.maximum(mons, pos.min() if pos.size else 1) # avoid log(0) + + vol = np.stack([iradon(np.log(m.max() / m).T, theta) + for m in mons.transpose(1, 0, 2)], axis=2) # (x, z, y) + if save: + np.save(os.path.join(datadir, 'vol.npy'), vol) + + title = 'slice from %s, Hann filtered' % geometry + c = [s // 2 for s in vol.shape] + fig, ax = plt.subplots(1, 3, figsize=(15, 5)) + for a, img, name in zip(ax, (vol[c[0]].T, vol[:, c[1]].T, vol[:, :, c[2]]), 'xyz'): + a.imshow(img, origin='lower') + a.set_title('Central %s %s' % (name, title)) + plt.show() + return vol + + +if __name__ == '__main__': + if len(sys.argv) < 2 or not os.path.isdir(sys.argv[1]): + sys.exit('usage: tomo_recon.py [--save]') + main(sys.argv[1], '--save' in sys.argv) diff --git a/mcstas-comps/examples/Tests_polarization/SE_example/SE_example.instr b/mcstas-comps/examples/Tests_polarization/SE_example/SE_example.instr index 1a9a1190d4..8d743001ac 100644 --- a/mcstas-comps/examples/Tests_polarization/SE_example/SE_example.instr +++ b/mcstas-comps/examples/Tests_polarization/SE_example/SE_example.instr @@ -19,6 +19,12 @@ * MeanPolLambda_montior is used to monitor beam polarisation. * * Example: mcrun SE_example.instr dBz=-0.0001,0.0001 -N41 -n1e5 +* %Scan: mcrun SE_example.instr dBz=-0.0001,0.0001 -N41 -n1e5 Detector: detector_I={ +* 86.04,101.892,108.731,94.7348,89.0012,85.6308,108.935,102.351,69.0694,92.3323, +* 136.083,64.8421,57.8734,158.608,99.3713,18.8129,133.474,164.129,20.9553,74.1501, +* 193.637,76.1058,20.3687,168.351,133.934,18.4337,110.268,161.186,59.7896,65.6851, +* 141.567,96.8757,64.2307,106.846,110.065,90.8057,100.937,92.2149,106.729,106.231, +* 85.406 } * * %Parameters * POL_ANGLE: [deg] Reflection angle of polarizer/analyzer diff --git a/mcstas-comps/examples/Tests_polarization/SE_example2/SE_example2.instr b/mcstas-comps/examples/Tests_polarization/SE_example2/SE_example2.instr index 6103e207ee..070ec58a08 100644 --- a/mcstas-comps/examples/Tests_polarization/SE_example2/SE_example2.instr +++ b/mcstas-comps/examples/Tests_polarization/SE_example2/SE_example2.instr @@ -18,7 +18,13 @@ * Pol_FieldBox is used to define guidefields and flippers. * MeanPolLambda_montior is used to monitor beam polarisation. * -* Example: mcrun SE_example2.instr dBz=-0.0001,0.0001 -N41 -n1e5 +* %Scan: mcrun SE_example2.instr dBz=-0.0001,0.0001 -N41 -n1e5 Detector: detector_I={ +* 131.576,155.634,159.227,142.21,139.368,130.51,165.979,155.359, +* 99.6514,142.272,209.645,95.7382,90.1825,237.694,146.859,28.2776, +* 199.631,244.785,32.2439,112.725,277.63,115.369,29.3693,243.049, +* 201.203,28.3267,163.182,240.944,90.531,95.7656,212.13,145.286, +* 93.5099,161.98,164.791,134.604,147.311,143.726,153.511,155.375, +* 129.977 } * * %Parameters * POL_ANGLE: [deg] Reflection angle of polarizer/analyzer diff --git a/mcstas-comps/examples/Tests_samples/Samples_Incoherent/Samples_Incoherent.instr b/mcstas-comps/examples/Tests_samples/Samples_Incoherent/Samples_Incoherent.instr index 78eeb42a07..8066109e48 100644 --- a/mcstas-comps/examples/Tests_samples/Samples_Incoherent/Samples_Incoherent.instr +++ b/mcstas-comps/examples/Tests_samples/Samples_Incoherent/Samples_Incoherent.instr @@ -59,6 +59,12 @@ * %Example: SAMPLE=8 STOP=1 -n 1e5 Detector: PSD_Sphere_4pi_I=1.3233e+06 * %Example: SAMPLE=9 STOP=1 -n 1e5 Detector: PSD_Sphere_4pi_I=1.3239e+06 * %Example: SAMPLE=10 STOP=1 -n 1e5 Detector: PSD_Sphere_4pi_I=1.28905e+06 +* %Scan: mcrun Samples_Incoherent.instr SAMPLE=1,5 -N5 STOP=1 Detector: PSD_Sphere_4pi_I={ +* 1326220.0,1245050.0,1322310.0,1321210.0,1326550.0 } +* %Scan: mcrun Samples_Incoherent.instr SAMPLE=2,5 -N4 STOP=0 Detector: PSD_Sphere_4pi_I={ +* 12415200.0,7154040.0,7759140.0,1327980.0 } +* %Scan: mcrun Samples_Incoherent.instr SAMPLE=2,5 -N4 STOP=0 DB=1 Detector: Dirbeam_I={ +* 11195400.0,5833780.0,6441470.0,5783.69 } * * %Parameters * L_min: [AA] Minimum wavelength of neutrons diff --git a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr index 05873c9444..46b4fe57ad 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr @@ -43,8 +43,6 @@ * along x, y, z for both samples), and (h, l) = (1, 0) is the magnetic zone * centre MnF2 (1 0 0). * -* Example: mcrun Test_Dispersion_relation_TAS.instr -M -N 21,21 -n 2e5 h=0.6,1.4 l=0.6,1.4 dE=5 sample_choice=0 -* * The dispersion files are written by the SHELL command below when the * instrument is compiled, if they do not exist yet, by * generate_dispersion_files.py (Python with numpy). @@ -54,6 +52,11 @@ * %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10 +* %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.6,1.4 l=0.6,1.4 dE=5 sample_choice=0 Detector: detector_I={ +* 0.0,0.0,0.0,0.0,0.0,0.0,1.2167e-12,3.52579e-11, +* 2.26406e-12,0.0,0.0,1.17383e-11,0.0,1.39525e-11,0.0,0.0, +* 6.10596e-12,1.00579e-11,8.66923e-12,0.0,0.0,0.0,0.0,0.0, +* 0.0 } * * %Parameters * h: [r.l.u.] Scattering vector component along a* (sample x axis) diff --git a/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr b/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr index d0e034c658..5e079410c0 100644 --- a/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr +++ b/mcstas-comps/examples/Tests_samples/Test_Magnon_bcc_TAS/Test_Magnon_bcc_TAS.instr @@ -13,7 +13,10 @@ * Generic TAS instrument for test of samples with dispersions. Modeled over the RITA-2 instrument at PSI (one analyzer only). The instrument is able to position in q-space at q=(h 0 0) and fixed Ei and Ef. This can be used to longitudinal constant-E scans or constant-q scans. In addition, transverse constant-E scans can be made by scanning the sample orientation A3. * * %Example: Test_Magnon_bcc_TAS.instr hw=0.5 Detector: e_monitor1_I=460000 -* Example: Test_Magnon_bcc_TAS.instr hw=0,1 -N21 +* %Scan: mcrun Test_Magnon_bcc_TAS.instr hw=0,1 -N21 -n1e6 Detector: e_monitor_I={ +* 5541440.0,5681320.0,4956300.0,4947280.0,4886770.0,4370310.0,3596940.0,3900100.0, +* 3656700.0,3712310.0,3128430.0,3254960.0,2708640.0,2961760.0,2536780.0,2524220.0, +* 2570960.0,2493820.0,2173830.0,2275870.0,2174740.0 } * * %P * hw: [meV] Neutron energy transfer diff --git a/mcstas-comps/examples/Tests_samples/Test_PowderN_Res/Test_PowderN_Res.instr b/mcstas-comps/examples/Tests_samples/Test_PowderN_Res/Test_PowderN_Res.instr index f1ec2cbe7e..a29e9f8c3c 100644 --- a/mcstas-comps/examples/Tests_samples/Test_PowderN_Res/Test_PowderN_Res.instr +++ b/mcstas-comps/examples/Tests_samples/Test_PowderN_Res/Test_PowderN_Res.instr @@ -18,7 +18,10 @@ * * By default, absorption in the sample is effectively suppressed. * -* Example: Lambda0=2.5667 dLambda=0.01 radius=0.001,0.021 TT=71.8 D_PHI=6 SPLITS=117 Distance=30 sig_abs=1e-09 -N21 +* %Scan: mcrun Test_PowderN_Res.instr Lambda0=2.5667 dLambda=0.01 radius=0.001,0.021 TT=71.8 D_PHI=6 SPLITS=117 Distance=30 sig_abs=1e-09 -N21 Detector: DetectorSmall_I={ +* 3.86698,15.1625,33.2541,57.6607,85.9154,119.618,158.827,201.183, +* 243.131,289.07,340.294,393.375,443.613,498.612,549.899,612.472, +* 671.715,723.508,778.461,837.07,895.16 } * %Example: Lambda0=2.5667 dLambda=0.01 radius=0.01 TT=71.8 D_PHI=6 SPLITS=117 Distance=30 sig_abs=1e-09 Detector: DetectorSmall_I=288.57 * * %Parameters diff --git a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr index 6eb10defa4..5ceec336f9 100644 --- a/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr +++ b/mcxtrace-comps/examples/BARC/Czerny_Turner/Czerny_Turner.instr @@ -20,6 +20,10 @@ * Example: Do a scan (scan parameter =1) from x_screw=20 to 103 mm. The calculated monitor "Wavelength(Ang) as a function of x_screw(mm)" will be a linear function. The monochromator functions properly. * * %Example: x_screw=60 Detector: w_monitor_I=1.10881e-06 +* %Scan: mxrun Czerny_Turner.instr -N21 x_screw=20,103 scan=1 -n1e6 Detector: Wavelength_I={ +* 2007.08,2406.85,2820.95,3242.17,3659.02,4075.98,4485.76,4901.23, +* 5324.84,5733.95,6156.24,6569.68,6979.66,7395.58,7805.61,8222.2, +* 8637.53,9053.14,9467.18,9882.52,10288.8 } * * %Parameters * x_screw: [mm] Displacement perpendicular to initial position of the lever (sine drive mechanism) diff --git a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr index 588dc6b6cc..e4918a1750 100644 --- a/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr +++ b/mcxtrace-comps/examples/SOLEIL/SOLEIL_DISCO/SOLEIL_DISCO.instr @@ -13,7 +13,10 @@ * %Description * DISCO beamline description here: https://hal.archives-ouvertes.fr/hal-01479318/document * You may scan the grating monochromator with e.g.: -* mxrun -n 1e5 SOLEIL_DISCO.instr -N84 x_screw=20,103 scan=1 +* %Scan: mxrun SOLEIL_DISCO.instr -N21 x_screw=20,103 scan=1 Detector: Wavelength_I={ +* 2019.39,2365.97,2837.41,3242.59,3655.92,4070.23,4487.58,4921.15, +* 5306.77,5738.39,6151.91,6565.99,6979.47,7395.13,7786.17,8199.52, +* 8631.53,9051.46,9483.21,9913.98,10299.5 } * * %Example: grazing_angle_one=22.5 -n 1e5 Detector: second_HFP_I=2.08292e+11 * diff --git a/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_not_white_src/SSRL_bl_11_2_not_white_src.instr b/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_not_white_src/SSRL_bl_11_2_not_white_src.instr index a3b1eef5ca..88fc158472 100644 --- a/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_not_white_src/SSRL_bl_11_2_not_white_src.instr +++ b/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_not_white_src/SSRL_bl_11_2_not_white_src.instr @@ -19,7 +19,11 @@ * Documentation and results of the .instr here: https://gitlab.synchrotron-soleil.fr/grades/beamlines/-/tree/main/glitches * A jupyter notebook will soon be available regrouping everything. * You may scan like so for example: -* mxrun -n 1e7 SSRL_bl_11_2_not_white_src.instr -N601 Etohit=6900,7500 +* %Scan: mxrun SSRL_bl_11_2_not_white_src.instr -N31 Etohit=6900,7500 -n1e6 Detector: EnergyMonitor_first_I={ +* 13328900000.0,13634200000.0,14541900000.0,13884900000.0,13898100000.0,13682000000.0,13694000000.0,14507100000.0, +* 13206900000.0,13610500000.0,13946900000.0,15010800000.0,14594000000.0,15530700000.0,14609000000.0,13841600000.0, +* 14696400000.0,14927300000.0,15831200000.0,15905500000.0,16147000000.0,16044900000.0,16097800000.0,15847300000.0, +* 15588300000.0,16880100000.0,16423800000.0,18149700000.0,17975600000.0,17924000000.0,18274500000.0 } * * Example: Etohit=6900 Detector: EnergyMonitor_first_I=1.3e+10 * diff --git a/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_white_src/SSRL_bl_11_2_white_src.instr b/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_white_src/SSRL_bl_11_2_white_src.instr index 9dbc2eacd2..7faea4ec76 100644 --- a/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_white_src/SSRL_bl_11_2_white_src.instr +++ b/mcxtrace-comps/examples/SSRL/SSRL_bl_11_2_white_src/SSRL_bl_11_2_white_src.instr @@ -18,7 +18,11 @@ * Documentation and results of the .instr here: https://gitlab.synchrotron-soleil.fr/grades/beamlines/-/tree/main/glitches * A jupyter notebook will soon be available regrouping everything. * You may scan like so for example: -* mxrun -n 1e7 SSRL_bl_11_2_white_src.instr -N601 Etohit=6900,7500 detuning_percentage=40 +* %Scan: mxrun SSRL_bl_11_2_white_src.instr -N31 Etohit=6900,7500 detuning_percentage=40 -n1e7 Detector: EnergyMonitor_first_I={ +* 9.6608e+06,9.36768e+06,8.8759e+06,8.45813e+06,8.60571e+06,2.19181e+07,2.5781e+07,2.47468e+07, +* 2.33347e+07,2.36432e+07,2.44606e+07,2.54374e+07,2.50395e+07,2.47904e+07,2.48351e+07,2.48054e+07, +* 2.25687e+07,2.16413e+07,2.15756e+07,2.08025e+07,2.13743e+07,2.14039e+07,2.09034e+07,1.55632e+07, +* 1.75271e+07,1.73956e+07,1.73963e+07,1.72547e+07,1.72151e+07,1.74949e+07,1.71321e+07 } * * %Example: Etohit=6900 Detector: EnergyMonitor_first_I=4.24366e+07 * diff --git a/tools/Python/mccodelib/mcplotdiffloader.py b/tools/Python/mccodelib/mcplotdiffloader.py index 283255e667..45ce7ff59f 100644 --- a/tools/Python/mccodelib/mcplotdiffloader.py +++ b/tools/Python/mccodelib/mcplotdiffloader.py @@ -203,6 +203,24 @@ def resolve_simfile(path): return simfile, simdir +def monitor_basename(data): + """ Unique per-monitor file name: the monitor's own output file (e.g. + 'PSD.dat'), or for a scan curve - which all share filename + 'mccode.dat' - the per-curve name the loader gives it (e.g. + 'detector_I'), which is also what mcplot-html names its pages by. """ + return os.path.basename(data.filepath) if data.filepath else data.filename + + +def _dat_name(data, prefix): + name = prefix + monitor_basename(data) + return name if os.path.splitext(name)[1] else name + '.dat' + + +def _total(data): + """ Total intensity; scan curves carry no 'values' triplet, so sum the curve """ + return data.values[0] if data.values else float(np.sum(data.yvals)) + + def load_monitors(path): """ Loads a simulation (folder or single monitor file) and returns (monitors, directory), where monitors is a dict mapping a monitor @@ -227,7 +245,7 @@ def load_monitors(path): if not isinstance(data, (Data1D, Data2D)): # skip 0D / event-list / unknown monitors: nothing sensible to subtract continue - key = data.filename if data.filename else data.component + key = monitor_basename(data) or data.component if key in monitors: # extremely unlikely (duplicate filenames), disambiguate by component key = '%s__%s' % (key, data.component) @@ -301,7 +319,7 @@ def diff_1d(key, a, b, label_a, label_b): N = 0 d.values = (I, Ierr, N) d.statistics = '%s: %s\n%s: %s' % (label_a, a.statistics, label_b, b.statistics) - pct_str = _pct_diff_str(a.values[0], b.values[0]) + pct_str = _pct_diff_str(_total(a), _total(b)) d.diff_pct_str = pct_str d.title = ' - Diff (A - B), with:\n \nA=%s\nB=%s\nDiff: %s\n \n%s' % (label_a, label_b, pct_str, a.title) @@ -536,8 +554,7 @@ def _common_header_lines(data): def _write_1d_dat(data, outdir, prefix): - filename = prefix + data.filename - filepath = os.path.join(outdir, filename) + filepath = os.path.join(outdir, _dat_name(data, prefix)) lines = _common_header_lines(data) # Required: mcplotloader.py's _load_monitor() reads this line to decide @@ -548,9 +565,13 @@ def _write_1d_dat(data, outdir, prefix): lines.append('# xlabel: %s' % _sanitize(data.xlabel)) lines.append('# ylabel: %s' % _sanitize(data.ylabel)) lines.append('# xvar: %s' % _sanitize(data.xvar)) - lines.append('# yvar: (%s,%s)' % (data.yvar[0], data.yvar[1])) + yvar = data.yvar + if isinstance(yvar, str): + # scan curve: single 'detector_I' name rather than an (I, ERR) pair + yvar = (yvar, re.sub(r'_I$', '', yvar) + '_ERR') + lines.append('# yvar: (%s,%s)' % (yvar[0], yvar[1])) lines.append('# xlimits: %s %s' % (_fmt(data.xlimits[0]), _fmt(data.xlimits[1]))) - lines.append('# variables: %s %s %s N' % (data.xvar, data.yvar[0], data.yvar[1])) + lines.append('# variables: %s %s %s N' % (data.xvar, yvar[0], yvar[1])) lines.append('# values: %s %s %s' % (_fmt(data.values[0]), _fmt(data.values[1]), _fmt(data.values[2]))) # The standard "# statistics: X0=...; dX=...;" field is a # weighted centroid/width, which assumes a non-negative intensity @@ -565,15 +586,15 @@ def _write_1d_dat(data, outdir, prefix): with open(filepath, 'w') as f: f.write('\n'.join(lines) + '\n') - for x, y, yerr, n in zip(data.xvals, data.yvals, data.y_err_vals, data.Nvals): + Nvals = data.Nvals or [0] * len(data.xvals) # scan curves have no N column + for x, y, yerr, n in zip(data.xvals, data.yvals, data.y_err_vals, Nvals): f.write('%s %s %s %s\n' % (_fmt(x), _fmt(y), _fmt(yerr), _fmt(n))) return filepath def _write_2d_dat(data, outdir, prefix): - filename = prefix + data.filename - filepath = os.path.join(outdir, filename) + filepath = os.path.join(outdir, _dat_name(data, prefix)) lines = _common_header_lines(data) zshape = np.shape(data.zvals) if data.zvals else (0, 0) @@ -715,7 +736,7 @@ def write_mccode_sim(diffs, outdir, label_a=None, label_b=None, instrument='diff # This is the one line mcplotloader.py's _get_filenames_from_mccodesim() # actually looks for - it must match the real file written by # write_mccode_dat()/write_all_mccode_dat() with the same prefix. - lines.append(' filename: %s' % (prefix + data.filename)) + lines.append(' filename: %s' % _dat_name(data, prefix)) lines.append('end data') lines.append('') diff --git a/tools/Python/mccodelib/mcplotloader.py b/tools/Python/mccodelib/mcplotloader.py index 3033eff03a..d5ffc9adb7 100644 --- a/tools/Python/mccodelib/mcplotloader.py +++ b/tools/Python/mccodelib/mcplotloader.py @@ -483,6 +483,10 @@ def _load_multiplot_1D_lst(f_dat): data.yvar = yvars[2*i] data.title = '%s' % (data.yvar) data.component = data.title + # no per-curve file exists, so give each curve a unique path in + # the scan dir (frontends like mcplot-html derive output names + # from it; empty would make all curves collide on '.html') + data.filepath = join(dirname(f_dat), data.yvar) data.y_err_vals = yvals_err_lst[i] data_handle_lst.append(DataHandle(load_fct=None, data=data)) diff --git a/tools/Python/mccodelib/utils.py b/tools/Python/mccodelib/utils.py index 8dccbd4344..0d00fdf29c 100644 --- a/tools/Python/mccodelib/utils.py +++ b/tools/Python/mccodelib/utils.py @@ -395,7 +395,9 @@ def parse_header(text): tag_E=2 tag_P=3 tag_L=4 - lst = [text.find('%I'), text.find('%D'), text.find('%Ex'), text.find('%P'), text.find('%L')] + # %Scan: lines are tests too, so the test section starts at the first %Example or %Scan + tag_E_pos = min([p for p in (text.find('%Ex'), text.find('%Scan:')) if p != -1], default=-1) + lst = [text.find('%I'), text.find('%D'), tag_E_pos, text.find('%P'), text.find('%L')] # missing %Example tag if lst[tag_E] == -1: lst[tag_E] = lst[tag_P] @@ -493,9 +495,9 @@ def format_examples(test_text): ''' Splits the raw %Example header text (as stored in InstrCompHeaderInfo.test) into an ordered list of (line, is_example) - tuples. Lines starting with the '%Example:' tag have the leading - '%Example:' tag replaced with 'Test:' (so they render as - 'Test: ...') and are flagged True so that doc writers can highlight + tuples. Lines starting with the '%Example:' (or '%Scan:') tag have the + leading tag replaced with 'Test:' (or 'Scan:') (so they render as + 'Test: ...' / 'Scan: ...') and are flagged True so that doc writers can highlight them (e.g. in bold); all other lines are passed through unchanged and flagged False. Order is preserved even when %Example: lines are interspersed with other text. @@ -503,10 +505,12 @@ def format_examples(test_text): lines = [] if not test_text: return lines + # a %Scan: {...} list of target values may span several lines, show it as "{N values}" + test_text = re.sub(r'\{([-0-9.+eE,\s]*)\}', lambda m: '{%d values}' % len([v for v in m.group(1).split(',') if v.strip()]), test_text) for l in test_text.splitlines(): - m = re.match(r'%Example:\s*(.*)', l) + m = re.match(r'%(Example|Scan):\s*(.*)', l) if m: - lines.append(('Test: ' + m.group(1), True)) + lines.append((('Test: ' if m.group(1) == 'Example' else 'Scan: ') + m.group(2), True)) else: lines.append((l, False)) return lines diff --git a/tools/Python/mcplot/html/mccoplot.py b/tools/Python/mcplot/html/mccoplot.py index 4c948cac87..164df71f3a 100644 --- a/tools/Python/mcplot/html/mccoplot.py +++ b/tools/Python/mcplot/html/mccoplot.py @@ -42,7 +42,7 @@ from mccodelib import mccode_config from mccodelib.mcplotdiffloader import ( path_base_name, resolve_labels, resolve_colours, dirsafe_name, - load_monitors, find_original_plot, match_monitors_multi, DEFAULT_PALETTE, + load_monitors, find_original_plot, match_monitors_multi, DEFAULT_PALETTE, monitor_basename, ) # The bare html-plotter (mcplot-html itself, or mxplot-html under McXtrace), @@ -195,7 +195,7 @@ def coplot_single(key, datas, outdir, use_logscale, colours, dat_links, identiti global logscale logscale = use_logscale - basename = 'coplot_' + path_base_name(datas[0].filename) + basename = 'coplot_' + path_base_name(monitor_basename(datas[0])) fname = basename + ('_log.html' if use_logscale else '.html') f = os.path.join(outdir, fname) @@ -447,7 +447,7 @@ def main(args): # mcplotdiff.py does dat_links = [] for d, data in zip(dirs, datas): - lin, log = find_original_plot(d, data.filename) + lin, log = find_original_plot(d, monitor_basename(data)) dat_links.append(_relhref(lin, outdir)) f = coplot_single(key, datas, outdir, False, colours, dat_links, identities, diff --git a/tools/Python/mcplot/html/mcplot.py b/tools/Python/mcplot/html/mcplot.py index 0f453a7053..fdbd5a6226 100644 --- a/tools/Python/mcplot/html/mcplot.py +++ b/tools/Python/mcplot/html/mcplot.py @@ -374,9 +374,6 @@ def plotfunc_single(data, f = None, use_logscale=False): else: f = data.filepath + ".html" - if os.path.exists(f): - os.remove(f) - # create 1D html if type(data) is Data1D: text = get_html('template_1d.html', get_params_str_1D(data), os.path.basename(data.filename), use_logscale) diff --git a/tools/Python/mcplot/html/mcplotdiff.py b/tools/Python/mcplot/html/mcplotdiff.py index 59ce193f13..3d9d3a0cf9 100644 --- a/tools/Python/mcplot/html/mcplotdiff.py +++ b/tools/Python/mcplot/html/mcplotdiff.py @@ -36,7 +36,7 @@ from mccodelib import mcplotdiffloader as diffloader from mccodelib.mcplotdiffloader import ( path_base_name, default_labels, dirsafe_name, load_monitors, compute_diffs, find_original_plot, - write_mccode_dat, write_mccode_sim, + write_mccode_dat, write_mccode_sim, monitor_basename, ) global WIDTH, HEIGHT @@ -280,7 +280,7 @@ def plot_diff_single(data, outdir, use_logscale, dat_basename=None): global logscale logscale = use_logscale - basename = 'diff_' + path_base_name(data.filename) + basename = 'diff_' + path_base_name(monitor_basename(data)) fname = basename + ('_log.html' if use_logscale else '.html') f = os.path.join(outdir, fname) @@ -522,8 +522,8 @@ def main(args): # locate any pre-existing mcplot-html plots for this monitor, so # we can link to the original a/b data alongside the diff plot - a_lin, a_log = find_original_plot(dir_a, data.filename) - b_lin, b_log = find_original_plot(dir_b, data.filename) + a_lin, a_log = find_original_plot(dir_a, monitor_basename(data)) + b_lin, b_log = find_original_plot(dir_b, monitor_basename(data)) if not (a_lin or a_log): print("Note: no existing mcplot-html output found for '%s' in '%s'" % (data.filename, args.a)) if not (b_lin or b_log): diff --git a/tools/Python/mcrun.md b/tools/Python/mcrun.md index fc72515ecc..96d0420acf 100644 --- a/tools/Python/mcrun.md +++ b/tools/Python/mcrun.md @@ -34,7 +34,7 @@ in help text, and the McStas-only `-g`/`--gravitation` flag, differ). | `-N NP`, `--numpoints NP` | Set number of scan points. Two input modes available: 1) A single integer applies the same point count to every scanned parameter (default, and only valid form without -M) 2) Together with -M/--multi, a comma-separated list (e.g. -N=5,10,20) gives each scanned parameter its own point count, in the order in which parameters are listed on the command line. If a parameter is given as par="min:delta:max" the point count is instead computed from the requested bin width. | | `-L`, `--list` | Use list-mode scanning. Multiple input modes available: 1) If multiple lists (of identical length) are given (and -M is not requested) the lists are scanned together in lockstep. 2) Combined with -M/--multi, the cartesian product of each parameter's own list is used to set up a multidimensional 'grid' scan (lists may have different lengths) 3) Any parameter given as "min:delta:max" is expanded into its own explicit list of equidistant points and may be freely mixed with other, explicitly-listed parameters (e.g. a list of filenames) under -L. | | `-M`, `--multi` | Run a multi-dimensional scan (cartesian product of every scanned parameter's points, rather than a co-linear scan). Combine with -L/--list or give -N as a comma-separated list (see -N/--numpoints). | -| `--scan_split scan_split` | Scan by parallelising steps as individual cpu threads. Initialise by number of wanted threads (e.g. your number of cores). | +| `--scan_split scan_split` | Scan by parallelising steps as individual cpu threads. Initialise by number of wanted threads (e.g. your number of cores), or "auto" (or 0) for the number of cores minus one. | | `--seeds SEEDS` | Set range of seeds to scan (each must be: SEED != 0) | | `--optimize` | Optimize instrument variable parameters to maximize monitors | | `--optimize-maxiter optimize_maxiter` | Maximum number of optimization iterations to perform. Default=1000 | @@ -60,7 +60,7 @@ in help text, and the McStas-only `-g`/`--gravitation` flag, differ). | `--vecsize VECSIZE` | vector length in OpenACC parallel scenarios | | `--numgangs NUMGANGS` | number of 'gangs' in OpenACC parallel scenarios | | `--gpu_innerloop INNERLOOP` | Maximum particles in an OpenACC kernel run. (If INNERLOOP is smaller than ncount we repeat) | -| `-s SEED`, `--seed SEED` | Set random seed (must be: SEED != 0) | +| `-s SEED`, `--seed SEED` | Set random seed (must be: SEED != 0). In a scan, point i uses SEED+i*1024 (without --seed, the base seed is taken from the current time and logged) | | `-n COUNT`, `--ncount COUNT` | Set number of neutrons to simulate | | `-t trace`, `--trace trace` | Enable trace of neutrons through instrument | | `--no-trace` | Disable trace of neutrons in instrument (combine with -c) | diff --git a/tools/Python/mcrun/mcrun.py b/tools/Python/mcrun/mcrun.py index cd5bb77583..a62cb694d4 100644 --- a/tools/Python/mcrun/mcrun.py +++ b/tools/Python/mcrun/mcrun.py @@ -143,9 +143,9 @@ def add_mcrun_adv_options(parser): '(see -N/--numpoints). ') add("--scan_split", - type=int, metavar="scan_split", - help='Scan by parallelising steps as individual cpu threads. Initialise by number of wanted threads (e.g. your number of cores).') + help='Scan by parallelising steps as individual cpu threads. Initialise by number of wanted threads (e.g. your number of cores), ' + 'or "auto" (or 0) for the number of cores minus one.') add('--seeds', metavar='SEEDS', @@ -323,7 +323,8 @@ def add_mcstas_options(parser): add('-s', '--seed', metavar='SEED', type=int, action='callback', callback=check_seed, - help='Set random seed (must be: SEED != 0)') + help='Set random seed (must be: SEED != 0). In a scan, point i uses SEED+i*1024 ' + '(without --seed, the base seed is taken from the current time and logged)') add('-n', '--ncount', metavar='COUNT', type=float, default=1000000, @@ -904,8 +905,12 @@ def main(): interval_points = LinearInterval.from_range(options.numpoints, intervals) - # Check that mpi and scan split are not both used. Default to mpi if they are - if options.scan_split is not None and options.mpi is not None: + # Check that mpi and scan split are not both used. Default to mpi if they are. + # --mpi defaults to 1 (compile with MPI, run single-process), so only an + # explicit multi-process --mpi request should override --scan_split + if options.scan_split is not None and options.mpi is not None and str(options.mpi) != '1': + LOG.warning('--mpi=%s given, ignoring --scan_split=%s (scan points run serially via MPI)', + options.mpi, options.scan_split) options.scan_split = None # Parameters for linear scanning present @@ -917,8 +922,9 @@ def main(): scanner.run() # in optimisation.py elif options.scan_split is not None: - if options.scan_split == 0: - options.scan_split = multiprocessing.cpu_count()-1 + if options.scan_split in ("auto", "0"): + options.scan_split = max(1, multiprocessing.cpu_count()-1) + options.scan_split = int(options.scan_split) split_scanner = Scanner_split(mcstas, intervals, options.scan_split) split_scanner.set_points(interval_points) if (not options.dir == ''): diff --git a/tools/Python/mcrun/optimisation.py b/tools/Python/mcrun/optimisation.py index b1556124da..0fc2fbea57 100644 --- a/tools/Python/mcrun/optimisation.py +++ b/tools/Python/mcrun/optimisation.py @@ -392,8 +392,26 @@ def from_list(intervals): class InvalidInterval(McRunException): pass +def point_seed(base, i): + """ Seed for scan point i: the base seed shifted by i*1024, so every point + gets its own random sequence, reproducibly (a seed of 0 is not allowed) """ + return (base + i*1024) or 1 + +def scan_base_seed(mcstas, intervals): + """ The base seed for point_seed(): --seed if given, else the current Unix + epoch (logged, so the scan can be reproduced). None when --seed is itself + scanned (--seeds), in which case the scanned seeds are used unchanged. """ + if '--seed' in intervals: + return None + if mcstas.options.seed is None: + mcstas.options.seed = int(datetime.now().timestamp()) + LOG.info('No incoming seed from cmdline, scan base seed set from current Unix epoch: %d' % mcstas.options.seed) + else: + LOG.info('Scan base seed: %d' % mcstas.options.seed) + return mcstas.options.seed + def _simulate_point(args): - i, point, intervals, mcstas_config, mcstas_dir = args + i, point, intervals, mcstas_config, mcstas_dir, base_seed = args from shutil import copyfile from os.path import join @@ -405,8 +423,8 @@ def _simulate_point(args): # Ensure we get a mccode.sim pr. thread subdir (e.g. for monitoring seed value mcstas.simfile = join(mcstas_dir, 'mccode.sim') - # Shift thread seed to avoid duplicate simulations / biasing - mcstas.options.seed = (i*1024)+mcstas.options.seed + if base_seed is not None: + mcstas.options.seed = point_seed(base_seed, i) for key in intervals: mcstas.set_parameter(key, point[key]) @@ -487,9 +505,12 @@ def run(self): points = list(self.points) header_written = False skipped = [] + base_seed = scan_base_seed(self.mcstas, self.intervals) with open(self.outfile, 'w') as outfile: for i, point in enumerate(points): + if base_seed is not None: + self.mcstas.options.seed = point_seed(base_seed, i) par_values = [] for key in self.intervals: self.mcstas.set_parameter(key, point[key]) @@ -637,14 +658,11 @@ def run(self): mcstas_dir = self.mcstas.options.dir or '.' - if self.mcstas.options.seed is None: - dt=datetime.now() - LOG.info('No incoming seed from cmdline, setting to current Unix epoch (%d)!' % dt.timestamp()) - self.mcstas.options.seed=dt.timestamp() + base_seed = scan_base_seed(self.mcstas, self.intervals) # Prepare data to pass into processes args_list = [ - (i, point, self.intervals, self.mcstas, mcstas_dir) + (i, point, self.intervals, self.mcstas, mcstas_dir, base_seed) for i, point in enumerate(self.points) ] diff --git a/tools/Python/mctest.md b/tools/Python/mctest.md index 6c8ff10c1f..c42da3ed7f 100644 --- a/tools/Python/mctest.md +++ b/tools/Python/mctest.md @@ -10,6 +10,17 @@ Runs every `%Example` test embedded in the instrument library, compiling, displaying (single-particle), and running each one, then compares the result against the target value recorded in the instrument header. +A `%Scan:` header line works the same way for a parameter scan. The +parameters use the `mcrun` scan syntax, followed by one target value per +scan point, e.g. + + %Scan: dBz=-0.0001,0.0001 -N41 -n1e5 Detector: detector_I=79.2,102.7,... + +`mcrun` and `*.instr` tokens on `%Example:`/`%Scan:` lines are ignored. A scan +uses its own `-n` if given, otherwise at most `1e5` per point, and runs its +points in parallel (`--scan_split=auto`) unless MPI is used. A scan passes when +every point is within 20% of its target. + | Option | Description | |---|---| | `--ncount N`, `-n N` | ncount sent to `mcrun` (default: `1e6`) | @@ -32,6 +43,7 @@ result against the target value recorded in the instrument header. | `--compilemax S` | max seconds allowed per compilation (default: 1800; x100 with `--lint`) | | `--runmax S` | max seconds allowed per test run (default: 3600) | | `--displaymax S` | max seconds allowed per test display run (default: 60) | +| `--noscans` | skip the `%Scan` tests, only run the `%Example` tests | | `--noplots` | do not generate plots (`01_overview.pdf`, `02_plots.html`) of the test output | | `--permissive` | exit 0 even if some tests fail | | `--strict` | let instruments without `%Example` line(s) fail immediately (cannot be combined with `--permissive`) | diff --git a/tools/Python/mctest/main.template b/tools/Python/mctest/main.template index 53c9555a8b..d893bd4052 100644 --- a/tools/Python/mctest/main.template +++ b/tools/Python/mctest/main.template @@ -47,7 +47,7 @@ a:visited { color: lightgrey; background-color: transparent; text-decoration: no {%- for c in row %} {%- if loop.index0 > 0 %} {%- if c.0 == 1 %} - C:{{ c.1 }}/R:{{ c.2 }}
I={{ c.3 }} ({{ c.4 }}) + C:{{ c.1 }}/R:{{ c.2 }}
I={{ c.3 }} ({{ c.4 }})
{{ c.6 }}
PLOTS: [Interactive HTML] [PDF overview] {%- if c.8 %} @@ -58,7 +58,7 @@ a:visited { color: lightgrey; background-color: transparent; text-decoration: no {%- endif %} {%- elif c.0 == 2 %} - C:{{ c.1 }}/R:{{ c.2 }}
I={{ c.3 }} ({{ c.4 }}) + C:{{ c.1 }}/R:{{ c.2 }}
I={{ c.3 }} ({{ c.4 }})
{{ c.6 }}
PLOTS: [Interactive HTML] [PDF overview] {%- if c.8 %} diff --git a/tools/Python/mctest/mctest.py b/tools/Python/mctest/mctest.py index e359a99595..6bc4c6f5af 100644 --- a/tools/Python/mctest/mctest.py +++ b/tools/Python/mctest/mctest.py @@ -42,30 +42,72 @@ def paths_overlap(a: pathlib.Path, b: pathlib.Path) -> bool: # Functionality # -def create_instr_test_objs(sourcefile, localfile, header): - ''' returns a list containing one initialized test object pr %Example within the instr file ''' +def split_test_line(line): + ''' Splits the parameter part of a %Example/%Scan line into (parvals, ncount). + "mcrun"/"mxrun" and "*.instr" tokens are dropped, since the test setup defines those, + and a -n/--ncount is taken out of parvals and returned separately (or None). ''' + toks = line.split() + parvals = [] + ncount = None + i = 0 + while i < len(toks): + t = toks[i] + m = re.match(r"(-n|--ncount)=?(.*)$", t) + if m: + ncount = m.group(2) + if not ncount and i + 1 < len(toks): + i += 1 + ncount = toks[i] + elif t not in ("mcrun", "mxrun") and not t.endswith(".instr"): + parvals.append(t) + i += 1 + return " ".join(parvals), ncount + +def scan_first_point(parvals): + ''' parvals for a single run at the first point of a %Scan, e.g. for mcdisplay: + scan options (-N, -L, -M, ...) are dropped and each "par=a,b,..." or "par=a:delta:b" + is reduced to "par=a" ''' + out = [] + for t in parvals.split(): + if t.startswith("-") or "=" not in t: + continue + key, value = t.split("=", 1) + numeric = re.fullmatch(r"[0-9.eE+:,-]+", value) + out.append(key + "=" + re.split(r"[:,]" if numeric else r"(? 0: - testnb = 1 - for m in ms: - parvals = m[0].strip() - detector = m[1].strip() - targetval = float(m[2].strip()) - tests.append(InstrExampleTest(sourcefile, localfile, parvals, detector, targetval, testnb)) - testnb = testnb + 1 - else: + for m in re.findall(r"\%Example:([^\n]*)Detector\:([^\n]*)_I=([0-9.+-e]+)", header): + parvals, ncount = split_test_line(m[0]) + tests.append(InstrExampleTest(sourcefile, localfile, parvals, m[1].strip(), float(m[2].strip()), len(tests) + 1, ncount=ncount)) + if not noscans: + # the target values are a {}-enclosed, comma-separated list that may span several header lines + for m in re.findall(r"\%Scan:([^\n]*)Detector\:([^\n]*)_I=\{([^}]*)\}", header): + parvals, ncount = split_test_line(m[0]) + targetvals = [float(v) for v in m[2].replace("*", " ").split(",") if v.strip()] + tests.append(InstrExampleTest(sourcefile, localfile, parvals, m[1].strip(), targetvals, len(tests) + 1, scan=True, ncount=ncount)) + if not tests: tests.append(InstrExampleTest(sourcefile, localfile)) return tests class InstrExampleTest: - ''' instruent test house keeping object ''' - def __init__(self, sourcefile, localfile, parvals=None, detector=None, targetval=None, testnb=0): + ''' instruent test house keeping object, for a %Scan targetval and testval are lists ''' + def __init__(self, sourcefile, localfile, parvals=None, detector=None, targetval=None, testnb=0, scan=False, ncount=None): self.sourcefile = sourcefile self.localfile = localfile self.instrname = splitext(basename(sourcefile))[0] self.testnb = testnb - + self.scan = scan + self.ncount = ncount + self.parvals = parvals self.detector = detector self.targetval = targetval @@ -87,6 +129,8 @@ def get_json_repr(self): "localfile" : self.localfile, "instrname" : self.instrname, "testnb" : self.testnb, + "scan" : self.scan, + "ncount" : self.ncount, "parvals" : self.parvals, "detector" : self.detector, @@ -116,6 +160,8 @@ def load(self,testnb=0): self.localfile=obj['localfile'] self.instrname=obj['instrname'] self.testnb=obj['testnb'] + self.scan=obj.get('scan', False) + self.ncount=obj.get('ncount') self.parvals=obj['parvals'] self.detector=obj['detector'] self.targetval=obj['targetval'] @@ -208,6 +254,28 @@ def extract_testvals(datafolder, monitorname): return (I, I_err, N) break +def extract_scanvals(datafolder, monitorname): + ''' + Extract the list of monitor I values pr. scan point from the mccode.dat that mcrun + writes for a scan (in both the McCode and NeXus formats). + + Returns the list, or None if mccode.dat or the monitor column is missing. + ''' + datfile = join(datafolder, "mccode.dat") + if not os.path.isfile(datfile): + return None + column = None + vals = [] + for l in open(datfile).read().splitlines(): + if l.startswith("# variables:"): + variables = l.split(":", 1)[1].split() + if monitorname + "_I" not in variables: + return None + column = variables.index(monitorname + "_I") + elif l.strip() and not l.startswith("#") and column is not None: + vals.append(float(l.split()[column])) + return vals if column is not None else None + def parse_detector_I_value(resfile_path, detector_name): """ Return (value_float, success_bool, raw_value_str_or_None). @@ -283,7 +351,7 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt text = open(f, encoding='utf-8').read() f_new=str(pathlib.Path(join(instrdir,os.path.basename(f))).as_posix()) # create a test object for every test defined in the instrument header - instrtests = create_instr_test_objs(sourcefile=f, localfile=f_new, header=text) + instrtests = create_instr_test_objs(sourcefile=f, localfile=f_new, header=text, noscans=noscans) try: shutil.copytree(os.path.dirname(f),instrdir) tests = tests + instrtests @@ -384,7 +452,8 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt # Run mcdisplay (single particle only) t1 = time.time() if test.testnb>0: - cmd = mccode_config.configuration["MCDISPLAY"]+'-classic %s --nobrowse %s %s -n0 -d display > displaylog.txt 2>&1' % (mpiswitch, test.instrname+'.instr', test.parvals if test.parvals else '-y') + dispvals = scan_first_point(test.parvals) if test.scan else test.parvals + cmd = mccode_config.configuration["MCDISPLAY"]+'-classic %s --nobrowse %s %s -n0 -d display > displaylog.txt 2>&1' % (mpiswitch, test.instrname+'.instr', dispvals if dispvals else '-y') else: cmd = mccode_config.configuration["MCDISPLAY"]+'-classic %s --nobrowse %s -y -n0 -d display > displaylog.txt 2>&1' % (mpiswitch, test.instrname+'.instr') retcode = utils.run_subtool_noread(cmd, cwd=join(testdir, test.instrname), timeout=displaymax) @@ -460,6 +529,15 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt # told to use the default values (-y can not be combined with # parameter values, since the instrument then ignores those): parvals = test.parvals if test.parvals else "-y" + # A %Scan uses its own -n if given, else at most 1e5 pr. point, and runs + # its points in parallel unless MPI already does that within each point + testncount = ncount + timeout = runmax + if test.scan: + testncount = test.ncount or "%g" % min(float(ncount), 1e5) + timeout = runmax * len(test.targetval) + if mpi is None: + parvals = parvals + " --scan_split=auto" # Did test run already? if not os.path.exists(join(testdir, test.instrname, str(test.testnb))): if nexus: @@ -468,23 +546,25 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt if openacc is True: if configdir: cmd = cmd + " --override-config=" + configdir - cmd = cmd + " -s %s %s %s -n%s --openacc --mpi=%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, ncount, mpi, test.testnb, test.testnb) + cmd = cmd + " -s %s %s %s -n%s --openacc --mpi=%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, testncount, mpi, test.testnb, test.testnb) else: if configdir: cmd = cmd + " --override-config=" + configdir - cmd = cmd + " -s %s %s %s -n%s --mpi=%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, ncount, mpi, test.testnb, test.testnb) + cmd = cmd + " -s %s %s %s -n%s --mpi=%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, testncount, mpi, test.testnb, test.testnb) else: if configdir: cmd = cmd + " --no-mpi --override-config=" + configdir - cmd = cmd + " --no-mpi -s %s %s %s -n%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, ncount, test.testnb, test.testnb) + cmd = cmd + " --no-mpi -s %s %s %s -n%s -d%d > run_stdout_%d.txt 2>&1" % (seed, test.instrname, parvals, testncount, test.testnb, test.testnb) - retcode = utils.run_subtool_noread(cmd, cwd=join(testdir, test.instrname),timeout=runmax) + retcode = utils.run_subtool_noread(cmd, cwd=join(testdir, test.instrname),timeout=timeout) t2 = time.time() didwrite = os.path.exists(join(testdir, test.instrname, str(test.testnb), "mccode.sim")) didwrite_nexus = os.path.exists(join(testdir, test.instrname, str(test.testnb), "mccode.h5")) # retcode is a tuple: (returncode, timed_out) - test.didrun = retcode[0] == 0 and not retcode[1] and (didwrite or didwrite_nexus) + # a scan always writes mccode.dat, but not always a top-level mccode.sim/mccode.h5 (NeXus + --scan_split) + didwrite_scan = test.scan and os.path.exists(join(testdir, test.instrname, str(test.testnb), "mccode.dat")) + test.didrun = retcode[0] == 0 and not retcode[1] and (didwrite or didwrite_nexus or didwrite_scan) test.runtime = t2 - t1 else: suffix=" (cached)" @@ -506,7 +586,12 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt resbase="(No file)" # test value extraction - if not didwrite_nexus: + if test.scan: + test.testval = extract_scanvals(join(testdir, test.instrname, str(test.testnb)), test.detector) + if test.testval is None: + runfailed=True + resbase ="run_stdout_%d.txt" % (test.testnb) + elif not didwrite_nexus: extraction = extract_testvals(join(testdir, test.instrname, str(test.testnb)), test.detector) if type(extraction) is tuple: test.testval = extraction[0] @@ -535,7 +620,18 @@ def mccode_test(branchdir, testdir, limitinstrs=None, instrfilter=None, compfilt suffix += " + !! RUNTIME FAILURE - see %s !! " % (resbase) formatstr = "%-" + "%ds: " % (maxnamelen+1) + \ "{:3d}.".format(math.floor(test.runtime)) + str(test.runtime-int(test.runtime)).split('.')[1][:2] - if test.targetval!=0: # Normal situation, non-zero target value + if test.scan: + testvals = test.testval or [] + percents = [round(percent_of(t, r)) for t, r in zip(testvals, test.targetval)] + numoff = len([p for p in percents if p<80 or p>120]) + if len(testvals) != len(test.targetval) or numoff > 0: + suffix += " <--- BIG DISCREPANCY??" + num_valfail = num_valfail + 1 + anyfailed=True + worst = max(percents, key=lambda p: abs(p-100)) if percents else 0 + logging.info(formatstr % test.get_display_name() + " [scan: %d/%d points within 20%%, worst %d %%, %d/%d points run]" + % (len(percents)-numoff, len(test.targetval), worst, len(testvals), len(test.targetval)) + suffix) + elif test.targetval!=0: # Normal situation, non-zero target value percent=round(100.0*test.testval/test.targetval) if percent<80 or percent>120: suffix += " <--- BIG DISCREPANCY??" @@ -792,6 +888,7 @@ def get_config_files(configfltr): compilemax = None displaymax = None noplots = None +noscans = None def main(args): configfilter = args.config # test only config matching this label (default: as installed) @@ -849,7 +946,7 @@ def main(args): quit(1) logging.debug("") - global ncount, no_mpi, mpi, skipnontest, openacc, nexus, lint, permissive, runLocal, compilemax, displaymax, runmax, seed, strict, noplots + global ncount, no_mpi, mpi, skipnontest, openacc, nexus, lint, permissive, runLocal, compilemax, displaymax, runmax, seed, strict, noplots, noscans ncount = "1e6" no_mpi = False if args.ncount: @@ -968,6 +1065,9 @@ def main(args): logging.info("Strict mode, tool will report failure for instruments without %Example") noplots = args.noplots + noscans = args.noscans + if noscans: + logging.info("%Scan tests are skipped") if noplots: logging.info("No plots of the test output will be generated") @@ -1003,6 +1103,7 @@ def main(args): parser.add_argument('--permissive', action='store_true', help='Use zero return-value even if some tests fail. Useful for full test con systems that are only partially functional. Can not be combined with --strict.') parser.add_argument('--strict', action='store_true', help='Let instruments without %%Example line(s) instantly fail. Can not be combined with --permissive.') parser.add_argument('--noplots', action='store_true', help='Do not generate plots (01_overview.pdf and 02_plots.html) of the test output. Useful e.g. in CI, where the plots are not looked at, and can take long for instruments with many monitors.') + parser.add_argument('--noscans', action='store_true', help='Skip the %%Scan tests, only run the %%Example tests.') parser.add_argument('--local', help='Instruments to test are NOT picked up from MCCODE installation, instead from --local=DIR. Local path and --testdir can not overlap!') args = parser.parse_args() diff --git a/tools/Python/mctest/mcviewtest.py b/tools/Python/mctest/mcviewtest.py index dfe010a913..5befec6573 100644 --- a/tools/Python/mctest/mcviewtest.py +++ b/tools/Python/mctest/mcviewtest.py @@ -17,6 +17,12 @@ sys.path.append(os.path.join(os.path.dirname(__file__), '..')) from mccodelib import utils, mccode_config +def percent_of(testval, targetval): + ''' testval in percent of targetval, a 0 target counts as 100% only for a 0 testval ''' + if targetval == 0: + return 100 if testval == 0 else 0 + return 100.0 * testval / targetval + def scantree(path): """Recursively yield DirEntry objects for given directory.""" for entry in os.scandir(path): @@ -249,6 +255,24 @@ def get_cell_tuple(cellobj, refval=None, refcellobj=None, row=None, col_idx=None compiletime = "" state = 2 return (state, compiletime, runtime, testval, "", url, display, displayurl) + elif cellobj.get("scan"): + # one cell pr. %Scan, per-point details in the hover tooltip + runtime = "%.2f s" % cellobj["runtime"] + compiletime = "%.2f s" % cellobj["compiletime"] if cellobj["testnb"] <= 1 else "" + testvals = cellobj["testval"] + refvals = cellobj["targetval"] + percents = [percent_of(t, r) for t, r in zip(testvals, refvals)] + numok = len([p for p in percents if abs(p-100) <= ERROR_PERCENT_THRESSHOLD_ACCEPT]) + state = 1 if numok == len(percents) == len(refvals) else 2 + testval = "scan, %d/%d pts" % (len(testvals), len(refvals)) + refp = "%d/%d pts OK" % (numok, len(refvals)) + tooltip = "\n".join("%d: %.4g / %.4g = %.0f%%" % (i, t, r, p) for i, (t, r, p) in enumerate(zip(testvals, refvals, percents))) + diffurl = None + coplotUrl = None + if diffall or state == 2: + diffurl = plan_diff_link(refcellobj, url, label, row, col_idx) + coplotUrl = plan_coplot_link(refcellobj, url, label, row, col_idx) + return (state, compiletime, runtime, testval, refp, url, display, displayurl, diffurl, coplotUrl, tooltip) else: testval = "%.2g" % float(cellobj["testval"]) runtime = "%.2f s" % cellobj["runtime"] @@ -259,14 +283,7 @@ def get_cell_tuple(cellobj, refval=None, refcellobj=None, row=None, col_idx=None # Always use embedded target value refval = float(cellobj["targetval"]) testval = float(cellobj["testval"]) - # Special case, target test value is 0 explicitly: - if refval==0: - if testval==0: - refp=100 - else: - refp=0 - else: # Standard case, target test value is non-zero - refp = abs(testval/refval*100) + refp = abs(percent_of(testval, refval)) if abs(refp-100) > ERROR_PERCENT_THRESSHOLD_ACCEPT: state = 2 else: