List of usage examples for org.apache.commons.math3.analysis UnivariateFunction value
double value(double x);
From source file:com.opengamma.strata.math.impl.util.CommonsMathWrapperTest.java
@Test public void test1DFunction() { UnivariateFunction commons = CommonsMathWrapper.wrapUnivariate(OG_FUNCTION_1D); for (int i = 0; i < 100; i++) { assertEquals(OG_FUNCTION_1D.apply((double) i), commons.value(i), 1e-15); }//from w ww . j av a 2 s . c o m }
From source file:au.gov.ga.conn4d.utils.BicubicSplineInterpolator.java
public BicubicSplineInterpolatingFunction interpolate(final double[] xval, final double[] yval, final double[][] fval) throws NoDataException, DimensionMismatchException, NonMonotonicSequenceException, NumberIsTooSmallException { if (xval.length == 0 || yval.length == 0 || fval.length == 0) { throw new NoDataException(); }/*from w w w . jav a 2s. c o m*/ if (xval.length != fval.length) { throw new DimensionMismatchException(xval.length, fval.length); } MathArrays.checkOrder(xval); MathArrays.checkOrder(yval); final int xLen = xval.length; final int yLen = yval.length; // Samples (first index is y-coordinate, i.e. subarray variable is x) // 0 <= i < xval.length // 0 <= j < yval.length // fX[j][i] = f(xval[i], yval[j]) final double[][] fX = new double[yLen][xLen]; for (int i = 0; i < xLen; i++) { if (fval[i].length != yLen) { throw new DimensionMismatchException(fval[i].length, yLen); } for (int j = 0; j < yLen; j++) { fX[j][i] = fval[i][j]; } } final SplineInterpolator spInterpolator = new SplineInterpolator(); // For each line y[j] (0 <= j < yLen), construct a 1D spline with // respect to variable x final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; for (int j = 0; j < yLen; j++) { ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); } // For each line x[i] (0 <= i < xLen), construct a 1D spline with // respect to variable y generated by array fY_1[i] final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; for (int i = 0; i < xLen; i++) { xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); } // Partial derivatives with respect to x at the grid knots final double[][] dFdX = new double[xLen][yLen]; for (int j = 0; j < yLen; j++) { final UnivariateFunction f = ySplineX[j].derivative(); for (int i = 0; i < xLen; i++) { dFdX[i][j] = f.value(xval[i]); } } // Partial derivatives with respect to y at the grid knots final double[][] dFdY = new double[xLen][yLen]; for (int i = 0; i < xLen; i++) { final UnivariateFunction f = xSplineY[i].derivative(); for (int j = 0; j < yLen; j++) { dFdY[i][j] = f.value(yval[j]); } } // Cross partial derivatives final double[][] d2FdXdY = new double[xLen][yLen]; for (int i = 0; i < xLen; i++) { final int nI = nextIndex(i, xLen); final int pI = previousIndex(i); for (int j = 0; j < yLen; j++) { final int nJ = nextIndex(j, yLen); final int pJ = previousIndex(j); d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); } } // Create the interpolating splines return new BicubicSplineInterpolatingFunction(xval, yval, fval, dFdX, dFdY, d2FdXdY); }
From source file:au.gov.ga.conn4d.utils.BicubicSplineInterpolator.java
/** * {@inheritDoc}//from ww w.j a v a 2s .c o m */ public BicubicSplineInterpolatingFunction interpolate(final double[] xval, final double[] yval, final float[][] fval) throws NoDataException, DimensionMismatchException, NonMonotonicSequenceException, NumberIsTooSmallException { if (xval.length == 0 || yval.length == 0 || fval.length == 0) { throw new NoDataException(); } if (xval.length != fval.length) { throw new DimensionMismatchException(xval.length, fval.length); } MathArrays.checkOrder(xval); MathArrays.checkOrder(yval); final int xLen = xval.length; final int yLen = yval.length; // Samples (first index is y-coordinate, i.e. subarray variable is x) // 0 <= i < xval.length // 0 <= j < yval.length // fX[j][i] = f(xval[i], yval[j]) final double[][] fX = new double[yLen][xLen]; for (int i = 0; i < xLen; i++) { if (fval[i].length != yLen) { throw new DimensionMismatchException(fval[i].length, yLen); } for (int j = 0; j < yLen; j++) { fX[j][i] = fval[i][j]; } } final SplineInterpolator spInterpolator = new SplineInterpolator(); // For each line y[j] (0 <= j < yLen), construct a 1D spline with // respect to variable x final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; for (int j = 0; j < yLen; j++) { ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); } // For each line x[i] (0 <= i < xLen), construct a 1D spline with // respect to variable y generated by array fY_1[i] final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; for (int i = 0; i < xLen; i++) { xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); } // Partial derivatives with respect to x at the grid knots final double[][] dFdX = new double[xLen][yLen]; for (int j = 0; j < yLen; j++) { final UnivariateFunction f = ySplineX[j].derivative(); for (int i = 0; i < xLen; i++) { dFdX[i][j] = f.value(xval[i]); } } // Partial derivatives with respect to y at the grid knots final double[][] dFdY = new double[xLen][yLen]; for (int i = 0; i < xLen; i++) { final UnivariateFunction f = xSplineY[i].derivative(); for (int j = 0; j < yLen; j++) { dFdY[i][j] = f.value(yval[j]); } } // Cross partial derivatives final double[][] d2FdXdY = new double[xLen][yLen]; for (int i = 0; i < xLen; i++) { final int nI = nextIndex(i, xLen); final int pI = previousIndex(i); for (int j = 0; j < yLen; j++) { final int nJ = nextIndex(j, yLen); final int pJ = previousIndex(j); d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); } } // Create the interpolating splines return new BicubicSplineInterpolatingFunction(xval, yval, fval, dFdX, dFdY, d2FdXdY); }
From source file:com.sciaps.utils.ImportExportSpectrumCSV.java
public void exportSpectrumFile(File saveFile, PiecewiseSpectrum spectrum) throws IOException { if (spectrum == null || saveFile == null) { logger_.warn("", "will not save spectrum csv file"); return;/*from w w w .ja v a2s . co m*/ } final UnivariateFunction intensity = spectrum.getIntensityFunction(); BufferedWriter bout = new BufferedWriter(new OutputStreamWriter(new FileOutputStream(saveFile))); try { bout.append("wavelength, intensity"); bout.newLine(); final DoubleRange range = spectrum.getValidRange(); for (double x = range.getMinimumDouble(); x <= range.getMaximumDouble(); x += 1.0 / EXPORT_SAMPLE_RATE) { double y = intensity.value(x); if (Double.isNaN(y)) { y = 0; } bout.append(Double.toString(x)); bout.append(", "); bout.append(Double.toString(y)); bout.newLine(); } } finally { bout.close(); } logger_.info("saved spectrum csv file to " + saveFile.getAbsolutePath()); }
From source file:de.biomedical_imaging.traJ.features.SplineCurveDynamicsFeature.java
@Override public double[] evaluate() { PolynomialSplineFunction spline = null; if (recalculate) { splinefit = new TrajectorySplineFit(t, nSegments); spline = splinefit.calculateSpline(); } else {/*w w w . jav a2s.c om*/ spline = splinefit.getSpline(); } if (!splinefit.wasSuccessfull()) { return new double[] { Double.NaN, Double.NaN, Double.NaN }; } Trajectory tr = splinefit.getRotatedTrajectory(); UnivariateFunction derivative = spline.derivative(); int N = 0; double sumParallel = 0; double sumPerpendicular = 0; //Split each step into replacment rependicular and parallel to spline tangent for (int i = timelag; i < t.size(); i += timelag) { Point2D.Double pRef = splinefit.minDistancePointSpline(new Point2D.Double(tr.get(i).x, tr.get(i).y), 50); Point2D.Double pTangend = new Point2D.Double(pRef.x + 1, derivative.value(pRef.x) * (pRef.x + 1 - pRef.x) + spline.value(pRef.x)); Point2D.Double pNormal = new Point2D.Double(-1 * pTangend.y, pTangend.x); Point2D.Double dp = new Point2D.Double(pRef.x + tr.get(i).x - tr.get(i - timelag).x, pRef.y + tr.get(i).y - tr.get(i - timelag).y); sumParallel += Math.pow(splinefit.distancePointLine(pRef, pNormal, dp), 2); sumPerpendicular += Math.pow(splinefit.distancePointLine(pRef, pTangend, dp), 2); N++; } //System.out.println("N: " +spline.getN()); double msdParallel = sumParallel / N; double msdPerpendicular = sumPerpendicular / N; result = new double[] { msdParallel, msdPerpendicular, N }; recalculate = false; return result; }
From source file:area.math.nm.ui.GUIInterpolationJPanel.java
private void plot_btnActionPerformed(java.awt.event.ActionEvent evt) {//GEN-FIRST:event_plot_btnActionPerformed double xmin = Double.parseDouble(xmin_txt.getText()); double xinc = Double.parseDouble(xinc_txt.getText()); double xmax = Double.parseDouble(xmax_txt.getText()); Arithmetic al = new Arithmetic(); List<Double> xi = new ArrayList<>(); List<Double> yi = new ArrayList<>(); //Extract from table for (int i = 0; i < jTable1.getRowCount(); i++) { if (jTable1.getValueAt(i, 0) != null && jTable1.getValueAt(i, 1) != null) { String cstr = (String) jTable1.getValueAt(i, 0); String estr = (String) jTable1.getValueAt(i, 1); double cf = Double.parseDouble(cstr); double ep = Double.parseDouble(estr); xi.add(cf);/* w w w . j a v a2 s .c o m*/ yi.add(ep); } } double[] Xin = new double[xi.size()]; double[] Yin = new double[xi.size()]; for (int i = 0; i < xi.size(); i++) { Xin[i] = xi.get(i); Yin[i] = yi.get(i); } Interpolation ip = new Interpolation(); UnivariateFunction function; function = ip.splinePlot(Xin, Yin); double x[] = al.linArray(xmin, xinc, xmax); int dim = (int) (((xmax - xmin) / xinc) + 1); double[] y = new double[dim]; for (int i = 0; i < dim; i++) y[i] = function.value(x[i]); Plot p = new Plot(graph_container); p.plot(x, y, "Cubic-Spline-Curve"); }
From source file:gdsc.smlm.function.PoissonGammaGaussianFunction.java
/** * Compute the likelihood//www . jav a 2s.c o m * * @param o * The observed count * @param e * The expected count * @return The likelihood */ public double likelihood(final double o, final double e) { // Use the same variables as the Python code final double cij = o; final double eta = alpha * e; // convert to photons if (sigma == 0) { // No convolution with a Gaussian. Simply evaluate for a Poisson-Gamma distribution. // Any observed count above zero if (cij > 0.0) { // The observed count converted to photons final double nij = alpha * cij; if (eta * nij > 10000) { // Approximate Bessel function i1(x) when using large x: // i1(x) ~ exp(x)/sqrt(2*pi*x) // However the entire equation is logged (creating transform), // evaluated then raised to e to prevent overflow error on // large exp(x) final double transform = 0.5 * Math.log(alpha * eta / cij) - nij - eta + 2 * Math.sqrt(eta * nij) - Math.log(twoSqrtPi * Math.pow(eta * nij, 0.25)); return FastMath.exp(transform); } else { // Second part of equation 135 return Math.sqrt(alpha * eta / cij) * FastMath.exp(-nij - eta) * Bessel.I1(2 * Math.sqrt(eta * nij)); } } else if (cij == 0.0) { return FastMath.exp(-eta); } else { return 0; } } else if (useApproximation) { return mortensenApproximation(cij, eta); } else { // This code is the full evaluation of equation 7 from the supplementary information // of the paper Chao, et al (2013) Nature Methods 10, 335-338. // It is the full evaluation of a Poisson-Gamma-Gaussian convolution PMF. final double sk = sigma; // Read noise final double g = 1.0 / alpha; // Gain final double z = o; // Observed pixel value final double vk = eta; // Expected number of photons // Compute the integral to infinity of: // exp( -((z-u)/(sqrt(2)*s)).^2 - u/g ) * sqrt(vk*u/g) .* besseli(1, 2 * sqrt(vk*u/g)) ./ u; final double vk_g = vk * alpha; // vk / g final double sqrt2sigma = Math.sqrt(2) * sk; // Specify the function to integrate UnivariateFunction f = new UnivariateFunction() { public double value(double u) { return eval(sqrt2sigma, z, vk_g, g, u); } }; // Integrate to infinity is not necessary. The convolution of the function with the // Gaussian should be adequately sampled using a nxSD around the maximum. // Find a bracket containing the maximum double lower, upper; double maxU = Math.max(1, o); double rLower = maxU; double rUpper = maxU + 1; double f1 = f.value(rLower); double f2 = f.value(rUpper); // Calculate the simple integral and the range double sum = f1 + f2; boolean searchUp = f2 > f1; if (searchUp) { while (f2 > f1) { f1 = f2; rUpper += 1; f2 = f.value(rUpper); sum += f2; } maxU = rUpper - 1; } else { // Ensure that u stays above zero while (f1 > f2 && rLower > 1) { f2 = f1; rLower -= 1; f1 = f.value(rLower); sum += f1; } maxU = (rLower > 1) ? rLower + 1 : rLower; } lower = Math.max(0, maxU - 5 * sk); upper = maxU + 5 * sk; if (useSimpleIntegration && lower > 0) { // If we are not at the zero boundary then we can use a simple integration by adding the // remaining points in the range for (double u = rLower - 1; u >= lower; u -= 1) { sum += f.value(u); } for (double u = rUpper + 1; u <= upper; u += 1) { sum += f.value(u); } } else { // Use Legendre-Gauss integrator try { final double relativeAccuracy = 1e-4; final double absoluteAccuracy = 1e-8; final int minimalIterationCount = 3; final int maximalIterationCount = 32; final int integrationPoints = 16; // Use an integrator that does not use the boundary since u=0 is undefined (divide by zero) UnivariateIntegrator i = new IterativeLegendreGaussIntegrator(integrationPoints, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); sum = i.integrate(2000, f, lower, upper); } catch (TooManyEvaluationsException ex) { return mortensenApproximation(cij, eta); } } // Compute the final probability //final double f1 = z / sqrt2sigma; return (FastMath.exp(-vk) / (sqrt2pi * sk)) * (FastMath.exp(-(f1 * f1)) + sum); } }
From source file:com.github.lindenb.jvarkit.tools.redon.CopyNumber01.java
private void normalizeCoverage() { final Median medianOp = new Median(); final Mean meanOp = new Mean(); if (medianOp.evaluate(new double[] { 20, 1000, 19 }) != 20) { throw new RuntimeException("boum"); }// w w w .j a va2s . co m int autosome_count = 0; Collections.sort(this.interval2row, CopyNumber01.sortOnXY); for (int j = 0; j < this.interval2row.size(); ++j) { GCAndDepth r = this.interval2row.get(j); if (isSexualChrom(r.getChrom())) continue; autosome_count++; } double x[] = new double[autosome_count]; double y[] = new double[autosome_count]; int i = 0; for (int j = 0; j < this.interval2row.size(); ++j) { GCAndDepth r = this.interval2row.get(j); if (isSexualChrom(r.getChrom())) continue; x[i] = r.getX(); y[i] = r.getY(); ++i; } final double min_x = x[0]; final double max_x = x[x.length - 1]; /* merge adjacent x having same values */ i = 0; int k = 0; while (i < x.length) { int j = i + 1; while (j < x.length && Double.compare(x[i], x[j]) == 0) { ++j; } x[k] = x[i]; y[k] = meanOp.evaluate(y, i, j - i); ++k; i = j; } /* reduce size of x et y */ if (k != x.length) { info("Compacting X from " + x.length + " to " + k); x = Arrays.copyOf(x, k); y = Arrays.copyOf(y, k); } //min depth cal double min_depth = Double.MAX_VALUE; UnivariateInterpolator interpolator = createInterpolator(); UnivariateFunction spline = interpolator.interpolate(x, y); int points_removed = 0; i = 0; while (i < this.interval2row.size()) { GCAndDepth r = this.interval2row.get(i); if (r.getX() < min_x || r.getX() > max_x) { this.interval2row.remove(i); ++points_removed; } else { double norm = spline.value(r.getX()); if (Double.isNaN(norm) || Double.isInfinite(norm)) { info("NAN " + r); this.interval2row.remove(i); ++points_removed; continue; } r.depth -= norm; min_depth = Math.min(min_depth, r.depth); ++i; } } info("Removed " + points_removed + " because GC% is too small (Sexual chrom)"); spline = null; //fit to min, fill new y for median calculation info("min:" + min_depth); y = new double[this.interval2row.size()]; for (i = 0; i < this.interval2row.size(); ++i) { GCAndDepth gc = this.interval2row.get(i); gc.depth -= min_depth; y[i] = gc.depth; } //normalize on median double median_depth = medianOp.evaluate(y, 0, y.length); info("median:" + median_depth); for (i = 0; i < this.interval2row.size(); ++i) { GCAndDepth gc = this.interval2row.get(i); gc.depth /= median_depth; } //restore genomic order Collections.sort(this.interval2row, CopyNumber01.sortOnPosition); /** smoothing values with neighbours */ final int SMOOTH_WINDOW = 5; y = new double[this.interval2row.size()]; for (i = 0; i < this.interval2row.size(); ++i) { y[i] = this.interval2row.get(i).getY(); } for (i = 0; i < this.interval2row.size(); ++i) { GCAndDepth gc = this.interval2row.get(i); int left = i; int right = i; while (left > 0 && i - left < SMOOTH_WINDOW && this.interval2row.get(left - 1).tid == gc.tid) { left--; } while (right + 1 < this.interval2row.size() && right - i < SMOOTH_WINDOW && this.interval2row.get(right + 1).tid == gc.tid) { right++; } gc.depth = medianOp.evaluate(y, left, (right - left) + 1); } }
From source file:imagingbook.pub.fd.FourierDescriptor.java
/** * Calculates the 'canonical' start point. This version uses * (a) a coarse search for a global maximum of fp() and subsequently * (b) a numerical optimization using Brent's method * (implemented with Apache Commons Math). *///ww w .j ava2s . c o m public double getStartPointPhase(int Mp) { Mp = Math.min(Mp, (G.length - 1) / 2); UnivariateFunction fp = new TargetFunction(Mp); // search for the global maximum in coarse steps double cmax = Double.NEGATIVE_INFINITY; int kmax = -1; int K = 25; // number of steps over 180 degrees for (int k = 0; k < K; k++) { final double phi = Math.PI * k / K; // phase to evaluate final double c = fp.value(phi); if (c > cmax) { cmax = c; kmax = k; } } // optimize using previous and next point as the bracket. double minPhi = Math.PI * (kmax - 1) / K; double maxPhi = Math.PI * (kmax + 1) / K; UnivariateOptimizer optimizer = new BrentOptimizer(1E-4, 1E-6); int maxIter = 20; UnivariatePointValuePair result = optimizer.optimize(new MaxEval(maxIter), new UnivariateObjectiveFunction(fp), GoalType.MAXIMIZE, new SearchInterval(minPhi, maxPhi)); double phi0 = result.getPoint(); return phi0; // the canonical start point phase }
From source file:eu.crisis_economics.abm.fund.FundTest.java
/** * This test is in two parts.//from w ww .j av a2 s . c o m * Part 1: a new fund is created with an initial deposit. This deposit is * reserved as private equity. A collection of household stubs are * created with ample starting deposits. For a number of business * cycles, the households all withdraw or contribute random sums to * the fund. The size of these sums is drawn from a Gaussian * distribution with custom mean and variance. After a fixed number * of business cycles have elapsed, the total amount of cash owned * by each household (be this in fund investments or bank deposits) * is counted. Since the fund has no investment sources and cannot * make a profit, it is expected that no household ever loses or * gains assets (cash plus fund account value). At every step, it is * further asserted that the equity of the fund is not significantly * different to the private (constant) equity. * Part 2: the fund and all households are destroyed and respawned. Households * begin with ample deposits. The test now continues as in Part 1, except * that the fund deposit account is directly injected with exogenous * profits during each business cycle. The size of these profits is * drawn from a Gaussian distribution with custom mean and variance. * At the conclusion of Part 2, the amount of cash in the system * is tallied. It is then asserted that the sum of all household deposits * and fund account values is equal to the initial cash in the system * plus total exogenous profits. * Part 3: as Part 2, except the fund now makes exogenous losses as well as * profits. In this test, the fund is deliberately taken to the edge of * its private equity due to market losses. */ public void testPrivateEquityConstancyAndConservationOfCash() { /* * Test parameters */ final int numberOfBusinessCycles = 1000; final double meanFundProfit = 1.e6, fundProfitVariance = 2.e5, meanHouseholdInvestment = 0., householdInvestmentVariances = 1000., householdInitialDeposits = 1.e7, fundInitialDeposits = 0.; Random dice = new Random(12345); myFund = new MutualFund(myBank, new HomogeneousPortfolioWeighting(), new FundamentalistStockReturnExpectationFunction(), clearingHouse, new NoSmoothingAlgorithm()); // Fresh fund myFund.debit(fundInitialDeposits); final Contract exogenousReturns = myFund.mDepositor.depositAccount; Assert.assertEquals(myFund.getEquity(), fundInitialDeposits); for (int i = 0; i < myHouseholds.length; ++i) { myHouseholds[i] = // Ample household liquidity new HouseholdStub(myFund, myBank, householdInitialDeposits); Assert.assertEquals(myHouseholds[i].getTotalAssets(), householdInitialDeposits, 1.e-5); } /* * Part 1 */ for (int i = 0; i < numberOfBusinessCycles; ++i) { // Random household investments and withdrawals for (int j = 0; j < N_HOUSEHOLDS; ++j) { final double investment = dice.nextGaussian() * householdInvestmentVariances + meanHouseholdInvestment; try { if (investment > 0.) myFund.invest(myHouseholds[j], investment); else if (investment < 0.) { final double maximumWithdrawal = myFund.getBalance(myHouseholds[j]); if (maximumWithdrawal == 0.) continue; myFund.requestWithdrawal(myHouseholds[j], Math.min(-investment, maximumWithdrawal)); } } catch (final InsufficientFundsException noFunds) { Assert.fail(); } } myFund.preClearingProcessing(); // No profits from fund investments. try { myFund.postClearingProcessing(); } catch (final InsufficientFundsException unexpectedException) { Assert.fail(); } System.out.printf("MutualFund assets: %16.10g\n", myFund.getTotalAssets()); System.out.printf("MutualFund liabilities: %16.10g\n", myFund.getTotalLiabilities()); System.out.printf("MutualFund equity: %16.10g\n", myFund.getEquity()); // Assertions final double postClearingEquity = myFund.getEquity(); Assert.assertEquals(postClearingEquity, fundInitialDeposits, 1.e-5); continue; } // Assert no cash has been created System.out.printf("household deposits [open] fund investments [open]" + " deposits [close] fund investments [close] sum [close]\n"); for (int i = 0; i < N_HOUSEHOLDS; ++i) { // Assert no money creation final double closingHouseholdDeposits = myHouseholds[i].bankAccount.getBalance(), closingFundAccountValue = myFund.getBalance(myHouseholds[i]); System.out.printf("%5d %16.10g %16.10g %16.10g %16.10g %16.10g\n", i, householdInitialDeposits, 0., closingHouseholdDeposits, closingFundAccountValue, closingHouseholdDeposits + closingFundAccountValue); Assert.assertEquals(closingHouseholdDeposits + closingFundAccountValue, householdInitialDeposits, 1.e-6); } /* * Part 2 */ // Rerun with exogenous fund profits double totalFundProfits = 0.; for (int i = 0; i < numberOfBusinessCycles; ++i) { // Random household investments and withdrawals for (int j = 0; j < N_HOUSEHOLDS; ++j) { final double investment = dice.nextGaussian() * householdInvestmentVariances + meanHouseholdInvestment; try { if (investment > 0.) myFund.invest(myHouseholds[j], investment); else if (investment < 0.) { final double maximumWithdrawal = myFund.getBalance(myHouseholds[j]); if (maximumWithdrawal == 0.) continue; myFund.requestWithdrawal(myHouseholds[j], Math.min(-investment, maximumWithdrawal)); } } catch (final InsufficientFundsException noFunds) { Assert.fail(); } } myFund.preClearingProcessing(); final double fundProfits = dice.nextGaussian() * fundProfitVariance + meanFundProfit; exogenousReturns.setValue(exogenousReturns.getValue() + fundProfits); totalFundProfits += fundProfits; try { myFund.postClearingProcessing(); } catch (final InsufficientFundsException unexpectedException) { Assert.fail(); } System.out.printf("MutualFund profits: %16.10g\n", fundProfits); System.out.printf("MutualFund assets: %16.10g\n", myFund.getTotalAssets()); System.out.printf("MutualFund liabilities: %16.10g\n", myFund.getTotalLiabilities()); System.out.printf("MutualFund equity: %16.10g\n", myFund.getEquity()); System.out.println("Number of shares: " + myFund.mInvestmentAccount.getNumberOfEmittedShares()); // Assertions final double postClearingEquity = myFund.getEquity(); Assert.assertEquals(postClearingEquity, fundInitialDeposits, 1.e-1); continue; } // Tally the total amount of cash in the system. double totalHouseholdProfits = 0.; System.out.printf("household deposits [open] fund investments [open]" + " deposits [close] fund investments [close] sum [close]\n"); for (int i = 0; i < N_HOUSEHOLDS; ++i) { // Assert no money creation final double closingHouseholdDeposits = myHouseholds[i].bankAccount.getBalance(), closingFundAccountValue = myFund.getBalance(myHouseholds[i]); System.out.printf("%5d %16.10g %16.10g %16.10g %16.10g %16.10g\n", i, householdInitialDeposits, 0., closingHouseholdDeposits, closingFundAccountValue, closingHouseholdDeposits + closingFundAccountValue); totalHouseholdProfits += (closingHouseholdDeposits + closingFundAccountValue) - householdInitialDeposits; } System.out.printf("Aggregate household profits: %16.10g. MutualFund profits: %16.10g\n", totalHouseholdProfits, totalFundProfits); Assert.assertEquals(totalHouseholdProfits / Math.pow(10., (int) Math.log10(totalHouseholdProfits)), totalFundProfits / Math.pow(10., (int) Math.log10(totalFundProfits)), 1.e-12); /* * Part 3 */ // Rerun with predetermined fund losses. UnivariateFunction forcedFundAssetPosition = new UnivariateFunction() { @Override public double value(double time) { return fundInitialDeposits + 1.e-5 + householdInvestmentVariances * (1. + Math.sin(2. * Math.PI * time / 50.)); } }; for (int i = 0; i < numberOfBusinessCycles; ++i) { // Random household investments and withdrawals for (int j = 0; j < N_HOUSEHOLDS; ++j) { final double investment = dice.nextGaussian() * householdInvestmentVariances + meanHouseholdInvestment; try { if (investment > 0.) { myFund.invest(myHouseholds[j], investment); } else if (investment < 0.) { final double maximumWithdrawal = myFund.getBalance(myHouseholds[j]); if (maximumWithdrawal == 0.) continue; double withdrawal = Math.min(-investment, maximumWithdrawal); myFund.requestWithdrawal(myHouseholds[j], withdrawal); } } catch (final InsufficientFundsException noFunds) { Assert.fail(); } } myFund.preClearingProcessing(); final double profitsNow = forcedFundAssetPosition.value(i) - exogenousReturns.getValue(); totalFundProfits += profitsNow; exogenousReturns.setValue(forcedFundAssetPosition.value(i)); try { myFund.postClearingProcessing(); } catch (final InsufficientFundsException unexpectedException) { Assert.fail(); } System.out.printf("MutualFund profits: %16.10g\n", profitsNow); System.out.printf("MutualFund assets: %16.10g\n", myFund.getTotalAssets()); System.out.printf("MutualFund liabilities: %16.10g\n", myFund.getTotalLiabilities()); System.out.printf("MutualFund equity: %16.10g\n", myFund.getEquity()); System.out.println("Number of shares: " + myFund.mInvestmentAccount.getNumberOfEmittedShares()); // Assertions final double postClearingEquity = myFund.getEquity(); Assert.assertEquals(postClearingEquity, fundInitialDeposits, 1.e-1); continue; } // Tally the total amount of cash in the system. totalHouseholdProfits = 0.; System.out.printf("household deposits [open] fund investments [open]" + " deposits [close] fund investments [close] sum [close]\n"); for (int i = 0; i < N_HOUSEHOLDS; ++i) { // Assert no money creation final double closingHouseholdDeposits = myHouseholds[i].bankAccount.getBalance(), closingFundAccountValue = myFund.getBalance(myHouseholds[i]); System.out.printf("%5d %16.10g %16.10g %16.10g %16.10g %16.10g\n", i, householdInitialDeposits, 0., closingHouseholdDeposits, closingFundAccountValue, closingHouseholdDeposits + closingFundAccountValue); totalHouseholdProfits += (closingHouseholdDeposits + closingFundAccountValue) - householdInitialDeposits; } System.out.printf("Aggregate household profits: %16.10g. MutualFund profits: %16.10g\n", totalHouseholdProfits, totalFundProfits); Assert.assertEquals( totalHouseholdProfits / Math.pow(10., (int) Math.log10(Math.abs(totalHouseholdProfits))), totalFundProfits / Math.pow(10., (int) Math.log10(Math.abs(totalFundProfits))), 1.e-12); return; }