/*
 * Supplementary Code S1
 * ---------------------------------------------------------------------------
 * "Allocating Public Resources for Maximum Societal Benefit:
 *  An Optimization Approach to County Budgeting"
 *
 * Reproduces every quantitative result reported in the paper:
 *   - Optimized allocations and award scores for all weighting scenarios
 *     (Tables 3-5, Figure 1)
 *   - The uniform-allocation comparison (Table 4)
 *   - Concentration measures, HHI and Gini (Table 6)
 *   - The Monte Carlo parameter-perturbation analysis (Results,
 *     "Optimization Dominance"), using a fixed seed for exact
 *     reproducibility.
 *
 * Run with:  node supplementary_code.js     (Node.js, no dependencies)
 *
 * NOTE: This file is distributed as .txt because email providers block .js
 * attachments. To run: save as supplementary_code.js (or run directly with
 *   node supplementary_code.txt
 * — Node.js does not require the .js extension). No dependencies needed.
 * ---------------------------------------------------------------------------
 */
"use strict";

const S = 23.17; // total surplus, $M (FY 2027 Advertised Budget Plan, Vol. 1)
const DOMAINS = ["Public Safety", "Health", "Social Services", "Schools (FCPS)", "Economic Dev."];
const BETA = [0.55, 0.65, 0.70, 0.75, 0.45]; // diminishing-returns parameters (held constant)

const SCENARIOS = {
  "Baseline":     [0.18, 0.22, 0.22, 0.33, 0.05],
  "Equal alpha":  [0.20, 0.20, 0.20, 0.20, 0.20],
  "Health-first": [0.15, 0.35, 0.25, 0.15, 0.10],
  "Equity-first": [0.18, 0.22, 0.35, 0.20, 0.05],
};

// Optimal share s_i/S for a given Lagrange multiplier lambda (KKT condition):
//   s_i/S = ( lambda * S / (alpha_i * beta_i) )^( 1 / (beta_i - 1) )
function shareForLambda(alpha, beta, lambda) {
  return Math.pow((lambda * S) / (alpha * beta), 1 / (beta - 1));
}

// Bisect on lambda until shares sum to 1 (tolerance ~ $0.001M equivalent).
function solve(alphas, betas) {
  let lo = 1e-12, hi = 1e9;
  const totalShare = (lam) =>
    alphas.reduce((acc, a, i) => acc + shareForLambda(a, betas[i], lam), 0);
  for (let it = 0; it < 500; it++) {
    const mid = (lo + hi) / 2;
    if (totalShare(mid) > 1) lo = mid; else hi = mid;
  }
  const lam = (lo + hi) / 2;
  const shares = alphas.map((a, i) => shareForLambda(a, betas[i], lam));
  const allocs = shares.map(sh => sh * S);
  const awards = alphas.map((a, i) => a * Math.pow(shares[i], betas[i]));
  const A = awards.reduce((x, y) => x + y, 0);
  const A_uniform = alphas.reduce((acc, a, i) => acc + a * Math.pow(0.2, betas[i]), 0);
  return { shares, allocs, awards, A, A_uniform };
}

function hhiGini(allocs) {
  const tot = allocs.reduce((x, y) => x + y, 0);
  const pct = allocs.map(v => 100 * v / tot);
  const hhi = pct.reduce((acc, x) => acc + x * x, 0);
  const xs = [...pct].sort((a, b) => a - b);
  const n = xs.length, mean = xs.reduce((x, y) => x + y, 0) / n;
  let cum = 0;
  for (const a of xs) for (const b of xs) cum += Math.abs(a - b);
  return { hhi, gini: cum / (2 * n * n * mean) };
}

const f2 = x => x.toFixed(2), f3 = x => x.toFixed(3), f1 = x => x.toFixed(1);

// Largest-remainder rounding to cents so each displayed row sums exactly to S,
// matching the convention used in Tables 3 and 5 of the paper.
function roundRow(vals, total) {
  const cents = vals.map(v => Math.floor(v * 100));
  let deficit = Math.round(total * 100) - cents.reduce((x, y) => x + y, 0);
  const order = vals.map((v, i) => [v * 100 - cents[i], i]).sort((a, b) => b[0] - a[0]);
  for (let k = 0; k < deficit; k++) cents[order[k][1]] += 1;
  return cents.map(c => c / 100);
}

