#! /usr/bin/python3
# Last edited on 2026-10-30 16:02:49 by stolfi

import os, sys, re
from sys import stderr as err
import numpy as np
from math import inf, isfinite, floor, pow
import datetime as dtt
import time

import bencao_starps_size_position_funcs as spf
import bimatching_eval_funcs as bef
import bitemplate_funcs as btf
import bimatching_funcs as bmf
# For tests:
import bitemplate_match_test_funcs as bmtestf

def match_bitemplate(text0, text1, bitemp, eval_one_gap, debug_level):
  # Matches two given strings {text0,text1} to a bitemplate {bitemp}. If
  # there is at least one possible bimatching, chooses the one with the
  # lowest badness score. Otherwise returns an infinite badness score
  # and no bimatching.
  #
  # A /bitemplate/ is a list {bitemp[0..nr]} of /bipatterns/.
  # Each bipattern is a list of /variants/.  Each variant is a
  # triple {(pena, pat0, pat1)} where {pena} is a non-negative float,
  # typically in {[0 _ 1]}; and {pat1,pat2} are two RE patterns,
  # neither of which matches the empty string.
  #
  # Let {nt0} and {nt1} be the lengths of {text0} and {text1}. A
  # /bimatching/ of {text0,text1} to {bitemp} is a list {bimatch} of
  # {nr} /rungs/. A rung {bimatch[ih]} is a quintuple
  # {(i0,j0,i1,j1,bpvar)} where {bpvar} is a variant of the bipattern
  # {bitemp[ih]}, as above; {i0,j0} are indices in {0..nt0} with {i0 <=
  # j0}, and the string {text0[i0:j0]} matches the RE pattern {pat0} in
  # {bpvar}; and ditto for {i1,j1,text1} with {pat1} of {bpvar}. Among
  # successive
  #
  # The procedure effectively enumerates all the possible bimatchings
  # {bimatch[0..nr-1]} of {text0,text1} against the bitemplate
  # {bitemp[0..nr-1]}. For each bimatching it computes a badness
  # {score}, and selects the bimatching that yielded the lowest {score}.
  #
  # Actually, this version uses a dynamic programming algorithm
  # to find the best match in time {O(nr*nt0^2*nt1^2)} -- instead of
  # enumerating all bimatching, which would take exponential time.
  #
  # The badness {score} of a bimatching is computed by converting it
  # implicitly to a pair of macro-parsings {segs0,segs1} of {text0} and
  # {text1}, and calling {eval_one_gap(len(gap0),len(gap1),ig,ng)} on
  # each pair of corresponding gaps. The {score} is the sum of all the
  # gap scores, plus the penalties of the variants used to match the
  # keywords.
  #
  # Returns the bimatching {opt_bimatch} that yielded the smallest
  # {score}, as well as its score.
  #
  # If there are no possible bimatchings, returns {None,+inf}.
  
  # debug_all = (text1[:15] == "pokararkeeyokee")
  # debug_all = (text1[:13] == "poarkeeodaiin")

  debug_top = debug_level >= 4
  debug_tableau = debug_level >= 8
  debug_hcands = debug_level >= 7
  debug_plane = debug_level >= 6
  debug_vertex = debug_level >= 10
  debug_cost = debug_level >= 12
  debug_edge = debug_level >= 20

  verbose = True

  nt0 = len(text0); nt1 = len(text1)  # Lengths of the texts.
  nr = len(bitemp)                    # Number of bipatterns.
  nh = nr; ng = nh + 1; ns = ng + nh  # Number of hits, gaps, segs.
  assert ns >= 1

  if debug_top:
    err.write(f"!~ enter match_bitemplate\n")
    err.write(f"!~ {nr = } {nh = } {ng = } {ns = }\n")
    err.write(f"!~ {nt0 = :3d} {text0 = !r}\n")
    err.write(f"!~ {nt1 = :3d} {text1 = !r}\n")

  if (debug_cost): plot_cost(eval_one_gap)

  tot = { }
  tot['evgap'] = 0  # Number of gap scores computed.
  tot['evhit'] = 0  # Number of attempts to match a hit pair.
  tot['nohit'] = 0  # Number of hit pairs that failed to match.
  tot['finel'] = 0  # Number of finite elements in whole tableau.

  # The dynamic programming tableau is a graph with the vertices organized into
  # {ns+1} /planes/, inclding a /source vertex/ (plane {-1}), 
  # {ns-1} ordinary planes indexed {0..ns-2}, and a /sink vertex/ (plane {ns-1}).
  # 
  # A /breakpoint/ on {text0} is an index in {0..nt0} (note: not {0..nt0-1}).
  # Likewise, a breakpoint on text1 is an index in {0..nt1}.
  #
  # In each plane, each vertex is identified by two indices {(i0,i1)}
  # which are breakpoints of {text0} and {text1}. The source vertex has
  # indices {(0,0)} and the sink vertex has indices {(nt0,nt1)}. In the
  # ordinary planes the indices {i0} and {i1} are in {0..nt0} and
  # {0..nt1}.  Thus each plane could have up to {(nt0+1)*(nt1+1)} vertices
  # But usually most of those potential vertices are never created.
  #
  # A /biparsing/ of size {ms <= ns} is a list of pairs {bm[ks] =
  # (b0[ks],b1[ks])} for {ks} in {0..ms-1}, such that {b0[0..ms-1]} and
  # {b1[0..ks-1]} are non-decreasing breakpoints on each text.
  # If {ms=ns}, {bm[ms-1]} must be {(nt0,nt1)}. By convention {bm[-1]}
  # is defined as {(0,0)}.  Thus a biparsing defines a path 
  # in the graph from the source vertex to a vertex with indices{bm[ms-1]}
  # on plane {ms-1}.
  #
  # A biparsing of size {ms} defines a partition of a prefix of {text0}
  # and a prefix of {text1} into {ms} pairs of /segments/
  # {(segs0[ks],segs1[ks])} for {ks} in {0..ms-1}; where {segs0[ks]} is
  # the substring of {text0} between breakpoints {b0[ks-1]} and
  # {b0[ks]}. The pair {bm[ms-1]} is said to be the /end/ of the
  # biparsing.
  #
  # If {ks} is even, segments {segs0[ks]} and {segs1[ks]} are called
  # /gaps/. If {ks} is odd they are called /hits/. The number of
  # segments {ns} is always odd, so in a full biparsing (with {ms=ns})
  # the last segment pairs are gaps.
  # 
  # The edges of the graph connect consecutive planes {ks-1} to {ks}.
  # Each edge represents a possible segment pair of a biparsing of
  # {text0} and {text1}.
  #
  # The tableau is represented here by a list {tb} with {ns} dicts. The
  # attributes of each vertex with indices {(i0,i1)} of plane {ks} in
  # {0..ns-1} (thus except the source) are represented by a tuple
  # {(k0,k1,score,ivar)} described below. These tuples are stored in the
  # dict {tb[ks]} with {(i0,i1)} as the key.

  tb = []
  for ks in range(ns): tb.append(dict())

  # Each edge between a vertex {(k0,k1) of plane {ks-1} and vertex
  # {(i0,i1) of plane {ks} has an /edge score/ that measures its badness. 
  #
  # If the edge is a gap edge, namely {ks = 2*ig} for some {ig}, the
  # score is a function of the sizes of the two segments, namely {i0-k0}
  # and {i1-k1}, weighted by a factor that depends on the index {ig} and
  # the total number of gaps {ng=ns//2+1}.
  #
  # If the edge is a hit edge, namely {ks = 2*ih+1} for some {ih}, the
  # two segments are supposed to match the two patterns of some variant
  # {(pena, pat0, pat1)} of the bipattern {bitemp[ih]}, and the score of
  # the edge is the penalty {pena} defined for that variant. If there
  # are multiple variants of {bitemp[ih]} that match those two segments,
  # the score of the edge is the minimum of the corresponding penalties.
  # 
  # The /score/ of a biparsing of any size {ms<=ns} is
  # the sum of the scores of its edges.
  #
  # Each vertex on a plane {ks} has a /vertex score/ that is the minimum score
  # of any biparsing of size {ks+1} that ends at that vertex of that plane.
  # Only vertices with finite score are stored in the tableau.
  # 
  # Each entry {tb[ks][(i0,i1)]} is a tuple {(k0,k1,score,ivar,nb)} where
  # {(k0,k1)} is the previous vertex on the minimum-cost bimatching that
  # ends at vertex {(i0,i1)} of plane {ks}, {score} is the vertex score
  # (the total score of that bimatching), and {nb} is the total number
  # of bimatchings with finite score that end at that vertex.
  #
  # If the last segment was a hit, namely {ks=2*ih+1}, then 
  # the field {ivar} is the index of the variant of the bipattern {bitemp[ih]}
  # that matches the two corresponding segments, and thus defined the 
  # cost of the edge.  If the last segment was a gap, the field {ivar}
  # is {-1}.
  #
  # The goal is to fill each dict {tb[ks]} for {ks} in {0..ns-1} in succession,

  def gap_score(ks, k0,i0, k1,i1):
    # Badness score of the gap segments {seg00[ks]} and {segs1[ks]},
    # between breakpoints {(k0,k1)} and {(i0,i1)}.
    # Requires {ks} even, {ks == 2*ig}.
    nonlocal tot
    
    debug_gap_score = False
    if debug_gap_score:
      err.write(f"!$   computing score for gaps {k0}..{i0} and {k1}..{i1}\n")
    assert 0 <= k0 and k0 <= i0 and i0 <= nt0
    assert 0 <= k1 and k1 <= i1 and i1 <= nt1
    assert ks % 2 == 0
    gsz0 = i0 - k0; gsz1 = i1 - k1
    ig = ks//2
    tot['evgap'] += 1
    score_gap = eval_one_gap(gsz0, gsz1, ig, ng)
    return score_gap
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def hit_score(ks, k0,i0, k1,i1):
    # Badness score of the hit segments {segs0[ks]} and {segs1[ks]},
    # between breakpoints {(k0,k1)} and {(i0,i1)}.
    # Requires {ks} odd, {ks == 2*ih + 1}.
    # If the segments match the bipattern {bitemp[ih]}, returns
    # the penalty {pena} of the variant that matched, and its index {ivar}.
    # Otherwise returns {+inf,-1}.
    nonlocal tot
    assert 0 <= k0 and k0 <= i0 and i0 <= nt0
    assert 0 <= k0 and k0 <= i0 and i0 <= nt0
    assert 0 <= k1 and k1 <= i1 and i1 <= nt1
    assert ks % 2 == 1
    ih =  ks//2
    bipat = bitemp[ih]
    nv = len(bipat)
    pena_prev = -1
    tot['evhit'] += 1
    seg0 = text0[k0:i0]
    seg1 = text1[k1:i1]
    for iv in range(nv):
      pena, pat0, pat1 = bipat[iv]
      assert pena > pena_prev, f"key penalties not increasing {pena_prev} {pena}"
      if re.fullmatch(pat0, seg0):
        if re.fullmatch(pat1, seg1):
          return pena, iv
      pena_prev = pena
    tot['nohit'] += 1
    return +inf, -1
    # :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def print_vertex(vert):
    # Prints the vertex tuple {vert}.
    k0, k1, score, ivar, nb = vert
    err.write(f"( {k0:3d} {k1:3d}  {score:9.5f}  {ivar:2d}  {nb:11d} )")
    return 
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  
  def print_tableau_plane(ks):
    # Prints plane {ks} of the tableau.
    # Only the score for now.
    nt_max = 30
    if nt0 <= nt_max and nt1 <= nt_max:
      rt0 = nt0; st0 = nt0 + 1
      rt1 = nt1; st1 = nt1 + 1
    else:
      ht = (nt_max - 2)//2
      rt0 = ht; st0 = nt0 - ht
      rt1 = ht; st1 = nt1 - ht

    # !!! Not much sense printin in matrix form... !!!
    for key, val in tb[ks].items():
      i0, i1 = key
      err.write(f"!o ( {i0:3d} {i1:3d} ) = ")
      print_vertex(val)
      err.write("\n")
    return
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def dump_tableau():
    # Dumps the tableau to file "/tmp/.tableau-py-{SS}.txt", where {SS}
    # is a random two-digit number, in a format that can be compared
    # with the C tableau.
    
    cdate = dtt.datetime.now(dtt.UTC).isoformat(timespec='seconds', sep='-')
    sec = re.sub(r"[+]00:00", "", cdate)
    sec = re.sub(r"^[-0-9:]+[:]", "", sec)
    fname = f"/tmp/.tableau-py-{sec}.txt"
    wr = open(fname, "w")
    err.write(f"dumping tableau to \"{fname}\" ... \n");
    wr.write(f"nt = {nt0} {nt1}\n");
    wr.write(f"np = {ns}\n");
    for ks in range(ns):
      vts = tb[ks].items()
      wr.write(f"plane {ks} - {len(vts)} vertices\n")
      for key, vert in vts:
        i0, i1 = key
        wr.write(f"[{ks:3d}][{i0:3d},{i1:3d}]");
        k0, k1, score, ivar, nb = vert
        iscore = int(floor(score*1000 + 0.5))
        wr.write(f"  {iscore:5d}")
        wr.write(f"  from [{ks-1:3d}][{k0:3d},{k1:3d}]");
        if (ks % 2) == 1: wr.write(f" [{ivar:3d}]")
        wr.write("\n")
    wr.close()
    time.sleep(1.0)
    return
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  
  def extract_bimatch():
    # Recover the optimum bimatching from {tb}, assuming the 
    # sink is the only vertex on plane {ns-1}.
    #
    # Returns the bimatching {bimatch}, its badness score {score_fin},
    # and the number {nb_fin} of distinct paths (bimatchings) that would lead
    # lead to that vertex.
    #
    # If there is no vable bimatching, returns {None,+inf,0}.

    assert ns >= 1 and (ns % 2) == 1
    if len(tb[ns-1]) == 0:
      return None, +inf, 0
    assert len(tb[ns-1]) == 1

    i0 = nt0; i1 = nt1
    bimatch = [ None ] * nr  # Rungs of partial bimatching.

    score_fin = None
    nb_fin = None
    ih = nr - 1; ks = ns - 1
    while ih >= 0:
      assert ks >= 1
      assert ks % 2 == 0
      assert i0 >= 0 and i1 >= 0

      # Skip the gap edge ending at {(i0,i1)}:
      ig = ks//2
      k0, k1, score, ivar, nb = tb[ks][(i0,i1)]
      assert isfinite(score) and score >= 0
      assert ivar == -1
      assert nb >= 1
      assert 0 <= k0 and k0 <= i0
      assert 0 <= k1 and k1 <= i1

      if ks == ns-1: 
        score_fin = score; nb_fin = nb

      if debug_plane: err.write(f"!B skipped gap {ig} segs {ks} {k0}..{i0} {k1}..{i1}\n")

      # Get rung {ih}.
      i0 = k0; i1 = k1
      ks = ks - 1;

      assert ks >= 0
      assert ks % 2 == 1
      assert ks == 2*ih + 1

      k0, k1, score, ivar, nb = tb[ks][(i0,i1)]
      assert isfinite(score) and score >= 0
      assert ivar >= 0 and ivar < len(bitemp[ih])
      assert nb >= 1
      assert 0 <= k0 and k0 <= i0
      assert 0 <= k1 and k1 <= i1

      bpvar = bitemp[ih][ivar]
      bimatch[ih] = (k0,i0, k1,i1, bpvar)
      pena, pat0, pat1 = bpvar

      if debug_plane: 
        err.write(f"!B took hits {ih} segs {ks} {k0}..{i0} {k1}..{i1} {ivar = } {pena = :.8f}\n")

      # Move back one more:
      ks = ks - 1
      ih = ih - 1
      i0 = k0; i1 = k1
    assert ks == 0
    assert ih == -1
    return bimatch, score_fin, nb_fin
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def find_best_extension(ks, i0, i1, ecand_list):
    # Finds the best biparsing with {ks} segments from {(0,0)} to {(i0,i0)},
    # by extending one of the optimal biparsings with {ks-1} segments
    # already found.
    #
    # The index {ks} must be in {0..ns}. If {ks} is {ns-1}, then
    # {(i0,i1)} must be {(nt0,nt1)}. If {ks > 0}, assumes that plane
    # {tb[ks-1]} has been filled.
    #
    # Also assumes that {ecand_list} is a list of /candidate
    # edges/, the triples {(k0,k1,score,ivar)} such that {k0 <= i0}, {k1 <= i1}, 
    # and the segments {k0..i0} and {k1..i1} are worth considering,
    # where {score} is the additional score amount that would be incurred
    # by adding that pair of segments to any biparsing that ends at
    # {(k0,k1)}.
    #
    # If {ks} is even (a gap segment), the field {ivar} of the edge must be
    # {-1}. If {ks} is 0, there must be at most one candidate
    # {(0,0,score,-1)} in {ecand_list} with finite score.
    #
    # If {ks} is odd, {ks=2*ih+1} (a hit segment) {ivar} should be the
    # index of the variant of bipattern {bitemp[ih]} that matched the
    # susbstrings {text0[k0:i0]} and {text1[k1:i1]).
    #
    # The procedure returns the score {score_best} of the extended
    # biparsing and the candidate edge {edge_best} used by it. Also
    # returns the number {nb} of biparsings with finite score that end
    # at {(i0,i1)}.
    #
    # If there is no viable extension, returns {+inf} as the
    # score, {None} as {edge_best}, and 0 as {nb}.

    if debug_vertex: 
      err.write(f"!X segment {ks}, trying {len(ecand_list)} edges to ({i0},{i1})\n")
    assert ks >= 0 and ks <= ns
    if ks == ns:
      assert i0 == nt0 and i1 == nt1
    else:
      assert i0 <= nt0 and i1 <= nt1
    score_best = +inf
    edge_best = None
    nb = 0
    for k0, k1, score_seg, ivar in ecand_list:
      assert isfinite(score_seg) and score_seg >= 0
      if debug_edge: 
        err.write(f"!X   {k0 = } {k1 = } {score_seg = :.8f}")
      if ks == 0:
        assert k0 == 0 and k1 == 0
        k0_old, k1_old, score_old, iv_old, nb_old = ( -1, -1, 0.000, -1, 1)
      else:
        assert k0 <= i0 and k1 <= i1
        k0_old, k1_old, score_old, iv_old, nb_old = tb[ks-1][(k0,k1)]
      assert isfinite(score_old)
      assert nb_old >= 1
      score_new = score_old + score_seg
      assert isfinite(score_new)
      # !!! Should clip at {max_score}
      if debug_edge: 
        err.write(f" {score_old = :.8f} {score_new = :.8f}\n")
      nb += nb_old
      if score_new < score_best:
        edge_best = (k0, k1, score_seg, ivar)
        score_best = score_new
    return score_best, edge_best, nb
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  hcand_lists = btf.get_hit_candidate_lists(bitemp, text0, text1) 
  if hcand_lists is None: return None, +inf

  def get_candidate_gap_vertices(ks):
    # Given the even index {ks} in {0..ns-1} of a gap segment, computes a
    # list of vertex indices {(i0,i1)} on plane {ks} that may have at
    # least one valid input edge.
    # 
    # These are all pairs of possible hit starts in the two sides; that
    # is, all pairs {(i0,i1)} such that {i0>=k0_min}, {i1>= k1_min}, and
    # there is a candidate hit in {hcands[ih]} of the form
    # {(k0c,i0c,k1c,i1c, ...)} with {i0=k0c} and {i1=k1c}; where {ih} is
    # the index of the hit following the gap, and {k0_min}, {k1_min}
    # are the min indices of any vertex in {tb[ks]}. However, if {ks==ns-1}
    # (meaning {ih==nh}), the only possible end is {(nt0,nt1)}.

    ig = ks//2 # Index of gap in {0..ng-1}
    assert 0 <= ig and ig < ng
    ih = ks//2 # Index of next rung, or {nh}.
    assert 0 <= ih and ih <= nh
    
    # Get the min indices on plane {ks-1}:
    if ks == 0:
      k0_min = 0; k1_min = 0
    else:
      if debug_plane: err.write(f"!v plane {ks-1} has {len(tb[ks-1])} vertices\n")
      k0_min = min(k0 for (k0,k1) in tb[ks-1].keys())
      k1_min = min(k1 for (k0,k1) in tb[ks-1].keys())

    # Get the relevant values {i0_list} and {i1_list} of {i0} and {i1}:
    i0_list = None; i1_list = None
    if ih >= nh:
      assert ks == ns-1
      i0_i1_list = [ (nt0,nt1,) ]
    else:
      i0_i1_list = [ \
        (k0c,k1c,) \
        for k0c, i0c, k1c,i1c, pena, ivc in hcand_lists[ih] \
          if k0c >= k0_min and k1c >= k1_min \
      ]
      # Remove duplicate vertices:
      i0_i1_list = sorted(set(i0_i1_list))
    return i0_i1_list
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def get_candidate_hit_vertices(ks):      
    # Given the odd index {ks} in {1..ns-2} of a hit segment, computes a list of
    # vertex indices {(i0,i1)} on plane {ks} that may have at least one
    # valid input edge.
    # 
    # The possible ends are all {(i0,i1)} such that there is a hit candidate
    # {(k0,i0,k1,i1,pena,ivar)} in {hcands[ih]} where {(k0,k1)} is a vertex
    # of plane {ks-1} with finite score.
    # 
    assert 1 <= ks and ks <= ns-2
    ih = (ks-1)//2 # Index of relevant rung.
    i0_i1_list = [ ]
    assert 0 <= ih and ih < nh
    for k0, i0, k1, i1, pena, ivar in hcand_lists[ih]:
      assert 0 <= k0 and k0 < i0 and i0 <= nt0
      assert 0 <= k1 and k1 < i1 and i1 <= nt1
      if (k0,k1) in tb[ks-1]:
        i0_i1_list.append((i0,i1,))
            
    # Remove duplicate vertices:
    i0_i1_list = sorted(list(set(i0_i1_list)))
    return i0_i1_list
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def get_candidate_gap_edges(ks, i0, i1):
    # Given the even index {ks} in {0..ns-1} of a gap segment,  
    # and the indices {(i0,i1)} of a vertex on plane {ks},
    # gets all the the viable candidate gap edges that end
    # at that vertex.
    #
    # If {ks < ns-1} (meaning {ig < ng-1}), assumes that {i0} is a
    # possible start of a hit of {bitemp[ih]} on side 0 and ditto
    # for{i1} on side 1, where {ih} is the index of the next hit (that
    # is, {ih = ig}). If {ks == ns-1} (meaning {ig == ng-1}, then
    # {(i0,i1)} must be {(nt0,nt1)}.
    # 
    # The list will include every triple {(k0,k1,score,-1,)} such that
    # {(k0,k1)} a vertex of plane {ks-1}, {k0} is in {0..i0}, {k1}
    # is in {0..i1}, and {score} is finite; where {score} is the the gap
    # score computed from the gap lengths {i0-k0} and {i1-k1}
    #
    # If {ks} is zero the plane {ks-1} is assumed to have only the
    # pair {(0,0)}.

    assert (ks % 2) == 0
    ig = ks//2
    
    k0_k1_list = iter(((0,0,),)) if ks == 0 else tb[ks-1].keys()
    
    ecand_list = []
    for k0, k1 in k0_k1_list:
      if ks == 0:
        assert k0 == 0 and k1 == 0
      else:
        assert k0 <= nt0
        assert k1 <= nt1
      if k0 <= i0 and k1 <= i1:
        score = gap_score(ks, k0,i0, k1,i1)
        if isfinite(score):
          ecand_list.append((k0,k1,score,-1,))
    return ecand_list
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  def get_candidate_hit_edges(ks, i0, i1):      
    # Given the odd index {ks} in {1..ns-2} of a hit segment,  
    # and the indices {(i0,i1)} of a vertex on plane {ks}, gets all
    # the the viable candidate hit edges that end at vertex {(i0,i1)} of
    # plane {ks}.
    #
    # Assumes that {i0} is a possible end of a hit of {bitemp[ih]} on
    # side 0 and ditto for{i1} on side 1. Assumes that hits are never
    # empty.
    # 
    # The list will include every triple {(k0,k1,pena,ivar)} such that
    # there is a hit candidate {(k0,i0,k1,i1,pena,ivar)} in the list
    # {hcands[ih]}, and vertex {(k0,k1)} from plane {ks-1} has finite
    # score. There shoud not be two entries in the {hcands[ih]} list
    # with the same {k0,i0,k1,i1}.

    assert 1 <= ks <= ns-2
    ih = (ks-1)//2
    ecand_list = [ ]
    for k0c, i0c, k1c, i1c, penac, ivc in hcand_lists[ih]:
      if debug_plane:
        err.write(f"!o   cand hit edge ( {k0c:3d} {i0c:3d}  {k1c:3d} {i1c:3d} {penac:6.3f} {ivc:2d} )\n")
      assert 0 <= k0c and k0c < i0c and i0c <= nt0
      assert 0 <= k1c and k1c < i1c and i1c <= nt1
      if i0c == i0 and i1c == i1 and (k0c,k1c) in tb[ks-1]:
        ecand_list.append((k0c,k1c,penac,ivc))
    return ecand_list
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    
  def print_candidate_hits(ks):
    assert (ks % 2) == 1
    ih = ks//2
    nch = len(hcand_lists[ih])
    err.write(f"there are {nch} candidate hits for hit {ih} segment {ks}:\n")
    for k0, i0, k1, i1, pena, ivar in hcand_lists[ih]:
      err.write(f"{k0:3d}:{i0:<3d} {k1:3d}:{i1:<3d}")
      err.write(f" {pena:5.3f}")
      err.write(f" {ivar:3d}\n")
    return
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  # Fill the tableau with partial biparsings:
  for ks in range(ns):
    if debug_plane: err.write(f"!o filling tableau plane {ks} ...\n")
    # Get a list of {i0_i1_list} of all vertices in plane {ks} that may have an input edge:
    i0_i1_list = None
    if (ks % 2) == 0:
      # Edges from plane {ks-1} to {ks} are gap edges:
      i0_i1_list = get_candidate_gap_vertices(ks)
    else:
      # Edges from plane {ks-1} to {ks} are hit edges:
      if debug_hcands: print_candidate_hits(ks)
      i0_i1_list = get_candidate_hit_vertices(ks)
      
    if len(i0_i1_list) == 0: return None, +inf

    # Now compute the entry {[ks][(i0,i1)]} of the tableau for all vertces {(i0,i1)}
    # of plane {ks} that may have an input edge:
    for i0,i1 in i0_i1_list:

      if debug_vertex:
        err.write(f"!~  looking for best hit extension to ({i0},{i1})\n")

      # Collect the edges [ecand_list} into {(i0,i1)}:
      ecand_list = None
      if (ks % 2) == 0:
        ecand_list = get_candidate_gap_edges(ks, i0, i1)
      else:
        ecand_list = get_candidate_hit_edges(ks, i0, i1)

      if debug_vertex:
        err.write(f"!~    got {len(ecand_list)} edges\n")

      score_best, edge_best, nb = find_best_extension(ks, i0,i1, ecand_list)
      if debug_vertex:
        err.write(f"!~    {score_best = :.8f} {edge_best = !r} {nb = }\n")
      if isfinite(score_best):
        assert score_best >= 0 # !!! Should check {max_score}.
        k0_best, k1_best, score_seg_best, iv_best = edge_best
        assert 0 <= k0_best and k0_best <= i0
        assert 0 <= k1_best and k1_best <= i1
        assert nb >= 1
        tb[ks][(i0,i1)] = ( k0_best, k1_best, score_best, iv_best, nb)
        tot['finel'] += 1

      else:
        assert edge_best is None
        assert nb == 0

    if debug_plane:
      err.write(f"!o tableau plane {ks} has {len(tb[ks])} finite elements.\n")
      print_tableau_plane(ks)

  if debug_tableau:
    dump_tableau()

  if len(tb[ns-1]) == 0:
    # No bimatch possible:
    bimatch_fin = None
    score_fin = +inf
  else:
    assert(len(tb[ns-1])) == 1
    bimatch_fin, score_fin, nb_fin = extract_bimatch()
    err.write(f"!~ FINAL {score_fin = :.8f}\n")
    assert isfinite(score_fin)
    bmf.check_bimatch(bimatch_fin, bitemp, text0, text1)

  if verbose:
    err.write(f"compared lengths of {tot['evgap']} gap pairs.\n")
    err.write(f"considered {tot['evhit']} possible hit pairs ({tot['nohit']} failed).\n")
    err.write(f"total {tot['finel']} finite entries in tableau (incl sink).\n")
    err.write(f"potentially {nb_fin} distinct bimatchings.\n")
    err.write(f"final score {score_fin:.8f}.\n")

  if debug_top: err.write("!~ exit match_bitemplate\n")

  return bimatch_fin, score_fin
  # ----------------------------------------------------------------------

