path: root/cad/src/dna/model/DnaLadder.py

# Copyright 2007-2008 Nanorex, Inc.  See LICENSE file for details. 
"""
DnaLadder.py - automatically recognized guide objects for understanding PAM DNA
duplexes, sticky end regions, and single strands.

@author: Bruce
@version: $Id$
@copyright: 2007-2008 Nanorex, Inc.  See LICENSE file for details.

Context: [warning: some of this info might be obs]

The PAM atoms in any PAM DNA structure will be internally divided into
"ladders", in which all the bonding strictly follows the pattern

 ...  -+-+-+-+-+-+-> ... (strand 1)
       | | | | | | 
 ...  -+-+-+-+-+-+-  ... (axis; missing for free-floating single strand)
       | | | | | | 
 ... <-+-+-+-+-+-+-  ... (strand 2; missing for any kind of single strand)

(with the strand bond directions antiparallel and standardized
to 1 on top (measured left to right == towards increasing baseatoms
index), -1 on bottom).

The dna updater will maintain these ladders, updating them
when their structure is changed. So a changed atom in a ladder will
either dissolve or fragment its ladder (or mark it for the updater to
do that when it next runs), and the updater, dealing with all changed
atoms, will scan all the atoms not in valid ladders to make them into
new ladders (merging new ladders end-to-end with old ones, when that's
valid, and whe they don't become too long to be handled efficiently).

A Dna Segment will therefore consist of a series of ladders (a new one
at every point along it where either strand has a nick or crossover,
or at other points to make it not too long). Each "ladder rail"
(one of the three horizontal atom-bond chains in the figure")
will probably be a separate Chunk, though the entire ladder should
have a single display list for efficiency (so its rung bonds are
in a display list) as soon as that's practical to implement.

So the roles for a ladder include:
- guide the dna updater in making new chunks
- have display code and a display list

A DnaLadder is not visible to copy or undo (i.e. it contains no
undoable state), and is not stored in the mmp file.
"""

from Numeric import dot
from geometry.VQT import cross

from model.elements import Pl5

from utilities.constants import noop
from utilities.constants import MODEL_PAM5

from utilities.Log import orangemsg, redmsg, quote_html

from utilities import debug_flags
from utilities.debug import print_compact_stack

from platform_dependent.PlatformDependent import fix_plurals

import foundation.env as env

from model.bond_constants import atoms_are_bonded

from dna.model.DnaChain import merged_chain

from dna.model.DnaLadderRailChunk import make_or_reuse_DnaAxisChunk
from dna.model.DnaLadderRailChunk import make_or_reuse_DnaStrandChunk
    # note: if these imports are an issue, they could be moved
    # to a controller class, since they are needed only by remake_chunks method

from dna.model.pam_conversion_mmp import DnaLadder_writemmp_mapping_memo

def _DEBUG_LADDER_FINISH_AND_MERGE():
    return debug_flags.DEBUG_DNA_UPDATER_VERBOSE # for now

def _DEBUG_REVERSE_STRANDS():
    # same as _DEBUG_LADDER_FINISH_AND_MERGE, for now
    return debug_flags.DEBUG_DNA_UPDATER_VERBOSE

from dna.model.dna_model_constants import LADDER_ENDS
from dna.model.dna_model_constants import LADDER_END0
from dna.model.dna_model_constants import LADDER_END1
from dna.model.dna_model_constants import LADDER_BOND_DIRECTION_TO_OTHER_AT_END_OF_STRAND1
from dna.model.dna_model_constants import LADDER_STRAND1_BOND_DIRECTION

from dna.model.dna_model_constants import MAX_LADDER_LENGTH

from dna.model.DnaLadder_pam_conversion import DnaLadder_pam_conversion_methods

from dna.model.PAM_atom_rules import PAM_atoms_allowed_in_same_ladder

from dna.updater.dna_updater_prefs import pref_per_ladder_colors
from dna.updater.dna_updater_prefs import pref_permit_bare_axis_atoms

from dna.updater.dna_updater_globals import _f_ladders_with_up_to_date_baseframes_at_ends
from dna.updater.dna_updater_globals import _f_atom_to_ladder_location_dict
from dna.updater.dna_updater_globals import _f_invalid_dna_ladders

from dna.updater.dna_updater_globals import get_dnaladder_inval_policy
from dna.updater.dna_updater_globals import DNALADDER_INVAL_IS_OK
from dna.updater.dna_updater_globals import DNALADDER_INVAL_IS_ERROR
from dna.updater.dna_updater_globals import DNALADDER_INVAL_IS_NOOP_BUT_OK

from dna.updater.dna_updater_globals import rail_end_atom_to_ladder

# ==

### REVIEW: should a DnaLadder contain any undoable state?
# (guess: yes... maybe it'll be a Group subclass, for sake of copy & undo?
#  for now, assume just an object (not a Node), and try to fit it into
#  the copy & undo code in some other way... or I might decide it's entirely
#  implicit re them.) If that gets hard, make it a Group. (Where in the internal MT?
#  whereever the chunks would have been, without it.)

