Make grow formulas faster and more accurate. (#1044)

This commit is contained in:
David Walker 2024-01-31 16:27:31 -08:00 committed by GitHub
parent b6b4788845
commit 55e21d1e19
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
5 changed files with 184 additions and 140 deletions

@ -24,8 +24,8 @@ export const CONSTANTS: {
NeuroFluxGovernorLevelMult: number;
NumNetscriptPorts: number;
HomeComputerMaxRam: number;
ServerBaseGrowthRate: number;
ServerMaxGrowthRate: number;
ServerBaseGrowthIncr: number;
ServerMaxGrowthLog: number;
ServerFortifyAmount: number;
ServerWeakenAmount: number;
PurchasedServerLimit: number;
@ -125,8 +125,8 @@ export const CONSTANTS: {
// Server-related constants
HomeComputerMaxRam: 1073741824, // 2 ^ 30
ServerBaseGrowthRate: 1.03, // Unadjusted Growth rate
ServerMaxGrowthRate: 1.0035, // Maximum possible growth rate (max rate accounting for server security)
ServerBaseGrowthIncr: 0.03, // Unadjusted growth increment (growth rate is this * adjustment + 1)
ServerMaxGrowthLog: 0.00349388925425578, // Maximum possible growth rate accounting for server security, precomputed as log1p(.0035)
ServerFortifyAmount: 0.002, // Amount by which server's security increases when its hacked/grown
ServerWeakenAmount: 0.05, // Amount by which server's security decreases when weakened

@ -1,9 +1,8 @@
import { GetServer, createUniqueRandomIp, ipExists } from "./AllServers";
import { Server, IConstructorParams } from "./Server";
import { BaseServer } from "./BaseServer";
import { calculateServerGrowth } from "./formulas/grow";
import { calculateServerGrowth, calculateServerGrowthLog } from "./formulas/grow";
import { currentNodeMults } from "../BitNode/BitNodeMultipliers";
import { CONSTANTS } from "../Constants";
import { Player } from "@player";
import { CompletedProgramName, LiteratureName } from "@enums";
@ -52,24 +51,7 @@ export function safelyCreateUniqueServer(params: IConstructorParams): Server {
*/
export function numCycleForGrowth(server: IServer, growth: number, cores = 1): number {
if (!server.serverGrowth) return Infinity;
const hackDifficulty = server.hackDifficulty ?? 100;
let ajdGrowthRate = 1 + (CONSTANTS.ServerBaseGrowthRate - 1) / hackDifficulty;
if (ajdGrowthRate > CONSTANTS.ServerMaxGrowthRate) {
ajdGrowthRate = CONSTANTS.ServerMaxGrowthRate;
}
const serverGrowthPercentage = server.serverGrowth / 100;
const coreBonus = getCoreBonus(cores);
const cycles =
Math.log(growth) /
(Math.log(ajdGrowthRate) *
Player.mults.hacking_grow *
serverGrowthPercentage *
currentNodeMults.ServerGrowthRate *
coreBonus);
return cycles;
return Math.log(growth) / calculateServerGrowthLog(server, 1, Player, cores);
}
/**
@ -91,130 +73,106 @@ export function numCycleForGrowthCorrected(
): number {
if (!server.serverGrowth) return Infinity;
const moneyMax = server.moneyMax ?? 1;
const hackDifficulty = server.hackDifficulty ?? 100;
if (startMoney < 0) startMoney = 0; // servers "can't" have less than 0 dollars on them
if (targetMoney > moneyMax) targetMoney = moneyMax; // can't grow a server to more than its moneyMax
if (targetMoney <= startMoney) return 0; // no growth --> no threads
// exponential base adjusted by security
const adjGrowthRate = 1 + (CONSTANTS.ServerBaseGrowthRate - 1) / hackDifficulty;
const exponentialBase = Math.min(adjGrowthRate, CONSTANTS.ServerMaxGrowthRate); // cap growth rate
// total of all grow thread multipliers
const serverGrowthPercentage = server.serverGrowth / 100.0;
const coreMultiplier = getCoreBonus(cores);
const threadMultiplier =
serverGrowthPercentage * person.mults.hacking_grow * coreMultiplier * currentNodeMults.ServerGrowthRate;
const k = calculateServerGrowthLog(server, 1, person, cores);
/* To understand what is done below we need to do some math. I hope the explanation is clear enough.
* First of, the names will be shortened for ease of manipulation:
* n:= targetMoney (n for new), o:= startMoney (o for old), b:= exponentialBase, t:= threadMultiplier, c:= cycles/threads
* c is what we are trying to compute.
* n:= targetMoney (n for new), o:= startMoney (o for old), k:= calculateServerGrowthLog, x:= threads
* x is what we are trying to compute.
*
* After growing, the money on a server is n = (o + c) * b^(c*t)
* c appears in an exponent and outside it, this is usually solved using the productLog/lambert's W special function
* this function will be noted W in the following
* The idea behind lambert's W function is W(x)*exp(W(x)) = x, or in other words, solving for y, y*exp(y) = x, as a function of x
* This function is provided in some advanced math library but we will compute it ourself here.
* After growing, the money on a server is n = (o + x) * exp(k*x)
* x appears in an exponent and outside it, this is usually solved using the productLog/lambert's W special function,
* but it turns out that due to floating-point range issues this approach is *useless* to us, so it will be ignored.
*
* Let's get back to solving the equation. It cannot be rewrote using W immediately because the base of the exponentiation is b
* b^(c*t) = exp(ln(b)*c*t) (this is how a^b is defined on reals, it matches the definition on integers)
* so n = (o + c) * exp(ln(b)*c*t) , W still cannot be used directly. We want to eliminate the other terms in 'o + c' and 'ln(b)*c*t'.
*
* A change of variable will do. The idea is to add an equation introducing a new variable (w here) in the form c = f(w) (for some f)
* With this equation we will eliminate all references to c, then solve for w and plug the result in the new equation to get c.
* The change of variable performed here should get rid of the unwanted terms mentioned above, c = w/(ln(b)*t) - o should help.
* This change of variable is allowed because whatever the value of c is, there is a value of w such that this equation holds:
* w = (c + o)*ln(b)*t (see how we used the terms we wanted to eliminate in order to build this variable change)
*
* We get n = (o + w/(ln(b)*t) - o) * exp(ln(b)*(w/(ln(b)*t) - o)*t) [ = w/(ln(b)*t) * exp(w - ln(b)*o*t) ]
* The change of variable exposed exp(w - o*ln(b)*t), we can rewrite that with exp(a - b) = exp(a)/exp(b) to isolate 'w*exp(w)'
* n = w/(ln(b)*t) * exp(w)/exp(ln(b)*o*t) [ = w*exp(w) / (ln(b) * t * b^(o*t)) ]
* Almost there, we just need to cancel the denominator on the right side of the equation:
* n * ln(b) * t * b^(o*t) = w*exp(w), Thus w = W(n * ln(b) * t * b^(o*t))
* Finally we invert the variable change: c = W(n * ln(b) * t * b^(o*t))/(ln(b)*t) - o
*
* There is still an issue left: b^(o*t) doesn't fit inside a double precision float
* because the typical amount of money on servers is around 10^6~10^9
* We need to get an approximation of W without computing the power when o is huge
* Thankfully an approximation giving ~30% error uses log immediately so we will use
* W(n * ln(b) * t * b^(o*t)) ~= log(n * ln(b) * t * b^(o*t)) = log(n * ln(b) * t) + log(exp(ln(b) * o * t))
* = log(n * ln(b) * t) + ln(b) * o * t
* (thanks to Drak for the grow formula, f4113nb34st and Wolfram Alpha for the rewrite, dwRchyngqxs for the explanation)
*/
const x = threadMultiplier * Math.log(exponentialBase);
const y = startMoney * x + Math.log(targetMoney * x);
/* Code for the approximation of lambert's W function is adapted from
* https://git.savannah.gnu.org/cgit/gsl.git/tree/specfunc/lambert.c
* using the articles [1] https://doi.org/10.1007/BF02124750 (algorithm above)
* and [2] https://doi.org/10.1145/361952.361970 (initial approximation when x < 2.5)
*/
let w;
if (y < Math.log(2.5)) {
/* exp(y) can be safely computed without overflow.
* The relative error on the result is better when exp(y) < 2.5
* using Padé rational fraction approximation [2](5)
*/
const ey = Math.exp(y);
w = (ey + (4 / 3) * ey * ey) / (1 + (7 / 3) * ey + (5 / 6) * ey * ey);
} else {
/* obtain initial approximation from rough asymptotic [1](4.18)
* w = y [- log y when 0 <= y]
*/
w = y;
if (y > 0) w -= Math.log(y);
}
let cycles = w / x - startMoney;
/* Iterative refinement, the goal is to correct c until |(o + c) * b^(c*t) - n| < 1
* or the correction on the approximation is less than 1
* The Newton-Raphson method will be used, this method is a classic to find roots of functions
* (given f, find c such that f(c) = 0).
* Instead, we proceed directly to Newton-Raphson iteration. We first rewrite the equation in
* log-form, since iterating it this way has faster convergence: log(n) = log(o+x) + k*x.
* Now our goal is to find the zero of f(x) = log((o+x)/n) + k*x.
* (Due to the shape of the function, there will be a single zero.)
*
* The idea of this method is to take the horizontal position at which the horizontal axis
* intersects with of the tangent of the function's curve as the next approximation.
* It is equivalent to treating the curve as a line (it is called a first order approximation)
* If the current approximation is c then the new approximated value is c - f(c)/f'(c)
* If the current approximation is x then the new approximated value is x - f(x)/f'(x)
* (where f' is the derivative of f).
*
* In our case f(c) = (o + c) * b^(c*t) - n, f'(c) = d((o + c) * b^(c*t) - n)/dc
* = (ln(b)*t * (c + o) + 1) * b^(c*t)
* And the update step is c[new] = c[old] - ((o + c) * b^(c*t) - n)/((ln(b)*t * (o + c) + 1) * b^(c*t))
* In our case f(x) = log((o+x)/n) + k*x, f'(x) = d(log((o+x)/n) + k*x)/dx
* = 1/(o + x) + k
* And the update step is x[new] = x - (log((o+x)/n) + k*x)/(1/(o+x) + k)
* We can simplify this by bringing the first term up into the fraction:
* = (x * (1/(o+x) + k) - log((o+x)/n) - k*x) / (1/(o+x) + k)
* = (x/(o+x) - log((o+x)/n)) / (1/(o+x) + k) [multiplying top and bottom by (o+x)]
* = (x - (o+x)*log((o+x)/n)) / (1 + (o+x)*k)
*
* The main question to ask when using this method is "does it converges?"
* The main question to ask when using this method is "does it converge?"
* (are the approximations getting better?), if it does then it does quickly.
* DOES IT CONVERGES? In the present case it does. The reason why doesn't help explaining the algorithm.
* If you are interested then check out the wikipedia page.
* Since the derivative is always positive but also strictly decreasing, convergence is guaranteed.
* This also provides the useful knowledge that any x which starts *greater* than the solution will
* undershoot across to the left, while values *smaller* than the zero will continue to find
* closer approximations that are still smaller than the final value.
*
* Of great importance for reducing the number of iterations is starting with a good initial
* guess. We use a very simple starting condition: x_0 = n - o. We *know* this will always overshot
* the target, usually by a vast amount. But we can run it manually through one Newton iteration
* to get a better start with nice properties:
* x_1 = ((n - o) - (n - o + o)*log((n-o+o)/n)) / (1 + (n-o+o)*k)
* = ((n - o) - n * log(n/n)) / (1 + n*k)
* = ((n - o) - n * 0) / (1 + n*k)
* = (n - o) / (1 + n*k)
* We can do the same procedure with the exponential form of Newton's method, starting from x_0 = 0.
* This gives x_1 = (n - o) / (1 + o*k), (full derivation omitted) which will be an overestimate.
* We use a weighted average of the denominators to get the final guess:
* x = (n - o) / (1 + (1/16*n + 15/16*o)*k)
* The reason for this particular weighting is subtle; it is exactly representable and holds up
* well under a wide variety of conditions, making it likely that the we start within 1 thread of
* correct. It particularly bounds the worst-case to 3 iterations, and gives a very wide swatch
* where 2 iterations is good enough.
*
* The accuracy of the initial guess is good for many inputs - often one iteration
* is sufficient. This means the overall cost is two logs (counting the one in calculateServerGrowthLog),
* possibly one exp, 5 divisions, and a handful of basic arithmetic.
*/
let bt = exponentialBase ** threadMultiplier;
if (bt == Infinity) bt = 1e300;
let corr = Infinity;
// Two sided error because we do not want to get stuck if the error stays on the wrong side
const guess = (targetMoney - startMoney) / (1 + (targetMoney * (1 / 16) + startMoney * (15 / 16)) * k);
let x = guess;
let diff;
do {
// c should be above 0 so Halley's method can't be used, we have to stick to Newton-Raphson
let bct = bt ** cycles;
if (bct == Infinity) bct = 1e300;
const opc = startMoney + cycles;
let diff = opc * bct - targetMoney;
if (diff == Infinity) diff = 1e300;
corr = diff / (opc * x + 1.0) / bct;
cycles -= corr;
} while (Math.abs(corr) >= 1);
/* c is now within +/- 1 of the exact result.
* We want the ceiling of the exact result, so the floor if the approximation is above,
* the ceiling if the approximation is in the same unit as the exact result,
* and the ceiling + 1 if the approximation is below.
const ox = startMoney + x;
// Have to use division instead of multiplication by inverse, because
// if targetMoney is MIN_VALUE then inverting gives Infinity
const newx = (x - ox * Math.log(ox / targetMoney)) / (1 + ox * k);
diff = newx - x;
x = newx;
} while (diff < -1 || diff > 1);
/* If we see a diff of 1 or less we know all future diffs will be smaller, and the rate of
* convergence means the *sum* of the diffs will be less than 1.
* In most cases, our result here will be ceil(x).
*/
const fca = Math.floor(cycles);
if (targetMoney <= (startMoney + fca) * Math.pow(exponentialBase, fca * threadMultiplier)) {
return fca;
const ccycle = Math.ceil(x);
if (ccycle - x > 0.999999) {
// Rounding-error path: It's possible that we slightly overshot the integer value due to
// rounding error, and more specifically precision issues with log and the size difference of
// startMoney vs. x. See if a smaller integer works. Most of the time, x was not close enough
// that we need to try.
const fcycle = ccycle - 1;
if (targetMoney <= (startMoney + fcycle) * Math.exp(k * fcycle)) {
return fcycle;
}
const cca = Math.ceil(cycles);
if (targetMoney <= (startMoney + cca) * Math.pow(exponentialBase, cca * threadMultiplier)) {
return cca;
}
return cca + 1;
if (ccycle >= x + ((diff <= 0 ? -diff : diff) + 0.000001)) {
// Fast-path: We know the true value is somewhere in the range [x, x + |diff|] but the next
// greatest integer is past this. Since we have to round up grows anyway, we can return this
// with no more calculation. We need some slop due to rounding errors - we can't fast-path
// a value that is too small.
return ccycle;
}
if (targetMoney <= (startMoney + ccycle) * Math.exp(k * ccycle)) {
return ccycle;
}
return ccycle + 1;
}
//Applied server growth for a single server. Returns the percentage growth
@ -226,7 +184,7 @@ export function processSingleServerGrowth(server: Server, threads: number, cores
}
const oldMoneyAvailable = server.moneyAvailable;
server.moneyAvailable += 1 * threads; // It can be grown even if it has no money
server.moneyAvailable += threads; // It can be grown even if it has no money
server.moneyAvailable *= serverGrowth;
// in case of data corruption

@ -2,24 +2,31 @@ import { CONSTANTS } from "../../Constants";
import { currentNodeMults } from "../../BitNode/BitNodeMultipliers";
import { Person as IPerson, Server as IServer } from "@nsdefs";
export function calculateServerGrowth(server: IServer, threads: number, p: IPerson, cores = 1): number {
if (!server.serverGrowth) return 0;
// Returns the log of the growth rate. When passing 1 for threads, this gives a useful constant.
export function calculateServerGrowthLog(server: IServer, threads: number, p: IPerson, cores = 1): number {
if (!server.serverGrowth) return -Infinity;
const hackDifficulty = server.hackDifficulty ?? 100;
const numServerGrowthCycles = Math.max(Math.floor(threads), 0);
const numServerGrowthCycles = Math.max(threads, 0);
//Get adjusted growth rate, which accounts for server security
const growthRate = CONSTANTS.ServerBaseGrowthRate;
let adjGrowthRate = 1 + (growthRate - 1) / hackDifficulty;
if (adjGrowthRate > CONSTANTS.ServerMaxGrowthRate) {
adjGrowthRate = CONSTANTS.ServerMaxGrowthRate;
//Get adjusted growth log, which accounts for server security
//log1p computes log(1+p), it is far more accurate for small values.
let adjGrowthLog = Math.log1p(CONSTANTS.ServerBaseGrowthIncr / hackDifficulty);
if (adjGrowthLog >= CONSTANTS.ServerMaxGrowthLog) {
adjGrowthLog = CONSTANTS.ServerMaxGrowthLog;
}
//Calculate adjusted server growth rate based on parameters
const serverGrowthPercentage = server.serverGrowth / 100;
const numServerGrowthCyclesAdjusted =
numServerGrowthCycles * serverGrowthPercentage * currentNodeMults.ServerGrowthRate;
const serverGrowthPercentageAdjusted = serverGrowthPercentage * currentNodeMults.ServerGrowthRate;
//Apply serverGrowth for the calculated number of growth cycles
const coreBonus = 1 + (cores - 1) / 16;
return Math.pow(adjGrowthRate, numServerGrowthCyclesAdjusted * p.mults.hacking_grow * coreBonus);
const coreBonus = 1 + (cores - 1) * (1 / 16);
// It is critical that numServerGrowthCycles (aka threads) is multiplied last,
// so that it rounds the same way as numCycleForGrowthCorrected.
return adjGrowthLog * serverGrowthPercentageAdjusted * p.mults.hacking_grow * coreBonus * numServerGrowthCycles;
}
export function calculateServerGrowth(server: IServer, threads: number, p: IPerson, cores = 1): number {
if (!server.serverGrowth) return 0;
return Math.exp(calculateServerGrowthLog(server, threads, p, cores));
}

76
test/jest/Grow.test.ts Normal file

@ -0,0 +1,76 @@
import { PlayerObject } from "../../src/PersonObjects/Player/PlayerObject";
import { Server } from "../../src/Server/Server";
import { calculateServerGrowth } from "../../src/Server/formulas/grow";
import { numCycleForGrowthCorrected } from "../../src/Server/ServerHelpers";
test("Grow is accurate", () => {
// Tests that certain special values work out to exactly what we'd expect,
// given the formulas. The growth multiplier maxes out at 1.0035, and the
// increment is .03 so with a difficulty of 10 it should be .003.
// These tests are *exact* because the whole point is that the math should
// get exactly the right value (or as close as is possible with floating-point).
const server = new Server({ hostname: "foo", hackDifficulty: 5, serverGrowth: 100 });
const player = new PlayerObject();
expect(calculateServerGrowth(server, 1, player)).toBe(1.0035);
expect(calculateServerGrowth(server, 2, player)).toBe(1.00701225);
server.hackDifficulty = 10;
expect(calculateServerGrowth(server, 1, player)).toBe(1.003);
expect(calculateServerGrowth(server, 2, player)).toBe(1.006009);
expect(calculateServerGrowth(server, 3, player)).toBe(1.009027027);
expect(calculateServerGrowth(server, 4, player)).toBe(1.012054108081);
});
describe("numCycleForGrowthCorrected reverses calculateServerGrowth", () => {
// Test that numCycleForGrowthCorrected() functions properly as the inverse
// of calculateServerGrowth().
// calculateServerGrowth() works for *any* given number of threads, but is
// usually passed an integer. numCycleForGrowthCorrected() *always* returns
// an integer, the ceiling of the floating-point value. When reversing an
// integer number of threads, it should *always* return the same integer.
// Similarly, if we pass the next-highest representable number as a target,
// it should *always* return the next largest integer, since the previous
// number is no longer sufficient.
// This is an arbitrary transcedental constant.
const multiplier = Math.exp(1.4);
const server = new Server({ hostname: "foo", hackDifficulty: 10 * multiplier, serverGrowth: 100 });
server.moneyMax = 1e308; // Not available as a constructor param
const player = new PlayerObject();
const tests = [];
while (server.moneyAvailable < 5e49) {
tests.push(server.moneyAvailable);
server.moneyAvailable = (server.moneyAvailable + 59) * calculateServerGrowth(server, 59, player);
}
test.each(tests)("startMoney: %f", (money: number) => {
// Do fewer threads to save on test time
for (let threads = 0; threads < 30; threads++) {
const newMoney = (money + threads) * calculateServerGrowth(server, threads, player);
const eps = newMoney ? 2 ** (Math.floor(Math.log2(newMoney)) - 52) : Number.MIN_VALUE;
expect(numCycleForGrowthCorrected(server, newMoney, money, 1, player)).toBe(threads);
const value = numCycleForGrowthCorrected(server, newMoney + eps, money, 1, player);
// Write our own check because Jest is a goblin that can't provide context
if (value !== threads + 1) {
console.log(
`newMoney: ${newMoney} eps: ${eps} newMoney+eps: ${newMoney + eps} value: ${value} threads: ${threads}`,
);
expect(value).toBe(threads + 1);
}
}
});
const tests2 = [];
for (let t = 0; t <= 900000; t += 2000) {
tests2.push([t, t * calculateServerGrowth(server, t, player)]);
}
test.each(tests2)("threads: %f newMoney: %f", (threads: number, newMoney: number) => {
const eps = newMoney ? 2 ** (Math.floor(Math.log2(newMoney)) - 52) : Number.MIN_VALUE;
expect(numCycleForGrowthCorrected(server, newMoney, 0, 1, player)).toBe(threads);
const value = numCycleForGrowthCorrected(server, newMoney + eps, 0, 1, player);
// Write our own check because Jest is a goblin that can't provide context
if (value !== threads + 1) {
console.log(
`newMoney: ${newMoney} eps: ${eps} newMoney+eps: ${newMoney + eps} value: ${value} threads: ${threads}`,
);
expect(value).toBe(threads + 1);
}
});
});

@ -477,7 +477,10 @@ describe("Stock Market Tests", function () {
const initValue = initialValues[stock.symbol];
expect(initValue.price).not.toEqual(stock.price);
if (initValue.otlkMag === stock.otlkMag && initValue.b === stock.b) {
throw new Error("expected either price or otlkMag to be different");
throw new Error(
"expected either price or otlkMag to be different: " +
`stock: ${stockName} otlkMag: ${stock.otlkMag} b: ${stock.b}`,
);
}
}
});