|
| 1 | +package rtss.data.curves; |
| 2 | + |
| 3 | +import rtss.util.Util; |
| 4 | + |
| 5 | +/* |
| 6 | + * Деформировать "носик" кривой l(x) с дневным разрешением в пределах первого года. |
| 7 | + * Использовать полином 4-й степени. |
| 8 | + */ |
| 9 | +public class SculptDailyLX |
| 10 | +{ |
| 11 | + public static double[] scultDailyLX(double[] f, double convexityControl) throws Exception |
| 12 | + { |
| 13 | + int x1 = 365; |
| 14 | + |
| 15 | + // Clone the original array |
| 16 | + double[] fClone = f.clone(); |
| 17 | + |
| 18 | + // Boundary conditions |
| 19 | + double f0 = f[0]; // f(0) |
| 20 | + double fx1 = f[x1]; // f(x1) |
| 21 | + double dfx1 = (f[x1 + 1] - f[x1]) / 1.0; // Approximate f'(x1) |
| 22 | + double d2fx1 = (f[x1 + 1] - 2 * f[x1] + f[x1 - 1]) / 1.0; // Approximate f''(x1) |
| 23 | + |
| 24 | + // Coefficients of the quartic polynomial f(x) = ax^4 + bx^3 + cx^2 + dx + e |
| 25 | + double e = f0; // From f(0) = f0 |
| 26 | + |
| 27 | + // Solve the system of equations for a, b, c, d: |
| 28 | + // 1. a*x1^4 + b*x1^3 + c*x1^2 + d*x1 + e = fx1 |
| 29 | + // 2. 4a*x1^3 + 3b*x1^2 + 2c*x1 + d = dfx1 |
| 30 | + // 3. 12a*x1^2 + 6b*x1 + 2c = d2fx1 |
| 31 | + // 4. a = convexityControl (extra condition to control convexity) |
| 32 | + |
| 33 | + double a = convexityControl; // Control parameter for convexity |
| 34 | + |
| 35 | + // Substitute a into the equations and solve for b, c, d |
| 36 | + double[][] A = { |
| 37 | + { x1 * x1 * x1, x1 * x1, x1 }, |
| 38 | + { 3 * x1 * x1, 2 * x1, 1 }, |
| 39 | + { 6 * x1, 2, 0 } |
| 40 | + }; |
| 41 | + |
| 42 | + double[] B = { |
| 43 | + fx1 - e - a * x1 * x1 * x1 * x1, |
| 44 | + dfx1 - 4 * a * x1 * x1 * x1, |
| 45 | + d2fx1 - 12 * a * x1 * x1 |
| 46 | + }; |
| 47 | + |
| 48 | + // Solve the system Ax = B for x = [b, c, d] |
| 49 | + double[] coefficients = solveLinearSystem(A, B); |
| 50 | + double b = coefficients[0]; |
| 51 | + double c = coefficients[1]; |
| 52 | + double d = coefficients[2]; |
| 53 | + |
| 54 | + // Fill in the missing values |
| 55 | + for (int x = 1; x < x1; x++) |
| 56 | + { |
| 57 | + fClone[x] = a * x * x * x * x + b * x * x * x + c * x * x + d * x + e; |
| 58 | + } |
| 59 | + |
| 60 | + validate(Util.splice(fClone, 0, x1 + 100)); |
| 61 | + |
| 62 | + return fClone; |
| 63 | + } |
| 64 | + |
| 65 | + // Helper method to solve a 3x3 linear system Ax = B |
| 66 | + private static double[] solveLinearSystem(double[][] A, double[] B) |
| 67 | + { |
| 68 | + int n = B.length; |
| 69 | + double[] x = new double[n]; |
| 70 | + |
| 71 | + // Using Gaussian elimination |
| 72 | + for (int i = 0; i < n; i++) |
| 73 | + { |
| 74 | + // Pivot for the current row |
| 75 | + double pivot = A[i][i]; |
| 76 | + for (int j = i + 1; j < n; j++) |
| 77 | + { |
| 78 | + double factor = A[j][i] / pivot; |
| 79 | + B[j] -= factor * B[i]; |
| 80 | + for (int k = i; k < n; k++) |
| 81 | + { |
| 82 | + A[j][k] -= factor * A[i][k]; |
| 83 | + } |
| 84 | + } |
| 85 | + } |
| 86 | + |
| 87 | + // Back substitution |
| 88 | + for (int i = n - 1; i >= 0; i--) |
| 89 | + { |
| 90 | + x[i] = B[i]; |
| 91 | + for (int j = i + 1; j < n; j++) |
| 92 | + { |
| 93 | + x[i] -= A[i][j] * x[j]; |
| 94 | + } |
| 95 | + x[i] /= A[i][i]; |
| 96 | + } |
| 97 | + |
| 98 | + return x; |
| 99 | + } |
| 100 | + |
| 101 | + private static void validate(double[] f) throws Exception |
| 102 | + { |
| 103 | + double[] d1 = derivative(f); |
| 104 | + double[] d2 = derivative(d1); |
| 105 | + |
| 106 | + if (!Util.isMonotonicallyDecreasing(f, true)) |
| 107 | + throw new Exception("Improperly scuplted lx (non-decreasing)"); |
| 108 | + |
| 109 | + if (!Util.isNegative(d1)) |
| 110 | + throw new Exception("Improperly scuplted lx (non-negative d1)"); |
| 111 | + |
| 112 | + if (!Util.isPositive(d2)) |
| 113 | + throw new Exception("Improperly scuplted lx (non-positive d2)"); |
| 114 | + } |
| 115 | + |
| 116 | + private static double[] derivative(double[] p) |
| 117 | + { |
| 118 | + if (p.length <= 1) |
| 119 | + return new double[0]; |
| 120 | + |
| 121 | + double[] d = new double[p.length - 1]; |
| 122 | + for (int i = 0; i <= p.length - 2; i++) |
| 123 | + d[i] = p[i + 1] - p[i]; |
| 124 | + return d; |
| 125 | + } |
| 126 | + |
| 127 | + private static double[] d2(double[] p) |
| 128 | + { |
| 129 | + return derivative(derivative(p)); |
| 130 | + } |
| 131 | +} |
0 commit comments