fn takens_embed(series: Tensor, tau: Int, dim: Int) -> Tensor {
    let n = tensor_len(series)
    let rows = n - (dim - 1) * tau
    let result = t_zeros(rows, dim)
    for i in 0..dim {
        let offset = i * tau
        for r in 0..rows {
            let val = t_index_get(series, r + offset)
            let result = t_index_set(result, r, i, val)
        }
    }
    return result
}

fn optimal_tau(series: Tensor) -> Int {
    let n = tensor_len(series)
    let mean = tensor_sum(series) / n
    let found = 0
    let best_tau = 1
    for lag in 1..50 {
        if lag < n {
            if found == 0 {
                let sum = 0.0
                for i in 0..(n - lag) {
                    let a = t_index_get(series, i) - mean
                    let b = t_index_get(series, i + lag) - mean
                    let sum = sum + a * b
                }
                if sum < 0.0 {
                    let best_tau = lag
                    let found = 1
                }
            }
        }
    }
    return best_tau
}

fn curvature_3d(traj: Tensor) -> Tensor {
    let dr = t_gradient(traj)
    let ddr = t_gradient(dr)
    let cross = t_cross(dr, ddr)
    let cross_mag = t_norm(cross)
    let dr_mag = t_norm(dr)
    let rows = tensor_len(cross_mag)
    let kappa = t_zeros(rows)
    for i in 0..rows {
        let cm = t_index_get(cross_mag, i)
        let dm = t_index_get(dr_mag, i)
        let dm3 = dm * dm * dm
        if dm3 > 0.0000001 {
            let kappa = t_index_set(kappa, i, cm / dm3)
        }
    }
    return kappa
}

fn willmore_energy(traj: Tensor) -> Float {
    let kappa = curvature_3d(traj)
    let dr = t_gradient(traj)
    let ds = t_norm(dr)
    let n = tensor_len(kappa)
    let energy = 0.0
    for i in 0..n {
        let k = t_index_get(kappa, i)
        let s = t_index_get(ds, i)
        let energy = energy + k * k * s
    }
    return energy
}

fn curvature_to_volatility_ratio(series: Tensor) -> Float {
    let tau = optimal_tau(series)
    let traj = takens_embed(series, tau, 3)
    let kappa = curvature_3d(traj)
    let n = tensor_len(kappa)
    let mean_kappa = tensor_sum(kappa) / n
    let mean_s = tensor_sum(series) / tensor_len(series)
    let var_sum = 0.0
    for i in 0..tensor_len(series) {
        let diff = t_index_get(series, i) - mean_s
        let var_sum = var_sum + diff * diff
    }
    let sigma = 0.0
    if tensor_len(series) > 1 {
        let sigma = (var_sum / (tensor_len(series) - 1))
    }
    if sigma > 0.0001 {
        return mean_kappa / sigma
    }
    return 0.0
}

fn fraying_tensor(series: Tensor, micro_kappa: Float, macro_kappa: Float, prev_crossings: Int, curr_crossings: Int, vol_current: Float, vol_avg: Float) -> Tensor {
    let f_stress = 0.0
    if macro_kappa > 0.0000001 {
        let f_stress = micro_kappa / macro_kappa
    }
    if f_stress > 10.0 {
        let f_stress = 10.0
    }
    let f_isolation = 0.0
    if prev_crossings > 0 {
        let f_isolation = (prev_crossings - curr_crossings) / prev_crossings
        if f_isolation < 0.0 {
            let f_isolation = 0.0
        }
    }
    let f_thinning = 0.0
    if vol_avg > 0.0001 {
        let f_thinning = 1.0 - (vol_current / vol_avg)
        if f_thinning < 0.0 {
            let f_thinning = 0.0
        }
        if f_thinning > 1.0 {
            let f_thinning = 1.0
        }
    }
    return tensor(f_stress, f_isolation, f_thinning)
}

fn loom_metric(traj: Tensor, jones: Float) -> Float {
    let w = willmore_energy(traj)
    let algebraic = jones - 1.0
    if algebraic < 0.0 {
        let algebraic = 0.0 - algebraic
    }
    return 0.3 * algebraic + 0.7 * w
}

fn segments_intersect(ax1: Float, ay1: Float, ax2: Float, ay2: Float, bx1: Float, by1: Float, bx2: Float, by2: Float) -> Int {
    let d1 = (bx2 - bx1) * (ay1 - by1) - (by2 - by1) * (ax1 - bx1)
    let d2 = (bx2 - bx1) * (ay2 - by1) - (by2 - by1) * (ax2 - bx1)
    let d3 = (ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1)
    let d4 = (ax2 - ax1) * (by2 - ay1) - (ay2 - ay1) * (bx2 - ax1)
    if d1 * d2 < 0.0 {
        if d3 * d4 < 0.0 {
            return 1
        }
    }
    return 0
}

