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