SL2-CTL KPS LENS/VISOR TRACKER

tic_tac_sl2_ctl_1.py

# tic_tac_sl2_ctl_1.py, September 16 2025

#_____________________________________________________________________________________
# SL2-CTL Gauge(KPS Track Lens) Notes:

# 1. Ensure coords are in meters if is physical wave-
#    lengths -> key depth values --> non-negative disc

# 2. For ray-optical corrections, consider local slope
#    deflection angle -> transverse shift @ detect dist
	
# 3. Concerning speed @ tic-tac, possible to replace
#    disp_scale with propagation distance: Physical
#    fidelity, calculate a local surface norm from a
#    @height surface; and then, translate the normal
#    tilt to the prime deflection angle°

# 4. Also possible to match wave optics rather than
#    ray deflection...consider applying complex field
#    multiplication(s), example exp(i*phase)) or else
#    and calculating PSFs instead of shifting KDE Cx

# 5. Synthetic sampling can be improved by joining
#    probabilistic kde mass layering to the rotation
#    cycles as height dispersion reduction; requires
#    howbeit more broadened covar weight structs
#_____________________________________________________________________________________
# Zechariah 14:9

from typing import List, Tuple
import random
import time
import math
import json
import os
#_____________________________________________________________________________________

def _tt_mat_mult_vec(m: tuple, v: tuple) -> tuple:
    
    return (m[0][0]*v[0]+m[0][1]*v[1], m[1][0]*v[0]+m[1][1]*v[1])
#_____________________________________________________________________________________

def _tt_rotation_matrix(th: float) -> tuple:
    
    c, s = math.cos(th), math.sin(th)
    return ((c,-s),(s,c))
#_____________________________________________________________________________________

def _tt_weighted_centroid(p: list, wt: list) -> tuple:
    
    ttl = sum(wt)
    if ttl == 0.0:
        raise Exception('<SL2-CTL ERR> a total weight sum=zero for centroid impossible')
    return (sum(z.real*w for z, w in zip(p, wt))/ttl, sum(z.imag*w for z, w in zip(p, wt))/ttl)
#_____________________________________________________________________________________

def _tt_weighted_mean_weight(wt: float) -> float:
    
    if wt:
        return sum(wt)/float(len(wt))
    return 0.0
#_____________________________________________________________________________________

def _tt_weighted_covariance(p: complex, wt: float, C: tuple) -> tuple:
    
    ttl = sum(wt)
    if ttl == 0.0:
        raise Exception(f'<SL2-CTL ERR> a total weight sum=zero for covariance impossible')
    cx, cy = C
    sxx = syy = sxy = 0.0
    for z, w in zip(p, wt):
        dx = z.real-cx
        dy = z.imag-cy
        sxx+=w*dx*dx
        syy+=w*dy*dy
        sxy+=w*dx*dy
    return ((sxx/ttl, sxy/ttl),(sxy/ttl, syy/ttl))
#_____________________________________________________________________________________

def _tt_eig2_symmetric(m: list) -> tuple:
    
    a, b, c = m[0][0], m[0][1], m[1][1]
    tr, d, = a+c, a*c-b*b
    disc = tr*tr/4.0-d
    if disc < 0.0 and disc > -1e-12: disc = 0.0
    if disc < 0.0:
        raise Exception(f'<SL2-CTL ERR> negative disc found; no root disc halt')
    sd = math.sqrt(disc)
    l1, l2 = tr/2.0+sd, tr/2.0-sd
    def vec(lam):
        vx, vy = b, lam-a
        if abs(vx) < 1e-12 and abs(vy) < 1e-12:
            return (1.0, 0.0) if a >= c else (0.0, 1.0)
        n = math.hypot(vx, vy)
        return (vx/n, vy/n)
    return (l1, vec(l1)), (l2, vec(l2))
#_____________________________________________________________________________________

def _tt_apply_transform_point_xy(p_xy: tuple, s: float, r: tuple, t: tuple) -> tuple:
    
    v = _tt_mat_mult_vec(r, p_xy)
    return (v[0]*s+t[0], v[1]*s+t[1])
#_____________________________________________________________________________________

def _tt_gaussian_kernel(d: float, h: float) -> float:
    
    return math.exp(-0.5*d/(h*h))/(2.0*math.pi*h*h)
#_____________________________________________________________________________________

def _tt_kde_mass_at_point(xy: tuple, c_xy: tuple, wt: tuple, h: float) -> float:
    
    ws = sum(wt)
    if ws == 0.0:
        return 0.0
    s=0.0
    x, y = xy
    for (cx, cy), w in zip(c_xy, wt):
        dx, dy = x-cx, y-cy
        s+=w*_tt_gaussian_kernel(dx*dx+dy*dy, h)
    return s/ws
