#! /usr/bin/python3
# Last edited on 2026-02-16 11:21:50 by stolfi

# Functions to simulate the delete-shift-insert process
# for mutation of value lists.

from sys import stderr as err
from math import sqrt, sin, cos, log, exp, floor, pi, inf
import random
     
def simulate_del_shift_add(LA, nK, nB, max_val, min_sep, pert_dev):
  # Simulates the effect of the delete-shif-insert process on the list {LA}.
  #
  # The list {LA} must have zero or more increasing float values in the range
  # {M = [0 _ max_val]}, with {min_sep} minimum difference between them.
  # 
  # The process consists in deleting a random susbet of the elements
  # of {LA}, leaving only {nK} of them; applying to each surviving element
  # a normal perturbation with deviation {pert_dev};  shifting the 
  # perturbed values by a random amount that will keep all of them inside
  # {M}; then inserting {nB-nK} values before, between and after those elements,
  # all within {M}, so that the resulting list {LB} still has min separation {min_sep}.
  #
  # Returns the list {KA} of the {nK} elements of {LA} that were
  # retained, the list {LB} od new values, the list {KB} of the indices
  # into {LB} of the element that were preserved.

  nA = len(LA)
  
  # Check randge and min separation, just incase:
  if nA > 0:
    assert LA[0] >= 0
    for iA in range(1,nA): assert LA[iA] >= LA[iA-1] + min_sep
    assert LA[nA-1] <= max_val
  
  # Choose the values in {LA} that will survive:
  assert nK >= 0 and nK <= nA
  KA = choose_surviving_values(nA, nK)
  assert len(KA) == nK
  
  # Create the list {LK} of values that survived, perturbed:
  LK = [ LA[k] for k in KA ]
  assert len(LK) == nK
  
  LP = perturb_values_monotonically(LK, max_val, min_sep, pert_dev)
  
  # Apply a random shift {sh}:
  min_sh = 0 - LP[0]
  max_sh = max_val - LP[nK-1]
  sh = random.uniform(min_sh, max_sh)
  LP = [ LP[iK] + sh for iK in range(nK) ]
  
  # Insert the new values, monotonically:
  LB, KB = insert_random_values(nB, max_val, min_sep, LP)
  return KA, LB, KB
  # ----------------------------------------------------------------------

def choose_surviving_values(nA, nK):
  # Chooses a random subset {KA} of {0..nA} with exactly {nK}
  # elements.  Surely there is a faster way but who cares;
  assert 0 <= nK and nK <= nA
  
  # Enumerate all subsets of size {nK}, save in {C}
  JA_bins = [] # Subsets packed in binary as ints.
  for JA_bin in range(2**nA):
    m = 0; 
    for i in range(nA): 
      if (JA_bin & (1 << i)) != 0: m += 1
    if m == nK: JA_bins.append(JA_bin)
  
  # Pick a susbset at random:
  KA_bin = JA_bins[random.randint(0,len(JA_bins)-1)]
  
  # Unpack it:
  KAA = [ i for i in range(nA) if (KA_bin & (1 << i)) != 0 ]
  assert len(KAA) == nK
  return KAA
  # ----------------------------------------------------------------------
  
def perturb_values_monotonically(LK, max_val, min_sep, pert_dev):
  # Applies random perturbations to the elements of {LK}.
  # The elements of {LK} must be in increasin order, in the range {[0 _ max_val]},
  # with min separation {min-sep} between them.
  # The perturbations are Gaussian with zero mean and deviation {pert-dev}.
  # Retries the perturbation until the perturbed elements 
  # satisfy the same conditions.
  
  nK = len(LK)
  trials = 0
  # Fudge if needed to keep the perturbed values sorted and less than {max_val}:
  ok = False
  LP = [ None, ]*nK
  while not ok:
    trials += 1; assert trials <= 5, "failed to perturb monotonically"
    ok = True
    for iK in range(nK):
      LP[iK] = LK[iK] + random.gauss(0, pert_dev)
      if iK > 0 and LP[iK] - LP[iK-1] < min_sep: ok = False; break
    if LP[nK-1] - LP[0] > max_val: ok = False
  return LP
  # ----------------------------------------------------------------------
 
def insert_random_values(nB, max_val, min_sep, LP):
  # Creates a sorted list {LB} of {nB} values in {[0_max_val]}
  # that includes the values in the list {LP} and keeps 
  # the minimum distance {min_sep} between all points.
  # Returns {LB} and the list of indices in {LB} of the copied elements.
  
  nP = len(LP)

  debug  = True
  if debug: err.write(f"~~~~ insert_random_values {nB = } {nP = } {max_val = } {min_sep = } ~~~~\n")
  if debug: err.write(f"~~   {LP = }\n")
  
  # Choose the list {LS} of values to be inserted, compressed:
  nS = nB - nP
  max_val_free = max_val - (nB-1)*min_sep
  assert max_val_free > 0, "the {nB} is too big for {max_val}"
  LS = [ random.uniform(0, max_val_free) for kB in range(nS) ]
  LS.sort()
  if debug: err.write(f"~~ {LS = }\n")
  
  # Insert the values of {LS} between those of {LP}, with proper spacing:
  LB = [ None ]*nB
  KB = []
  iP = 0
  iS = 0
  used_valS = 0 # How much of {LS} we used so far.
  for iB in range(nB):
    min_valB = 0 if iB == 0 else LB[iB-1] + min_sep
    if debug:
      # Sanity check:
      err.write(f"~~  {iB = } {min_valB = } {used_valS = }\n")
      sum_used = 0
      for jB in range(0,iB):
        sum_used += LB[jB] - (0 if jB == 0 else LB[jB-1] + min_sep)
      assert abs(sum_used - used_valS) < 0.001

    # Next elem of {LB}, if we insert from {LS}:
    next_valS = +inf if iS >= nS else min_valB + (LS[iS] - used_valS)
    assert next_valS >= min_valB - 0.000001*min_sep
    
    # Next elem of {LB}, if we copy from {LP}:
    next_valP = +inf if iP >= nP else LP[iP]
    assert next_valP >= min_valB - 0.000001*min_sep
    
    if debug: err.write(f"~~    {next_valS = } {next_valP = } \n")
    
    if next_valS < next_valP - min_sep:
      # Insert from {LS}:
      LB[iB] = next_valS; iS += 1
      # Consume it:
    else:
      # Copy from {LP}:
      LB[iB] = next_valP; iP += 1; KB.append(iB)
      # That will consume some of {SS}:
    used_valS += LB[iB] - min_valB
    if debug: err.write(f"~~    {LB[iB] = } {used_valS = } \n")
    assert LB[iB] >= min_valB - 0.000001*min_sep
  assert LB[nB-1] <= max_val, f"overflow {LB[nB-1] = :24.16e} {max_val = :24.16e}"
  assert iS >= nS and iP >= nP
  assert len(KB) == nP
  return LB, KB
  # ----------------------------------------------------------------------