diff --git a/tools/Python/mccodelib/mcplotloader.py b/tools/Python/mccodelib/mcplotloader.py index 1828e7577..3033eff03 100644 --- a/tools/Python/mccodelib/mcplotloader.py +++ b/tools/Python/mccodelib/mcplotloader.py @@ -502,7 +502,7 @@ def walkfunc(arg, dirname, fnames): dirsignature = (dirname, mnames) for f in fnames: # NOTE: this will attempt to load all files except for mccode.sim - if f not in mnames and f != 'mccode.sim' and f != 'mcstas.sim': + if f not in mnames and f != 'mccode.sim': mnames.append(f) arg.append(dirsignature) @@ -512,16 +512,35 @@ def walkfunc(arg, dirname, fnames): walkfunc(subdirtuple, root, files) del subdirtuple[0] # remove root dir subdirs = [t[0] for t in subdirtuple] - # get the right order of subdirs by recreating them a little bit - subdirs = [join(dirname(subdirs[i]), str(i)) for i in range(len(subdirs))] # sortalpha(subdirs) + # Sort by the REAL, actual numeric folder name (e.g. "0", "1", "3", + # "4" - numerically, not alphabetically, since alphabetical sort would + # put "10" before "2") - rather than the previous approach of + # renaming every subdir BY ITS POSITION in os.walk()'s traversal order + # (`join(dirname(subdirs[i]), str(i))`). That renaming assumed + # scan-step subfolders are always named consecutively with no gaps, + # which no longer holds now that a failed scan step's subfolder is + # deliberately never created (mcrun's Scanner.run()/Scanner_split now + # tolerate individual failed points rather than aborting the whole + # scan - see tools/Python/mcrun/optimisation.py). A gap used to + # silently rename every subfolder from that point on to a DIFFERENT, + # wrong path, misattributing each remaining step's data to the wrong + # index - surfacing as anything from wrong data to the "list index + # out of range" this caused (a later step's own mccode.sim declaring + # fewer monitors than the one it got misaligned with expected). + try: + subdirs = sorted(subdirs, key=lambda p: int(basename(p))) + except ValueError: + # Subfolder names aren't purely numeric for some reason - fall + # back to a plain alphabetical ordering rather than crashing + # outright (matches this function's previous fallback comment, + # "sortalpha(subdirs)", which was never actually reachable before). + subdirs = sorted(subdirs) # get the monitor ordering right by snooping the ' filename:' labels out of the scan point file 0/mccode.sim def get_subdir_monitors(subdir): mons = [] if exists(join(subdir, 'mccode.sim')): indexfile='mccode.sim' - elif exists(join(subdir, 'mccode.sim')): - indexfile='mcstas.sim' else: return @@ -544,7 +563,30 @@ def get_subdir_monitors(subdir): monitors_by_subdir = [] for s in subdirs: - monitors_by_subdir.append(get_subdir_monitors(s)) + mons = get_subdir_monitors(s) + if mons: + monitors_by_subdir.append(mons) + else: + # A discovered subfolder without usable monitor data - either + # no mccode.sim at all (get_subdir_monitors() returns None), + # or a mccode.sim that exists but has zero "begin data" blocks + # (returns [] - e.g. the simulation crashed after writing its + # header but before writing monitor results + # Either way, skip it here rather than letting an empty (or + # None) entry propagate into monitors_by_subdir and crash the + # indexing below; load_sweep()'s own root/secondary + # length-mismatch check further down already warns (rather + # than crashing) if this ever leaves the secondary monitor + # count out of sync with mccode.dat's row count. + print("_load_sweep_monitors: skipping subdir %s (no usable monitor data found)" % s) + + if not monitors_by_subdir: + # No subfolder had usable data at all (e.g. every scan step + # failed) - return no secondary monitors rather than crashing on + # monitors_by_subdir[0] below. load_sweep() still has the root + # sweep curves from mccode.dat either way; it just won't have a + # per-step drill-down to offer. + return [] # notice that columns and rows are swapped, so we get to use a # list-of-lists data structure, with rows the same monitor @@ -599,7 +641,7 @@ def has_filename(args): def is_mccodesim_or_mccodedat(args): f = args['simfile'] f_name = basename(f) - return (f_name == 'mccode.sim' or f_name == 'mcstas.sim' or f_name == 'mccode.dat') and isfile(f) + return (f_name == 'mccode.sim' or f_name == 'mccode.dat') and isfile(f) def is_monitorfile(args): @@ -623,13 +665,24 @@ def is_sweepfolder(args): def is_broken_sweepfolder(args): - ''' not implemented (returns trivial answer) ''' - return False - - -def is_sweep_data_present(args): - ''' not implemented ''' - raise Exception('is_sweep_data_present has not been implemented.') + ''' A sweep/scan directory that has mccode.dat (the combined scan-curve + file written by mcrun's Scanner/Scanner_split/Optimizer - see + tools/Python/mcrun/optimisation.py) but is missing mccode.sim - + e.g. because only mccode.dat was kept or shared on its own, or the + per-scan-step subfolders and their individual monitor files were + cleaned up/never transferred. is_sweepfolder() (checked just + before this in the flowchart) already requires BOTH files to be + present for the full, richer sweep view (load_sweep(), with its + per-step secondary drill-down) - reaching this function means that + check already failed, so finding mccode.dat here specifically + means mccode.sim must be the one missing. mccode.dat is + self-describing via its own '#'-prefixed header, though (see + _load_multiplot_1D_lst()), so the overlaid sweep curves themselves + can still be recovered and plotted even without mccode.sim or any + of the underlying per-monitor detector files - see + load_sweep_dat_only(). ''' + d = args['directory'] + return isfile(join(d, 'mccode.dat')) def is_mccodesim_w_monitors(args): @@ -638,8 +691,6 @@ def is_mccodesim_w_monitors(args): # checks mccode.sim existence if isfile(join(d, 'mccode.sim')): indexfile='mccode.sim' - elif isfile(join(d, 'mcstas.sim')): - indexfile='mcstas.sim' else: return False @@ -650,8 +701,6 @@ def is_mccodesim_w_monitors(args): datfiles = glob.glob(join(d, '*')) if 'mccode.sim' in datfiles: datfiles.remove('mccode.sim') - if 'mcstas.sim' in datfiles: - datfiles.remove('mcstas.sim') if 'mccode.dat' in datfiles: datfiles.remove('mccode.dat') return len(datfiles) > 0 @@ -664,8 +713,6 @@ def has_datfile(args): datfiles = glob.glob(join(d, '*')) if 'mccode.sim' in datfiles: datfiles.remove('mccode.sim') - if 'mcstas.sim' in datfiles: - datfiles.remove('mcstas.sim') if 'mccode.dat' in datfiles: datfiles.remove('mccode.dat') if len(datfiles) > 0: @@ -685,8 +732,6 @@ def has_multiple_datfiles(args): datfiles = glob.glob(join(d, '*')) if 'mccode.sim' in datfiles: datfiles.remove('mccode.sim') - if 'mcstas.sim' in datfiles: - datfiles.remove('mcstas.sim') if 'mccode.dat' in datfiles: datfiles.remove('mccode.dat') for f in datfiles: @@ -705,7 +750,6 @@ def test_decfuncs(simfile): print('is_monitorfile: %s' % str(is_monitorfile(args))) print('is_sweepfolder: %s' % str(is_sweepfolder(args))) print('is_broken_sweepfolder: %s' % str(is_broken_sweepfolder(args))) - #print('is_sweep_data_present: %s' % str(is_sweep_data_present(args))) # should not be called until implemented print('is_mccodesim_w_monitors: %s' % str(is_mccodesim_w_monitors(args))) print('has_datfile: %s' % str(has_datfile(args))) print('has_multiple_datfiles: %s' % str(has_multiple_datfiles(args))) @@ -733,8 +777,6 @@ def load_simulation(args): # load monitor data handles if isfile(join(d, 'mccode.sim')): indexfile='mccode.sim' - elif isfile(join(d, 'mcstas.sim')): - indexfile='mcstas.sim' else: indexfile='' data_lst = _load_data_from_mcfiles(_get_filenames_from_mccodesim(join(d, indexfile))) @@ -754,8 +796,6 @@ def load_simulation(args): def load_sweep(args): d = args['directory'] f_dat = join(d, 'mccode.dat') - if isfile(join(d, 'mcstas.sim')): - f_dat = join(d, 'mcstas.sim') # load primary data_handle, 1D sweep values data_handle_lst_sweep1D = _load_multiplot_1D_lst(f_dat) @@ -853,6 +893,45 @@ def load_sweep_c(args): raise Exception('load_sweep_c is not implemented.') +def load_sweep_dat_only(args): + ''' Fallback for a sweep/scan directory that has mccode.dat (the + combined scan-curve file) but not mccode.sim and/or the + underlying per-step subfolders and monitor files - see + is_broken_sweepfolder(). mccode.dat is self-describing via its own + '#'-prefixed header (component/filename/title/xvars/xlimits/ + variables/yvars - the same fields build_header() in + tools/Python/mcrun/optimisation.py writes), so + _load_multiplot_1D_lst() can parse it directly with no other + input at all. + + Only builds the root+primary levels (mirroring load_simulation()/ + load_monitor_folder()'s two-level graphs): an overview showing + every monitor's sweep curve overlaid, and a primary per-monitor + drill-down to see just that one curve. There deliberately are no + secondaries here - load_sweep()'s secondary level lets you drill + further into an individual scan step's own raw monitor data + (loaded from that step's own subfolder+mccode.sim), which simply + doesn't exist in this fallback - only the combined sweep curves + do. Leaving secondaries as the default empty list is a normal, + supported plot-graph state (matching how a PNSingle leaf node, or + load_simulation()'s single-level case, has none either), not an + incomplete/degraded one - frontends already handle it as "nothing + further to click into" rather than an error. ''' + d = args['directory'] + f_dat = join(d, 'mccode.dat') + + data_handle_lst_sweep1D = _load_multiplot_1D_lst(f_dat) + root = PNMultiple(data_handle_lst_sweep1D) + + primnodes_lst = [] + for data_handle in data_handle_lst_sweep1D: + primnode = PNSingle(data_handle) + primnodes_lst.append(primnode) + root.set_primaries(primnodes_lst) + + return root + + def load_monitor_folder(args): # assume simfile is folder with multiple dat files d = args['directory'] @@ -892,8 +971,7 @@ def load(self): exit_term_case1 = FCNTerminal(key = "case1", fct = load_monitor) exit_term_case2 = FCNTerminal(key = "case2", fct = load_simulation) exit_term_case3 = FCNTerminal(key = "case3", fct = load_sweep) - exit_term_case3b = FCNTerminal(key = "case3b", fct = throw_error) - exit_term_case3c = FCNTerminal(key = "case3c", fct = throw_error) + exit_term_case3fallback = FCNTerminal(key = "case3-fallback", fct = load_sweep_dat_only) exit_term_case4 = FCNTerminal(key = "case4", fct = load_monitor_folder) # decision nodes (assembled in backwards order) @@ -906,11 +984,8 @@ def load(self): dec_ismccodesimwmonitors = FCNDecisionBool(fct = is_mccodesim_w_monitors, node_T = exit_term_case2, node_F = dec_hasdatfile) - dec_datafolderspresent = FCNDecisionBool(fct = is_sweep_data_present, - node_T = exit_term_case3b, - node_F = exit_term_case3c) dec_isbrokensweep = FCNDecisionBool(fct = is_broken_sweepfolder, - node_T = dec_datafolderspresent, + node_T = exit_term_case3fallback, node_F = dec_ismccodesimwmonitors) dec_issweepfolder = FCNDecisionBool(fct = is_sweepfolder, node_T = exit_term_case3, diff --git a/tools/Python/mcplot/matplotlib/mcplot.py b/tools/Python/mcplot/matplotlib/mcplot.py index 3e5e9a914..05a23c18b 100644 --- a/tools/Python/mcplot/matplotlib/mcplot.py +++ b/tools/Python/mcplot/matplotlib/mcplot.py @@ -45,6 +45,19 @@ def main(args): if (h5file): if not os.path.isabs(h5file): h5file = os.path.join(os.getcwd(),h5file) + # A NeXus-format scan (as opposed to a single simulation) + # also writes a scan-summary mccode.dat alongside + # mccode.h5. If it's there, also spawn a separate, + # independent instance of this same mcplot variant pointed + # directly at mccode.dat, running in the background + # alongside the HDFVIEW. + datfile = os.path.join(os.path.dirname(h5file), 'mccode.dat') + if os.path.isfile(datfile): + try: + print('Also spawning %s on %s' % (os.path.basename(__file__), datfile)) + subprocess.Popen([sys.executable, os.path.abspath(__file__), datfile]) + except Exception as e: + print('Could not launch a second mcplot instance on ' + datfile + ': ' + e.__str__()) try: cmd = mccode_config.configuration['HDFVIEW'] + ' ' + h5file print('Spawning ' + mccode_config.configuration['HDFVIEW']) diff --git a/tools/Python/mcplot/pyqtgraph/mcplot.py b/tools/Python/mcplot/pyqtgraph/mcplot.py index bda1ee06b..9dbc63673 100644 --- a/tools/Python/mcplot/pyqtgraph/mcplot.py +++ b/tools/Python/mcplot/pyqtgraph/mcplot.py @@ -43,6 +43,19 @@ def main(args): if (h5file): if not os.path.isabs(h5file): h5file = os.path.join(os.getcwd(),h5file) + # A NeXus-format scan (as opposed to a single simulation) + # also writes a scan-summary mccode.dat alongside + # mccode.h5. If it's there, also spawn a separate, + # independent instance of this same mcplot variant pointed + # directly at mccode.dat, running in the background + # alongside the HDFVIEW. + datfile = os.path.join(os.path.dirname(h5file), 'mccode.dat') + if os.path.isfile(datfile): + try: + print('Also spawning %s on %s' % (os.path.basename(__file__), datfile)) + subprocess.Popen([sys.executable, os.path.abspath(__file__), datfile]) + except Exception as e: + print('Could not launch a second mcplot instance on ' + datfile + ': ' + e.__str__()) try: cmd = mccode_config.configuration['HDFVIEW'] + ' ' + h5file print('Spawning ' + mccode_config.configuration['HDFVIEW']) diff --git a/tools/Python/mcrun/mcrun.py b/tools/Python/mcrun/mcrun.py index 0e987a192..04ffda612 100644 --- a/tools/Python/mcrun/mcrun.py +++ b/tools/Python/mcrun/mcrun.py @@ -92,8 +92,13 @@ def add_mcrun_options(parser): help='Read parameters from file FILE') add('-N', '--numpoints', - type=int, metavar='NP', - help='Set number of scan points') + metavar='NP', + help='Set number of scan points. A single integer applies the same ' + 'point count to every scanned parameter (the default, and the ' + 'only valid form without -M). With -M/--multi, a comma-separated ' + 'list (e.g. -N=5,10,20) instead gives each scanned parameter its ' + 'own point count, in the same order the parameters are listed ' + 'on the command line.') add('--seeds', metavar='SEEDS', @@ -101,11 +106,20 @@ def add_mcrun_options(parser): add('-L', '--list', action='store_true', - help='Use a fixed list of points for linear scanning') + help='Use a fixed list of points for scanning, walking every scanned ' + 'parameter\'s list together in lockstep (all lists must then be ' + 'the same length). Combine with -M/--multi instead to take the ' + 'cartesian product of each parameter\'s own list (lists may then ' + 'have different lengths).') add('-M', '--multi', action='store_true', - help='Run a multi-dimensional scan') + help='Run a multi-dimensional scan (the cartesian product of every ' + 'scanned parameter\'s own points, rather than walking them all ' + 'in lockstep). Combine with -L/--list (each parameter\'s ' + 'explicit list can then have a different length) or give -N ' + 'a comma-separated list (see -N/--numpoints) for per-parameter ' + 'point counts.') add("--scan_split", type=int, @@ -460,9 +474,17 @@ def get_parameters(options): if '=' in param: key, value = param.split('=', 1) interval = value.split(',') + # Protect against trailing (or doubled) commas (empty-string entries + n_before = len(interval) + interval = [v for v in interval if v != ''] + if len(interval) != n_before: + LOG.warning('Ignoring %d empty value(s) in parameter "%s" ' + '(check for a trailing or doubled comma)', n_before - len(interval), key) # When just one point is present, fix as constant if len(interval) == 1: - fixed_params[key] = value + fixed_params[key] = interval[0] + elif len(interval) == 0: + LOG.warning('Ignoring parameter "%s": no values left after removing empty entries', key) else: LOG.debug('interval[%s]: %s', key, interval) intervals[key] = interval @@ -595,31 +617,89 @@ def main(): if options.list and options.seeds: raise OptionValueError('--seeds cannot be used with --list') + # Parse -N/--numpoints (a plain string now, not auto-int'd by optparse - + # see add_mcrun_options()): with -M/--multi it may be a comma-separated + # list of integers, one per scanned parameter in the same order the + # parameters were given on the command line, rather than a single + # integer applied uniformly to every dimension. A list form without -M + # is rejected outright: a plain (co-linear) scan walks every parameter + # in lockstep over the same number of steps, so per-dimension point + # counts don't apply there. Unreachable when --list was also given, + # thanks to the check just above. + numpoints_list = None + if options.numpoints is not None: + numpoints_parts = str(options.numpoints).split(',') + if len(numpoints_parts) > 1: + if not options.multi: + raise OptionValueError( + 'A comma-separated list for -N/--numpoints (e.g. -N=5,10,20) is only valid ' + 'together with -M/--multi.') + try: + numpoints_list = [int(p) for p in numpoints_parts] + except ValueError: + raise OptionValueError('-N/--numpoints list must contain only integers: "%s"' % options.numpoints) + if any(n < 2 for n in numpoints_list): + raise OptionValueError( + 'Cannot scan using only one data point - every entry in -N/--numpoints must be at least 2.') + options.numpoints = None # resolved into numpoints_list/numpoints_dict instead, below + else: + try: + options.numpoints = int(numpoints_parts[0]) + except ValueError: + raise OptionValueError( + '-N/--numpoints must be an integer (or, with -M, a comma-separated list of integers): "%s"' + % options.numpoints) + if options.list: if len(intervals) == 0: raise OptionValueError( '--list was chosen but no lists was presented.') - pointlist = list(intervals.values()) - points = len(pointlist[0]) - if not (all(map(lambda i: len(i) == points, intervals.values()))): + if options.multi: + # -L + -M: cartesian product across each parameter's own + # explicit list of points - unlike plain -L (which walks every + # list together in lockstep, requiring them all to be the same + # length), each dimension is independent here, so the lists + # may have different lengths. + interval_points = MultiInterval.from_list(intervals) + options.numpoints = 1 + for values in intervals.values(): + options.numpoints *= len(values) + else: + pointlist = list(intervals.values()) + points = len(pointlist[0]) + if not (all(map(lambda i: len(i) == points, intervals.values()))): + raise OptionValueError( + 'All variables must have an equal amount of points.') + interval_points = LinearInterval.from_list( + points, intervals) + options.numpoints = points + + elif numpoints_list is not None: + # -M + -N=a,b,c,...: per-dimension point counts, no explicit lists + if len(numpoints_list) != len(intervals): raise OptionValueError( - 'All variables must have an equal amount of points.') - interval_points = LinearInterval.from_list( - points, intervals) + '-N/--numpoints list has %d entr%s but %d parameter%s being scanned (%s); ' + 'provide exactly one point-count per scanned parameter, in the same order.' % ( + len(numpoints_list), 'y' if len(numpoints_list) == 1 else 'ies', + len(intervals), '' if len(intervals) == 1 else 's are', + ', '.join(intervals))) + numpoints_dict = dict(zip(intervals.keys(), numpoints_list)) + interval_points = MultiInterval.from_range(numpoints_dict, intervals) + total = 1 + for n in numpoints_list: + total *= n + options.numpoints = total - scan = options.multi or options.numpoints - if (options.numpoints is not None and options.numpoints < 2) or (scan and options.numpoints is None): - raise OptionValueError((f'Cannot scan variable(s) {", ".join(intervals)} using only one data point. ' - 'Please use -N to specify the number of points.')) - ## ## This *was* unreachable due to its indentation. Should it be removed entirely? - # # Check that input is valid decimals - # if not all(map(lambda i: len(i) == 2 and all(map(is_decimal, i)), intervals.values())): - # raise OptionValueError(f'Could not parse intervals -- result: {intervals}') + else: + scan = options.multi or options.numpoints + if (options.numpoints is not None and options.numpoints < 2) or (scan and options.numpoints is None): + raise OptionValueError((f'Cannot scan variable(s) {", ".join(intervals)} using only one data point. ' + 'Please use -N to specify the number of points.')) - if options.multi is not None: - interval_points = MultiInterval.from_range(options.numpoints, intervals) - elif options.numpoints is not None: - interval_points = LinearInterval.from_range(options.numpoints, intervals) + if options.multi is not None: + interval_points = MultiInterval.from_range(options.numpoints, intervals) + elif options.numpoints is not None: + interval_points = LinearInterval.from_range(options.numpoints, intervals) # Check that mpi and scan split are not both used. Default to mpi if they are @@ -628,9 +708,6 @@ def main(): # Parameters for linear scanning present if interval_points and (options.scan_split is None): - # In case of list, update with number of list points - if options.list: - options.numpoints=len(pointlist[0]) scanner = Scanner(mcstas, intervals) scanner.set_points(interval_points) if (not options.dir == ''): diff --git a/tools/Python/mcrun/optimisation.py b/tools/Python/mcrun/optimisation.py index 1543dbaaa..670cb8d83 100644 --- a/tools/Python/mcrun/optimisation.py +++ b/tools/Python/mcrun/optimisation.py @@ -5,6 +5,7 @@ from os.path import join from multiprocessing import Pool import copy +import re try: from scipy.optimize import minimize @@ -16,6 +17,37 @@ LOG = getLogger('optimisation') +def _list_scan_xlimits(lst): + """ Computes the (xmin, xmax) header hint for an -L/--list scan's + first scanned parameter, matching whatever will actually end up + plotted as that parameter's x-values. + + A genuinely numeric list (the common case, e.g. -L lambda=2,3) + uses its own real min/max - this MUST match resolve_scan_value()'s + numeric passthrough for the actual per-point column written into + mccode.dat, since the matplotlib frontend's plot_single_data() + uses this value directly via pylab.xlim(xmin, xmax) to set the + visible axis range. + + A non-numeric list (e.g. -L filename=Na2Ca3Al2F14.laz,...) uses + the 0-based index range (0..N-1) instead, matching + resolve_scan_value()'s own index-substitution fallback for that + case - a literal min()/max() of the raw strings would be + lexicographic and meaningless there anyway. """ + try: + numeric_vals = [float(v) for v in lst] + return min(numeric_vals), max(numeric_vals) + except (TypeError, ValueError): + # Non-numeric: resolve_scan_value() substitutes each value with + # its own 0-based index within intervals[key] + # (list(intervals[key]).index(value)), so the matching range is + # (0, len(lst)-1) - NOT (1, len(lst)), which would itself clip the + # first data point (plotted at x=0) outside the visible axis, the + # same class of bug this function exists to avoid for the numeric + # case above. + return 0, len(lst) - 1 + + def build_header(options, params, intervals, detectors): template = """ # Instrument-source: '%(instr)s' @@ -43,8 +75,7 @@ def build_header(options, params, intervals, detectors): xvars = ', '.join(hdrparams) lst = intervals[list(params)[0]] if options.list: - xmin=1 - xmax=len(lst) + xmin, xmax = _list_scan_xlimits(lst) else: xmin = min(lst) xmax = max(lst) @@ -131,6 +162,19 @@ def build_mccodesim_header(options, intervals: dict, detectors: list, version: s # TODO: figure out correct scan type numpoints = 1 if options.optimize else options.numpoints + # -L list scan: use the shared helper, which keeps the real min/max + # for a numeric list (matching what actually gets plotted - see + # _list_scan_xlimits()'s own docstring for why this matters), not just + # for a non-numeric one (e.g. filenames), where a literal min()/max() + # of the raw strings would be lexicographic and meaningless, so a + # 0-based index range (matching resolve_scan_value()'s own index + # substitution) is used instead. Equidistant (-N/-M, non-list) scans + # are untouched, keeping their existing min()/max() behaviour. + if options.list: + xmin, xmax = _list_scan_xlimits(first_key_interval) + else: + xmin, xmax = min(first_key_interval), max(first_key_interval) + values = { 'instr': options.instr, 'date': datetime.strftime(datetime.now(), '%a %b %d %H %M %Y'), @@ -145,8 +189,8 @@ def build_mccodesim_header(options, intervals: dict, detectors: list, version: s 'xvars': interval_names, 'yvars': ' '.join(f'({d}_I,{d}_ERR' for d in detectors), - 'xmin': min(first_key_interval), - 'xmax': max(first_key_interval), + 'xmin': xmin, + 'xmax': xmax, 'filename': basename(options.optimise_file) or 'mccode.dat', 'variables': ' '.join(intervals.keys()) + ' '.join(f'{d}_I {d}_ERR' for d in detectors), @@ -183,12 +227,102 @@ def mcsimdetectors(directory_name: str): return [Detector(d['component'], *d['values'].split(), d['filename'], d['statistics']) for d in blocks] +# Matches one of a simulation binary's own "Detector: ..." summary lines, +# e.g.: +# Detector: PSDbefore_guides_I=2.34581e+09 PSDbefore_guides_ERR=2.34585e+06 PSDbefore_guides_N=999991 "PSDbefore_guides.dat" +# The detector name itself can contain underscores (as in the example +# above), so a plain \w+ before "_I=" isn't reliable - a backreference +# instead requires the SAME name to reappear before "_ERR=" and "_N=", +# which correctly anchors the split point regardless of what characters +# the name itself contains. +DETECTOR_STDOUT_RE = re.compile( + r'Detector:\s*(.+?)_I=(\S+)\s+\1_ERR=(\S+)\s+\1_N=(\S+)\s+"([^"]*)"' +) + + +def parse_detectors_from_stdout(stdout_text): + """ Parses a simulation's own "Detector: NAME_I=... NAME_ERR=... + NAME_N=... "file.dat"" summary lines directly out of its stdout, + and returns them as the same list of Detector objects + mcsimdetectors() builds from a per-step mccode.sim file. + + Needed specifically for --format=NeXus scans: the default McCode + output format writes one mccode.sim/mccode.dat pair per scan step, + each in its own "dir/0", "dir/1", ... subfolder, which + mcsimdetectors() reads back after each step. NeXus format instead + (intentionally) accumulates every step into a single shared .h5 + file (see Scanner.run()'s options.append=True for the NeXus + branch) - so there is no per-step mccode.sim to read detector + values back from at all; mcsimdetectors() finds only a .h5 file + there and returns nothing. The underlying simulation binary still + prints its normal per-run "Detector: ..." summary to stdout + regardless of output format, though, so that's used as the source + of per-step detector values in the NeXus case instead. """ + from mccode import Detector + if not stdout_text: + return [] + detectors = [] + for match in DETECTOR_STDOUT_RE.finditer(stdout_text): + name, intensity, error, count, path = match.groups() + # Detector()'s "statistics" argument is normally the ';'-separated + # X0=...;dX=...; block that also appears in mccode.sim's per-monitor + # header block - stdout's one-line summary doesn't carry that, so + # fall back to Detector's own defaults (X0=0, dX=1, Y0=0, dY=1) by + # passing an empty string. + detectors.append(Detector(name, intensity, error, count, path, '')) + return detectors + + def point_at(N, key, minmax, step): """ Helper to compute the point for key at step """ low, high = map(Decimal, minmax) return step * (high - low) / Decimal(N - 1) + low +def resolve_scan_value(key, value, intervals): + """ Returns a numeric representation of one scanned parameter's value + for one scan point, for writing into mccode.dat's per-point data + row (mccode.dat's format is a matrix of numbers - see module + docstring/build_header() - so every column needs one, regardless + of what kind of value the parameter itself actually is). + + A genuinely numeric value (the overwhelming majority of scans, and + the only kind LinearInterval/MultiInterval.from_range() ever + produce) passes straight through unchanged - this function has no + effect at all outside of an -L/--list scan with a non-numeric + list. + + A non-numeric value (e.g. a -L scan like + filename=Na2Ca3Al2F14.laz,YBaCuO.lau,Fe.laz,Cu.laz) is replaced by + its own *index* within intervals[key] - the position it appears + at in the original -L list, e.g. that list gives indices 0,1,2,3 + respectively - so mccode.dat keeps a properly numeric column for + this parameter too, and remains plottable against it (as a + categorical/index axis) rather than needing the actual string + embedded in a number matrix. + + Each scanned parameter is resolved independently, unlike the + previous behaviour of collapsing the ENTIRE row down to a single + step-index the moment ANY ONE scanned parameter was non-numeric - + which silently discarded every OTHER parameter's real value too + (numeric ones included), and produced only one parameter column + regardless of how many were actually being scanned - a mismatch + against the header's declared xvars/variables count that broke + every downstream plotting tool, since they parse a fixed number + of parameter columns based on that count. """ + try: + return float(value) + except (TypeError, ValueError): + pass + try: + return float(list(intervals[key]).index(value)) + except (KeyError, ValueError): + # value isn't literally in intervals[key] (shouldn't normally + # happen, since scan points are always built FROM intervals[key] - + # but fall back to a stable value rather than crashing outright) + return float(abs(hash(value)) % 1000000) + + class LinearInterval: """ Intervals for linear scanning """ @@ -211,6 +345,13 @@ class MultiInterval: @staticmethod def from_range(N, intervals): + """ N is either a single int (the same point count applied to + every scanned dimension - the original behaviour) or a dict + mapping each interval key to its own point count (mcrun's + -N=a,b,c,... list-form, only valid together with -M, letting + different parameters be sampled at different resolutions - + e.g. a coarse 3-point sweep on one axis against a fine + 20-point sweep on another). """ print(f"MultiInterval from {N=} and {intervals=}") # base case: no intervals yields empty dict if len(intervals) == 0: @@ -219,12 +360,34 @@ def from_range(N, intervals): # recursively generate the multi dict intervals = intervals.copy() key, minmax = intervals.popitem() - for step in range(N): - point = point_at(N, key, minmax, step) + n_here = N[key] if isinstance(N, dict) else N + for step in range(n_here): + point = point_at(n_here, key, minmax, step) for dic in MultiInterval.from_range(N, intervals): dic[key] = point yield dic + @staticmethod + def from_list(intervals): + """ Cartesian product across each key's own explicit list of + points (mcrun's -L/--list combined with -M/--multi). Unlike + LinearInterval.from_list() (co-linear: every key's list is + walked together in lockstep, so all lists must be the same + length), each key here is varied independently, so the lists + may have different lengths - which is also how different + parameters naturally end up with different numbers of scan + points in this mode, without needing a separate -N. """ + print(f"MultiInterval from_list {intervals=}") + if len(intervals) == 0: + yield {} + return + intervals = intervals.copy() + key, values = intervals.popitem() + for value in values: + for dic in MultiInterval.from_list(intervals): + dic[key] = value + yield dic + class InvalidInterval(McRunException): pass @@ -247,11 +410,46 @@ def _simulate_point(args): for key in intervals: mcstas.set_parameter(key, point[key]) - par_values.append(point[key]) + # set_parameter() above needs the real value (a genuine instrument + # filename parameter needs the actual string, not an index) - only + # what goes into the OUTPUT ROW (par_values, eventually written to + # mccode.dat) needs the numeric-or-index resolution. Unlike + # Scanner.run() (only reachable for -L scans), this path is shared + # with plain equidistant multi-dim scans too, but + # resolve_scan_value() is a no-op for those - every value is + # already numeric there. + par_values.append(resolve_scan_value(key, point[key], intervals)) current_dir = f'{mcstas_dir}/{i}' - mcstas.run(pipe=False, extra_opts={'dir': current_dir}) - detectors = mcsimdetectors(current_dir) + is_nexus = mcstas.options.format.lower() == 'nexus' + # See Scanner.run()'s matching NeXus branch: there is no per-step + # mccode.sim to read detector values back from in NeXus mode, so + # capture stdout (pipe=True) and parse its "Detector: ..." summary + # lines directly instead of calling mcsimdetectors(). + try: + stdout_text = mcstas.run(pipe=is_nexus, extra_opts={'dir': current_dir}) + if is_nexus: + detectors = parse_detectors_from_stdout(stdout_text) + else: + detectors = mcsimdetectors(current_dir) + if not detectors: + # No exception, but nothing usable either (e.g. a NeXus step + # whose stdout didn't contain any "Detector: ..." lines at + # all) - treated the same as a runtime failure below: skip + # this point rather than writing an empty/malformed row. + LOG.warning('Scan step %d produced no detector data - skipping this point. Parameters were: %s', + i, ', '.join(f'{k}={v}' for k, v in point.items())) + detectors = None + except Exception as e: + # A single failed scan point (simulation crash, non-zero exit, + # unreadable output, ...) shouldn't take down the whole scan - + # log it and report no detectors for this point; Scanner_split.run() + # already skips any result with detectors=None rather than writing + # a row for it, so the scan carries on to the remaining points and + # mccode.dat simply omits this one. + LOG.warning('Scan step %d failed (%s: %s) - skipping this point and continuing with the rest of the scan. ' + 'Parameters were: %s', i, type(e).__name__, e, ', '.join(f'{k}={v}' for k, v in point.items())) + detectors = None result = { 'index': i, @@ -286,55 +484,132 @@ def run(self): if mcstas_dir == '': mcstas_dir = '.' + points = list(self.points) + header_written = False + skipped = [] + with open(self.outfile, 'w') as outfile: - for i, point in enumerate(self.points): + for i, point in enumerate(points): par_values = [] for key in self.intervals: self.mcstas.set_parameter(key, point[key]) LOG.debug("%s: %s", key, point[key]) par_values.append(point[key]) - if not self.mcstas.options.format.lower() == 'nexus': - LOG.info(', '.join(f'{name}: {value}' for name, value in point.items())) - # Change subdirectory as an extra option (dir/1 -> dir/2) - current_dir = f'{mcstas_dir}/{i}' - LOG.info(f"Output step into scan directory {current_dir}") - self.mcstas.run(pipe=False, extra_opts={'dir': current_dir}) - else: - current_dir = mcstas_dir - LOG.info(f"NeXus output step into scan directory {current_dir}") - self.mcstas.options.append=True - self.mcstas.run(pipe=False, extra_opts={'dir': current_dir}) - - LOG.info("Finish running step, get detectors") - detectors = mcsimdetectors(current_dir) - if detectors is not None: - LOG.info("Got detectors") - if i == 0: - LOG.info("Write headers") - names = [det.name for det in detectors] - outfile.write(build_header(self.mcstas.options, self.intervals.keys(), self.intervals, names)) + try: + if not self.mcstas.options.format.lower() == 'nexus': + LOG.info(', '.join(f'{name}: {value}' for name, value in point.items())) + # Change subdirectory as an extra option (dir/1 -> dir/2) + current_dir = f'{mcstas_dir}/{i}' + LOG.info(f"Output step into scan directory {current_dir}") + self.mcstas.run(pipe=False, extra_opts={'dir': current_dir}) + LOG.info("Finish running step, get detectors") + detectors = mcsimdetectors(current_dir) + else: + current_dir = mcstas_dir + LOG.info(f"NeXus output step into scan directory {current_dir}") + self.mcstas.options.append=True + # NeXus (intentionally) accumulates every step into one + # shared .h5 file rather than a per-step mccode.sim/ + # mccode.dat, so there is no per-step results file to + # read detector values back from at all - + # mcsimdetectors() would just find a .h5 file here and + # return nothing. Capture the simulation's own stdout + # instead (pipe=True) and parse its "Detector: ..." + # summary lines directly (see + # parse_detectors_from_stdout()) - the underlying + # binary always prints that per-run summary regardless + # of output format. + stdout_text = self.mcstas.run(pipe=True, extra_opts={'dir': current_dir}) + if stdout_text: + # pipe=True suppresses the simulation's live + # console output in favour of capturing it for + # parsing - echo it back so nothing is silently + # lost, just delayed until the step completes + # rather than streamed in real time. + print(stdout_text, end='' if stdout_text.endswith('\n') else '\n') + LOG.info("Finish running step, get detectors from stdout") + detectors = parse_detectors_from_stdout(stdout_text) + except Exception as e: + # A single failed scan point (simulation crash, + # non-zero exit, unreadable output, a bad parameter + # combination the instrument itself rejects, ...) + # shouldn't take down the whole scan - log it clearly + # (which parameters were in play) and move on to the + # next point rather than writing anything for this one. + LOG.warning( + 'Scan step %d/%d failed (%s: %s) - skipping this point and continuing with the rest of ' + 'the scan. Parameters were: %s', i + 1, len(points), type(e).__name__, e, + ', '.join(f'{k}={v}' for k, v in point.items())) + skipped.append(i) + continue + + if not detectors: + # No exception, but nothing usable either (e.g. a NeXus + # step whose stdout didn't contain any + # "Detector: ..." lines at all) - skip this point too, + # rather than writing an empty/malformed row that would + # desync the column count from the header. + LOG.warning('Scan step %d/%d produced no detector data - skipping this point. Parameters were: %s', + i + 1, len(points), ', '.join(f'{k}={v}' for k, v in point.items())) + skipped.append(i) + continue + + LOG.info("Got detectors") + if not header_written: + # Written on the first SUCCESSFUL point, not + # unconditionally at index 0 - point 0 might itself be + # the one that failed above. + LOG.info("Write headers") + names = [det.name for det in detectors] + outfile.write(build_header(self.mcstas.options, self.intervals.keys(), self.intervals, names)) + # NeXus format writes every scan step's data into + # its own combined mccode.h5 rather than per-step + # mccode.sim/detector files - a scan-level + # mccode.sim here would describe mccode.dat + # correctly on its own, but would misleadingly + # look like the usual pairing with per-monitor + # detector files that don't actually exist in + # NeXus mode, so it's skipped. mccode.dat itself is + # still written either way, and stays directly + # plottable via its own embedded header alone (see + # mcplotloader.py's load_sweep_dat_only()). + if self.mcstas.options.format.lower() != 'nexus': # Opening a file inside of this loop seems like a bad idea ... oh well with open(self.simfile, 'w') as simfile: simfile.write(build_mccodesim_header(self.mcstas.options, self.intervals, names, version=self.mcstas.version)) - LOG.info("Wrote headers") - LOG.info(f"Write step detectors line into {self.outfile}") - values = ['%s %s' % (d.intensity, d.error) for d in detectors] + LOG.info("Wrote headers") + header_written = True + LOG.info(f"Write step detectors line into {self.outfile}") + values = ['%s %s' % (d.intensity, d.error) for d in detectors] + + if not self.mcstas.options.list: + # Normal equidistant scan: LinearInterval/MultiInterval + # .from_range() only ever produce numeric values, so + # this is unchanged. + line = '%s %s\n' % (' '.join(map(str, par_values)), ' '.join(values)) + else: + # -L list scan: resolve each scanned parameter's + # value independently (see resolve_scan_value()) - + # a genuinely numeric value passes straight + # through, and only a non-numeric one (e.g. a + # filename) becomes its own index within that + # parameter's own list, keeping one proper numeric + # column per scanned parameter either way. + resolved = [resolve_scan_value(key, val, self.intervals) + for key, val in zip(self.intervals.keys(), par_values)] + line = '%s %s\n' % (' '.join(map(str, resolved)), ' '.join(values)) + outfile.write(line) + outfile.flush() - # Normal equidistant scan - if not self.mcstas.options.list: - line = '%s %s\n' % (' '.join(map(str, par_values)), ' '.join(values)) - else: - try: - # Check if parameters are numeric/float - par_floats = [float(x) for x in par_values] - line = '%s %s\n' % (' '.join(map(str, par_floats)), ' '.join(values)) - except: - # otherwise use simple 'index' (may be scanning e.g. a filename) - line = '%s %s\n' % (str(i), ' '.join(values)) - outfile.write(line) - outfile.flush() + if skipped: + LOG.warning('%d of %d scan point(s) failed or produced no data and were skipped ' + '(step indices: %s). %s contains only the %d successful point(s).', + len(skipped), len(points), ', '.join(str(s) for s in skipped), + self.outfile, len(points) - len(skipped)) + else: + LOG.info('Scan complete: all %d point(s) succeeded.', len(points)) class Scanner_split: @@ -379,22 +654,29 @@ def run(self): # Sort results to preserve order results.sort(key=lambda r: r['index']) + skipped = [r['index'] for r in results if not r['detectors']] + with open(self.outfile, 'w') as outfile: wrote_headers = False for result in results: - if result['detectors'] is None: + if not result['detectors']: continue if not wrote_headers: names = [d.name for d in result['detectors']] outfile.write(build_header(self.mcstas.options, self.intervals.keys(), self.intervals, names)) - with open(self.simfile, 'w') as simfile: - simfile.write(build_mccodesim_header( - self.mcstas.options, - self.intervals, - names, - version=self.mcstas.version - )) + # See Scanner.run()'s matching NeXus branch: skip the + # scan-level mccode.sim for NeXus format, for the same + # reason - mccode.dat itself is still written and + # stays plottable on its own. + if self.mcstas.options.format.lower() != 'nexus': + with open(self.simfile, 'w') as simfile: + simfile.write(build_mccodesim_header( + self.mcstas.options, + self.intervals, + names, + version=self.mcstas.version + )) wrote_headers = True values = ['%s %s' % (d.intensity, d.error) for d in result['detectors']] @@ -402,6 +684,14 @@ def run(self): outfile.write(line) outfile.flush() + if skipped: + LOG.warning('%d of %d scan point(s) failed or produced no data and were skipped ' + '(step indices: %s). %s contains only the %d successful point(s).', + len(skipped), len(results), ', '.join(str(s) for s in skipped), + self.outfile, len(results) - len(skipped)) + else: + LOG.info('Scan complete: all %d point(s) succeeded.', len(results)) + class Optimizer: """ Optimize monitors by varying the parameters within interval """