class DnaLadder(object, DnaLadder_pam_conversion_methods):
    """
    [see module docstring]
        
    @note: a valid ladder is not always in
    correspondence with atom chunks for its rails
    (since the chunks are only made or fixed after ladder merging),
    so the methods in self can't find an atom's ladder via atom.molecule.
    Instead we assume that no atom is in more than one valid ladder
    at a time, and set a private atom attribute to point to that ladder.
    """
    #TODO: split this class into a common superclass with DnaSingleStrandDomain
    # (see comment there) and a specific subclass for when we have axis_rail.
    valid = False # public for read, private for set; whether structure is ok and we own our atoms
    error = "" # ditto; whether num_strands or strand bond directions are wrong (parallel), etc # todo: use more?
        # on error, should be set to a short string (suitable for tooltip)
        # without instance-specific data (i.e. one of a small fixed set of
        # possible strings, so also suitable as part of a summary_format string)

    _class_for_writemmp_mapping_memo = DnaLadder_writemmp_mapping_memo
        # subclass can override if needed (presumably with a subclass of this value)
        ### REVIEW: what about for single strands? (our subclass) @@@@@@

    def __init__(self, axis_rail):
        self.axis_rail = axis_rail
        self.assy = axis_rail.baseatoms[0].molecule.assy #k
        self.strand_rails = []
        assert self.assy is not None, "%r.__init__: assy is None" % self

    def _f_make_writemmp_mapping_memo(self, mapping):
        """
        [friend method for class writemmp_mapping.get_memo_for(self)]
        """
        # note: same code in some other classes that implement this method
        return self._class_for_writemmp_mapping_memo(mapping, self)

    def baselength(self):
        return len(self.axis_rail)
    
    def __len__(self):
        return self.baselength()
    
    def __nonzero__(self): # 080311
        # avoid Python calling __len__ for this [review: need __eq__ as well?]
        return True
    
    def add_strand_rail(self, strand_rail): # review: _f_, since requires a lot of calling code?
        """
        This is called while constructing a dna ladder (self)
        which already has an axis rail and 0 or 1 strand rails.
        
        Due to how the calling code works, it's guaranteed
        that the atomlist in the new strand rail corresponds to that
        in the existing axis rail, either directly or in reverse order --
        even if one or both of those rails is a ring rather than a chain.
        (That is, it's never necessary to "rotate a ring" to make this
        alignment possible. The caller splits rings into chains if
        necessary to avoid this need.)
        We depend on this and assert it, in self.finished().
        """
        assert strand_rail.baselength() == self.axis_rail.baselength(), \
               "baselengths in %r don't match: %r (%d) != %r (%d)" % \
               (self,
                strand_rail,
                strand_rail.baselength(),
                self.axis_rail,
                self.axis_rail.baselength())
        self.strand_rails.append(strand_rail)
        if _DEBUG_REVERSE_STRANDS():
            strand_rail.debug_check_bond_direction("in %r.add_strand_rail" % self)
        return
    
    def finished(self): # Q. rename to 'finish'?
        """
        This is called once to signify that construction of self is done
        as far as the caller is concerned (i.e. it's called add_strand_rail
        all the times it's going to), and it should be finished internally.
        This involves reversing rails as needed to make their baseatoms
        correspond (so each base pair's 3 PAM atoms (or each lone base's 2 PAM atoms)
        has one atom per rail, at the same relative baseindex),
        and perhaps(?) storing pointers at each rail-end to neighboring rails.
        See add_strand_rail docstring for why optional reversing is sufficient,
        even when some rails are rings.
        """
        # note: this method is only for the subclass with axis_rail present;
        # at the end it calls a common method.
        assert not self.valid # not required, just a useful check on the current caller algorithm
        ## assert len(self.strand_rails) in (1, 2)
        # happens in mmkit - leave it as just a print at least until we implem "delete bare atoms" -
        from dna.updater.dna_updater_prefs import legal_numbers_of_strand_neighbors_on_axis
        legal_nums = legal_numbers_of_strand_neighbors_on_axis()
        if not ( len(self.strand_rails) in legal_nums ):
            print "error: axis atom %r has %d strand_neighbors (should be %s)" % \
                  (self, len(self.strand_rails), " or ".join(map(str, legal_nums)))
            self.error = "bug: illegal number (%d) of strand rails" % len(self.strand_rails)
        axis_rail = self.axis_rail
        # make sure rungs are aligned between the strand and axis rails
        # (note: this is unrelated to their index_direction,
        #  and is not directly related to strand bond direction)
        # (note: it's trivially automatically true if our length is 1;
        #  the following alg does nothing then, which we assert)
        axis_left_end_baseatom = axis_rail.end_baseatoms()[LADDER_END0]
        for strand_rail in self.strand_rails:
            if strand_rail.end_baseatoms()[LADDER_END0].axis_neighbor() is not axis_left_end_baseatom:
                # we need to reverse that strand
                # (note: we might re-reverse it below, if bond direction wrong)
                assert strand_rail.end_baseatoms()[LADDER_END1].axis_neighbor() is axis_left_end_baseatom
                assert self.baselength() > 1 # shouldn't happen otherwise
                strand_rail.reverse_baseatoms()
                assert strand_rail.end_baseatoms()[LADDER_END0].axis_neighbor() is axis_left_end_baseatom
            continue
        del axis_left_end_baseatom # would be wrong after the following code
        self._finish_strand_rails()
        return
    
    def _finish_strand_rails(self): # also used in DnaSingleStrandDomain
        """
        verify strand bond directions are antiparallel, and standardize them;
        then self._ladder_set_valid()
        (note there might be only the first strand_rail;
         this code works even if there are none)
        """
        if _DEBUG_REVERSE_STRANDS():
            print "\n *** _DEBUG_REVERSE_STRANDS start for %r" % self
        for strand_rail in self.strand_rails:
            if _DEBUG_REVERSE_STRANDS():
                print "deciding whether to reverse:", strand_rail
            if strand_rail is self.strand_rails[0]:
                desired_dir = 1
                reverse = True # if dir is wrong, reverse all three rails
            else:
                desired_dir = -1
                reverse = False # if dir is wrong, error
                    # review: how to handle it in later steps? mark ladder error, don't merge it?
            have_dir = strand_rail.bond_direction() # 1 = right, -1 = left, 0 = inconsistent or unknown
                # @@@ IMPLEM note - this is implemented except for merged ladders; some bugs for length==1 chains.
                # strand_rail.bond_direction must check consistency of bond
                # directions not only throughout the rail, but just after the
                # ends (thru Pl too), so we don't need to recheck it for the
                # joining bonds as we merge ladders. BUT as an optim, it ought
                # to cache the result and use it when merging the rails,
                # otherwise we'll keep rescanning rails as we merge them. #e
            if have_dir == 0:
                # this should never happen now that fix_bond_directions is active
                # (since it should mark involved atoms as errors using _dna_updater__error
                #  and thereby keep them out of ladders), so always print it [080201]
                msg = "bug: %r strand %r has unknown or inconsistent bond " \
                      "direction - should have been caught by fix_bond_directions" % \
                      (self, strand_rail)
                print "\n*** " + msg
                ## env.history.redmsg(quote_html(msg))
                self.error = "bug: strand with unknown or inconsistent bond direction"
                reverse = True # might as well fix the other strand, if we didn't get to it yet
            else:
                if have_dir != desired_dir:
                    # need to reverse all rails in this ladder (if we're first strand,
                    #  or if dir is arb since len is 1) 
                    # or report error (if we're second strand)
                    if strand_rail.bond_direction_is_arbitrary():
                        if _DEBUG_REVERSE_STRANDS():
                            print "reversing %r (bond_direction_is_arbitrary)" % strand_rail
                        strand_rail._f_reverse_arbitrary_bond_direction()
                        have_dir = strand_rail.bond_direction() # only needed for assert
                        assert have_dir == desired_dir
                    elif reverse:
                        if _DEBUG_REVERSE_STRANDS():
                            print "reversing all rails in %r:" % (self,)
                        for rail in self.all_rails(): # works for ladder or single strand domain
                            if _DEBUG_REVERSE_STRANDS():
                                print "   including: %r" % (rail,)
                            ### review: this reverses even the ones this for loop didn't get to, is that desired? @@@
                            rail.reverse_baseatoms()
                    else:
                        # this can happen for bad input (not a bug);
                        # don't print details unless verbose,
                        # use a summary message instead
                        if debug_flags.DEBUG_DNA_UPDATER_VERBOSE:
                            print "\ndna updater input data error: %r's strands have parallel bond directions." % self
                            print "data about that:"
                            print "this rail %r in list of rails %r" % (strand_rail, self.strand_rails)
                            print "desired_dir = %r, reverse = %r, have_dir = %r" % (desired_dir, reverse, have_dir)
                            print "self.error now = %r (is about to be set)" % (self.error,)
                        if strand_rail.bond_direction_is_arbitrary():
                            print " *** bug: parallel strands but strand_rail.bond_direction_is_arbitrary() is true, should be false"
                        self.error = "parallel strand bond directions"
                        # should we just reverse them? no, unless we use minor/major groove to decide which is right.
                    pass
                else:
                    # strand bond direction is already correct
                    if _DEBUG_REVERSE_STRANDS():
                        print "no need to reverse %r" % strand_rail
            continue
        if _DEBUG_REVERSE_STRANDS():
            print "done reversing if needed; current rail baseatom order is:"
            for rail in self.all_rails():
                print "  %r" % rail.baseatoms # see if it gets it right after the first step
            for rail in self.strand_rails:
                strand_rail.debug_check_bond_direction("end of %r._finish_strand_rails" % self)
            print " *** _DEBUG_REVERSE_STRANDS end for %r\n" % self
        self._ladder_set_valid(True) # desired even if self.error was set above, or will be set below
        self._duplex_geometric_checks() # might set or append to self.error
        if self.error:
            # We want to display the error, but not necessarily by setting
            # atom._dna_updater__error for every atom (or at least not right here,
            # since we won't know all the atoms until we remake_chunks).
            # For now, just report it: [#e todo: fix inconsistent case in dna updater warnings @@@@]
            # Later display it in remake_chunks or in some ladder-specific display code. (And in tooltip too, somehow.)
            # Ideally we draw an unmissable graphic which points out the error.
            # (And same for the per-atom _dna_updater__error, which might be too missable
            #  as mere pink wireframe.) @@@@
            error_string = self.error
            if error_string.startswith("bug:"):
                prefix = "Bug"
                colorfunc = redmsg
                # remove "bug: " from error_string used in summary_format
                error_string = error_string[len("bug:"):].strip()
            else:
                prefix = "Warning"
                colorfunc = orangemsg
            summary_format = "%s: DNA updater: found [N] dna ladder(s) with %s" % \
                             ( prefix, error_string )
            env.history.deferred_summary_message( colorfunc( summary_format))
                # note: error prevents ladder merge, so no issue of reporting an error more than once
                # (as there would be otherwise, since same error would show up in the merged ladder)
        return

    def set_error(self, error_string): #bruce 080401 REVIEW: use this above?
        """
        """
        self.error = error_string # REVIEW: more? relation to valid?
        return
    
    def _duplex_geometric_checks(self):
        """
        Set (or append to) self.error, if self has any duplex geometric errors.
        """
        if len(self.strand_rails) == 2 and len(self) > 1:
            # use the middle two base pairs; check both strands
            i1 = (len(self) / 2) - 1
            i2 = i1 + 1
            axis_atoms = self.axis_rail.baseatoms[i1:i1+2] # note: for PAM5 these might be Gv, so formulas may differ... @@@
            strand1_atoms = self.strand_rails[0].baseatoms[i1:i1+2]
            strand2_atoms = self.strand_rails[1].baseatoms[i1:i1+2]
            error_from_this = self._check_geom( axis_atoms, strand1_atoms, strand2_atoms)
            if not error_from_this:
                axis_atoms.reverse()
                strand1_atoms.reverse()
                strand2_atoms.reverse()
                self._check_geom( axis_atoms, strand2_atoms, strand1_atoms)
        return
    
    def _check_geom(self, axis_atoms, strand1_atoms, strand2_atoms ):
        """
        Check the first 4 atoms (out of 6 passed to us, as 3 lists of 2)
        for geometric duplex errors.
        If strand2 is being passed in place of strand1 (and thus vice versa),
        all 3 atom lists should be reversed.
        If we set or add to self.error, return something true
        (which will be the error string we add, not including whatever
         separator we use to add it to an existing self.error).
        """
        error_here = "" # sometimes modified below
        if axis_atoms[0].element.pam == MODEL_PAM5:
            # These tests don't yet work properly for PAM5.
            # This should be fixed -- but the error responses
            # also do more harm than good in some cases,
            # so those need fixing too. In the meantime it's
            # best to disable it for PAM5. [bruce 080516]
            return error_here
        old_self_error = self.error # for assert only
        axis_posns = [atom.posn() for atom in axis_atoms]
        strand1_posns = [atom.posn() for atom in strand1_atoms]
        rung0_vec = strand1_posns[0] - axis_posns[0]
        rung1_vec = strand1_posns[1] - axis_posns[1]
        axis_vec = axis_posns[1] - axis_posns[0]
        righthandedness = dot( cross(rung0_vec, rung1_vec), axis_vec )
        # todo: make sure it's roughly the right magnitude, not merely positive
        if not self.error or not self.error.startswith("bug:"):
            # check strand1 for reverse twist (append to error string if we already have one)
            if righthandedness <= 0:
                error_here = "reverse twist"
                if not self.error:
                    self.error = error_here
                else:
                    self.error += "\n" + error_here # is '\n' what we want?
        if not self.error:
            assert not error_here
            # check for wrong strand directions re major groove
            # (since no self.error, we can assume strands are antiparallel and correctly aligned,
            #  i.e. that bond directions point from strand1_atoms[0] to strand1_atoms[1])
            strand2_posns = [atom.posn() for atom in strand2_atoms]
            rung0_vec_2 = strand2_posns[0] - axis_posns[0]
            minor_groove_where_expected_ness = dot( cross( rung0_vec, rung0_vec_2 ), axis_vec )
            if minor_groove_where_expected_ness <= 0:
                error_here = "bond directions both wrong re major groove"
                if 0: # turn this on to know why it reported that error [080210]
                    print "data for %s:" % error_here
                    print axis_atoms, strand1_atoms, strand2_atoms
                    print rung0_vec, rung0_vec_2, axis_vec
                self.error = error_here
        if old_self_error or error_here:
            assert self.error
        assert self.error.startswith(old_self_error)
        assert self.error.endswith(error_here)
        return error_here

    def num_strands(self):
        return len(self.strand_rails)
    
    def _ladder_set_valid(self, val):
        if val != self.valid:
            self.valid = val
            if val:
                # tell the rail end atoms their ladder;
                # see rail_end_atom_to_ladder for how this hint is interpreted
                for atom in self.rail_end_baseatoms():
                    # (could debug-warn if already set to a valid ladder)
                    if atom._DnaLadder__ladder is not None and atom._DnaLadder__ladder.valid:
                        print "\n*** likely bug: %r is owned by %r as %r takes it over" % \
                              (atom, atom._DnaLadder__ladder, self)
                    atom._DnaLadder__ladder = self
                    if debug_flags.DEBUG_DNA_UPDATER_VERBOSE:
                        print "%r owning %r" % (self, atom)
            else:
                # un-tell them, but first, warn if this might cause or indicate
                # a bug in pam conversion
                ladders_dict = _f_ladders_with_up_to_date_baseframes_at_ends
                if self in ladders_dict:
                    msg = "\n*** likely bug: %r is in ladders_dict (value %r) but is becoming invalid" % \
                          (self, ladders_dict[self])
                    print_compact_stack(msg + ": ")
                if self.strand_rails:
                    locator = _f_atom_to_ladder_location_dict
                    sample_atom = self.strand_rails[0].baseatoms[0]
                    if sample_atom.key in locator:
                        msg = "\n*** likely bug: %r's atom %r is in locator (value %r) but self is becoming invalid" % \
                          (self, sample_atom, locator[sample_atom.key])
                        print_compact_stack(msg + ": ")
                    pass
                # now un-tell them (even if we warned)
                for atom in self.rail_end_baseatoms():
                    # (could debug-warn if not set to self)
                    assert atom._DnaLadder__ladder is self # 080411
                    # don't do this:
                    ## atom._DnaLadder__ladder = None
                    ## del atom._DnaLadder__ladder
                    # since it makes it harder to debug things;
                    # uses of atom._DnaLadder__ladder check self.valid anyway
                    if debug_flags.DEBUG_DNA_UPDATER_VERBOSE:
                        print "%r de-owning %r" % (self, atom)
                # tell the next run of the dna updater we're invalid
                _f_invalid_dna_ladders[id(self)] = self
        return
    
    def rail_end_baseatoms(self):
        # rename? end_baseatoms is logical, but occurs too much in other code.
        # ladder_end_baseatoms is good too, but is too confusable
        # with ladder_end_atoms (which also needs renaming).
        # end_atoms occurs only as a localvar, but is misleading
        # re Pl or X at ends.
        """
        yield the 3 or 6 atoms which are end-atoms for one of our 3 rails,
        in arbitrary order #k
        """
        for rail in self.all_rails():
            # note: rail is a DnaChain, not a DnaLadderRailChunk!
            atom1, atom2 = rail.end_baseatoms()
            yield atom1
            if atom2 is not atom1:
                yield atom2
        return

    # external invalidation methods, which check dnaladder_inval_policy
    # [bruce 080413 separated existing invalidate method calls into calls
    #  of these two variants, also giving them text-searchable methodnames]

        # note: one of these is called from dna updater and from
        # DnaLadderRailChunk.delatom, so changed atoms inval their
        # entire ladders
    
    def ladder_invalidate_and_assert_permitted(self): #bruce 080413
        dnaladder_inval_policy = get_dnaladder_inval_policy()
        if dnaladder_inval_policy != DNALADDER_INVAL_IS_OK:
            print "\n*** BUG: dnaladder_inval_policy is %r, but we'll inval %r anyway" % \
                  (dnaladder_inval_policy, self)
        self._ladder_set_valid(False)
        return

    def ladder_invalidate_if_not_disabled(self): #bruce 080413
        dnaladder_inval_policy = get_dnaladder_inval_policy()
        if dnaladder_inval_policy == DNALADDER_INVAL_IS_OK:
            # usual case should be fastest
            self._ladder_set_valid(False)            
        elif dnaladder_inval_policy == DNALADDER_INVAL_IS_NOOP_BUT_OK:
            print "fyi: disabled inval of %r" % self ### remove when works, or put under debug_flags
            pass # don't invalidate it
        elif dnaladder_inval_policy == DNALADDER_INVAL_IS_ERROR:
            # this will happen for true errors, but also for places where we
            # need to disable but permit the inval. To show those clearly,
            # also useful for true errors, we always print the stack, even
            # though this will be verbose when there are bugs.
            msg = "\n*** BUG: dnaladder inval is an error now, but we'll inval %r anyway" % \
                  (self, )
            print_compact_stack( msg + ": " )
            self._ladder_set_valid(False) # an error, but fewer bugs predicted from doing it than not
        else:
            assert 0, "unrecognized dnaladder_inval_policy %r" % dnaladder_inval_policy
        return

    # ==
    
    def get_ladder_end(self, endBaseAtom): # by Ninad # todo: rename to indicate it gets it from an atom, not an end-index
        """
        Return the end (as a value in LADDER_ENDS) of the ladder self
        which has the given end base atom, or None if the given atom
        is not one of self's end_baseatoms. (If self has length 1, either
        end might be returned.)
        
        @param endBaseAtom: One of the end atoms of the ladder (or any
                            other value, for which we will return None,
                            though passing a non-Atom is an error).
        """
        for ladderEnd in LADDER_ENDS:
            if endBaseAtom in self.ladder_end_atoms(ladderEnd):
                return ladderEnd
        
        return None
    
    def get_endBaseAtoms_containing_atom(self, baseAtom): # by Ninad
        """
        Returns a list of end base atoms that contain <baseAtom> (including that
        atom.

        @see: other methods that return these atoms, sometimes in reverse order
              depending on which end they are from, e.g. ladder_end_atoms.
        """
        endBaseAtomList = []
        sortedEndBaseAtomList = []
        end_strand_base_atoms = ()
        for ladderEnd in LADDER_ENDS:
            if baseAtom in self.ladder_end_atoms(ladderEnd):
                endBaseAtomList = self.ladder_end_atoms(ladderEnd)
                break
        if len(endBaseAtomList) > 2:
            end_strand_base_atoms = (endBaseAtomList[0], endBaseAtomList[2])
        else:
            #@TODO: single stranded dna. Need to figure out if it has 
            #strand1 or strand2 as its strand.
            # [always strand1, it guarantees this. --bruce]
            pass
        
        strand1Atom = None
        axisAtom =  endBaseAtomList [1]
        strand2Atom = None
        for atom in end_strand_base_atoms:     
            if atom is not None:
                temp_strand = atom.molecule        
                rail = temp_strand.get_ladder_rail()
                #rail goes from left to right. (increasing chain index order)                
                
                if rail.bond_direction() == 1:
                    strand1Atom = atom
                elif rail.bond_direction() == -1:
                    strand2Atom = atom
     
        endBaseAtomList = [strand1Atom, axisAtom, strand2Atom]
        
        return endBaseAtomList                
     
    def is_ladder_end_baseatom(self, baseAtom):
        """
        Returns true if the given atom is an end base atom of a ladder (either
        at ladder end0 or ladder end1)
        """
        for ladderEnd in LADDER_ENDS:
            if endBaseAtom in self.ladder_end_atoms(ladderEnd):
                return True
        
        return False

    def whichrail_and_index_of_baseatom(self, baseatom): # bruce 080402
        """
        Given one of self's baseatoms, possibly on any rail at any baseindex,
        determine where it is among all self's baseatoms,
        and return this in the form described below. Error if baseatom
        is not found. basetom.element.role may be used to narrow down
        which rails of self to examine.

        @warning: slow for non-end atoms of very long ladders, due to a
                  linear search for self within the ladder. Could be optimized
                  by passing an index hint if necessary. Current code is
                  optimized for baseatom being an end-atom, or a strand atom.

        @return: (whichrail, index), where whichrail is in LADDER_RAIL_INDICES
                 [i.e. 0, 1, or 2, maybe not yet formally defined]
                 and index is from 0 to len(self) - 1, such that
                 self.rail_at_index(whichrail).baseatoms[index] is baseatom.

        @see: related method, _f_store_locator_data
        """
        if not self.valid: #bruce 080411
            # maybe make this an assertion failure?
            print "likely bug: invalid ladder in whichrail_and_index_of_baseatom(%r, %r)" % \
                  (self, baseatom)
        # TODO: pass index hint to optimize?
        look_at_rails = self.rail_indices_and_rails(baseatom)
        for index in range(-1, len(self) - 1):
            # ends first (indices -1 and (if long enough) 0)
            for whichrail, rail in look_at_rails:
                if baseatom is rail.baseatoms[index]:
                    # found it
                    if index < 0:
                        index += len(self)
                    assert self.rail_at_index(whichrail).baseatoms[index] is baseatom # slow, remove when works ###
                    return whichrail, index
        assert 0, "can't find %r in %r, looking at %r" % \
                  (baseatom, self, look_at_rails)
        pass

    def rail_indices_and_rails(self, baseatom = None): # bruce 080402
        """
        #doc

        @see: LADDER_RAIL_INDICES
        """
        del baseatom # todo: optim: use .element.role to filter rails
        # note: to optim some callers, return strand rails first
        res = []
        for i in range(len(self.strand_rails)):
            # assert i < 2
            res.append( (i * 2, self.strand_rails[i]) )
        if self.axis_rail is not None:
            res.append( (1, self.axis_rail) )
        return res

    def rail_at_index(self, whichrail): # bruce 080402
        """
        #doc

        @note: call only on a value of whichrail that is known to
               correspond to a rail present in self,
               e.g. a value which was returned from a method in self,
               e.g. rail_indices_and_rails or whichrail_and_index_of_baseatom.

        @see: LADDER_RAIL_INDICES
        """
        if whichrail == 1:
            return self.axis_rail
        return self.strand_rails[whichrail / 2]

    # == ladder-merge methods
    
    def can_merge(self):
        """
        Is self valid, and mergeable (end-end) with another valid ladder?

        @return: If no merge is possible, return None; otherwise, for one
                 possible merge (chosen arbitrarily if more than one is
                 possible), return the tuple (other_ladder, merge_info)
                 where merge_info is private info to describe the merge.
        """
        if not self.valid or self.error:
            return None # not an error
        for end in LADDER_ENDS:
            other_ladder_and_merge_info = self._can_merge_at_end(end) # might be None
            if other_ladder_and_merge_info:
                return other_ladder_and_merge_info
        return None
    
    def do_merge(self, other_ladder_and_merge_info):
        """
        Caller promises that other_ladder_and_merge_info was returned by
        self.can_merge() (and is not None).
        Do the specified permitted merge (with another ladder, end to end,
        both valid).
        
        Invalidate the merged ladders; return the new valid ladder
        resulting from the merge.
        """
        other_ladder, merge_info = other_ladder_and_merge_info
        end, other_end = merge_info
        return self._do_merge_with_other_at_ends(other_ladder, end, other_end)
    
    def _can_merge_at_end(self, end): # TODO: update & clean up docstring
        """
        Is the same valid other ladder (with no error) attached to each rail of self
        at the given end (0 or 1), and if so, can we merge to it
        at that end (based on which ends of which of its rails
        our rails attach to)? (If so, return other_ladder_and_merge_info, otherwise None.)

        ### REVIEW (what's done now, what;s best): Note that self or the other ladder might
        be length==1 (in which case only end 0 is mergeable, as an
        arbitrary decision to make only one case mergeable), ### OR we might use bond_dir, then use both cases again
        and that the other ladder might only merge when flipped.
        
        Also note that bond directions (on strands) need to match. ### really? are they always set?
        
        TODO: Also worry about single-strand regions. [guess: just try them out here, axis ends are None, must match]
        """
        # algorithm:
        # use strand1 bond direction to find only possible other ladder
        # (since that works even if we're of baselength 1);
        # then if other ladder qualifies,
        # check both its orientations to see if all necessary bonds
        # needed for a merge are present. (Only one orientation can make sense,
        # given the bond we used to find it, but checking both is easier than
        # checking which one fits with that bond, especially if other or self
        # len is 1. But to avoid confusion when other len is 1, make sure only
        # one orientation of other can match, by using the convention
        # that the ladder-end-atoms go top to bottom on the right,
        # but bottom-to-top on the left.)

        if not self.strand_rails:
            # simplest just to never merge bare axis, for now,
            # especially since len 1 case would be nontrivial [bruce 080407]
            return False

        strand1 = self.strand_rails[0]
        assert 1 == LADDER_STRAND1_BOND_DIRECTION # since it's probably also hardcoded implicitly
        assert strand1.bond_direction() == LADDER_STRAND1_BOND_DIRECTION, \
               "wrong bond direction %r in strand1 of %r" % \
               (strand1.bond_direction, self)
        end_atom = strand1.end_baseatoms()[end]
        assert not end_atom._dna_updater__error # otherwise we should not get this far with it
        assert self is rail_end_atom_to_ladder(end_atom) # sanity check
        bond_direction_to_other = LADDER_BOND_DIRECTION_TO_OTHER_AT_END_OF_STRAND1[end]
        next_atom = end_atom.strand_next_baseatom(bond_direction = bond_direction_to_other)
            # (note: strand_next_baseatom returns None if end_atom or the atom it
            #  might return has ._dna_updater__error set, or if it reaches a non-Ss atom.)
        if next_atom is None:
            # end of the chain (since bondpoints are not baseatoms), or
            # inconsistent bond directions at or near end_atom
            # (report error in that case??)
            return None
        assert not next_atom._dna_updater__error
        other = rail_end_atom_to_ladder(next_atom)
            # other ladder (might be self in case of ring)
        if other.error:
            return None
        # other.valid was checked in rail_end_atom_to_ladder
        if other is self:
            return None
        if len(self) + len(other) > MAX_LADDER_LENGTH:
            return None
        # Now a merge is possible and allowed, if all ladder-ends are bonded
        # in either orientation of other, AND if the ladders are compatible
        # for merging in attributes which we don't wish to mix (e.g. chunk.display_as_pam).
        # In checking those attrs, note that the atoms are still in their old chunks.

        if not self.compatible_for_merge_with(other):
            return None

        # optim: reject early if ladders have different numbers of rails --
        # enable this optim once i verify it works without it: @@
        ## if self.num_strands() != other.num_strands():
        ##    return None # optimization
        ## #e could also check axis present or not, same in both
        
        # try each orientation for other ladder;
        # first collect the atoms to test for bonding to other
        self_end_atoms = self.ladder_end_atoms(end)
        if debug_flags.DEBUG_DNA_UPDATER:
            assert end_atom in self_end_atoms
        for other_end in LADDER_ENDS:
            other_end_atoms = other.ladder_end_atoms(other_end, reverse = True)
                # note: .ladder_end_atoms returns a sequence of length 3,
                # with each element being a rail's end-atom or None for a missing rail;
                # the rail order is strand-axis-strand, which on the right is in
                # top to bottom order, but on the left is in bottom to top order.
                # strand1 is always present; strand2 is missing for either kind of
                # single strands; axis_rail is missing for free-floating single strands.
                #
                # bugfix 080122 (untested): reverse order for other_end to facilitate
                # comparison if the ladders abut (with standardized strand directions).
                # Not doing this presumably prevented all merges except those involving
                # length 1 ladders, and made those happen incorrectly when len(other) > 1.
            still_ok = True
            for atom1, atom2, strandQ in zip(self_end_atoms, other_end_atoms, (True, False, True)):
                ok = _end_to_end_bonded(atom1, atom2, strandQ)
                if not ok:
                    still_ok = False
                    break
            if still_ok:
                # success
                if debug_flags.DEBUG_DNA_UPDATER:
                    # end_atom and next_atom are bonded (perhaps via Pl),
                    # so if we succeeded then they ought to be in the same
                    # positions in these lists:
                    assert next_atom in other_end_atoms
                    assert self_end_atoms.index(end_atom) == other_end_atoms.index(next_atom)
                merge_info = (end, other_end)
                other_ladder_and_merge_info = other, merge_info
                return other_ladder_and_merge_info
        return None

    def compatible_for_merge_with(self, other): #bruce 080401
        our_atom = self.arbitrary_baseatom()
        other_atom = other.arbitrary_baseatom()
        return PAM_atoms_allowed_in_same_ladder( our_atom, other_atom)
        # note: above doesn't reject Pl5 bonded to Ss3...
        # review whether anything else will catch that.
        # Guess: yes, the ladder-forming code will catch it
        # (not sure what it would do with Ss3-Pl-Ss3, though -- REVIEW).
        # But note that Ss3-Pl-Ss5 is allowed, at the bond between a PAM3 ad PAM5 chunk.
        # (In which case the Pl should prefer to stay in the same chunk
        #  as the Ss5 atom, not the Ss3 atom -- this is now implemented in Pl_preferred_Ss_neighbor.)

    def arbitrary_baseatom(self): #bruce 080401, revised 080407
        # (this might also make it correct for ladders with different pam models per rail,
        #  in its current use in this file, which it wasn't until now,
        #  since now it returns an atom from a rail that's unchanged if self is flipped)
        """
        [overridden in DnaSingleStrandDomain]
        """
        return self.axis_rail.baseatoms[0]

    def ladder_end_atoms(self, end, reverse = False): # rename, and/or merge better with rail_end_baseatoms?
        """
        Return a list of our 3 rail-end atoms at the specified end,
        using None for a missing atom due to a missing strand2 or axis_rail,
        in the order strand1, axis, strand2 for the right end,
        or the reverse order for the left end (since flipping self
        also reverses that rail order in order to preserve strand bond
        directions as a function of position in the returned list,
        and for other reasons mentioned in some caller comments).

        @param end: LADDER_END0 or LADDER_END1

        @param reverse: caller wants result in reverse order compared to usual
        (it will still differ for the two ends)

        @warning: the order being flipped for different ends
        means that a caller looking for adjacent baseatoms
        from abutting ends of two ladders needs
        to reverse the result for one of them
        (always, given the strand direction convention).
        To facilitate and (nim)optimize this the caller can pass
        the reverse flag if it knows ahead of time it wants
        to reverse the result.
        """
        # note: top to bottom on right, bottom to top on left
        # (when reverse option is false);
        # None in place of missing atoms for strand2 or axis_rail
        assert len(self.strand_rails) in (1,2)
        strand_rails = self.strand_rails + [None] # + [None] needed in case len is 1
        rails = [strand_rails[0], self.axis_rail, strand_rails[1]] # order matters
        def atom0(rail):
            if rail is None:
                return None
            return rail.end_baseatoms()[end]
        res0 = map( atom0, rails)
        if end != LADDER_END1:
            assert end == LADDER_END0
            res0.reverse()
        if reverse:
            res0.reverse() #e could optim by xoring those booleans
        return res0
    
    def __repr__(self):
        ns = self.num_strands()
        if ns == 2:
            extra = ""
        elif ns == 1:
            extra = " (single strand)"
        elif ns == 0:
            # don't say it's an error if 0, since it might be still being made
            # (or be a bare axis ladder)
            extra = " (0 strands)"
        else:
            extra = " (error: %d strands)" % ns
        if not self.valid:
            extra += " (invalid)"
        if self.error:
            extra += " (.error set)"
        return "<%s at %#x, axis len %d%s>" % (self.__class__.__name__, id(self), self.axis_rail.baselength(), extra)

    def _do_merge_with_other_at_ends(self, other, end, other_end): # works in either class; when works, could optim for single strand
        """
        Merge two ladders, at specified ends of each one:
        self, at end, with other, at other_end.
        
        Invalidate the merged ladders; return the new valid ladder
        resulting from the merge.
        """
        # algorithm: (might differ for subclasses)
        # - make merged rails (noting possible flip of each ladder)
        #   (just append lists of baseatoms)
        # - put them into a new ladder
        # - when debugging, could check all the same things that
        #   finished does, assert no change
        flip_self = (end != LADDER_END1)
        flip_other = (other_end != LADDER_END0)

        if _DEBUG_LADDER_FINISH_AND_MERGE():
            print "dna updater: fyi: _do_merge_with_other_at_ends(self = %r, other_ladder = %r, end = %r, other_end = %r)" % \
              (self, other, end, other_end)
        if _DEBUG_LADDER_FINISH_AND_MERGE():
            # following was useful for a bug: 080122 noon
            print self.ladder_string("self", mark_end = end, flipQ = flip_self)
            print other.ladder_string("other", mark_end = other_end, flipQ = flip_other)
            print

        self_baseatom_lists = [rail_to_baseatoms(rail, flip_self)
                               for rail in self.all_rail_slots_from_top_to_bottom()]
        other_baseatom_lists = [rail_to_baseatoms(rail, flip_other)
                               for rail in other.all_rail_slots_from_top_to_bottom()]
        if flip_self:
            # the individual lists were already reversed in rail_to_baseatoms
            self_baseatom_lists.reverse()
        if flip_other:
            other_baseatom_lists.reverse()
        # now (after flipping) we need to attach self END1 to other END0,
        # i.e. "self rail + other rail" for each rail.
        def concatenate_rail_atomlists((r1, r2)):
            assert (r1 and r2) or (not r1 and not r2)
            res = r1 + r2
            assert res is not r1
            assert res is not r2
            return res
        new_baseatom_lists = map( concatenate_rail_atomlists,
                                  zip(self_baseatom_lists, other_baseatom_lists))
        # flip the result if we flipped self, so we know it defines a legal
        # set of rails for a ladder even if some rails are missing
        # (todo: could optim by never flipping self, instead swapping
        #  ladders and negating both flips)
        if flip_self:
            if _DEBUG_LADDER_FINISH_AND_MERGE():
                print "dna updater: fyi: new_baseatom_lists before re-flip == %r" % (new_baseatom_lists,)
            new_baseatom_lists.reverse() # bugfix 080122 circa 2pm: .reverse -> .reverse() [confirmed]
            for listi in new_baseatom_lists:
                listi.reverse()
        if _DEBUG_LADDER_FINISH_AND_MERGE():
            print "dna updater: fyi: new_baseatom_lists (after flip or no flip) == %r" % (new_baseatom_lists,)
        # invalidate old ladders before making new rails or new ladder
        self.ladder_invalidate_and_assert_permitted()
        other.ladder_invalidate_and_assert_permitted()
            # ok at this stage? should be -- in fact, caller
            # explicitly permits it to work, via global dnaladder_inval_policy
            # [new feature in caller, bruce 080413]
        # now make new rails and new ladder
        new_rails = [ _new_rail( baseatom_list, strandQ, bond_direction)
                      for baseatom_list, strandQ, bond_direction in
                      zip( new_baseatom_lists,
                           (True, False, True),
                           (1, 0, -1 ) # TODO: grab from a named constant?
                          ) ]
        new_ladder = _new_ladder(new_rails)
        if debug_flags.DEBUG_DNA_UPDATER:
            print "dna updater: fyi: merged %r and %r to produce %r" % (self, other, new_ladder,)
        return new_ladder

    def ladder_string(self, name = None, flipQ = False, mark_end = None):
        """
        Return a multi-line string listing all baseatoms in self,
        labelling self as name, marking one end of self if requested by mark_end,
        flipping self if requested (mark_end is interpreted as referring
        to the unflipped state). Even empty rails are shown.

        @type name: string or None
        @type flipQ: boolean
        @type mark_end: None or LADDER_END0 or LADDER_END1
        """
        if name is None:
            name = "%r" % self
        assert LADDER_END0 == 0 and LADDER_END1 == 1
            # make it easier to compute left/right mark
        classname = self.__class__.__name__.split('.')[-1]
        flipped_string = flipQ and " (flipped)" or ""
        label = "%s %r%s:" % (classname, name, flipped_string)
        baseatom_lists = [rail_to_baseatoms(rail, flipQ)
                               for rail in self.all_rail_slots_from_top_to_bottom()]
        if flipQ:
            # the individual lists were already reversed in rail_to_baseatoms
            baseatom_lists.reverse()
        if mark_end is not None:
            if flipQ:
                # make it relative to the printout
                mark_end = (LADDER_END0 + LADDER_END1) - mark_end
        left_mark = (mark_end == LADDER_END0) and "* " or "  "
        right_mark = (mark_end == LADDER_END1) and " *" or "  "
        strings = [left_mark +
                   (baseatoms and ("%r" % baseatoms) or "------") +
                   right_mark
                   for baseatoms in baseatom_lists]
        return "\n".join([label] + strings)
    
    # ==

    def all_rails(self):
        """
        Return a list of all our rails (axis and 1 or 2 strands), in arbitrary order,
        not including missing rails.
        [implem is subclass-specific]
        """
        return [self.axis_rail] + self.strand_rails

    # TODO: make self.strand_rails private; define a method of the same name for external use,
    # to be uniform with axis_rails and all_rails methods [bruce 080224 comment]

    def axis_rails(self):
        """
        Return a list of all our axis rails (currently, 1 in this class,
        0 in some subclasses), in arbitrary order,
        not including missing rails.
        [implem is subclass-specific]
        """
        return [self.axis_rail]
    
    def all_rail_slots_from_top_to_bottom(self):
        """
        Return a list of all our rails (axis and 1 or 2 strands), in top to bottom order,
        using None for missing rails.
        
        [implem is subclass-specific]
        """
        if len(self.strand_rails) == 2:
            return [self.strand_rails[0], self.axis_rail, self.strand_rails[1]]
        elif len(self.strand_rails) == 1:
            return [self.strand_rails[0], self.axis_rail, None]
        elif len(self.strand_rails) == 0 and pref_permit_bare_axis_atoms():
            return [None, self.axis_rail, None] # not sure this will work in callers;
                # note: this works even if self.axis_rail is None due to bugs
        else:
            # this happens, after certain bugs, in debug code with self.error already set;
            # make sure we notice if it happens without self.error being set: [bruce 080219]
            if not self.error:
                self.error = "late-detected wrong number of strands, %d" % len(self.strand_rails)
                print "\n***BUG: %r: %s" % (self, self.error)
            return [None, self.axis_rail, None] # same as above, see comment there
        pass

    # ==
    
    def remake_chunks(self, rails = None):
        """
        Make new chunks for self, as close to the old chunks as we can
        (definitely in the same Part; in the same Group if they were in one).

        @param rails: list of rails of self to make them for, or None (default)
                      to make them for all rails of self
        
        @return: list of all newly made (by this method) DnaLadderRailChunks
                 (or modified ones, if that's ever possible)
        """
        assert self.valid
        # but don't assert not self.error
        res = []
        ladder_color = None # might be reassigned
        if pref_per_ladder_colors():
            from dna.model.Dna_Constants import getNextStrandColor
            # todo: apply this in Chunk.drawing_color (extend that method),
            # so it doesn't erase the saved color
            ladder_color = getNextStrandColor()
        # note: if self.error, Chunk.draw will notice this and use black
        # (within its call of DnaLadderRailChunk.modify_color_for_error);
        # we won't reassign ladder_color, since we want old chunk colors
        # to last, so they can be preserved until errors are removed.
        if rails is None:
            rails = self.all_rails()
        for rail in rails:
            if rail is self.axis_rail:
                constructor = make_or_reuse_DnaAxisChunk
            else:
                constructor = make_or_reuse_DnaStrandChunk
            old_chunk = rail.baseatoms[0].molecule # from arbitrary baseatom in rail
            part = old_chunk.part
            if not part:
                # will this happen in MMKit? get it from assy
                assy = old_chunk.assy
                print "\nbug: %r in assy %r has no .part" % \
                      (old_chunk, assy)
                assert assy is self.assy
                # note: when some bugs happen, assy is None here [080227]
                assert assy is not None, "%r.remake_chunks notices self.assy is None" % self
                part = assy.part
                assert part, "%r.remake_chunks(): %r.part is None" % (self, assy)
                    # will this work in MMKit?
            # todo: assert rail atoms are in this part
            part.ensure_toplevel_group()
            group = old_chunk.dad # put new chunk in same group as old
                # (but don't worry about where inside it, for now)
            if group is None:
                # this happens for MMKit chunks with "dummy" in their names;
                # can it happen for anything else? should fix in mmkit code.
                if "dummy" not in old_chunk.name:
                    print "dna updater: why is %r.dad == None? (old_chunk.assy = %r)" % (old_chunk, old_chunk.assy) ###
                group = part.topnode
            assert group.is_group()
            ## assert group.part is part, \
            if not (group.part is part):
                print "\n*** BUG: " \
                   "group.part %r is not part %r, for old_chunk %r, .dad %r, part.topnode %r" % \
                   (group.part, part, old_chunk, old_chunk.dad, part.topnode)
                # this is failing in Undo of creating or bonding lone Ss3, or Undo/Redo/Undo of same,
                # with group.part None, group probably a DnaStrand which doesn't exist in new state(?).
                # [bruce 080405 comment, might relate to "tom's undo bug" reported by email yesterday]
            chunk = constructor(self.assy, rail, self)
                # Note: these constructors need to be passed self,
                # in case they reuse an old chunk, which self takes over.
                # (Even though in theory they could figure out self somehow.)
                #
                # this pulls in all atoms belonging to rail
                # (including certain kinds of attached atoms);
                # it also may copy in certain kinds of info
                # from their old chunk (e.g. its .hidden state?)
            if debug_flags.DEBUG_DNA_UPDATER_VERBOSE:
                print "%r.remake_chunks made %r" % (self, chunk)
            chunk.color = old_chunk.color # works, at least if a color was set
            #e put it into the model in the right place [stub - puts it in the same group]
            # (also - might be wrong, we might want to hold off and put it in a better place a bit later during dna updater)
            # also works: part.addnode(chunk), which includes its own ensure_toplevel_group
            if ladder_color is not None:
                chunk.color = ladder_color
            if self.error:
                #bruce 080401 kluge, not yet fully correct, explained in _f_convert_pam_if_desired
                chunk.display_as_pam = None
            group.addchild(chunk)
            # todo: if you don't want this location for the added node chunk,
            # just call chunk.moveto when you're done,
            # or if you know a group you want it in, call group.addchild(chunk) instead of this.
            res.append(chunk)
        return res

    def _f_reposition_baggage(self): #bruce 080404
        """
        [friend method for dna updater;
         unsafe to call until after self.remake_chunks()]
        """
        DEBUG_BONDPOINTS = False # DO NOT COMMIT WITH TRUE
        # only for atoms with _f_dna_updater_should_reposition_baggage set
        # (and optimizes by knowing where those might be inside the ladder)
        if DEBUG_BONDPOINTS:
            print "_f_reposition_baggage called on", self
        for atom in self._atoms_needing_reposition_baggage_check():
            # (for single strands, we need to scan all the baseatoms)
            if atom._f_dna_updater_should_reposition_baggage:
                if DEBUG_BONDPOINTS:
                    print " calling %r.reposition_baggage_using_DnaLadder()" % atom
                atom.reposition_baggage_using_DnaLadder()
                    # fyi: as of 080404, the only case this handles is giving
                    # an Ax with one axis bond a parallel open bond.
                    # Note: only safe to call this after our chunks were remade!
                    # and, if it turns out this ever needs to look inside
                    # neighbor atom dnaladders (not in the same base pair),
                    # we'll need to revise the caller to call it in a separate
                    # loop over ladders than the one that calls remake_chunks
                    # on all of them.
                atom._f_dna_updater_should_reposition_baggage = False
                    # probably redundant with that method
        #e assert no other atoms have that flag? Not true! interior Ax for single strand might have it.
        return

    def _atoms_needing_reposition_baggage_check(self): #bruce 080405
        """
        [private helper for _f_reposition_baggage]

        [subclasses with interior bondpoints must override this implem]
        """
        return self.rail_end_baseatoms() # for duplex ladders
    
    # == chunk access methods

    def strand_chunks(self):
        """
        Return a list of all the strand chunks (DnaStrandChunk objects) of this DnaLadder.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        ###FIX: self.strand_rails is an attr, self.axis_rails is a method.
        return [rail.baseatoms[0].molecule for rail in self.strand_rails]

    def axis_chunks(self):
        """
        Return a list of all the axis chunks (DnaAxisChunk objects) of this DnaLadder.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        return [rail.baseatoms[0].molecule for rail in self.axis_rails()]

    def all_chunks(self):
        """
        Return a list of all the strand and axis chunks (DnaLadderRailChunk objects) of this DnaLadder.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        return [rail.baseatoms[0].molecule for rail in self.all_rails()]

    # == convenience methods
    
    def kill_strand_chunks(self):
        """
        Kill all the strand chunks (DnaStrandChunk objects) of this DnaLadder.

        @note: doesn't invalidate the ladder -- that's up to the dna updater
               to do later if it wants to.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        for chunk in self.strand_chunks():
            chunk.kill()

    def kill_axis_chunks(self):
        """
        Kill all the axis chunks (DnaAxisChunk objects) of this DnaLadder.

        @note: doesn't invalidate the ladder -- that's up to the dna updater
               to do later if it wants to.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        for chunk in self.axis_chunks():
            chunk.kill()
        return
    
    def kill_all_chunks(self):
        """
        Kill all the strand and axis chunks (DnaLadderRailChunk objects) of this DnaLadder.

        @note: doesn't invalidate the ladder -- that's up to the dna updater
               to do later if it wants to.

        @see: our methods strand_chunks, all_chunks, axis_chunks,
              kill_ versions of each of those.
        """
        for chunk in self.all_chunks():
            chunk.kill()
        return
    
    # ==

    def dnaladder_menu_spec(self, selobj): #bruce 080401, renamed 080407, split 080409
        """
        [called by menu_spec methods for various selobj classes]
        
        Return a menu_spec list of context menu entries specific to selobj being
        a PAM atom or chunk whose DnaLadder is self. (If we have no entries,
        return an empty menu_spec list, namely [].)

        (Other kinds of selobj might be permitted later.)
        """
        if isinstance(selobj, self.assy.Chunk):
            # don't add these entries to Chunk's cmenu.
            # [bruce 080723 refactoring of recent Mark change in Chunk.py]
            return
        res = []
        nwhats = self._n_bases_or_basepairs_string()
        if self.error or not self.valid:
            # (this prevents reoffering a failed pam conversion
            #  until you modify the ladder atoms enough for the
            #  dna updater to remake the ladder)
            # If there's a problem with this ladder, report the status
            # and don't offer any other menu entries for it.
            if not self.valid:
                text = "invalid DnaLadder (bug?) (%s)" % nwhats
                res.append( (text, noop, 'disabled') )
            if self.error:
                text = "%s (%s)" % (self.error, nwhats)
                res.append( (text, noop, 'disabled') )
            if self.error and self.valid:
                text = "rescan %s for errors" % nwhats
                cmd = self.cmd_rescan_for_errors
                res.append( (text, cmd) )
            pass
        else:
            # valid, no error
            res.extend( self._pam_conversion_menu_spec(selobj) )
            text = "scan %s for errors" % nwhats
            cmd = self.cmd_rescan_for_errors
            res.append( (text, cmd) )
        return res

    def _n_bases_or_basepairs_string(self):
        nstrands = len(self.strand_rails)
        if nstrands == 2:
            what = "basepair"
        elif nstrands == 1:
            what = "base"
        elif nstrands == 0:
            what = "bare axis atom"
        else:
            what = "wrongly structured base"
        text = "%d %s(s)" % (len(self), what)
        text = fix_plurals(text)
            # note: works here, wouldn't work if text had another ')' after it
            # (due to trivial bug in fix_plurals)
        return text

    def cmd_rescan_for_errors(self): #bruce 080412
        self._dna_updater_rescan_all_atoms()
        # review: does this miss anything?
        # the following (or something like it) seems to be needed,
        # due to bugs in handling duplex errors I guess
        # (e.g. for relative motion of rails w/o distorting any rail)
        for chunk in self.all_chunks():
            chunk.changeapp(0)
        # also gl_update? doesn't seem needed by test
        return

    def _dna_updater_rescan_all_atoms(self): #bruce 080412 split this out
        """
        tells dna updater to remake self from scratch (when it next runs)
        """
        self.ladder_invalidate_and_assert_permitted()
        # BUG (but will probably become officially ok, with the following
        #  workaround officially required and documented; see comments in
        #  initial code called by full_dna_update [comment revised, bruce 080413]):
        # the above, alone, doesn't cause dna updater to actually run.
        # KLUGE WORKAROUND -- also do this:
        self.arbitrary_baseatom()._changed_structure()
        return

    def strand_neighbor_ladders(self): #bruce 080408
        # (mostly superseded by _f_ladders_with_up_to_date_baseframes_at_ends;
        #  but still used for bug safety)
        """
        Return a sequence of 0 to 4 ladders (or Nones) which includes all
        ladders which are connected to self at a "corner" (by a strand that
         leaves self and immediately enters the other ladder, possibly by
         crossing a "bridging Pl atom").

        Always return them in the order strand1-left,
        strand1-right, strand2-left, strand2-right, even if this list
        contains None (for an end there), self, or the same ladder twice.
        
        The length of the returned list is always exactly twice the number of
        strand rails of self.

        @note: this is safe to call on freshly made valid ladders with no
               errors set which didn't yet remake their chunks.
        """
        if not self.valid:
            # was assert, but happened for convert to pam3 of ladder
            # whose self & neighbor were converted to PAM5 [bruce 080413 336pm]
            print "possible bug: strand_neighbor_ladders called on %r " \
                  " which is not valid" % (self,)
        if self.error:
            print "possible bug: strand_neighbor_ladders called on %r " \
                  " which has error = %r" % (self, self.error)
        res = []
        for rail in self.strand_rails:
            for end in LADDER_ENDS:
                res += [self.neighbor_ladder_at_corner(rail, end)]
        return res

    def neighbor_ladder_at_corner(self, rail, end):
        """
        @param rail: an element of self.strand_rails (not checked)

        @end: an element of the constant sequence LADDER_ENDS

        @return: the neighbor ladder joined to self by a strand at the
                 specified corner of self (perhaps self), or None
                 if that strand ends at that corner.

        @note: can be wrong (or raise an exception) if self.error or
               not self.valid (not checked).

        @note: this is safe to call on freshly made valid ladders which
               didn't yet remake their chunks.
        """
        bond_direction_to_other = LADDER_BOND_DIRECTION_TO_OTHER_AT_END_OF_STRAND1[end]
        if rail is not self.strand_rails[0]:
            bond_direction_to_other = - bond_direction_to_other
                # only correct if not self.error
        # this is not set yet: next_atom = rail.neighbor_baseatoms[end]
        end_atom = rail.end_baseatoms()[end]
        next_atom = end_atom.strand_next_baseatom(bond_direction = bond_direction_to_other)
            # (note: strand_next_baseatom returns None if end_atom or the atom it
            #  might return has ._dna_updater__error set, or if it reaches a non-Ss atom.)
        if next_atom is not None:            
            return rail_end_atom_to_ladder(next_atom)
                # note: if self.error is set due to wrong strand directions,
                # then we may have looked in the wrong bond direction for
                # next_atom, and probably found an interior atom of self,
                # in which case this will raise an exception.
                # [bruce 080523 comment]
        return None
        
    pass # end of class DnaLadder

# ==

# @@ review points:
# - this subclass should be a sibling class -- these are both kinds of
#   "ladder-like DnaDomains", and this has less not more code than DnaLadder;
#   maybe about half of DnaLadder's code looks like it's common code
#   (what's the right classname for the interface?
#    single or duplex domain? domain with some strands?)
# - the ladder api from the rest of the code (e.g. chunk.ladder)
#   is really a "ladder-like dna domain" api, since it work for single strand
#   ladders too (i.e. DnaLadder is not a perfect name for the common superclass)

class DnaSingleStrandDomain(DnaLadder):
    """
    Represent a free floating single strand domain;
    act exactly like a ladder with no axis rail or second strand rail
    in terms of ability to interface with other code.
    """
    valid = False # same as in super; set to True in _finish_strand_rails
    
    def __init__(self, strand_rail):        
        self.axis_rail = None
        self.assy = strand_rail.baseatoms[0].molecule.assy
        self.strand_rails = [strand_rail]
        assert self.assy is not None, "%r.__init__: assy is None" % self
        self._finish_strand_rails()
            # check bond direction, reverse if needed
        return
    
    def baselength(self):
        return len(self.strand_rails[0])
    
    def finished(self):
        assert 0, "illegal in this class" # or could be a noop
        
    def add_strand_rail(self, strand_rail):
        assert 0, "illegal in this class"
        
    def all_rails(self):
        """
        Return a list of all our rails (1 strand).
        [implem is subclass-specific]
        """
        return self.strand_rails
    
    def axis_rails(self):
        """
        Return a list of all our axis rails (none, in this subclass).
        [implem is subclass-specific]
        """
        return []
    
    def all_rail_slots_from_top_to_bottom(self):
        """
        Return a list of all our rails (1 strand), in top to bottom order,
        using None for missing rails (so length of return value is 3).
        
        [implem is subclass-specific]
        """
        return [self.strand_rails[0], None, None]

    def arbitrary_baseatom(self): #bruce 080407
        """
        [overrides DnaLadder implem]
        """
        return self.strand_rails[0].baseatoms[0]

    def _atoms_needing_reposition_baggage_check(self): #bruce 080405
        """
        [private helper for _f_reposition_baggage]

        [overrides duplex implem in superclass]
        """
        return self.strand_rails[0].baseatoms
        
    def __repr__(self):
        ns = self.num_strands()
        if ns == 1:
            extra = ""
        else:
            extra = " (error: %d strands)" % ns
        if self.axis_rail:
            extra += " (error: has axis_rail)"
        classname = self.__class__.__name__.split('.')[-1]
        return "<%s at %#x, len %d%s>" % (classname, id(self), self.baselength(), extra)
    pass

# ==

def _end_to_end_bonded( atom1, atom2, strandQ):
    """
    Are the expected end-to-end bonds present between end atoms
    from different ladders? (The atoms are None if they are from
    missing rails, i.e. strand2 of any single-stranded ladder
    or axis_rail of a free-floating single-stranded ladder;
    we expect no bonds between two missing atoms, but we do need
    bonds between None and a real end-atom, so we return False then.)

    param atom1: an end atom of a ladder rail, or None (for a missing rail in a ladder)

    param atom2: like atom1, for a possibly-adjacent ladder
    
    param strandQ: whether these atoms are part of strand rails (vs axis rails).
    type strandQ: boolean
    """
    if atom1 is None and atom2 is None:
        # note: this can happen for strandQ or not
        return True
    if atom1 is None or atom2 is None:
        # note: this can happen for strandQ or not
        return False
    if atoms_are_bonded(atom1, atom2):
        # note: no check for strand bond direction (assume caller did that)
        # (we were not passed enough info to check for it, anyway)
        return True
    if strandQ:
        # also check for an intervening Pl in case of PAM5
        return strand_atoms_are_bonded_by_Pl(atom1, atom2)
    else:
        return False

def strand_atoms_are_bonded_by_Pl(atom1, atom2): #e refile 
    """
    are these two strand base atoms bonded (indirectly) by means of
    a single intervening Pl atom (PAM5)?
    """
    if atom1 is atom2:
        # assert 0??
        return False # otherwise following code could be fooled!
# one way to support PAM5:
##    set1 = atom1.Pl_neighbors() # [this is defined as of 080122]
##    set2 = atom2.Pl_neighbors()
##    for Pl in set1:
##        if Pl in set2:
##            return True
##    return False
    # another way, easier for now:
    set1 = atom1.neighbors()
    set2 = atom2.neighbors()
    for Pl in set1:
        if Pl in set2 and Pl.element is Pl5:
            return True
    return False

# ==

# helpers for _do_merge_with_other_at_ends

def rail_to_baseatoms(rail, flipQ):
    if not rail:
        return []
    elif not flipQ:
        return rail.baseatoms
    else:
        res = list(rail.baseatoms)
        res.reverse()
        return res

def _new_rail(baseatoms, strandQ, bond_direction):
    """
    """
    if baseatoms == []:
        # special case
        return None
    assert len(baseatoms) > 1 # since result of merge of nonempty rails
        # (might not be needed, but good to sanity-check)
    if strandQ:
        assert bond_direction
    else:
        assert not bond_direction
    if _DEBUG_LADDER_FINISH_AND_MERGE():
        # print all calls
        print "\n_new_rail(%r, %r, %r)" % (baseatoms, strandQ, bond_direction)
    if debug_flags.DEBUG_DNA_UPDATER:
        # check bondedness of adjacent baseatoms
        for i in range(len(baseatoms) - 1):
            atom1, atom2 = baseatoms[i], baseatoms[i+1]
            if strandQ:
                assert atoms_are_bonded(atom1, atom2) or \
                       strand_atoms_are_bonded_by_Pl(atom1, atom2), \
                       "*** missing bond between adjacent baseatoms [%d,%d of %d]: %r and %r" % \
                       (i, i+1, len(baseatoms), atom1, atom2)
            else:
                assert atoms_are_bonded(atom1, atom2), \
                       "*** missing bond between adjacent baseatoms [%d,%d of %d]: %r and %r" % \
                       (i, i+1, len(baseatoms), atom1, atom2)
    return merged_chain( baseatoms, strandQ, bond_direction)

def _new_ladder(new_rails):
    """
    Make a new DnaLadder of the appropriate subclass,
    given 3 new rails-or-None in top to bottom order.
    
    @param new_rails: a sequence of length 3; each element is a new rail,
                      or None; top to bottom order.
    """
    strand1, axis, strand2 = new_rails
    assert strand1
    if not axis:
        assert not strand2
        return _new_single_strand_domain(strand1)
    else:
        return _new_axis_ladder(strand1, axis, strand2)
    pass

def _new_single_strand_domain(strand1):
    """
    [private helper for _new_ladder]
    """
    return DnaSingleStrandDomain(strand1)

def _new_axis_ladder(strand1, axis, strand2):
    """
    [private helper for _new_ladder]
    """
    ladder = DnaLadder(axis)
    ladder.add_strand_rail(strand1)
    if strand2:
        ladder.add_strand_rail(strand2)
    ladder.finished()
    return ladder

# end