#_____________________________________________________________________________________

def _tt_lens_key_height_map_xy(xy, k_a=0.0, k_o=0.0, k_w=0.5, k_d=0.02, b_amp=0.0, b_sigma=0.5) -> float:
    
    x, y = xy
    c, s = math.cos(k_a), math.sin(k_a)
    crd = x*c+y*s
    if k_w <= 0:
        # ..key width projection is otherwise normalize scaling
        sp = 0.5*(1.0 if crd > k_o else -1.0)
    else:
        # Pilot's visor feedback same, tanh maps to (-0.5, 0.5)
        # when is mapped/logged scaled --> center to key offset
        sp = 0.5*math.tanh((crd-k_o)*(2.0/k_w))
    k_h, r = k_d*sp, x*x+y*y
    b_rtr = b_amp*math.exp(-0.5*r/(b_sigma*b_sigma)) if b_sigma > 0 else 0.0
    return k_h+b_rtr
#_____________________________________________________________________________________

def _tt_calculate_phase_from_height(h: float, wl: float) -> float:
    
    return 2.0*math.pi*h/wl
#_____________________________________________________________________________________

def _tt_phase_to_transverse_shift_at_xy(xy: tuple, h_func, wl: float, d=1e-4, dsp_scale=1.0) -> tuple:
    
    # Calculate tiny transverse displace vector, being proportional @ phase gradient;
    # --> will return in same xy units, meters @ height func, delta and displacement
    x, y = xy
    hxp, hxm, hyp, hym = h_func((x+d, y)), h_func((x-d, y)), h_func((x, y+d)), h_func((x, y-d))
    # Negative signals; local phase gradient deflects rays reverse of phase increase
    dx = -dsp_scale*(_tt_calculate_phase_from_height(hxp, wl)-_tt_calculate_phase_from_height(hxm, wl)/(2.0*d))
    dy = -dsp_scale*(_tt_calculate_phase_from_height(hyp, wl)-_tt_calculate_phase_from_height(hym, wl)/(2.0*d))
    return (dx, dy)
#_____________________________________________________________________________________

def _tt_apply_lens_key_to_centers(c_xy: tuple, wl=550e-9, k_a=0.0, k_o=0.0, k_w=0.5, k_d=1e-5, b_amp=0.0, b_sigma=0.5, d=1e-5, dsp_scale=1e-9) -> list:

    # XY units must be consistent here, meters @wavelength/height: tune
    # --> to converted phase gradient --> transverse & then shift @mass
    def h_fn(pt):
        return _tt_lens_key_height_map_xy(pt, k_a, k_o, k_w, k_d, b_amp, b_sigma)
    return [_tt_phase_to_transverse_shift_at_xy(p, h_fn, wl, d, dsp_scale) for p in c_xy]
#_____________________________________________________________________________________

def _tt_build_gauss_samples(cx: float, cy: float, sgm: float, b_span: float, e_uni: int, e_eml: int, e_dir: int, p_xy: list, wt: list, prnc_axis: list, rng: object) -> list:
    
    samples = [(rng.uniform(cx-b_span, cx+b_span), rng.uniform(cy-b_span, cy+b_span)) for _ in range(e_uni)]
    cm, ttl_w, s = [], sum(wt), 0.0
    for w in wt:
        s+=w
        cm.append(s)
    for _ in range(e_eml):
        r, idx = rng.uniform(0, ttl_w), 0
        while idx + 1 < len(cm) and cm[idx] < r: idx+=1
        cxp, cyp = p_xy[idx]
        jtr = 0.2*sgm if sgm > 0 else 0.05
        samples.append((cxp+rng.gauss(0.0, jtr), cyp+rng.gauss(0.0, jtr)))
    # Set directional samples along all principal axes:
    radii = [0.25*sgm, 0.5*sgm, 1.0*sgm, 1.5*sgm]
    for a in prnc_axis:
        ux, uy = a
        for r in radii:
            samples.append((cx+ux*r, cy+uy*r))
            samples.append((cx-ux*r, cy-uy*r))
    for _ in range(e_dir):
        a, r = rng.uniform(0, 2*math.pi), rng.uniform(0.0, b_span)
        samples.append((cx+r*math.cos(a), cy+r*math.sin(a)))
    return samples
#_____________________________________________________________________________________

