Unity Cup junior year simulation

unity.js

function factorial(n){
  let r=1;
  for(let i=2;i<=n;i++) r*=i;
  return r;
}

const BINOM = []; // BINOM[n][k] is the chance for n team members to have k arrows
const CUBE_DIST = [];

for (let n=0;n<=5;n++){
  BINOM[n]=[];
  for (let k=0;k<=n;k++){
    let comb =
      factorial(n) /
      (factorial(k)*factorial(n-k));
    BINOM[n][k] =
      comb *
      Math.pow(0.5,k) *
      Math.pow(0.5,n-k);
  }

  const cubeDist = new Map();
  for (let k = 0; k <= n; k++) {
    cubeDist.set((n + k) * 2, BINOM[n][k]);
  }
  CUBE_DIST[n] = cubeDist;
}

const BALL_MEMO = new Map();

// In case of BALL_MEMO, counts is caps
// In case of fillMemo, counts is fills
function ballKey(balls, counts, adjustCaps) {
  let key = 0;
  for (let i = 4; i >= 0; i--) {
    key = key * 6 + (adjustCaps ? Math.min(balls, counts[i]) : counts[i]);
  }
  return key * 6 + balls;
}

function decodeBallKey(key) {
  const balls = key % 6;
  const counts = new Array(5);
  for (let i = 0; i < 5; i++) {
    key = Math.trunc(key / 6);
    counts[i] = key % 6;
  }
  return [balls, counts];
}

function sortCaps(caps) {
  const first4 = caps.slice(0, 4).sort((a,b) => a - b);
  for (let i = 0; i < 4; i++) {
    caps[i] = first4[i];
  }
}

function getPrizeDist(balls, caps, prizeDist, baseP) {
  // We expect caps to be sorted
  const key = ballKey(balls, caps, true);
  if (!BALL_MEMO.has(key)) {
    populateBallMemo(balls, caps);
  }
  const fillDist = BALL_MEMO.get(key);

  for (let [fillKey, fillProb] of fillDist) {
    const fills = decodeBallKey(fillKey)[1];
    let heavyDist = new Map();
    heavyDist.set(0, 1);
    for (let i = 0; i < 4; i++) {
      const cubeDist = CUBE_DIST[5 - caps[i]];
      const bucketDist = new Map();
      for (let k = 0; k <= fills[i]; k++) {
        const prob = BINOM[fills[i]][k];
        for (let [s, p] of cubeDist) {
          const newScore = k + s;
          bucketDist.set(newScore, (bucketDist.get(newScore) || 0) + (prob * p));
        }
      }
      const newHeavyDist = new Map();
      for (let [score, prob] of bucketDist) {
        for (let [s, p] of heavyDist) {
          const newScore = Math.max(score, s);
          newHeavyDist.set(newScore, (newHeavyDist.get(newScore) || 0) + (prob * p));
        }
      }
      heavyDist = newHeavyDist;
    }
    const witCubeDist = CUBE_DIST[5 - caps[4]];
    const witDist = new Map();
    for (let k = 0; k <= fills[4]; k++) {
      const prob = BINOM[fills[4]][k] * fillProb * baseP;
      for (let [s, p] of witCubeDist) {
        const newScore = k + s;
        witDist.set(newScore, (witDist.get(newScore) || 0) + (prob * p));
      }
    }
    for (let [score, prob] of witDist) {
      for (let [s, p] of heavyDist) {
        prizeDist[s * 21 + score] += prob * p;
      }
    }
  }

  return prizeDist;
}

function populateBallMemo(balls, caps) {
  // We expect caps to be sorted
  const key = ballKey(balls, caps, true);
  const fillDist = new Map();
  if (balls === 0) {
    fillDist.set(0, 1);
    BALL_MEMO.set(key, fillDist);
    return;
  }
  const prevKey = ballKey(balls - 1, caps, true);
  if (!BALL_MEMO.has(prevKey)) {
    populateBallMemo(balls - 1, caps);
  }
  const prevDist = BALL_MEMO.get(prevKey);
  for (let [prevFillKey, prevProb] of prevDist) {
    const fills = decodeBallKey(prevFillKey)[1];
    const nonFull = [];
    for (let i = 0; i < 5; i++) {
      if (fills[i] < caps[i]) {
        nonFull.push(i);
      }
    }
    const x = nonFull.length;
    const pIncrement = prevProb / (x * 2 + 1);
    const missKey = ballKey(0, fills, false);
    fillDist.set(missKey, (fillDist.get(missKey) || 0) + pIncrement);
    for (let i of nonFull) {
      const newFills = fills.slice();
      newFills[i]++;
      const newKey = ballKey(0, newFills, false);
      fillDist.set(newKey, (fillDist.get(newKey) || 0) + pIncrement + pIncrement);
    }
  }
  BALL_MEMO.set(key, fillDist);
}

function visitDist(specialty, priority) {
  const denom = 550 + priority;
  const miss = 50 / denom;
  const train = miss + miss;
  const ret = [train, train, train, train, train, miss];
  ret[specialty] += priority / denom;
  return ret;
}

const PRIORITY = [
  visitDist(0, 100),
  visitDist(0, 50),
  visitDist(2, 35),
  visitDist(4, 17),
  visitDist(4, 25),
  visitDist(0, 0),
];

const capsDist = new Map();

