def compute_gaps_score(gsizes0, gsizes1, un1_per_un0):
  # Computes the gap sizes component of the badness score of a bimatching
  # given the actual sizes {gsizes0[0..ng-1],gsizes1[0..ng-1]}
  # of the gaps determined by the bimatching.
  debug = False
 
  ng = len(gsizes0); assert len(gsizes1) == ng
  gaps_score = 0
  for ig in range(ng):
    gsci = compute_single_gap_score(gsizes0[ig], gsizes1[ig], un1_per_un0, ig, ng)
    gaps_score += gsci
    if debug: err.write(f"!G gap score = {gsci:.6f} accum = {gaps_score:.6f}\n")
  if debug: err.write(f"!G gap score final = {gaps_score:.6f}\n")
  return gaps_score
  # ----------------------------------------------------------------------
   
def compute_single_gap_score(gsize0, gsize1, un1_per_un0, ig, ng):
  # Computes the contribution to the the badness score of a bimatching
  # due to a single gap, whose size is {gsize0} on a text of some type {utype0} and
  # the size {gsize1} of a supposedly matchng text of some type {uptype1},
  # given the number {un1_per_un0} of units of {utype1} per unit of {utype0}.
  #
  # The contribution is properly weighted based on the index {ig}
  # of the gap and the number of {ng} of gaps.
 
  debug = False
 
  if debug: err.write(f"\n!g {gsize0 = :3d} {gsize1 = :3d} {ig = :2d} {ng = :2d}")
  exp_gsize1 = expected_size1_from_size0(gsize0, un1_per_un0)
  wt = gap_score_weight(ig, ng)
  score = wt * compute_single_size_score(gsize1, exp_gsize1)
  if debug: err.write(f" {wt = :.6f} {score = :7.3f}\n")
  return score
  # ----------------------------------------------------------------------
def gap_score_weight(ig, ng):
  # Weight of gap with index {ig} in {ng} gaps (including
  # those before and after all hits).
  wt = 0.25 if ig == 0 else 0.50 if ig == ng-1 else 1.00
  wt /= ng
  return wt
  # ----------------------------------------------------------------------
def compute_single_size_score(sz, esz):
  # Mismatch score between a single size {sz} and the expected size range {esz}.
  debug = False
  if debug: err.write(f"!s {sz = :4.1f} esz = {esz[0]:3d}..{esz[1]:3d}")
 
  frac = compute_frac_error(sz, esz)
  score = 0.20 * frac**2
  if debug: err.write(f" {frac = :.6f} {score = :.6f}\n")
  return score
  # ----------------------------------------------------------------------
def expected_size1_from_size0(size0, un1_per_un0):
  # Computes the expected size {exp_size1} of some text of some type
  # {utype1} given the size of {size0} of the supposedly corresponding
  # text of type {utype0} and the number {un1_per_un0} of units of
  # {utype1} that correspond to one unit of {utype0} The result is a
  # narrow range {lo,hi}.
  eps = 0.0001
  exp_size1 = spf.scale_size_fuzzy(size0, un1_per_un0, eps)
  return exp_size1
  # ----------------------------------------------------------------------
def compute_dev_of_log(es):
  # Assumes that {es} is the expected value of a size variable with a
  # log-normal distribution with deviation proportional to {sqrt(es)}.
  # Returns the approximate deviation of the distributon of the log of
  # the size.
 
  debug = False
  if debug: err.write(f"!@ {es = :3d}")
   
  eps = 1
  es0 = hypot(es, eps) # In case {es} is zero.
  es0_log = log(es0)
  if debug: err.write(f" {es0 = :5.1f} {es0_log = :+7.3f}")
 
  var_abs = 0.10 * es0 # Rough guess - to be improved.
  dev_abs = sqrt(var_abs)
  if debug: err.write(f" {dev_abs = :7.3f}")
 
  es1 = hypot(es0 + dev_abs, eps)
  es1_log = log(es1)
  if debug: err.write(f" {es1 = :5.1f} {es1_log = :+7.3f}")
  dev_log = es1_log - es0_log
  if debug: err.write(f" {dev_log = :+7.3f}\n")
  return dev_log
  # ----------------------------------------------------------------------
def compute_frac_error(s, es):
  # Fractional error between a size {s} and the expected size range {es},
  # relative to the estimated deviation of the latter, in log scale.
  debug = False
  if debug: err.write(f"!~ {s = :3d} es = {es[0]:3d}..{es[1]:3d}")
  assert len(es) == 2 and es[0] <= es[1]
  avg = es[0] if s < es[0] else es[1] if s > es[1] else s
  if debug: err.write(f" {avg = :3d}")
  if s == avg:
    frac = 0
  else:
    eps = 1
    s0 = hypot(s, eps)
    s0_log = log(s0)
    if debug: err.write(f" {s0 = :5.1f} {s0_log = :+7.3f}")
    avg0 = hypot(avg, eps)
    avg0_log = log(avg0)
    if debug: err.write(f" {avg0 = :5.1f} {avg0_log = :+7.3f}")
    dev_log = compute_dev_of_log(avg)
    if debug: err.write(f" {dev_log = :7.3f}")
    frac = (s0_log - avg0_log)/dev_log
  if debug: err.write(f" {frac = :7.3f}\n")
  return frac
  # ----------------------------------------------------------------------
def compute_rorre_carf(frac, es):
  # Inverse of fractional error
  eps = 1
  if frac < 0:
    dev_log = compute_dev_of_log(es[0])
    y = frac*dev_log + log(hypot(es[0], eps))
    z = exp(y)
    s = 0 if z <= eps else sqrt(z*z - eps*eps)
  elif frac > 0:
    dev_log = compute_dev_of_log(es[1])
    y = frac*dev_log + log(hypot(es[1], eps))
    z = exp(y)
    s = 0 if z <= eps else sqrt(z*z - eps*eps)
  else:
    # Geometric mean of the range end values:
    s = sqrt(hypot(es[0], eps)*hypot(es[1], eps))
    assert s >= eps
    s = sqrt(s*s - eps*eps)
  return s
  # ----------------------------------------------------------------------