List of usage examples for java.util.logging Logger getGlobal
public static final Logger getGlobal()
From source file:fr.inria.atlanmod.instantiator.SpecimenGenerator.java
/** * @param string */ private void log(String string) { Logger.getGlobal().log(Level.INFO, string); }
From source file:fri.cbw.ThermodynamicSimulationEngine.DDE23.java
/** * * @param ddes delay differential equation * @param lags array of time delay values, does not have to be sorted * @param history values of variables for t < 0 * @param t0 time to start the integration, it should be 0 * @param tfinal the end time/*from www.j a v a 2 s . c o m*/ * @return Class representing the result of the integration */ public static IntegrationResult integrate(FirstOrderDelayDifferentialEquations ddes, double[] lags, double[] history, double t0, double tfinal) { IntegrationResult sol = new IntegrationResult(); //% Initialize statistics. int nsteps = 0; int nfailed = 0; int nfevals = 0; if (tfinal <= t0) { throw new IllegalArgumentException("Must have t0 < tfinal."); } double[] temp; temp = new double[history.length]; System.arraycopy(history, 0, temp, 0, history.length); double[] y0; y0 = new double[temp.length]; System.arraycopy(temp, 0, y0, 0, temp.length); int maxlevel = 4; double t = t0; double[] y; y = new double[y0.length]; System.arraycopy(y0, 0, y, 0, y0.length); int neq = ddes.getDimension(); //% If solving a DDE, locate potential discontinuities. We need to step to //% each of the points of potential lack of smoothness. Because we start at //% t0, we can remove it from discont. The solver always steps to tfinal, //% so it is convenient to add it to discont. List<Double> discont = new ArrayList<Double>(); double minlag = Double.POSITIVE_INFINITY; if (lags.length == 0) { discont.add(tfinal); } else { for (int i = 0; i < lags.length; i++) { if (lags[i] < minlag) { minlag = lags[i]; } } if (minlag <= 0) { throw new IllegalArgumentException("The lags must all be positive."); } List<Double> vl = new ArrayList<Double>(); vl.add(t0); for (int level = 1; level < maxlevel; level++) { // it is not <=, comparing discont with matlab List<Double> vlp1 = new ArrayList<Double>(); for (int i = 0; i < vl.size(); i++) { for (int x = 0; x < lags.length; x++) { // should probably be from 1 // probably not vlp1.add(vl.get(i) + lags[x]); } } //% Restrict to tspan. vl.clear(); // FIXME: this looks like a mistake for (double x : vlp1) { if (x <= tfinal) { vl.add(x); } if (vl.isEmpty()) { break; } int nvl = vl.size(); if (nvl > 1) //% Purge duplicates in vl. { vl = purgeDuplicates(vl); } } discont.addAll(vl); // FIXME: make sure this is where it should be } if (discont.size() > 1) { //% Purge duplicates. discont = purgeDuplicates(discont); } } //% Add tfinal to the list of discontinuities if it is not already included. int end = discont.size() - 1; if (Math.abs(tfinal - discont.get(end)) <= 10 * eps * Math.abs(tfinal)) { discont.set(end, tfinal); } else { discont.add(tfinal); } sol.setDiscont(discont); //% Discard t0 and discontinuities in the history. List<Integer> indices = new ArrayList<Integer>(); for (int i = 0; i < discont.size(); i++) { if (discont.get(i) <= t0) { indices.add(i); } } discont.removeAll(indices); int nextdsc = 0; // zero based List<Integer> ndxdsc = new ArrayList<Integer>(); ndxdsc.add(1); //% Get options, and set defaults. double rtol = 1e-3; // if (rtol < 100 * eps) { // rtol = 100 * eps; // warning(['RelTol has been increased to ' num2str(rtol) '.']); // } double atol = 1e-6; double threshold = atol / rtol; //% By default, hmax is 1/10 of the interval of integration. double hmax = Math.min((tfinal - t), Math.abs(0.1 * (tfinal - t))); if (hmax <= 0) { throw new IllegalArgumentException("MaxStep must be greater than zero."); } double minstep = 0; boolean rept_minh = false; boolean printstats = true; // FIXME: for debugging //% Allocate storage for output arrays and initialize them. int chunk = Math.min(100, (int) Math.floor((2 ^ 13) / neq)); double[] tout = new double[chunk]; double[][] yout = new double[chunk][neq]; // flip to C style double[][] ypout = new double[chunk][neq]; int nout = 1; // FIXME: should be 0-based? // probably not tout[nout] = t; for (int k = 0; k < neq; k++) { for (int l = 0; l < chunk; l++) { yout[l][k] = y[k]; // FIXME: think about this } } //% Initialize method parameters. double pow = 1.0 / 3; double[][] B = { // transposed { 1.0 / 2, 0, 0, 0 }, { 0, 3.0 / 4, 0, 0, 0 }, { 2.0 / 9, 1.0 / 3, 4.0 / 9.0, 0 } }; double[] E = { -5.0 / 72, 1.0 / 12, 1.0 / 9, -1.0 / 8 }; double[][] f = new double[4][neq]; // flipped to C indexing //% Evaluate initial history at t0 - lags. double[][] emptyarr = new double[0][]; double[] fakeX = { t }; // FIXME: try to understand what is happening here with t/X double[][] fakeY = { y }; // FIXME: DTTO double[][] Z = lagvals(t, lags, history, fakeX, fakeY, emptyarr); double[] f0 = new double[neq]; ddes.computeDerivatives(t, y, Z, f0); for (int k = 0; k < neq; k++) { ypout[nout][k] = f0[k]; // FIXME: think about this // should be correct now } nfevals = nfevals + 1; //% stats int m = f0.length; if (m != neq) { // FIXME: find better class to throw? throw new IllegalArgumentException("returned derivative and history vectors of different lengths."); } double hmin = 16 * eps * t; //% Compute an initial step size h using y'(t). double h = Math.min(hmax, tfinal - t0); double n = Double.NEGATIVE_INFINITY; // compute norm(xk, inf) for (int k = 0; k < f0.length; k++) { double xk = Math.abs(f0[k] / Math.max(Math.abs(y[k]), threshold)); if (xk > n) { n = xk; } } double rh = n / (0.8 * Math.pow(rtol, pow)); if (h * rh > 1) { h = 1.0 / rh; // NOTE: used to have 1 / rh there, did not work :( } h = Math.max(h, hmin); //% Make sure that the first step is explicit so that the code can //% properly initialize the interpolant. h = Math.min(h, 0.5 * minlag); System.arraycopy(f0, 0, f[0], 0, neq); //% THE MAIN LOOP double[] last_y = new double[neq]; boolean done = false; while (!done) { //% By default, hmin is a small number such that t+hmin is only slightly //% different than t. It might be 0 if t is 0. hmin = 16 * eps * t; h = Math.min(hmax, Math.max(hmin, h)); //% couldn't limit h until new hmin //% Adjust step size to hit discontinuity. tfinal = discont(end). boolean hitdsc = false; double distance = discont.get(nextdsc) - t; if (Math.min(1.1 * h, hmax) >= distance) { //% stretch h = distance; hitdsc = true; } else if (2 * h >= distance) { //% look-ahead h = distance / 2; } if (!hitdsc && (minlag < h) && (h < 2 * minlag)) { h = minlag; } //% LOOP FOR ADVANCING ONE STEP. double tnew; double ynew[] = new double[neq]; boolean itfail; boolean nofailed = true; //% no failed attempts double err; while (true) { double[][] hB = new double[3][4]; // dimensions of B for (int k = 0; k < 3; k++) { for (int l = 0; l < 4; l++) { hB[k][l] = h * B[k][l]; } } double t1 = t + 0.5 * h; double t2 = t + 0.75 * h; tnew = t + h; //% If a lagged argument falls in the current step, we evaluate the //% formula by iteration. Extrapolation is used for the evaluation //% of the history terms in the first iteration and the tnew,ynew, //% ypnew of the current iteration are used in the evaluation of //% these terms in the next iteration. int maxit; if (minlag < h) { maxit = 5; } else { maxit = 1; } // FIXME: maybe it should be implemented as appending to a List double[] X = new double[nout + 1]; // looks like +1 is unavoidable. nout is the largest index, it is 1 based System.arraycopy(tout, 0, X, 0, nout + 1); // so need nout+1 length double[][] Y = new double[nout + 1][neq]; System.arraycopy(yout, 0, Y, 0, nout + 1); double[][] YP = new double[nout + 1][neq]; System.arraycopy(ypout, 0, YP, 0, nout + 1); itfail = false; for (int iter = 1; iter <= maxit; iter++) { //FIXME: 0-based? // no, from 1, maybe <= // try <= double[] vecmul1 = new double[neq]; for (int k = 0; k < neq; k++) { // FIXME: merge the loops, it should be possible // it is not for (int l = 0; l < 4; l++) { vecmul1[k] += f[l][k] * hB[0][l]; } } double yt1[] = new double[neq]; for (int k = 0; k < f[0].length; k++) { // changed from f.length yt1[k] = y[k] + vecmul1[k]; // FIXME: indices? } Z = lagvals(t1, lags, history, X, Y, YP); ddes.computeDerivatives(t1, yt1, Z, f[1]); double[] vecmul2 = new double[neq]; for (int k = 0; k < neq; k++) { // FIXME: merge the loops, it should be possible // it is not for (int l = 0; l < 4; l++) { vecmul2[k] += f[l][k] * hB[1][l]; } } double yt2[] = new double[neq]; for (int k = 0; k < f[0].length; k++) { // changed from f.length yt2[k] = y[k] + vecmul2[k]; } Z = lagvals(t2, lags, history, X, Y, YP); ddes.computeDerivatives(t2, yt2, Z, f[2]); double[] vecmul3 = new double[neq]; for (int k = 0; k < neq; k++) { for (int l = 0; l < 4; l++) { vecmul3[k] += f[l][k] * hB[2][l]; } } for (int k = 0; k < f[0].length; k++) { ynew[k] = y[k] + vecmul3[k]; } Z = lagvals(tnew, lags, history, X, Y, YP); ddes.computeDerivatives(tnew, ynew, Z, f[3]); nfevals = nfevals + 3; if (maxit > 1) { if (iter > 1) { double errit = Double.NEGATIVE_INFINITY; for (int k = 0; k < neq; k++) { // another norm(erritk, inf) // missed this norm the first time double erritk = Math.abs((ynew[k] - last_y[k]) / Math.max(Math.max(Math.abs(y[k]), Math.abs(ynew[k])), threshold)); if (erritk > errit) { errit = erritk; } } if (errit <= 0.1 * rtol) { break; } } //% Use the tentative solution at tnew in the evaluation of the //% history terms of the next iteration. X = Arrays.copyOf(X, X.length + 1); X[X.length - 1] = tnew; Y = Arrays.copyOf(Y, Y.length + 1); Y[Y.length - 1] = ynew; YP = Arrays.copyOf(YP, YP.length + 1); YP[YP.length - 1] = f[3]; System.arraycopy(ynew, 0, last_y, 0, ynew.length); // had switched last_y and ynew itfail = (iter == maxit); } } // FIXME: matrix multiplication? // had bug: first compute n, then do h*n, not the other way double[] vecmul = new double[neq]; for (int l = 0; l < neq; l++) { for (int o = 0; o < E.length; o++) { vecmul[l] += (f[o][l] * E[o]); } } n = Double.NEGATIVE_INFINITY; for (int k = 0; k < neq; k++) { // had bug: infty norm is in abs value // another ona, had multiplication instead of division. took me two hours to realize :( double nk = Math .abs(vecmul[k] / Math.max(Math.max(Math.abs(y[k]), Math.abs(ynew[k])), threshold)); // FIXME: indexing of f if (nk > n) { n = nk; } } //% Estimate the error. err = h * n; //% If h <= minstep, adjust err so that the step will be accepted. //% Note that minstep < minlag, so maxit = 1 and itfail = false. Report //% once that a step of minimum size was taken. if (h <= minstep) { err = Math.min(err, rtol); if (rept_minh) { Logger.getGlobal().log(Level.INFO, "Steps of size MinStep were taken."); rept_minh = false; } } //% Accept the solution only if the weighted error is no more than the //% tolerance rtol. Estimate an h that will yield an error of rtol on //% the next step or the next try at taking this step, as the case may be, //% and use 0.8 of this value to avoid failures. if ((err > rtol) || itfail) { //% Failed step nfailed = nfailed + 1; Logger logger = Logger.getGlobal(); if (h <= hmin) { String msg = String.format("Failure at t=%e. Unable to meet integration " + "tolerances without reducing the step size below " + "the smallest value allowed (%e) at time t.\n", t, hmin); logger.log(Level.WARNING, msg); if (printstats) { logger.log(Level.INFO, "%g successful steps", nsteps); logger.log(Level.INFO, "%g failed attempts", nfailed); logger.log(Level.INFO, "%g function evaluations", nfevals); } //% Trim output arrays, place in solution structure, and return. sol = trim(sol, history, nout, tout, yout, ypout, ndxdsc); return sol; } if (itfail) { h = 0.5 * h; if (h < 2 * minlag) { h = minlag; } } else if (nofailed) { nofailed = false; h = Math.max(hmin, h * Math.max(0.5, 0.8 * Math.pow(rtol / err, pow))); } else { h = Math.max(hmin, 0.5 * h); } hitdsc = false; } else { //% Successful step break; } } nsteps = nsteps + 1; //% stats //% Advance the integration one step. t = tnew; y = ynew; System.arraycopy(f[3], 0, f[0], 0, neq); //% BS(2,3) is FSAL. nout = nout + 1; // reallocate arrays if (nout > tout.length - 1) { tout = Arrays.copyOf(tout, tout.length + chunk); yout = Arrays.copyOf(yout, yout.length + chunk); ypout = Arrays.copyOf(ypout, ypout.length + chunk); // allocate the second dimensions int upto = yout.length; for (int k = yout.length - chunk; k < upto; k++) { yout[k] = new double[neq]; } upto = ypout.length; for (int k = ypout.length - chunk; k < upto; k++) { ypout[k] = new double[neq]; } } tout[nout] = t; for (int k = 0; k < neq; k++) { yout[nout][k] = y[k]; ypout[nout][k] = f[0][k]; } //% If there were no failures, compute a new h. if (nofailed && !itfail) { //% Note that h may shrink by 0.8, and that err may be 0. double temp2 = 1.25 * Math.pow(err / rtol, pow); if (temp2 > 0.2) { h = h / temp2; } else { h = 5 * h; } h = Math.max(h, minstep); //FIXME: NetBeans says value never used } //% Have we hit tfinal = discont(end)? if (hitdsc) { nextdsc = nextdsc + 1; done = nextdsc > discont.size() - 1; if (!done) { ndxdsc.add(nout); } } } //% Successful integration: if (printstats) { if (printstats) { Logger logger = Logger.getGlobal(); logger.log(Level.INFO, "%g successful steps", nsteps); logger.log(Level.INFO, "%g failed attempts", nfailed); logger.log(Level.INFO, "%g function evaluations", nfevals); } } //% Trim output arrays, place in solution structure, and return. sol = trim(sol, history, nout, tout, yout, ypout, ndxdsc); return sol; //% End of function dde23. }