console.log("=== Optimized allocations by scenario (Tables 3 & 5) ===");
for (const [name, alphas] of Object.entries(SCENARIOS)) {
  const r = solve(alphas, BETA);
  const { hhi, gini } = hhiGini(r.allocs);
  const disp = roundRow(r.allocs, S);
  console.log(`\n${name}`);
  DOMAINS.forEach((d, i) =>
    console.log(`  ${d.padEnd(16)} $${f2(disp[i]).padStart(6)}M  ` +
      `share ${f1(100 * r.shares[i]).padStart(5)}%  award ${f3(r.awards[i])}`));
  console.log(`  Total award A = ${f3(r.A)}   Equal-dollar A = ${f3(r.A_uniform)}   gain = ${f3(r.A - r.A_uniform)}`);
  console.log(`  HHI (exact shares) = ${hhi.toFixed(0)}   Gini = ${f3(gini)}`);
}

console.log("\n=== Table 4: baseline optimized vs uniform, per-domain awards ===");
{
  const alphas = SCENARIOS["Baseline"];
  const r = solve(alphas, BETA);
  DOMAINS.forEach((d, i) => {
    const uni = alphas[i] * Math.pow(0.2, BETA[i]);
    console.log(`  ${d.padEnd(16)} optimized ${f3(r.awards[i])}   uniform ${f3(uni)}`);
  });
}

console.log("\n=== Table 6: HHI / Gini from printed Table 5 rows ===");
const printedRows = {
  "Uniform (flat)": [4.63, 4.63, 4.63, 4.63, 4.63],
  "Baseline":       [2.29, 3.39, 3.15, 14.10, 0.24],
  "Equal alpha":    [4.71, 4.81, 4.74, 4.55, 4.36],
  "Health-first":   [1.67, 14.36, 5.53, 0.71, 0.90],
  "Equity-first":   [2.33, 3.46, 15.18, 1.96, 0.24],
};
for (const [name, row] of Object.entries(printedRows)) {
  const { hhi, gini } = hhiGini(row);
  console.log(`  ${name.padEnd(14)} HHI ${Math.round(hhi).toLocaleString("en-US")}   Gini ${f3(gini)}`);
}

/* ---------------- Monte Carlo (fixed seed => exact reproducibility) -------
 * 10,000 draws. Each beta ~ U(assigned-0.10, assigned+0.10); each alpha ~
 * U(assigned-0.05, assigned+0.05), renormalized to sum to 1.00 per draw.
 * Deterministic PRNG: mulberry32 with seed 20260505 (FY27 adoption date).
 */
function mulberry32(seed) {
  return function () {
    let t = (seed += 0x6D2B79F5);
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}
const rand = mulberry32(20260505);
const A0 = SCENARIOS["Baseline"], N = 10000;
let schoolsLargest = 0, allPositive = 0;
const schoolShares = [], gains = [];
for (let d = 0; d < N; d++) {
  let a = A0.map(x => x + (rand() * 0.10 - 0.05));
  const sa = a.reduce((x, y) => x + y, 0);
  a = a.map(x => x / sa);
  const b = BETA.map(x => x + (rand() * 0.20 - 0.10));
  const r = solve(a, b);
  const shares = r.shares;
  const iSch = 3;
  if (shares.indexOf(Math.max(...shares)) === iSch) schoolsLargest++;
  schoolShares.push(shares[iSch]);
  const g = r.A - a.reduce((acc, x, i) => acc + x * Math.pow(0.2, b[i]), 0);
  gains.push(g);
  if (g > 0) allPositive++;
}
schoolShares.sort((x, y) => x - y);
gains.sort((x, y) => x - y);
const q = (arr, p) => arr[Math.min(arr.length - 1, Math.floor(p * arr.length))];
console.log("\n=== Monte Carlo (10,000 draws, seed 20260505) ===");
console.log(`  Schools allocation largest in ${(100 * schoolsLargest / N).toFixed(1)}% of draws`);
console.log(`  Schools share: median ${(100 * q(schoolShares, 0.5)).toFixed(1)}%  ` +
  `90% interval ${(100 * q(schoolShares, 0.05)).toFixed(1)}% to ${(100 * q(schoolShares, 0.95)).toFixed(1)}%`);
console.log(`  Gain over own equal-dollar baseline positive in ${(100 * allPositive / N).toFixed(1)}% of draws`);
console.log(`  Gain: median ${q(gains, 0.5).toFixed(3)}   minimum ${gains[0].toFixed(3)}`);