def _tt_fit_rigid_scale_basis(pt: list, wt: list, tc: tuple, t_axis=None) -> tuple:
    
    src_ctrd = _tt_weighted_centroid(pt, wt)
    (lam1, v1), (lam2, v2) = _tt_eig2_symmetric(_tt_weighted_covariance(pt, wt, src_ctrd))
    src_axis = v1
    if t_axis is None: a = 0.0
    else:
        tx, ty = t_axis
        t_nrm = math.hypot(tx, ty)
        if t_nrm == 0.0:
            raise Exception(f'<SL2-CTL ERR> zero target axis error, fit rigid scale halted')
        tu = (tx/t_nrm, ty/t_nrm)
        a = math.atan2(src_axis[0]*tu[1]-src_axis[1]*tu[0], src_axis[0]*tu[0]+src_axis[1]*tu[1])
    r = _tt_rotation_matrix(a)
    t_ctrd = _tt_apply_transform_point_xy(src_ctrd, 1.0, r, (0.0,0.0))
    return (1.0, a, r, (tc[0]-t_ctrd[0], tc[1]-t_ctrd[1]), src_ctrd, (lam1, lam2), src_axis)
#_____________________________________________________________________________________

def _tt_l2_kde_error_samples(cntr_A: tuple, cntr_B: tuple, wt: list, h: float, samples: list) -> float:
    
    es = []
    for xy in samples:
        d = _tt_kde_mass_at_point(xy, cntr_A, wt, h)-_tt_kde_mass_at_point(xy, cntr_B, wt, h)
        es.append(d*d)
    return sum(es)/max(1, len(es))
#_____________________________________________________________________________________
#
#                   ////////// SL2-CTL SIM TRACK RELAY \\\\\\\\\\
#_____________________________________________________________________________________
#.....................................................................................
#.....................................................................................
#.....................................................................................

def _tt_evaluate_scales_sim(points, weights, target_centroid, target_axis,
                            device_scales, iterations, alpha,
                            eval_uniform, eval_empirical, eval_directional,
                            batches, seed, progress_interval=10, out_tempfile=None,
                            lens_enabled=True, wavelength=550e-9,
                            key_angle=0.0, key_offset=0.0, key_width=0.2, key_depth=2e-6,
                            bump_amp=1e-6, bump_sigma=0.5, delta=1e-6, disp_scale=5e-9) -> dict:
    rng, n, pts_xy = random.Random(seed), len(points), [(z.real, z.imag) for z in points]
    s0_base, angle, R, t, src_centroid, lambdas, src_axis = _tt_fit_rigid_scale_basis(points, weights, target_centroid, target_axis)
    sigma = math.sqrt(max(0.0, lambdas[0]+lambdas[1]))
    h, bbox_span = 1.06*sigma*(n**(-0.2)) if sigma > 0 else 0.1, max(0.5, 2.0*sigma)
    mean_weight, total_tasks, task_count, best = _tt_weighted_mean_weight(weights), len(device_scales)*iterations, 0, None
    start_time = time.time()
    for d_idx, s0 in enumerate(device_scales):
        for k in range(iterations):
            task_count+=1
            decay = max(0.0, 1.0-alpha*k)
            s_k = s0*decay
            transformed_centers = [_tt_apply_transform_point_xy(pt, s_k, R, t) for pt in pts_xy]
            # *Applies lens-key perturbation only to the transformed centers:
            if lens_enabled:
                transformed_centers = _tt_apply_lens_key_to_centers(transformed_centers, wavelength,
                                      key_angle, key_offset, key_width, key_depth, bump_amp,
                                      bump_sigma if bump_sigma is not None else sigma, delta, disp_scale)
            batch_scores = []
            for b in range(batches):
                bs = seed+d_idx*1000+k*100+b
                brng = random.Random(bs)
                px, py = src_axis
                principal_axes = [(px, py), (-py, px)]
                samples = _tt_build_gauss_samples(src_centroid[0], src_centroid[1], sigma,
                bbox_span, eval_uniform, eval_empirical, eval_directional, pts_xy,
                weights, principal_axes, brng)
                err = _tt_l2_kde_error_samples(transformed_centers, pts_xy, weights, h, samples)
                batch_scores.append(err)
            mean_err = sum(batch_scores)/len(batch_scores)
            std_err = math.sqrt(sum((x-mean_err)**2 for x in batch_scores)/max(1, len(batch_scores)))
            score = mean_err
            if best is None or score < best['score']:
                best = {'score': score, 'device_index': d_idx, 'iter': k, 's_k': s_k, 'R_angle': angle, 't': t,
                        'h': h, 'mean_err': mean_err, 'std_err': std_err,
                        'lens': {'enabled': lens_enabled, 'wavelength': wavelength, 'key_angle': key_angle,
                                 'key_offset': key_offset, 'key_width': key_width, 'key_depth': key_depth,
                                 'bump_amp': bump_amp, 'bump_sigma': bump_sigma, 'delta': delta,
                                 'disp_scale': disp_scale}}
                print('[NEW BEST] device={} k={} s_k={:.5f} mean_err={:.6e} std_err={:.6e}'.format(
                      d_idx, k, s_k, mean_err, std_err))
            if task_count%progress_interval == 0 or task_count == total_tasks:
                elapsed = time.time()-start_time
                pct = 100.0*task_count/total_tasks
                span_ratio = bbox_span/(mean_weight if mean_weight != 0.0 else 1.0)
                print('[PROGRESS] {:.1f}% ({}/{}) pan={:.1f}s c_d={} k={} s_k={:.5f} s_r={:.4f} viz={:.6e}'.format(
                      pct, task_count, total_tasks, elapsed, d_idx, k, s_k, span_ratio, best['score']))
    if out_tempfile:
        with open(out_tempfile, 'w') as f:
            f.write(json.dumps(best))
    return best