def plot_cost(eval_one_gap):
  size_scale = 5.000
  iscale = int(floor(size_scale*10000 + 0.5))
  fname = f"/tmp/.cost-py-{iscale:06d}.txt"
  err.write(f"plotting the gap cost to file \"{fname}\"\n");
  wr = open(fname, "w")
  err.write(f"size_scale = {size_scale:.6f}\n")
  H = 2; N = 2*H+1
  gsz1 = [ None ]*N
  for kr in range(N):
    r = pow(2.0, kr-H)
    gsz1[kr] = int(r*50*size_scale + 0.5);
    err.write(f"gsz1[{kr:3d}] = {gsz1[kr]:3d}\n")

  for gsz0 in range(100):
    wr.write(f"{gsz0:3d} ");
    for kr in range(N):
      cost = eval_one_gap(gsz0, gsz1[kr], 1, 3)
      icost = int(floor(1000*cost + 0.5))
      wr.write(f" {icost:6d}")

    wr.write(f"\n");

  wr.close();
  return
  # ----------------------------------------------------------------------

def match_test_func(text0, text1, bitemp, un1_per_un0, debug_level):

  def eval_one_gap(gsize0, gsize1, ig, ng):
    debug = False
    if debug: err.write("!e {gsize0 = !r} {gsize1 = !r}\n")
    gap_score = bef.compute_single_gap_score(gsize0, gsize1, un1_per_un0, ig, ng)
    if debug: err.write(f"!e {gap_score = :.3f}\n")
    return gap_score
    # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  bimatch, score = match_bitemplate(text0, text1, bitemp, eval_one_gap, debug_level)
  return bimatch, score
  # ----------------------------------------------------------------------

if len(sys.argv) > 1 and sys.argv[1] == "BTMPF.TEST":
  bmtestf.all_tests(match_test_func, 20)