fn main() {
    println("============================================================")
    println("  LOOM-2: Topological Market Analysis")
    println("  Pure Scyla - No Python. Native Temporal Math.")
    println("  Aislinn McCubbins, February -> April 2026")
    println("============================================================")
    println("")

    let t = t_linspace(0.0, 3.0, 100)
    let sol = t_zeros(100)
    let btc = t_zeros(100)
    for i in 0..100 {
        let tv = t_index_get(t, i)
        let sol_v = 117.0 - 50.0 * tv + 20.0 * tv * tv
        let btc_v = 89000.0 - 4000.0 * tv
        let sol = t_index_set(sol, i, sol_v)
        let btc = t_index_set(btc, i, btc_v)
    }

    println("=== SOL (V-reversal) ===")
    let sol_tau = optimal_tau(sol)
    println("  optimal tau: " + sol_tau)
    let sol_traj = takens_embed(sol, sol_tau, 3)
    let sol_kappa = curvature_3d(sol_traj)
    let sol_n = tensor_len(sol_kappa)
    let sol_mean_k = tensor_sum(sol_kappa) / sol_n
    println("  mean curvature: " + sol_mean_k)
    let sol_W = willmore_energy(sol_traj)
    println("  willmore energy: " + sol_W)
    let sol_R = curvature_to_volatility_ratio(sol)
    println("  CVR: " + sol_R)
    if sol_R > 3.0 {
        println("  SIGNAL!")
    }
    if sol_R < 3.0 {
        println("  (noise)")
    }

    println("")
    println("=== BTC (Linear decline) ===")
    let btc_tau = optimal_tau(btc)
    println("  optimal tau: " + btc_tau)
    let btc_traj = takens_embed(btc, btc_tau, 3)
    let btc_kappa = curvature_3d(btc_traj)
    let btc_n = tensor_len(btc_kappa)
    let btc_mean_k = tensor_sum(btc_kappa) / btc_n
    println("  mean curvature: " + btc_mean_k)
    let btc_W = willmore_energy(btc_traj)
    println("  willmore energy: " + btc_W)
    let btc_R = curvature_to_volatility_ratio(btc)
    println("  CVR: " + btc_R)

    println("")
    println("=== Curvature Ratio ===")
    if btc_mean_k > 0.0000001 {
        let ratio = sol_mean_k / btc_mean_k
        println("  SOL/BTC curvature: " + ratio + "x")
    }

    println("")
    println("=== Loom Metric ===")
    let sol_L = loom_metric(sol_traj, 1.0)
    let btc_L = loom_metric(btc_traj, 1.0)
    println("  SOL: " + sol_L)
    println("  BTC: " + btc_L)

    println("")
    println("=== Crossing Detection ===")
    let crossings = 0
    for i in 1..sol_n {
        let ax1 = t_index_get(sol_traj, i - 1, 0)
        let ay1 = t_index_get(sol_traj, i - 1, 1)
        let ax2 = t_index_get(sol_traj, i, 0)
        let ay2 = t_index_get(sol_traj, i, 1)
        let bx1 = t_index_get(btc_traj, i - 1, 0)
        let by1 = t_index_get(btc_traj, i - 1, 1)
        let bx2 = t_index_get(btc_traj, i, 0)
        let by2 = t_index_get(btc_traj, i, 1)
        let hit = segments_intersect(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2)
        if hit == 1 {
            let crossings = crossings + 1
        }
    }
    println("  SOL x BTC crossings: " + crossings)

    println("")
    println("=== Fraying Tensor (SOL) ===")
    let fray = fraying_tensor(sol, sol_mean_k, btc_mean_k, 10, 8, 900000.0, 1000000.0)
    let f_stress = t_index_get(fray, 0)
    let f_iso = t_index_get(fray, 1)
    let f_thin = t_index_get(fray, 2)
    let f_total = f_stress * f_stress + f_iso * f_iso + f_thin * f_thin
    println("  stress:    " + f_stress)
    println("  isolation: " + f_iso)
    println("  thinning:  " + f_thin)
    println("  |F| total: " + f_total)
    if f_total < 0.25 {
        println("  status: HEALTHY")
    }
    if f_total > 0.25 {
        if f_total < 1.0 {
            println("  status: STRESSED")
        }
    }
    if f_total > 1.0 {
        if f_total < 2.25 {
            println("  status: FRAYING")
        }
    }
    if f_total > 2.25 {
        println("  status: SNAP IMMINENT")
    }

    println("")
    println("============================================================")
    println("  Pure Scyla. No Python. No numpy. No scipy.")
    println("  Takens embedding, curvature, crossings, fraying,")
    println("  Willmore energy, Loom metric - all native.")
    println("  The threads are alive.")
    println("============================================================")
}