for (let a = 0; a < 6; a++) {
  for (let b = 0; b < 6; b++) {
    for (let c = 0; c < 6; c++) {
      for (let d = 0; d < 6; d++) {
        for (let e = 0; e < 6; e++) {
          for (let f = 0; f < 6; f++) {
            const caps = [5, 5, 5, 5, 5, 5];
            caps[a]--;
            caps[b]--;
            caps[c]--;
            caps[d]--;
            caps[e]--;
//             caps[f]--;
            if (a != 5 && caps[a] === -1) {
              continue;
            }
            caps.length = 5;
            sortCaps(caps);
            const key = ballKey(0, caps, false);
            const p = PRIORITY[0][a] * PRIORITY[1][b] * PRIORITY[2][c] * PRIORITY[3][d] * PRIORITY[4][e] * PRIORITY[5][f];
            capsDist.set(key, (capsDist.get(key) || 0) + p);
          }
        }
      }
    }
  }
}

function prizeDistToCsv(pd) {
  const a = [];
  for (let i = 0; i < 21; i++) {
    a.push(pd.subarray(i * 21, i * 21 + 21).join(','));
  }
  return a.join('\n');
}

const totalPrizeDist = new Float64Array(21 * 21);
for (let [key, baseP] of capsDist) {
  getPrizeDist(3, decodeBallKey(key)[1], totalPrizeDist, baseP);
}

let checksum = 0;
for (let i = 0; i < 441; i++) {
  checksum += totalPrizeDist[i];
}
checksum = 1/checksum;
for (let i = 0; i < 441; i++) {
  totalPrizeDist[i] *= checksum;
}
checksum = 0;
for (let i = 0; i < 441; i++) {
  checksum += totalPrizeDist[i];
}
console.log(checksum);



/* ===============================
 *  Failure probabilities
 * =============================== */

const heavyFailTable = [
  1, // 0
  1, // 1
  0.964, // 2
  0.804, // 3
  0.654, // 4
  0.514, // 5
  0.394, // 6
  0.284, // 7
  0.184, // 8
  0.084  // 9
];

const lightFailTable = [
  0.844, // 0
  0.684, // 1
  0.524, // 2
  0.384, // 3
  0.264, // 4
  0.144, // 5
  0.044  // 6
];

function pFailHeavy(E) {
  return E > 9 ? 0 : heavyFailTable[E];
}

function pFailLight(E) {
  return E > 6 ? 0 : lightFailTable[E];
}

function heavyPenalty(pFail) {
  if (pFail < 0.20) return -5;
  if (pFail < 0.80) return -7;
  return -15;
}

const MAX_E = 20;
const NUM_PRIZES = 21;
const tol = 1e-10;

/* ===============================
 *  Relative Value Iteration
 * =============================== */

function solveMDP(prizeDist, maxIter) {

  const V = new Float64Array(MAX_E + 1);
  const Vnew = new Float64Array(MAX_E + 1);

  // for (let i = 0; i <= MAX_E; i++) {
  //   V[i] = i * 5;
  // }

  let rho = 0;

  for (let iter = 0; iter < maxIter; iter++) {

    /* ================================
     *  Compute in descending order
     * ================================ */

    for (let E = MAX_E; E >= 0; E--) {
      // Helper to read "current iteration if available"
      function Vmix(x) {
        // If x > E, it was already computed in this sweep
        return V[x];//return x > E ? Vnew[x] : V[x];
      }

      let expected = 0;

      for (let h = 0; h < NUM_PRIZES; h++) {
        for (let l = 0; l < NUM_PRIZES; l++) {

          const prob = prizeDist[h * 21 + l];
          if (prob < tol) continue;

          // REST
          const qRest = -0.1
            + 0.125 * Vmix(Math.min(20, E + 6))
            + 0.625 * Vmix(Math.min(20, E + 10))
            + 0.25  * Vmix(Math.min(20, E + 14));

          // HEAVY
          const pfHeavy = pFailHeavy(E);
          const penHeavy = heavyPenalty(pfHeavy);
          const EheavySucc = Math.max(0, E - 4);

          const qHeavy = (1 - pfHeavy) * (
            h
            + 0.9 * Vmix(EheavySucc)
            + 0.1 * Vmix(EheavySucc + 1)
          ) + pfHeavy * (penHeavy + Vmix(E));

          // LIGHT
          const pfLight = pFailLight(E);

          const qLight = (1 - pfLight) * (
            l
            + 0.9 * Vmix(Math.min(20, E + 1))
            + 0.1 * Vmix(Math.min(20, E + 2))
          ) + pfLight * Vmix(Math.min(20, E + 1));

          expected += Math.max(qRest, qLight, qHeavy) * prob;
        }
      }

      if (E === MAX_E) {
        Vnew[E] = expected;
      } else {
        Vnew[E] = 0.9 * expected + 0.1 * Vnew[Math.min(20, E + 2)];
      }
    }

    /* ================================
     *  Relative normalization
     * ================================ */

    const delta = Vnew[0];
    for (let E = 0; E <= MAX_E; E++) {
      Vnew[E] -= delta;
    }

    rho = delta;

    /* ================================
     *  Convergence check
     * ================================ */

    let maxDiff = 0;
    for (let E = 0; E <= MAX_E; E++) {
      const diff = Math.abs(Vnew[E] - V[E]);
      if (diff > maxDiff) maxDiff = diff;
      V[E] = Vnew[E];
    }

    if (maxDiff < tol) {
      console.log("Converged in", iter, "iterations");
      break;
    }
  }

  return { V, rho };
}

let result = solveMDP(totalPrizeDist, 10000);