From fd2f82fa8691e5ab3f7a21ffe39280cd0c17f45c Mon Sep 17 00:00:00 2001 From: Remis Date: Thu, 12 Nov 2026 20:00:07 +0000 Subject: [PATCH] Make the .mec reader say what is wrong instead of crashing or hanging A survey of one collection of 29 .mec files found every way mec_load and mec_get_list failed. Each is now handled: file is not a .mec file qHETGLY15.mec is an HJCFIT text printout with a .mec name. mec_get_list raises ValueError quoting how the file begins. index points at no record Bgtvlt30.mec: the bytes where its only record should start give version 0 and over a billion states, and mec_load then read a "drawing" of over a billion cells, one empty read past the end of the file at a time - for ever. The record header is now checked, and so is the drawing's size against what is left of the file. header contradicts itself qmechsDC1.mec @42458: 5 states, but 2+1+3+0 of the four kinds. ValueError, naming both. rates use undeclared states one record shared by qmechs.mec, qmechdem102.mec, qmechdemw.mec and qmechsDC1.mec declares 3 states; its rates connect 6. ValueError, naming the state. no state names, or names CKSIM3B.MEC stores none; independent-subunit for subunits only models store their subunits' state names (8 for 16 states) and their basic rate constants' names (12 for 48 rates). Stored names are now used only when there is one per state or rate; otherwise states are S1, S2, ... and rates q(i,j). Putting subunit names on the first states would label the mechanism wrongly. These records now load. fields read past the record six trailing integers, never used, are not in every record: in demo.mec each record ends at its state names, so reading them took 24 bytes of the next record and ran off the end of the file on the last. They are no longer read. Every failure is now a ValueError naming the file, the record and what is wrong - never an IndexError, an EOFError or a hang. Compared record by record with master across the whole collection: all 352 records that loaded before load identically (states, rate names, cycles, layout, Q to 12 decimals); 8 more now load; 5 are refused as contradictory. Not done, and documented in mec_load: the constraints an independent-subunit model stores are still read but not applied. Read as "rate = factor x rate" they hold in some records and contradict the stored rates by factors up to a million in others, so they are not reliable enough to fill in rates from. qmechsDC1.mec @68054, which stores only its 12 basic rate constants, loads with its other 36 rates at zero, as ten other records already load with zero rates. Also removes the hand-over to dcio.formats.mec. The module docstring said .mec reading had moved there; no such module exists in dcio, so the reader here was always the one used. Had one appeared, every .mec file would have been read by it instead, silently, without any of the above. The rules are tested without files; the tests that load the actual records need SCALCS_MEC_DIR and skip without it. Co-Authored-By: Claude Opus 5.5 --- scalcs/scalcsio.py | 267 +++++++++++++++++++++++++-------------- tests/test_mec_reader.py | 233 ++++++++++++++++++++++++++++++++++ 2 files changed, 402 insertions(+), 98 deletions(-) create mode 100644 tests/test_mec_reader.py diff --git a/scalcs/scalcsio.py b/scalcs/scalcsio.py index 7a851ac..8dfd1b0 100644 --- a/scalcs/scalcsio.py +++ b/scalcs/scalcsio.py @@ -1,17 +1,15 @@ """Reading MEC, SCN, SSD, ABF files and kinetic mechanisms from DCPROGS files. -.. deprecated:: - MEC reading is now provided by :mod:`dcio.formats.mec`. - The functions :func:`mec_get_list` and :func:`mec_load` below are thin - wrappers that delegate to ``dcio`` and convert the result to - :class:`scalcs.mechanism.Mechanism`. Import directly from ``dcio`` for - new code:: - - from dcio.formats.mec import read_list, read - +The .mec reader is the one in this module, :func:`mec_get_list` and +:func:`mec_load`. They used to hand over to ``dcio.formats.mec`` whenever that +could be imported, but no such module exists in dcio, so they never did. The +hand-over was removed rather than left waiting: had such a module appeared, +every .mec file would have been read by it instead, silently, and without the +checks below. """ import json +import os import struct import time from array import array @@ -20,68 +18,112 @@ from scalcs import mechanism # --------------------------------------------------------------------------- -# MEC reading — delegates to dcio.formats.mec, converts to Mechanism +# Checking a .mec file before trusting it +# +# Without these the reader trusted every count it read, so a damaged file or +# a file that is not a .mec file at all failed with an IndexError deep in the +# reader, or - when a count came out in the billions - read empty bytes past +# the end of the file indefinitely. Each check says what is wrong instead. # --------------------------------------------------------------------------- -try: - from dcio.formats.mec import ( - read_list as _dcio_mec_read_list, - read as _dcio_mec_read, - MECEntry as _MECEntry, - ) - _DCIO_AVAILABLE = True -except ImportError: - _DCIO_AVAILABLE = False +# DC's programs give their files and records three-digit version numbers - +# 102 and 200 in every file seen. Anything else is not a .mec header. +_MEC_VERSIONS = range(100, 1000) +# Far more states than any mechanism these programs could hold; a count above +# it comes from bytes that are not a header. +_MAX_STATES = 1000 -def _mec_record_to_mechanism(rec): - """Convert a dcio MECRecord to a scalcs Mechanism.""" - StateList = [] - j = 0 - for _ in range(rec.kA): - gamma = rec.conductances[j] if j < len(rec.conductances) else 0.0 - StateList.append(mechanism.State('A', rec.state_names[j], gamma)) - j += 1 - for _ in range(rec.kB): - StateList.append(mechanism.State('B', rec.state_names[j], 0)) - j += 1 - for _ in range(rec.kC): - StateList.append(mechanism.State('C', rec.state_names[j], 0)) - j += 1 - for _ in range(rec.kD): - StateList.append(mechanism.State('D', rec.state_names[j], 0)) - j += 1 - RateList = [] - for n in range(len(rec.irate)): - bound = None - for m in range(len(rec.cdep_from)): - if rec.cdep_from[m] == rec.irate[n] and rec.cdep_to[m] == rec.jrate[n]: - bound = 'c' - break - name = rec.rate_names[n] if n < len(rec.rate_names) else '' - RateList.append(mechanism.Rate( - rec.rate_values[n], - StateList[rec.irate[n] - 1], - StateList[rec.jrate[n] - 1], - name=name, - eff=bound, - )) +def _looks_like_text(head): + printable = sum(32 <= byte < 127 or byte in (9, 10, 13) for byte in head) + return len(head) > 0 and printable >= 0.95 * len(head) - CycleList = [] - for ci in range(rec.n_cycles): - states_in_cycle = [rec.state_names[idx - 1] - for idx in rec.cycle_state_indices[ci]] - CycleList.append(mechanism.Cycle(states_in_cycle)) - mec = mechanism.Mechanism( - RateList, CycleList, - mtitle=rec.mec_title, - rtitle=rec.rate_title, - ) - # a record that carries the stored drawing passes it straight through - mec.layout = getattr(rec, 'layout', None) - return mec +def _mec_file_problem(data): + """What makes `data`, the bytes of a whole file, unreadable as a .mec + file - or None if nothing does. + + The file begins with its version, the number of records, two words the + reader does not use, and then the byte at which each record starts. + """ + if len(data) < 16: + return 'it is not a .mec file: it is only %d bytes long' % len(data) + version, nrecs = struct.unpack_from('= k: + return list(stored[:k]) + return ['S{0:d}'.format(n) for n in range(1, k + 1)] + + +def _names_for_rates(stored, irate, jrate): + """A name for each rate, given its from- and to-states. + + As with states, the stored names are used only when there is one for + every rate. An independent-subunit model stores names for its basic rate + constants only - npar of them, fewer than its rates - with the rest tied + to those by constraints, and which rates the names belong to is not + recorded plainly enough to rely on. Every rate is then named by its + transition instead: q(13,9) is the rate from state 13 to state 9. + """ + if len(stored) >= len(irate): + return list(stored[:len(irate)]) + return ['q(%d,%d)' % (i, j) for i, j in zip(irate, jrate)] + + +def _rates_problem(irate, jrate, k): + """None if every rate joins two of the `k` declared states.""" + used = list(irate) + list(jrate) + bad = [state for state in used if not 1 <= state <= k] + if bad: + state = max(bad) if max(bad) > k else min(bad) + return ('its rates connect state %d, but the header declares %d ' + 'states' % (state, k)) + return None # --------------------------------------------------------------------------- @@ -288,15 +330,11 @@ def mec_get_list(mecfile): Number of different mechanisms in mec file. """ - if _DCIO_AVAILABLE: - version, entries, max_mecnum = _dcio_mec_read_list(mecfile) - meclist = [ - [e.start_byte, e.mec_num, e.mec_title, e.rate_title] - for e in entries - ] - return version, meclist, max_mecnum + with open(mecfile, 'rb') as whole: + problem = _mec_file_problem(whole.read()) + if problem: + raise ValueError('%s: %s' % (mecfile, problem)) - # --- legacy fallback (dcio not installed) --- f = open(mecfile, 'rb') ints = array('i') @@ -407,18 +445,35 @@ def mec_load(mecfile, start): Returns ------- mec.Mechanism(RateList, StateList, ncyc) : instance of Mechanism class. - """ - if _DCIO_AVAILABLE: - # start may be an MECEntry or a raw int (back-compat) - rec = _dcio_mec_read(mecfile, start) - return _mec_record_to_mechanism(rec) + Raises + ------ + ValueError + If there is no valid record at `start`, if the record contradicts + itself, or if its stored drawing is larger than the file. The message + names the file, the record and what is wrong. + + Notes + ----- + States and rates take their stored names only when the file stores one + for each; otherwise every state is named S1, S2, ... and every rate by + its transition, q(i,j). See _names_for_states and _names_for_rates. + + Rate constants are taken exactly as stored. An independent-subunit model + also stores constraints (ieq, jeq, ifq, jfq, efacq) tying its rates to + its basic rate constants; they are read but not applied. Read as "rate = + factor x rate", they hold in some records but contradict the stored + rates, by factors of up to a million, in others, so their meaning is not + certain enough to fill in rates from. A record that stores only its + basic constants therefore loads with its other rates at zero. + """ - # --- legacy fallback (dcio not installed) --- doubles = array('d') floats = array ('f') ints = array('i') + filesize = os.path.getsize(mecfile) + where = '%s, record at byte %d' % (mecfile, start) f=open(mecfile, 'rb') # open the .mec file as read only f.seek(start - 1); ints.fromfile(f, 1) @@ -439,6 +494,11 @@ def mec_load(mecfile, start): ints.fromfile(f,1) kD = ints.pop() + problem = _record_header_problem(version1, mecnum, k, kA, kB, kC, kD) + if problem: + f.close() + raise ValueError('%s: %s' % (where, problem)) + # In mec files version=103 all shut states are of type 'C'. # Check and leave just one in state 'C', others go as 'B'. if kB == 0: @@ -517,6 +577,16 @@ def mec_load(mecfile, start): # two-character cells, stored column by column. Kept, row by row, so the # positions of the states can be recovered once their names are known; # see _layout_from_board. + # + # Its size is checked against what is left of the file first. f.read() + # past the end returns empty bytes rather than failing, so a size read + # from damaged bytes - over a billion cells, in one file - made this loop + # run, in effect, for ever. + if ilast < 0 or jlast < 0 or f.tell() + 2 * ilast * jlast > filesize: + f.close() + raise ValueError( + '%s: the stored drawing of the mechanism would be %d x %d cells, ' + 'more than is left of the file' % (where, ilast, jlast)) Jch = [] for j in range(0, jlast): Ich = [] @@ -540,6 +610,10 @@ def mec_load(mecfile, start): for i in range(nrateq): ints.fromfile(f,1) jrate.append(ints.pop()) + problem = _rates_problem(irate, jrate, k) + if problem: + f.close() + raise ValueError('%s: %s' % (where, problem)) QT = np.zeros((k, k), 'float64') for i in range(nrateq): doubles.fromfile(f, 1) @@ -673,34 +747,32 @@ def mec_load(mecfile, start): #print "State name:", statename #print ("\n") - ints.fromfile(f,1) - nsub = ints.pop() - ints.fromfile(f,1) - kstat0 = ints.pop() - ints.fromfile(f,1) - npar0 = ints.pop() - ints.fromfile(f,1) - kcon = ints.pop() - ints.fromfile(f,1) - npar1 = ints.pop() - ints.fromfile(f,1) - ncyc0 = ints.pop() + # Some records go on with six more integers - nsub, kstat0, npar0, kcon, + # npar1 and ncyc0 - which nothing here uses. Not every record has them: + # in demo.mec each record ends right after its state names, so reading + # them took 24 bytes of the next record, and ran off the end of the file + # on the last one. They are no longer read. f.close() + # one name per state and per rate; see _names_for_states and + # _names_for_rates for when the stored ones cannot be used + names = _names_for_states(statenames, k) + ratename = _names_for_rates(ratename, irate, jrate) + StateList = [] j = 0 for i in range(kA): - StateList.append(mechanism.State('A', statenames[j], dgamma[j])) + StateList.append(mechanism.State('A', names[j], dgamma[j])) j += 1 for i in range(kB): - StateList.append(mechanism.State('B', statenames[j], 0)) + StateList.append(mechanism.State('B', names[j], 0)) j += 1 for i in range(kC): - StateList.append(mechanism.State('C', statenames[j], 0)) + StateList.append(mechanism.State('C', names[j], 0)) j += 1 for i in range(kD): - StateList.append(mechanism.State('D', statenames[j], 0)) + StateList.append(mechanism.State('D', names[j], 0)) j += 1 RateList = [] @@ -721,14 +793,13 @@ def mec_load(mecfile, start): # mrconstrained = False CycleStates = [] for j in range(nsc[i]): - CycleStates.append(statenames[im[i, j]-1]) + CycleStates.append(names[im[i, j]-1]) CycleList.append(mechanism.Cycle(CycleStates)) mec = mechanism.Mechanism(RateList, CycleList, mtitle=mectitle, rtitle=ratetitle) - # the file can list more names than the mechanism has states (kstat > k); - # only the first k are states, and only they have to be on the board - mec.layout = _layout_from_board(board, statenames[:len(StateList)]) + # names holds exactly one name per state, numbered as the drawing is + mec.layout = _layout_from_board(board, names) return mec def load_mec_from_prt(filename, verbose=False): diff --git a/tests/test_mec_reader.py b/tests/test_mec_reader.py new file mode 100644 index 0000000..024d6a2 --- /dev/null +++ b/tests/test_mec_reader.py @@ -0,0 +1,233 @@ +"""Tests for the .mec reader's handling of damaged and unusual files. + +A survey of one collection of 29 .mec files (374 records) found every way the +reader failed: a file that is not a .mec file at all; an index entry pointing at +no record, which made mec_load hang; records whose header contradicts itself or +whose rates use states the header does not declare; records that store no state +names, or only the names of the subunits of an independent-subunit model; and +six trailing fields, never used, that the reader read past the end of records +that do not have them. + +The rules are tested on their own, with no file. The tests at the end load the +actual records, and need SCALCS_MEC_DIR to point at that collection. +""" +import contextlib +import glob +import io +import os +import struct + +import numpy as np +import pytest + +from scalcs import scalcsio + + +def mec_header(version=102, starts=(20000,), size=22000): + """The first bytes of a .mec file, as mec_get_list reads them.""" + data = struct.pack('<4i', version, len(starts), 22000, starts[-1]) + data += struct.pack('<%di' % len(starts), *starts) + return data + b'\x00' * (size - len(data)) + + +class TestIsItAMecFile: + + def test_a_well_formed_header(self): + assert scalcsio._mec_file_problem(mec_header()) is None + + def test_a_text_file_says_what_it_starts_with(self): + """qHETGLY15.mec is an HJCFIT printout with a .mec extension. Read as + binary, its first words gave version 1344277005 and 1919381362 + records.""" + text = (b'\r\n Program HJCFIT Windows Version\r\n\r\nCopyright ' + b'D. Colquhoun, I. Vais and A.G. Hawkes\r\n' + b' ' * 200) + problem = scalcsio._mec_file_problem(text) + assert 'not a .mec file' in problem + assert 'Program HJCFIT Windows Version' in problem + + def test_a_record_that_would_start_beyond_the_end(self): + problem = scalcsio._mec_file_problem(mec_header(starts=(20000, 90000))) + assert '90000' in problem + + def test_an_impossible_number_of_records(self): + data = struct.pack('<4i', 102, 1919381362, 0, 0) + b'\x00' * 100 + assert scalcsio._mec_file_problem(data) is not None + + +class TestRecordHeader: + + def test_a_well_formed_record(self): + assert scalcsio._record_header_problem(102, 1, 10, 3, 6, 1, 0) is None + + def test_no_record_where_the_index_points(self): + """Bgtvlt30.mec: the index says byte 20000, where the bytes give + version 0, mechanism 0 and over a billion states. The reader went on + to read a drawing of over a billion cells, one empty read at a time.""" + problem = scalcsio._record_header_problem( + 0, 0, 1852793632, 1952671086, 1936617321, 1952801312, 1852138871) + assert 'no valid record' in problem + + def test_a_header_that_contradicts_itself(self): + """qmechsDC1.mec @42458: 5 states, but 2 + 1 + 3 + 0 of the four + kinds.""" + problem = scalcsio._record_header_problem(102, 23, 5, 2, 1, 3, 0) + assert '5 states' in problem and '6' in problem + + +class TestStateNames: + """A name is used only when the file stores one for every state.""" + + def test_one_name_per_state(self): + assert scalcsio._names_for_states(['AR*', 'AR', 'R'], 3) == \ + ['AR*', 'AR', 'R'] + + def test_no_names_stored(self): + """CKSIM3B.MEC stores none: kstat is 0.""" + assert scalcsio._names_for_states([], 3) == ['S1', 'S2', 'S3'] + + def test_fewer_names_than_states_are_not_used_at_all(self): + """ioana.mec, an independent-subunit model, has 16 states and stores + 8 names - the states of its two 4-state subunits, not states 1 to 8. + Putting those names on the first 8 states would label the mechanism + wrongly, so every state gets its positional name instead.""" + subunits = ['A2R*', 'A2R', 'AR', 'R', 'G2R*', 'G2R', 'GR', 'R'] + assert scalcsio._names_for_states(subunits, 16) == \ + ['S%d' % n for n in range(1, 17)] + + def test_more_names_than_states(self): + assert scalcsio._names_for_states(['A', 'B', 'C', 'D'], 3) == \ + ['A', 'B', 'C'] + + +class TestRateNames: + """As with states: a stored name is used only if every rate has one.""" + + def test_one_name_per_rate(self): + assert scalcsio._names_for_rates(['beta1', 'alpha1'], [1, 2], [2, 1]) \ + == ['beta1', 'alpha1'] + + def test_fewer_names_than_rates_name_every_rate_by_its_transition(self): + """An independent-subunit model stores names for its basic rate + constants only - 12 for ioana.mec's 48 rates - and which rates they + belong to is not something to guess.""" + names = scalcsio._names_for_rates(['k+'], [1, 2, 2, 3], [2, 1, 3, 2]) + assert names == ['q(1,2)', 'q(2,1)', 'q(2,3)', 'q(3,2)'] + + +class TestRatesUseDeclaredStates: + + def test_rates_within_the_states(self): + assert scalcsio._rates_problem([1, 2, 2, 3], [2, 1, 3, 2], 3) is None + + def test_a_rate_to_a_state_the_header_does_not_declare(self): + """qmechs.mec @22844 declares 3 states; its rates connect 6.""" + problem = scalcsio._rates_problem([1, 3, 1, 5, 2, 3, 3, 4, 4, 6], + [3, 1, 5, 1, 3, 2, 4, 3, 6, 4], 3) + assert '3 states' in problem and 'state 6' in problem + + +# --------------------------------------------------------------------------- +# The real files +# --------------------------------------------------------------------------- + +MEC_DIR = os.environ.get('SCALCS_MEC_DIR', '') + + +def _path(name): + path = os.path.join(MEC_DIR, name) + if not os.path.isfile(path): + pytest.skip('set SCALCS_MEC_DIR to a folder holding %s' % name) + return path + + +def _load(name, start): + with contextlib.redirect_stdout(io.StringIO()): + return scalcsio.mec_load(_path(name), start) + + +def _records(path): + with contextlib.redirect_stdout(io.StringIO()): + _, entries, _ = scalcsio.mec_get_list(path) + return sorted({entry[0] for entry in entries}) + + +def _q_is_proper(mec): + Q = np.asarray(mec.Q, dtype=float) + return np.abs(Q.sum(axis=1)).max() < 1e-6 * np.abs(Q).max() + + +class TestRealFiles: + + def test_not_a_mec_file(self): + with pytest.raises(ValueError, match='HJCFIT'): + with contextlib.redirect_stdout(io.StringIO()): + scalcsio.mec_get_list(_path('qHETGLY15.mec')) + + def test_an_index_pointing_at_no_record_is_an_error_not_a_hang(self): + with pytest.raises(ValueError, match='20000'): + _load('Bgtvlt30.mec', 20000) + + def test_a_drawing_larger_than_the_file_is_an_error_not_a_hang(self, tmp_path): + """bgtv2a.mec with the drawing's height overwritten with a huge + number: the header is otherwise sound, so only the size check on the + drawing stands between the reader and a hang.""" + data = bytearray(open(_path('bgtv2a.mec'), 'rb').read()) + struct.pack_into('= 360 + assert len(refused) <= 6, refused