#.....................................................................................
#.....................................................................................
#.....................................................................................

def _tt_run_scan(points, weights, target_centroid, target_axis=None,
                 devices=7, iterations=30, alpha=0.02,
                 eval_uniform=400, eval_empirical=200, eval_directional=40,
                 batches=3, processes_per_worker=100,
                 progress_interval=20,
                 # LENS-KEY OPTIONS/SETTINGS:
                 lens_enabled=True,
                 wavelength=550e-9,
                 key_angle=0.0, key_offset=0.0, key_width=0.2, key_depth=2e-6,
                 bump_amp=1e-6, bump_sigma=None,
                 delta=1e-6, disp_scale=5e-9):
    device_scales = [1.0+0.01*d for d in range(devices)]
    workers, temp_files = [], []
    nworkers = min(len(device_scales), processes_per_worker)
    chunk_size = int(math.ceil(len(device_scales)/float(nworkers)))
    for i in range(nworkers):
        chunk = device_scales[i*chunk_size:(i+1)*chunk_size]
        tf = "temp_worker_{}.json".format(i)
        temp_files.append(tf)
        from multiprocessing import Process
        p = Process(target=_tt_evaluate_scales_sim, args=(points, weights, target_centroid, target_axis,
                    chunk, iterations, alpha, eval_uniform, eval_empirical, eval_directional, batches,
                    12345+i*997, 5, tf, lens_enabled, wavelength, key_angle, key_offset, key_width,
                    key_depth, bump_amp, (bump_sigma if bump_sigma is not None else 0.5), delta, disp_scale))
        p.start()
        workers.append(p)
    for p in workers:
        p.join()
    best = None
    for tf in temp_files:
        if not os.path.exists(tf):
            continue
        try:
            with open(tf, 'r') as f:
                d = json.load(f)
        except Exception:
            continue
        try:
            os.remove(tf)
        except Exception:
            pass
        if best is None or d.get('score', 1e99) < best.get('score', 1e99):
            best = d
    return best
#.....................................................................................
#.....................................................................................
#.....................................................................................
#_____________________________________________________________________________________

def SL2_CTL_TEST():
    
    samples = [complex(1.0e-3, 0.5e-3), complex(0.8e-3, 1.2e-3),
               complex(1.5e-3, 0.9e-3), complex(0.9e-3, 0.2e-3)]
               
    weights = [0.4, 0.25, 0.2, 0.15]
    
    target_centroid = (1.1e-3, 0.8e-3)
    target_axis = (1.0, 0.0)

    best = _tt_run_scan(samples, weights, target_centroid, target_axis,
                        devices=6,
                        iterations=40,
                        alpha=0.018,
                        eval_uniform=300,
                        eval_empirical=150,
                        eval_directional=30,
                        batches=3,
                        progress_interval=10,
                        lens_enabled=True,
                        wavelength=550e-9,
                        key_angle=0.2,      # rad
                        key_offset=0.0,     # meters
                        key_width=0.2e-3,   # meters(smooth transition width)
                        key_depth=2e-6,     # meters(surface steps)
                        bump_amp=1e-6,      # meters
                        bump_sigma=0.5e-3,  # meters
                        delta=1e-6,         # meters(FD step for gradients)
                        disp_scale=5e-7)    # meters(per rad/m, tunable converts phase-grad --> transverse shift)
    print(f'\nBest Tracker:\n{best}')
SL2_CTL_TEST()