io.fair_acc.math.MathBase/Math in chartfx-math; necessary to have source code of io.fair_acc.math.MathBase/Math to write customized Renderer? #645
-
|
Background: for #641, try to implement In order to avoid compile error like class io.fair_acc.chartfx.renderer.spi.financial.CandleStickWithMALocalMaxMinRenderer cannot be cast to class io.fair_acc.chartfx.renderer.spi.financial.AbstractFinancialRenderer (io.fair_acc.chartfx.renderer.spi.financial.CandleStickWithMALocalMaxMinRenderer and io.fair_acc.chartfx.renderer.spi.financial.AbstractFinancialRenderer are in unnamed module of loader org.codehaus.mojo.exec.URLClassLoaderBuilder$ExecJavaClassLoader @17192ed3)
class io.fair_acc.chartfx.renderer.CandleStickWithMALocalMaxMinRenderer cannot be cast to class io.fair_acc.chartfx.renderer.spi.financial.AbstractFinancialRenderer (io.fair_acc.chartfx.renderer.spi.financial.CandleStickWithMALocalMaxMinRenderer and io.fair_acc.chartfx.renderer.spi.financial.AbstractFinancialRenderer are in unnamed module of loader org.codehaus.mojo.exec.URLClassLoaderBuilder$ExecJavaClassLoader @17192ed3)
git clone -b 11.3.0 git@github.com:fair-acc/chart-fx.git
chartfx-dataset/src/main/java/io/fair_acc/dataset
chartfx-chart/src/main/java/io/fair_acc/chartfx
chartfx-bench/src/main/java/io/fair_acc/bench
chartfx-math/src/main/java/io/fair_acc/math
<dependency>
<groupId>io.fair-acc</groupId>
<artifactId>chartfx</artifactId>
<version>11.3.0</version>
</dependency>Then confront Question: is it necessary to have source code of io.fair_acc.math.MathBase to write customized Renderer? Can recover from Math.class, MathBase.class as following but still debugging. |
Beta Was this translation helpful? Give feedback.
Answered by
yezhengli-Mr9
Nov 21, 2023
Replies: 3 comments
-
package io.fair_acc.math;
//
// Source code recreated from a .class file by IntelliJ IDEA
// (powered by FernFlower decompiler)
//
import java.util.Arrays;
import java.util.Locale;
public class Math extends MathBase {
private static final int kWorkMax = 100;
public Math() {
}
public double kolmogorovProb(double z) {
double[] fj = new double[]{-2.0, -8.0, -18.0, -32.0};
double[] r = new double[4];
double w = 2.50662827;
double c1 = -1.2337005501361697;
double c2 = -11.103304951225528;
double c3 = -30.842513753404244;
double u = abs(z);
double p;
if (u < 0.2) {
p = 1.0;
} else {
double v;
if (u < 0.755) {
v = 1.0 / (u * u);
p = 1.0 - 2.50662827 * (exp(-1.2337005501361697 * v) + exp(-11.103304951225528 * v) + exp(-30.842513753404244 * v)) / u;
} else if (u < 6.8116) {
r[1] = 0.0;
r[2] = 0.0;
r[3] = 0.0;
v = u * u;
int maxj = max(1, nInt(3.0 / u));
for(int j = 0; j < maxj; ++j) {
r[j] = exp(fj[j] * v);
}
p = 2.0 * (r[0] - r[1] + r[2] - r[3]);
} else {
p = 0.0;
}
}
return p;
}
public double kolmogorovTest(int na, double[] a, int nb, double[] b, String option) {
double prob = -1.0;
if (a != null && b != null && na > 2 && nb > 2) {
double rna = (double)na;
double rnb = (double)nb;
double sa = 1.0 / rna;
double sb = 1.0 / rnb;
double rdiff;
int ia;
int ib;
if (a[0] < b[0]) {
rdiff = -sa;
ia = 2;
ib = 1;
} else {
rdiff = sb;
ib = 2;
ia = 1;
}
double rdmax = abs(rdiff);
boolean ok = false;
for(int i = 0; i < na + nb; ++i) {
if (a[ia - 1] < b[ib - 1]) {
rdiff -= sa;
++ia;
if (ia > na) {
ok = true;
break;
}
} else if (a[ia - 1] > b[ib - 1]) {
rdiff += sb;
++ib;
if (ib > nb) {
ok = true;
break;
}
} else {
double x;
for(x = a[ia - 1]; a[ia - 1] == x && ia <= na; ++ia) {
rdiff -= sa;
}
while(b[ib - 1] == x && ib <= nb) {
rdiff += sb;
++ib;
}
if (ia > na) {
ok = true;
break;
}
if (ib > nb) {
ok = true;
break;
}
}
rdmax = max(rdmax, abs(rdiff));
}
if (ok) {
rdmax = max(rdmax, abs(rdiff));
double z = rdmax * sqrt(rna * rnb / (rna + rnb));
prob = this.kolmogorovProb(z);
}
String opt = option.toUpperCase(Locale.UK);
if (opt.contains("D")) {
System.out.println(String.format(" Kolmogorov Probability = %g, Max Dist = %g", prob, rdmax));
}
return opt.contains("M") ? rdmax : prob;
} else {
System.err.println("KolmogorovTest(): Sets must have more than 2 points");
return prob;
}
}
double kOrdStat(int n, double[] a, int k, int[] work) {
boolean isAllocated = false;
int[] workLocal = new int[100];
int[] ind;
if (work != null) {
ind = work;
} else {
ind = workLocal;
if (n > 100) {
isAllocated = true;
ind = new int[n];
}
}
int rk;
for(rk = 0; rk < n; ind[rk] = rk++) {
}
rk = k;
int l = 0;
int ir = n - 1;
int temp;
while(ir > l + 1) {
int mid = l + ir >> 1;
temp = ind[mid];
ind[mid] = ind[l + 1];
ind[l + 1] = temp;
if (a[ind[l]] > a[ind[ir]]) {
temp = ind[l];
ind[l] = ind[ir];
ind[ir] = temp;
}
if (a[ind[l + 1]] > a[ind[ir]]) {
temp = ind[l + 1];
ind[l + 1] = ind[ir];
ind[ir] = temp;
}
if (a[ind[l]] > a[ind[l + 1]]) {
temp = ind[l];
ind[l] = ind[l + 1];
ind[l + 1] = temp;
}
int i = l + 1;
int j = ir;
int arr = ind[l + 1];
while(true) {
do {
++i;
} while(a[ind[i]] < a[arr]);
do {
--j;
} while(a[ind[j]] > a[arr]);
if (j < i) {
ind[l + 1] = ind[j];
ind[j] = arr;
if (j >= rk) {
ir = j - 1;
}
if (j <= rk) {
l = i;
}
break;
}
temp = ind[i];
ind[i] = ind[j];
ind[j] = temp;
}
}
if (ir == l + 1 && a[ind[ir]] < a[ind[l]]) {
temp = ind[l];
ind[l] = ind[ir];
ind[ir] = temp;
}
double tmp = a[ind[rk]];
if (isAllocated) {
// int[] ind = null;
ind = null;
}
return tmp;
}
void quantiles(int n, int nprob, double[] x, double[] quantiles, double[] prob, boolean isSorted, int[] index, int type) {
if (type >= 1 && type <= 9) {
int[] ind = null;
if (!isSorted) {
if (index == null) {
ind = index;
} else {
ind = new int[n];
}
}
double npm = 0.0;
double xj;
double xjj;
int j;
int i;
if (type < 4) {
for(i = 0; i < nprob; ++i) {
npm = (double)n * prob[i];
if (npm < 1.0) {
if (isSorted) {
quantiles[i] = x[0];
} else {
quantiles[i] = this.kOrdStat(n, x, 0, ind);
}
} else {
j = max(floorNint(npm) - 1, 0);
if (npm - (double)j - 1.0 > 1.0E-14) {
if (isSorted) {
quantiles[i] = x[j + 1];
} else {
quantiles[i] = this.kOrdStat(n, x, j + 1, ind);
}
} else {
if (isSorted) {
xj = x[j];
} else {
xj = this.kOrdStat(n, x, j, ind);
}
if (type == 1) {
quantiles[i] = xj;
}
if (type == 2) {
if (isSorted) {
xjj = x[j + 1];
} else {
xjj = this.kOrdStat(n, x, j + 1, ind);
}
quantiles[i] = 0.5 * (xj + xjj);
}
if (type == 3) {
if (!even((long)(j - 1))) {
if (isSorted) {
xjj = x[j + 1];
} else {
xjj = this.kOrdStat(n, x, j + 1, ind);
}
quantiles[i] = xjj;
} else {
quantiles[i] = xj;
}
}
}
}
}
}
if (type > 3) {
for(i = 0; i < nprob; ++i) {
double np = (double)n * prob[i];
if (np < 1.0 && type != 7 && type != 4) {
quantiles[i] = this.kOrdStat(n, x, 0, ind);
} else {
if (type == 4) {
npm = np;
}
if (type == 5) {
npm = np + 0.5;
}
if (type == 6) {
npm = np + prob[i];
}
if (type == 7) {
npm = np - prob[i] + 1.0;
}
if (type == 8) {
npm = np + 0.3333333333333333 * (1.0 + prob[i]);
}
if (type == 9) {
npm = np + 0.25 * prob[i] + 0.375;
}
int intnpm = floorNint(npm);
j = max(intnpm - 1, 0);
double g = npm - (double)intnpm;
if (isSorted) {
xj = x[j];
xjj = x[j + 1];
} else {
xj = this.kOrdStat(n, x, j, ind);
xjj = this.kOrdStat(n, x, j + 1, ind);
}
quantiles[i] = (1.0 - g) * xj + g * xjj;
}
}
}
} else {
System.err.println("illegal value of type");
}
}
public boolean rootsCubic(double[] coef, double[] roots) {
if (coef[3] == 0.0) {
return false;
} else {
boolean complex = false;
double r = coef[2] / coef[3];
double s = coef[1] / coef[3];
double t = coef[0] / coef[3];
double p = s - r * r / 3.0;
double ps3 = p / 3.0;
double q = 2.0 * r * r * r / 27.0 - r * s / 3.0 + t;
double qs2 = q / 2.0;
double ps33 = ps3 * ps3 * ps3;
double d = ps33 + qs2 * qs2;
double tmp;
double y1;
double y2;
double y3;
double a;
double b;
double c;
if (d >= 0.0) {
complex = true;
d = sqrt(d);
double u = -qs2 + d;
double v = -qs2 - d;
tmp = 0.3333333333333333;
double lnu = log(abs(u));
double lnv = log(abs(v));
double su = sign(1.0, u);
double sv = sign(1.0, v);
u = su * exp(tmp * lnu);
v = sv * exp(tmp * lnv);
y1 = u + v;
y2 = -y1 / 2.0;
y3 = (u - v) * sqrt(3.0) / 2.0;
tmp = r / 3.0;
a = y1 - tmp;
b = y2 - tmp;
c = y3;
} else {
ps3 = -ps3;
ps33 = -ps33;
double cphi = -qs2 / sqrt(ps33);
double phi = aCos(cphi);
double phis3 = phi / 3.0;
double pis3 = 1.0471975511965976;
double c1 = cos(phis3);
double c2 = cos(pis3 + phis3);
double c3 = cos(pis3 - phis3);
tmp = sqrt(ps3);
y1 = 2.0 * tmp * c1;
y2 = -2.0 * tmp * c2;
y3 = -2.0 * tmp * c3;
tmp = r / 3.0;
a = y1 - tmp;
b = y2 - tmp;
c = y3 - tmp;
}
roots[0] = a;
roots[1] = b;
roots[2] = c;
return complex;
}
}
public double voigt(double xx, double sigma, double lg, int r) {
if (sigma < 0.0 || lg < 0.0 || sigma == 0.0 && lg == 0.0) {
return 0.0;
} else if (sigma == 0.0) {
return lg * 0.159154943 / (xx * xx + lg * lg / 4.0);
} else if (lg == 0.0) {
return 0.39894228 / sigma * exp(-xx * xx / (2.0 * sigma * sigma));
} else {
double x = xx / sigma / 1.41421356;
double y = lg / 2.0 / sigma / 1.41421356;
int rClamped = r < 2 ? 2 : (r > 5 ? 5 : r);
double r0 = 1.51 * exp(1.144 * (double)rClamped);
double r1 = 1.6 * exp(0.554 * (double)rClamped);
double rrtpi = 0.56418958;
double y0 = 1.5;
double y0py0 = 3.0;
double y0q = 2.25;
double[] c = new double[]{1.0117281, -0.75197147, 0.012557727, 0.010022008, -2.4206814E-4, 5.0084806E-7};
double[] s = new double[]{1.393237, 0.23115241, -0.15535147, 0.0062183662, 9.1908299E-5, -6.2752596E-7};
double[] t = new double[]{0.31424038, 0.94778839, 1.5976826, 2.2795071, 3.020637, 3.8897249};
double a0 = 0.0;
double d0 = 0.0;
double d2 = 0.0;
double e0 = 0.0;
double e2 = 0.0;
double e4 = 0.0;
double h0 = 0.0;
double h2 = 0.0;
double h4 = 0.0;
double h6 = 0.0;
double p0 = 0.0;
double p2 = 0.0;
double p4 = 0.0;
double p6 = 0.0;
double p8 = 0.0;
double z0 = 0.0;
double z2 = 0.0;
double z4 = 0.0;
double z6 = 0.0;
double z8 = 0.0;
double[] xp = new double[6];
double[] xm = new double[6];
double[] yp = new double[6];
double[] ym = new double[6];
double[] mq = new double[6];
double[] pq = new double[6];
double[] mf = new double[6];
double[] pf = new double[6];
boolean rg1 = true;
boolean rg2 = true;
boolean rg3 = true;
double yq = y * y;
double yrrtpi = y * 0.56418958;
double xlim0 = r0 - y;
double xlim1 = r1 - y;
double xlim3 = 3.097 * y - 0.45;
double xlim2 = 6.8 - y;
double xlim4 = 18.1 * y + 1.65;
if (y <= 1.0E-6) {
xlim1 = xlim0;
xlim2 = xlim0;
}
double abx = abs(x);
double xq = abx * abx;
double k;
if (abx > xlim0) {
k = yrrtpi / (xq + yq);
} else {
double d;
if (abx > xlim1) {
if (rg1) {
rg1 = false;
a0 = yq + 0.5;
d0 = a0 * a0;
d2 = yq + yq - 1.0;
}
d = 0.56418958 / (d0 + xq * (d2 + xq));
k = d * y * (a0 + xq);
} else if (abx > xlim2) {
if (rg2) {
h0 = 0.5625 + yq * (4.5 + yq * (10.5 + yq * (6.0 + yq)));
h2 = -4.5 + yq * (9.0 + yq * (6.0 + yq * 4.0));
h4 = 10.5 - yq * (6.0 - yq * 6.0);
h6 = -6.0 + yq * 4.0;
e0 = 1.875 + yq * (8.25 + yq * (5.5 + yq));
e2 = 5.25 + yq * (1.0 + yq * 3.0);
e4 = 0.75 * h6;
}
d = 0.56418958 / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
k = d * y * (e0 + xq * (e2 + xq * (e4 + xq)));
} else if (abx < xlim3) {
if (rg3) {
z0 = 272.1014 + y * (1280.829 + y * (2802.87 + y * (3764.966 + y * (3447.629 + y * (2256.981 + y * (1074.409 + y * (369.1989 + y * (88.26741 + y * (13.3988 + y)))))))));
z2 = 211.678 + y * (902.3066 + y * (1758.336 + y * (2037.31 + y * (1549.675 + y * (793.4273 + y * (266.2987 + y * (53.59518 + y * 5.0)))))));
z4 = 78.86585 + y * (308.1852 + y * (497.3014 + y * (479.2576 + y * (269.2916 + y * (80.39278 + y * 10.0)))));
z6 = 22.03523 + y * (55.02933 + y * (92.75679 + y * (53.59518 + y * 10.0)));
z8 = 1.49646 + y * (13.3988 + y * 5.0);
p0 = 153.5168 + y * (549.3954 + y * (919.4955 + y * (946.897 + y * (662.8097 + y * (328.2151 + y * (115.3772 + y * (27.93941 + y * (4.264678 + y * 0.3183291))))))));
p2 = -34.16955 + y * (-1.322256 + y * (124.5975 + y * (189.773 + y * (139.4665 + y * (56.81652 + y * (12.79458 + y * 1.2733163))))));
p4 = 2.584042 + y * (10.46332 + y * (24.01655 + y * (29.81482 + y * (12.79568 + y * 1.9099744))));
p6 = -0.07272979 + y * (0.9377051 + y * (4.266322 + y * 1.273316));
p8 = 5.480304E-4 + y * 0.3183291;
}
d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
k = d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
} else {
double ypy0 = y + 1.5;
double ypy0q = ypy0 * ypy0;
k = 0.0;
int j;
for(j = 0; j <= 5; ++j) {
d = x - t[j];
mq[j] = d * d;
mf[j] = 1.0 / (mq[j] + ypy0q);
xm[j] = mf[j] * d;
ym[j] = mf[j] * ypy0;
d = x + t[j];
pq[j] = d * d;
pf[j] = 1.0 / (pq[j] + ypy0q);
xp[j] = pf[j] * d;
yp[j] = pf[j] * ypy0;
}
if (abx <= xlim4) {
for(j = 0; j <= 5; ++j) {
k = k + c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
}
} else {
double yf = y + 3.0;
for(j = 0; j <= 5; ++j) {
k = k + (c[j] * (mq[j] * mf[j] - 1.5 * ym[j]) + s[j] * yf * xm[j]) / (mq[j] + 2.25) + (c[j] * (pq[j] * pf[j] - 1.5 * yp[j]) - s[j] * yf * xp[j]) / (pq[j] + 2.25);
}
k = y * k + exp(-xq);
}
}
}
return k / 2.506628 / sigma;
}
}
public static double besselI(int n, double x) {
if (n < 0) {
System.err.println("BesselI(): *I* Invalid argument(s) (n,x) = (" + n + ", " + x + ")");
return 0.0;
} else if (n == 0) {
return besselI0(x);
} else if (n == 1) {
return besselI1(x);
} else if (x == 0.0) {
return 0.0;
} else {
double kBigPositive = 1.0E10;
if (abs(x) > 1.0E10) {
return 0.0;
} else {
boolean iacc = true;
double kBigNegative = 1.0E-10;
double tox = 2.0 / abs(x);
double bip = 0.0;
double bim = 0.0;
double bi = 1.0;
double result = 0.0;
int m = 2 * (n + (int)sqrt((double)(40 * n)));
for(int j = m; j >= 1; --j) {
bim = bip + (double)j * tox * bi;
bip = bi;
bi = bim;
if (abs(bim) > 1.0E10) {
result *= 1.0E-10;
bi = bim * 1.0E-10;
bip *= 1.0E-10;
}
if (j == n) {
result = bip;
}
}
result *= besselI0(x) / bi;
if (x < 0.0 && n % 2 == 1) {
result = -result;
}
return result;
}
}
}
public static double besselI0(double x) {
double p1 = 1.0;
double p2 = 3.5156229;
double p3 = 3.0899424;
double p4 = 1.2067492;
double p5 = 0.2659732;
double p6 = 0.0360768;
double p7 = 0.0045813;
double q1 = 0.39894228;
double q2 = 0.01328592;
double q3 = 0.00225319;
double q4 = -0.00157565;
double q5 = 0.00916281;
double q6 = -0.02057706;
double q7 = 0.02635537;
double q8 = -0.01647633;
double q9 = 0.00392377;
double k1 = 3.75;
double ax = abs(x);
double y = 0.0;
double result = 0.0;
if (ax < 3.75) {
double xx = x / 3.75;
y = xx * xx;
result = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.0360768 + y * 0.0045813)))));
} else {
y = 3.75 / ax;
result = exp(ax) / sqrt(ax) * (0.39894228 + y * (0.01328592 + y * (0.00225319 + y * (-0.00157565 + y * (0.00916281 + y * (-0.02057706 + y * (0.02635537 + y * (-0.01647633 + y * 0.00392377))))))));
}
return result;
}
public static double besselI1(double x) {
double p1 = 0.5;
double p2 = 0.87890594;
double p3 = 0.51498869;
double p4 = 0.15084934;
double p5 = 0.02658733;
double p6 = 0.00301532;
double p7 = 3.2411E-4;
double q1 = 0.39894228;
double q2 = -0.03988024;
double q3 = -0.00362018;
double q4 = 0.00163801;
double q5 = -0.01031555;
double q6 = 0.02282967;
double q7 = -0.02895312;
double q8 = 0.01787654;
double q9 = -0.00420059;
double k1 = 3.75;
double ax = abs(x);
double y = 0.0;
double result = 0.0;
if (ax < 3.75) {
double xx = x / 3.75;
y = xx * xx;
result = x * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934 + y * (0.02658733 + y * (0.00301532 + y * 3.2411E-4))))));
} else {
y = 3.75 / ax;
result = exp(ax) / sqrt(ax) * (0.39894228 + y * (-0.03988024 + y * (-0.00362018 + y * (0.00163801 + y * (-0.01031555 + y * (0.02282967 + y * (-0.02895312 + y * (0.01787654 + y * -0.00420059))))))));
if (x < 0.0) {
result = -result;
}
}
return result;
}
public static double besselJ0(double x) {
double p1 = 5.7568490574E10;
double p2 = -1.3362590354E10;
double p3 = 6.516196407E8;
double p4 = -1.121442418E7;
double p5 = 77392.33017;
double p6 = -184.9052456;
double p7 = 5.7568490411E10;
double p8 = 1.029532985E9;
double p9 = 9494680.718;
double p10 = 59272.64853;
double p11 = 267.8532712;
double q1 = 0.785398164;
double q2 = -0.001098628627;
double q3 = 2.734510407E-5;
double q4 = -2.073370639E-6;
double q5 = 2.093887211E-7;
double q6 = -0.01562499995;
double q7 = 1.430488765E-4;
double q8 = -6.911147651E-6;
double q9 = 7.621095161E-7;
double q10 = 9.34935152E-8;
double q11 = 0.636619772;
double ax;
double y;
double result;
double result1;
double result2;
if ((ax = abs(x)) < 8.0) {
y = x * x;
result1 = 5.7568490574E10 + y * (-1.3362590354E10 + y * (6.516196407E8 + y * (-1.121442418E7 + y * (77392.33017 + y * -184.9052456))));
result2 = 5.7568490411E10 + y * (1.029532985E9 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y))));
result = result1 / result2;
} else {
double z = 8.0 / ax;
y = z * z;
double xx = ax - 0.785398164;
result1 = 1.0 + y * (-0.001098628627 + y * (2.734510407E-5 + y * (-2.073370639E-6 + y * 2.093887211E-7)));
result2 = -0.01562499995 + y * (1.430488765E-4 + y * (-6.911147651E-6 + y * (7.621095161E-7 - y * 9.34935152E-8)));
result = sqrt(0.636619772 / ax) * (cos(xx) * result1 - z * sin(xx) * result2);
}
return result;
}
public static double besselJ1(double x) {
double p1 = 7.2362614232E10;
double p2 = -7.895059235E9;
double p3 = 2.423968531E8;
double p4 = -2972611.439;
double p5 = 15704.4826;
double p6 = -30.16036606;
double p7 = 1.44725228442E11;
double p8 = 2.300535178E9;
double p9 = 1.858330474E7;
double p10 = 99447.43394;
double p11 = 376.9991397;
double q1 = 2.356194491;
double q2 = 0.00183105;
double q3 = -3.516396496E-5;
double q4 = 2.457520174E-6;
double q5 = -2.40337019E-7;
double q6 = 0.04687499995;
double q7 = -2.002690873E-4;
double q8 = 8.449199096E-6;
double q9 = -8.8228987E-7;
double q10 = 1.05787412E-7;
double q11 = 0.636619772;
double ax;
double y;
double result;
double result1;
double result2;
if ((ax = abs(x)) < 8.0) {
y = x * x;
result1 = x * (7.2362614232E10 + y * (-7.895059235E9 + y * (2.423968531E8 + y * (-2972611.439 + y * (15704.4826 + y * -30.16036606)))));
result2 = 1.44725228442E11 + y * (2.300535178E9 + y * (1.858330474E7 + y * (99447.43394 + y * (376.9991397 + y))));
result = result1 / result2;
} else {
double z = 8.0 / ax;
y = z * z;
double xx = ax - 2.356194491;
result1 = 1.0 + y * (0.00183105 + y * (-3.516396496E-5 + y * (2.457520174E-6 + y * -2.40337019E-7)));
result2 = 0.04687499995 + y * (-2.002690873E-4 + y * (8.449199096E-6 + y * (-8.8228987E-7 + y * 1.05787412E-7)));
result = sqrt(0.636619772 / ax) * (cos(xx) * result1 - z * sin(xx) * result2);
if (x < 0.0) {
result = -result;
}
}
return result;
}
public static double besselK(int n, double x) {
if (!(x <= 0.0) && n >= 0) {
if (n == 0) {
return besselK0(x);
} else if (n == 1) {
return besselK1(x);
} else {
double tox = 2.0 / x;
double bkm = besselK0(x);
double bk = besselK1(x);
double bkp = 0.0;
for(int j = 1; j < n; ++j) {
bkp = bkm + (double)j * tox * bk;
bkm = bk;
bk = bkp;
}
return bk;
}
} else {
System.err.println("BesselK(): *K* Invalid argument(s) (n,x) = (" + n + ", " + x + ")");
return 0.0;
}
}
public static double besselK0(double x) {
if (x <= 0.0) {
System.err.println("BesselK0(): *K0* Invalid argument x = " + x);
return 0.0;
} else {
double p1 = -0.57721566;
double p2 = 0.4227842;
double p3 = 0.23069756;
double p4 = 0.0348859;
double p5 = 0.00262698;
double p6 = 1.075E-4;
double p7 = 7.4E-6;
double q1 = 1.25331414;
double q2 = -0.07832358;
double q3 = 0.02189568;
double q4 = -0.01062446;
double q5 = 0.00587872;
double q6 = -0.0025154;
double q7 = 5.3208E-4;
double y = 0.0;
double result = 0.0;
if (x <= 2.0) {
y = x * x / 4.0;
result = -log(x / 2.0) * besselI0(x) + -0.57721566 + y * (0.4227842 + y * (0.23069756 + y * (0.0348859 + y * (0.00262698 + y * (1.075E-4 + y * 7.4E-6)))));
} else {
y = 2.0 / x;
result = exp(-x) / sqrt(x) * (1.25331414 + y * (-0.07832358 + y * (0.02189568 + y * (-0.01062446 + y * (0.00587872 + y * (-0.0025154 + y * 5.3208E-4))))));
}
return result;
}
}
public static double besselK1(double x) {
if (x <= 0.0) {
System.err.println("BesselK1(): *K1* Invalid argument x = " + x);
return 0.0;
} else {
double p1 = 1.0;
double p2 = 0.15443144;
double p3 = -0.67278579;
double p4 = -0.18156897;
double p5 = -0.01919402;
double p6 = -0.00110404;
double p7 = -4.686E-5;
double q1 = 1.25331414;
double q2 = 0.23498619;
double q3 = -0.0365562;
double q4 = 0.01504268;
double q5 = -0.00780353;
double q6 = 0.00325614;
double q7 = -6.8245E-4;
double y = 0.0;
double result = 0.0;
if (x <= 2.0) {
y = x * x / 4.0;
result = log(x / 2.0) * besselI1(x) + 1.0 / x * (1.0 + y * (0.15443144 + y * (-0.67278579 + y * (-0.18156897 + y * (-0.01919402 + y * (-0.00110404 + y * -4.686E-5))))));
} else {
y = 2.0 / x;
result = exp(-x) / sqrt(x) * (1.25331414 + y * (0.23498619 + y * (-0.0365562 + y * (0.01504268 + y * (-0.00780353 + y * (0.00325614 + y * -6.8245E-4))))));
}
return result;
}
}
public static double besselY0(double x) {
double p1 = -2.957821389E9;
double p2 = 7.062834065E9;
double p3 = -5.123598036E8;
double p4 = 1.087988129E7;
double p5 = -86327.92757;
double p6 = 228.4622733;
double p7 = 4.0076544269E10;
double p8 = 7.452499648E8;
double p9 = 7189466.438;
double p10 = 47447.2647;
double p11 = 226.1030244;
double p12 = 0.636619772;
double q1 = 0.785398164;
double q2 = -0.001098628627;
double q3 = 2.734510407E-5;
double q4 = -2.073370639E-6;
double q5 = 2.093887211E-7;
double q6 = -0.01562499995;
double q7 = 1.430488765E-4;
double q8 = -6.911147651E-6;
double q9 = 7.621095161E-7;
double q10 = -9.34945152E-8;
double q11 = 0.636619772;
double y;
double result;
double result1;
double result2;
if (x < 8.0) {
y = x * x;
result1 = -2.957821389E9 + y * (7.062834065E9 + y * (-5.123598036E8 + y * (1.087988129E7 + y * (-86327.92757 + y * 228.4622733))));
result2 = 4.0076544269E10 + y * (7.452499648E8 + y * (7189466.438 + y * (47447.2647 + y * (226.1030244 + y))));
result = result1 / result2 + 0.636619772 * besselJ0(x) * log(x);
} else {
double z = 8.0 / x;
y = z * z;
double xx = x - 0.785398164;
result1 = 1.0 + y * (-0.001098628627 + y * (2.734510407E-5 + y * (-2.073370639E-6 + y * 2.093887211E-7)));
result2 = -0.01562499995 + y * (1.430488765E-4 + y * (-6.911147651E-6 + y * (7.621095161E-7 + y * -9.34945152E-8)));
result = sqrt(0.636619772 / x) * (sin(xx) * result1 + z * cos(xx) * result2);
}
return result;
}
public static double besselY1(double x) {
double p1 = -4.900604943E12;
double p2 = 1.27527439E12;
double p3 = -5.153438139E10;
double p4 = 7.349264551E8;
double p5 = -4237922.726;
double p6 = 8511.937935;
double p7 = 2.49958057E13;
double p8 = 4.244419664E11;
double p9 = 3.733650367E9;
double p10 = 2.245904002E7;
double p11 = 102042.605;
double p12 = 354.9632885;
double p13 = 0.636619772;
double q1 = 2.356194491;
double q2 = 0.00183105;
double q3 = -3.516396496E-5;
double q4 = 2.457520174E-6;
double q5 = -2.40337019E-7;
double q6 = 0.04687499995;
double q7 = -2.002690873E-4;
double q8 = 8.449199096E-6;
double q9 = -8.8228987E-7;
double q10 = 1.05787412E-7;
double q11 = 0.636619772;
double y;
double result;
double result1;
double result2;
if (x < 8.0) {
y = x * x;
result1 = x * (-4.900604943E12 + y * (1.27527439E12 + y * (-5.153438139E10 + y * (7.349264551E8 + y * (-4237922.726 + y * 8511.937935)))));
result2 = 2.49958057E13 + y * (4.244419664E11 + y * (3.733650367E9 + y * (2.245904002E7 + y * (102042.605 + y * (354.9632885 + y)))));
result = result1 / result2 + 0.636619772 * (besselJ1(x) * log(x) - 1.0 / x);
} else {
double z = 8.0 / x;
y = z * z;
double xx = x - 2.356194491;
result1 = 1.0 + y * (0.00183105 + y * (-3.516396496E-5 + y * (2.457520174E-6 + y * -2.40337019E-7)));
result2 = 0.04687499995 + y * (-2.002690873E-4 + y * (8.449199096E-6 + y * (-8.8228987E-7 + y * 1.05787412E-7)));
result = sqrt(0.636619772 / x) * (sin(xx) * result1 + z * cos(xx) * result2);
}
return result;
}
public static double beta(double p, double q) {
return exp(lnGamma(p) + lnGamma(q) - lnGamma(p + q));
}
public static double betaCf(double x, double a, double b) {
boolean itmax = true;
double eps = 3.0E-14;
double fpmin = 1.0E-30;
double qab = a + b;
double qap = a + 1.0;
double qam = a - 1.0;
double c = 1.0;
double d = 1.0 - qab * x / qap;
if (abs(d) < 1.0E-30) {
d = 1.0E-30;
}
d = 1.0 / d;
double h = d;
int m;
for(m = 1; m <= 500; ++m) {
int m2 = m * 2;
double aa = (double)m * (b - (double)m) * x / ((qam + (double)m2) * (a + (double)m2));
d = 1.0 + aa * d;
if (abs(d) < 1.0E-30) {
d = 1.0E-30;
}
c = 1.0 + aa / c;
if (abs(c) < 1.0E-30) {
c = 1.0E-30;
}
d = 1.0 / d;
h *= d * c;
aa = -(a + (double)m) * (qab + (double)m) * x / ((a + (double)m2) * (qap + (double)m2));
d = 1.0 + aa * d;
if (abs(d) < 1.0E-30) {
d = 1.0E-30;
}
c = 1.0 + aa / c;
if (abs(c) < 1.0E-30) {
c = 1.0E-30;
}
d = 1.0 / d;
double del = d * c;
h *= del;
if (abs(del - 1.0) <= 3.0E-14) {
break;
}
}
if (m > 500) {
System.err.printf("BetaCf: a or b too big, or itmax too small, a=%e, b=%e, x=%e, h=%e, itmax=%e", a, b, x, h, 500);
}
return h;
}
public static double betaDist(double x, double p, double q) {
if (!(x < 0.0) && !(x > 1.0) && !(p <= 0.0) && !(q <= 0.0)) {
double beta = beta(p, q);
double r = pow(x, p - 1.0) * pow(1.0 - x, q - 1.0) / beta;
return r;
} else {
System.err.println("BetaDist(): - parameter value outside allowed range");
return 0.0;
}
}
public static double betaDistI(double x, double p, double q) {
if (!(x < 0.0) && !(x > 1.0) && !(p <= 0.0) && !(q <= 0.0)) {
double betai = betaIncomplete(x, p, q);
return betai;
} else {
System.err.println("BetaDistI(): parameter value outside allowed range");
return 0.0;
}
}
public static double betaIncomplete(double x, double a, double b) {
if (!(x < 0.0) && !(x > 1.0)) {
double bt;
if (x != 0.0 && x != 1.0) {
bt = pow(x, a) * pow(1.0 - x, b) / beta(a, b);
} else {
bt = 0.0;
}
return x < (a + 1.0) / (a + b + 2.0) ? bt * betaCf(x, a, b) / a : 1.0 - bt * betaCf(1.0 - x, b, a) / b;
} else {
System.err.println("BetaIncomplete(): X must between 0 and 1");
return 0.0;
}
}
public static long binarySearch(double[] array, double value) {
return binarySearch(array, 0, array.length, value);
}
public static long binarySearch(double[] array, int offset, int length, double value) {
if (array != null && array.length > 0) {
int nabove = length + offset + 1;
int nbelow = offset;
while(nabove - nbelow > 1) {
int middle = (nabove + nbelow) / 2;
if (value == array[middle - 1]) {
return (long)(middle - 1);
}
if (value < array[middle - 1]) {
nabove = middle;
} else {
nbelow = middle;
}
}
return (long)(nbelow - 1);
} else {
return -1L;
}
}
public static long binarySearch(float[] array, float value) {
return binarySearch(array, 0, array.length, value);
}
public static long binarySearch(float[] array, int offset, int length, float value) {
if (array != null && array.length > 0) {
int nabove = length + offset + 1;
int nbelow = offset;
while(nabove - nbelow > 1) {
int middle = (nabove + nbelow) / 2;
if (value == array[middle - 1]) {
return (long)(middle - 1);
}
if (value < array[middle - 1]) {
nabove = middle;
} else {
nbelow = middle;
}
}
return (long)(nbelow - 1);
} else {
return -1L;
}
}
public static long binarySearch(int[] array, int value) {
return binarySearch((int[])array, 0, array.length, (int)value);
}
public static long binarySearch(int[] array, int offset, int length, int value) {
if (array != null && array.length > 0) {
int nabove = length + offset + 1;
int nbelow = offset;
while(nabove - nbelow > 1) {
int middle = (nabove + nbelow) / 2;
if (value == array[middle - 1]) {
return (long)(middle - 1);
}
if (value < array[middle - 1]) {
nabove = middle;
} else {
nbelow = middle;
}
}
return (long)(nbelow - 1);
} else {
return -1L;
}
}
public static long binarySearch(long[] array, long value) {
return binarySearch(array, 0, array.length, value);
}
public static long binarySearch(long[] array, int offset, int length, long value) {
if (array != null && array.length > 0) {
int nabove = length + offset + 1;
int nbelow = offset;
while(nabove - nbelow > 1) {
int middle = (nabove + nbelow) / 2;
if (value == array[middle - 1]) {
return (long)(middle - 1);
}
if (value < array[middle - 1]) {
nabove = middle;
} else {
nbelow = middle;
}
}
return (long)(nbelow - 1);
} else {
return -1L;
}
}
public static long binarySearch(short[] array, short value) {
return binarySearch((short[])array, 0, array.length, (short)value);
}
public static long binarySearch(short[] array, int offset, int length, short value) {
if (array != null && array.length > 0) {
int nabove = length + offset + 1;
int nbelow = offset;
while(nabove - nbelow > 1) {
int middle = (nabove + nbelow) / 2;
if (value == array[middle - 1]) {
return (long)(middle - 1);
}
if (value < array[middle - 1]) {
nabove = middle;
} else {
nbelow = middle;
}
}
return (long)(nbelow - 1);
} else {
return -1L;
}
}
public static double binomial(int n, int k) {
if (k != 0 && n != k) {
if (n > 0 && k >= 0 && n >= k) {
int k1 = min(k, n - k);
int k2 = n - k1;
double fact = (double)(k2 + 1);
for(int i = k1; i > 1; --i) {
fact *= (double)(k2 + i) / (double)i;
}
return fact;
} else {
return 0.0;
}
} else {
return 1.0;
}
}
public static double binomialI(double p, int n, int k) {
if (k <= 0) {
return 1.0;
} else if (k > n) {
return 0.0;
} else {
return k == n ? pow(p, (double)n) : betaIncomplete(p, (double)k, (double)(n - k + 1));
}
}
public static double breitWigner(double x, double mean, double gamma) {
double bw = gamma / ((x - mean) * (x - mean) + gamma * gamma / 4.0);
return bw / 6.283185307179586;
}
public static double cauchyDist(double x, double t, double s) {
double temp = (x - t) * (x - t) / (s * s);
double result = 1.0 / (s * java.lang.Math.PI * (1.0 + temp));
return result;
}
public static double chisquareQuantile(double p, double ndf) {
if (ndf <= 0.0) {
return 0.0;
} else {
double[] c = new double[]{0.0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2, 4.67, 6.66, 6.73, 13.32, 60.0, 70.0, 84.0, 105.0, 120.0, 127.0, 140.0, 175.0, 210.0, 252.0, 264.0, 294.0, 346.0, 420.0, 462.0, 606.0, 672.0, 707.0, 735.0, 889.0, 932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0, 2520.0, 5040.0};
double e = 5.0E-7;
double aa = 0.6931471806;
double g = lnGamma(0.5 * ndf);
double xx = 0.5 * ndf;
double cp = xx - 1.0;
double ch;
double p1;
double p2;
double q;
double t;
double a;
if (ndf >= log(p) * -c[5]) {
if (ndf > c[3]) {
double x = normQuantile(p);
p1 = c[2] / ndf;
ch = ndf * pow(x * sqrt(p1) + 1.0 - p1, 3.0);
if (ch > c[6] * ndf + 6.0) {
ch = -2.0 * (log(1.0 - p) - cp * log(0.5 * ch) + g);
}
} else {
ch = c[4];
a = log(1.0 - p);
do {
q = ch;
p1 = 1.0 + ch * (c[7] + ch);
p2 = ch * (c[9] + ch * (c[8] + ch));
t = -0.5 + (c[7] + 2.0 * ch) / p1 - (c[9] + ch * (c[10] + 3.0 * ch)) / p2;
ch -= (1.0 - exp(a + g + 0.5 * ch + cp * 0.6931471806) * p2 / p1) / t;
} while(abs(q / ch - 1.0) > c[1]);
}
} else {
ch = pow(p * xx * exp(g + xx * 0.6931471806), 1.0 / xx);
if (ch < 5.0E-7) {
return ch;
}
}
boolean maxit = true;
for(int i = 0; i < 20; ++i) {
q = ch;
p1 = 0.5 * ch;
p2 = p - gamma(xx, p1);
t = p2 * exp(xx * 0.6931471806 + g + p1 - cp * log(ch));
double b = t / ch;
a = 0.5 * t - b * cp;
double s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] + c[11] * a))))) / c[24];
double s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
double s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
double s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
double s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
double s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
ch += t * (1.0 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
if (abs(q / ch - 1.0) > 5.0E-7) {
break;
}
}
return ch;
}
}
public static double correlationCoefficient(double[] x, double[] y, double[] w) {
int n = x.length;
if (y.length != n) {
throw new IllegalArgumentException("x and y array lengths must be equal");
} else if (w.length != n) {
throw new IllegalArgumentException("x and weight array lengths must be equal");
} else {
double sxy = covariance(x, y, w);
double sx = variance(x, w);
double sy = variance(y, w);
return sxy / sqrt(sx * sy);
}
}
public static double covariance(double[] xx, double[] yy, double[] ww) {
int n = xx.length;
if (n != yy.length) {
throw new IllegalArgumentException("length of x variable array, " + n + " and length of y array, " + yy.length + " are different");
} else if (n != ww.length) {
throw new IllegalArgumentException("length of x variable array, " + n + " and length of weight array, " + yy.length + " are different");
} else {
double nn = effectiveSampleNumber(ww);
double nterm = nn / (nn - 1.0);
double sumx = 0.0;
double sumy = 0.0;
double sumw = 0.0;
double[] weight = invertAndSquare(ww);
for(int i = 0; i < n; ++i) {
sumx += xx[i] * weight[i];
sumy += yy[i] * weight[i];
sumw += weight[i];
}
double meanx = sumx / sumw;
double meany = sumy / sumw;
double sum = 0.0;
for(int i = 0; i < n; ++i) {
sum += weight[i] * (xx[i] - meanx) * (yy[i] - meany);
}
return sum * nterm / sumw;
}
}
public static double[] cross(double[] v1, double[] v2, double[] out) {
out[0] = v1[1] * v2[2] - v1[2] * v2[1];
out[1] = v1[2] * v2[0] - v1[0] * v2[2];
out[2] = v1[0] * v2[1] - v1[1] * v2[0];
return out;
}
public static float[] cross(float[] v1, float[] v2, float[] out) {
out[0] = v1[1] * v2[2] - v1[2] * v2[1];
out[1] = v1[2] * v2[0] - v1[0] * v2[2];
out[2] = v1[0] * v2[1] - v1[1] * v2[0];
return out;
}
public static double diLog(double x) {
double hf = 0.5;
double pi = java.lang.Math.PI;
double pi2 = 9.869604401089358;
double pi3 = 3.289868133696453;
double pi6 = 1.6449340668482264;
double pi12 = 0.8224670334241132;
double[] c = new double[]{0.429966935608137, 0.4097598753307711, -0.01858843665014592, 0.00145751084062268, -1.430418444234E-4, 1.58841554188E-5, -1.90784959387E-6, 2.4195180854E-7, -3.193341274E-8, 4.34545063E-9, -6.057848E-10, 8.612098E-11, -1.244332E-11, 1.82256E-12, -2.7007E-13, 4.042E-14, -6.1E-15, 9.3E-16, -1.4E-16, 2.0E-17};
double b0 = 0.0;
double h;
if (x == 1.0) {
h = 1.6449340668482264;
} else if (x == -1.0) {
h = -0.8224670334241132;
} else {
double t = -x;
double y;
double s;
double a;
double b1;
double b2;
if (t <= -2.0) {
y = -1.0 / (1.0 + t);
s = 1.0;
b1 = log(-t);
b2 = log(1.0 + 1.0 / t);
a = -3.289868133696453 + 0.5 * (b1 * b1 - b2 * b2);
} else if (t < -1.0) {
y = -1.0 - t;
s = -1.0;
a = log(-t);
a = -1.6449340668482264 + a * (a + log(1.0 + 1.0 / t));
} else if (t <= -0.5) {
y = -(1.0 + t) / t;
s = 1.0;
a = log(-t);
a = -1.6449340668482264 + a * (-0.5 * a + log(1.0 + t));
} else if (t < 0.0) {
y = -t / (1.0 + t);
s = -1.0;
b1 = log(1.0 + t);
a = 0.5 * b1 * b1;
} else if (t <= 1.0) {
y = t;
s = 1.0;
a = 0.0;
} else {
y = 1.0 / t;
s = -1.0;
b1 = log(t);
a = 1.6449340668482264 + 0.5 * b1 * b1;
}
h = y + y - 1.0;
double alfa = h + h;
b1 = 0.0;
b2 = 0.0;
for(int i = 19; i >= 0; --i) {
b0 = c[i] + alfa * b1 - b2;
b2 = b1;
b1 = b0;
}
h = -(s * (b0 - h * b2) + a);
}
return h;
}
public static double effectiveSampleNumber(double[] ww) {
double[] weight = (double[])ww.clone();
int n;
for(n = 0; n < weight.length; ++n) {
weight[n] = 1.0 / MathBase.sqr(weight[n]);
}
n = weight.length;
double sum2w = 0.0;
double sumw2 = 0.0;
for(int i = 0; i < n; ++i) {
sum2w += weight[i];
sumw2 += weight[i] * weight[i];
}
sum2w *= sum2w;
double nEff = sum2w / sumw2;
return nEff;
}
public static double erf(double x) {
return 1.0 - erfc(x);
}
public static double erfc(double x) {
double v = 1.0;
double z = abs(x);
if (z <= 0.0) {
return v;
} else {
double a1 = -1.26551223;
double a2 = 1.00002368;
double a3 = 0.37409196;
double a4 = 0.09678418;
double a5 = -0.18628806;
double a6 = 0.27886807;
double a7 = -1.13520398;
double a8 = 1.48851587;
double a9 = -0.82215223;
double a10 = 0.17087277;
double t = 1.0 / (1.0 + 0.5 * z);
v = t * exp(-z * z + -1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277)))))))));
if (x < 0.0) {
v = 2.0 - v;
}
return v;
}
}
public static double erfInverse(double x) {
double kEps = 1.0E-14;
double kConst = 0.8862269254527579;
if (abs(x) <= 1.0E-14) {
return 0.8862269254527579 * x;
} else {
boolean kMaxit = true;
if (abs(x) < 1.0) {
double erfi = 0.8862269254527579 * abs(x);
double y0 = erf(0.9 * erfi);
double derfi = 0.1 * erfi;
for(int iter = 0; iter < 50; ++iter) {
double y1 = 1.0 - erfc(erfi);
double dy1 = abs(x) - y1;
if (abs(dy1) < 1.0E-14) {
if (x < 0.0) {
return -erfi;
}
return erfi;
}
double dy0 = y1 - y0;
derfi *= dy1 / dy0;
y0 = y1;
erfi += derfi;
if (abs(derfi / erfi) < 1.0E-14) {
if (x < 0.0) {
return -erfi;
}
return erfi;
}
}
}
return Double.NaN;
}
}
public static double factorial(int n) {
if (n <= 0) {
return 1.0;
} else {
double x = 1.0;
int b = 0;
do {
++b;
x *= (double)b;
} while(b != n);
return x;
}
}
public static double fDist(double F, double N, double M) {
if (!(F < 0.0) && !(N < 1.0) && !(M < 1.0)) {
double denom = gamma(N / 2.0) * gamma(M / 2.0) * pow(M + N * F, (N + M) / 2.0);
double div = gamma((N + M) / 2.0) * pow(N, N / 2.0) * pow(M, M / 2.0) * pow(F, 0.5 * N - 1.0);
return div / denom;
} else {
return 0.0;
}
}
public static double fDistI(double F, double N, double M) {
double fi = 1.0 - betaIncomplete(M / (M + N * F), M * 0.5, N * 0.5);
return fi;
}
public static double freq(double x) {
double c1 = 0.5641895835477563;
double w2 = 1.4142135623730951;
double p10 = 242.66795523053176;
double q10 = 215.0588758698612;
double p11 = 21.979261618294153;
double q11 = 91.1649054045149;
double p12 = 6.996383488619135;
double q12 = 15.082797630407788;
double p13 = -0.035609843701815386;
double q13 = 1.0;
double p20 = 300.4592610201616;
double q20 = 300.4592609569833;
double p21 = 451.9189537118729;
double q21 = 790.9509253278981;
double p22 = 339.3208167343437;
double q22 = 931.3540948506096;
double p23 = 152.9892850469404;
double q23 = 638.9802644656312;
double p24 = 43.162227222056735;
double q24 = 277.58544474398764;
double p25 = 7.2117582508830935;
double q25 = 77.00015293522948;
double p26 = 0.564195517478974;
double q26 = 12.782727319629423;
double p27 = -1.368648573827167E-7;
double q27 = 1.0;
double p30 = -0.002996107077035422;
double q30 = 0.010620923052846792;
double p31 = -0.04947309106232507;
double q31 = 0.19130892610782985;
double p32 = -0.22695659353968692;
double q32 = 1.051675107067932;
double p33 = -0.2786613086096478;
double q33 = 1.9873320181713525;
double p34 = -0.02231924597341847;
double q34 = 1.0;
double v = abs(x) / 1.4142135623730951;
double vv = v * v;
double ap;
double aq;
double h;
double hc;
if (v < 0.5) {
ap = -0.035609843701815386;
aq = 1.0;
ap = 6.996383488619135 + vv * ap;
ap = 21.979261618294153 + vv * ap;
ap = 242.66795523053176 + vv * ap;
aq = 15.082797630407788 + vv * aq;
aq = 91.1649054045149 + vv * aq;
aq = 215.0588758698612 + vv * aq;
h = v * ap / aq;
hc = 1.0 - h;
} else if (v < 4.0) {
ap = -1.368648573827167E-7;
aq = 1.0;
ap = 0.564195517478974 + v * ap;
ap = 7.2117582508830935 + v * ap;
ap = 43.162227222056735 + v * ap;
ap = 152.9892850469404 + v * ap;
ap = 339.3208167343437 + v * ap;
ap = 451.9189537118729 + v * ap;
ap = 300.4592610201616 + v * ap;
aq = 12.782727319629423 + v * aq;
aq = 77.00015293522948 + v * aq;
aq = 277.58544474398764 + v * aq;
aq = 638.9802644656312 + v * aq;
aq = 931.3540948506096 + v * aq;
aq = 790.9509253278981 + v * aq;
aq = 300.4592609569833 + v * aq;
hc = exp(-vv) * ap / aq;
h = 1.0 - hc;
} else {
double y = 1.0 / vv;
ap = -0.02231924597341847;
aq = 1.0;
ap = -0.2786613086096478 + y * ap;
ap = -0.22695659353968692 + y * ap;
ap = -0.04947309106232507 + y * ap;
ap = -0.002996107077035422 + y * ap;
aq = 1.9873320181713525 + y * aq;
aq = 1.051675107067932 + y * aq;
aq = 0.19130892610782985 + y * aq;
aq = 0.010620923052846792 + y * aq;
hc = exp(-vv) * (0.5641895835477563 + y * ap / aq) / v;
h = 1.0 - hc;
}
return x > 0.0 ? 0.5 + 0.5 * h : 0.5 * hc;
}
public static double gamCf(double a, double x) {
if (!(a <= 0.0) && !(x <= 0.0)) {
boolean itmax = true;
double eps = 3.0E-14;
double fpmin = 1.0E-30;
double gln = lnGamma(a);
double b = x + 1.0 - a;
double c = 9.999999999999999E29;
double d = 1.0 / b;
double h = d;
for(int i = 1; i <= 100; ++i) {
double an = (double)(-i) * ((double)i - a);
b += 2.0;
d = an * d + b;
if (abs(d) < 1.0E-30) {
d = 1.0E-30;
}
c = b + an / c;
if (abs(c) < 1.0E-30) {
c = 1.0E-30;
}
d = 1.0 / d;
double del = d * c;
h *= del;
if (abs(del - 1.0) < 3.0E-14) {
break;
}
}
double v = exp(-x + a * log(x) - gln) * h;
return 1.0 - v;
} else {
return 0.0;
}
}
public static double gamma(double z) {
if (z <= 0.0) {
return 0.0;
} else {
double v = lnGamma(z);
return exp(v);
}
}
public static double gamma(double a, double x) {
if (!(a <= 0.0) && !(x <= 0.0)) {
return x < a + 1.0 ? gamSer(a, x) : gamCf(a, x);
} else {
return 0.0;
}
}
public static double gammaDist(double x, double gamma, double mu, double beta) {
if (!(x < mu) && !(gamma <= 0.0) && !(beta <= 0.0)) {
double temp = (x - mu) / beta;
double temp2 = beta * gamma(gamma);
double result = pow(temp, gamma - 1.0) * exp(-temp) / temp2;
return result;
} else {
System.err.println("GammaDist(): illegal parameter values");
return 0.0;
}
}
public static double gamSer(double a, double x) {
if (!(a <= 0.0) && !(x <= 0.0)) {
boolean itmax = true;
double eps = 3.0E-14;
double gln = lnGamma(a);
double ap = a;
double sum = 1.0 / a;
double del = sum;
for(int n = 1; n <= 100; ++n) {
++ap;
del = del * x / ap;
sum += del;
if (abs(del) < abs(sum * 3.0E-14)) {
break;
}
}
double v = sum * exp(-x + a * log(x) - gln);
return v;
} else {
return 0.0;
}
}
public static double gauss(double x, double mean, double sigma, boolean norm) {
if (sigma == 0.0) {
return 1.0E30;
} else {
double arg = (x - mean) / sigma;
double res = exp(-0.5 * arg * arg);
return !norm ? res : res / (2.5066282746310002 * sigma);
}
}
public static double geometricMean(double[] a, int offset, int length) {
if (a != null && a.length > 0) {
double logsum = 0.0;
for(int i = offset; i < length + offset; ++i) {
if (a[i] == 0.0) {
return 0.0;
}
double absa = abs(a[i]);
logsum += log(absa);
}
return exp(logsum / (double)length);
} else {
return -1.0;
}
}
public static double geometricMean(double[] data) {
return geometricMean((double[])data, 0, data.length);
}
public static float geometricMean(float[] a, int offset, int length) {
if (a != null && a.length > 0) {
double logsum = 0.0;
for(int i = offset; i < length + offset; ++i) {
if (a[i] == 0.0F) {
return 0.0F;
}
double absa = (double)abs(a[i]);
logsum += log(absa);
}
return (float)exp(logsum / (double)length);
} else {
return -1.0F;
}
}
public static float geometricMean(float[] data) {
return geometricMean((float[])data, 0, data.length);
}
public static int geometricMean(int[] a, int offset, int length) {
if (a != null && a.length > 0) {
double logsum = 0.0;
for(int i = offset; i < length + offset; ++i) {
if (a[i] == 0) {
return 0;
}
double absa = (double)abs(a[i]);
logsum += log(absa);
}
return (int)exp(logsum / (double)length);
} else {
return -1;
}
}
public static int geometricMean(int[] data) {
return geometricMean((int[])data, 0, data.length);
}
public static long geometricMean(long[] a, int offset, int length) {
if (a != null && a.length > 0) {
double logsum = 0.0;
for(int i = offset; i < length + offset; ++i) {
if (a[i] == 0L) {
return 0L;
}
double absa = (double)abs(a[i]);
logsum += log(absa);
}
return (long)exp(logsum / (double)length);
} else {
return -1L;
}
}
public static long geometricMean(long[] data) {
return geometricMean((long[])data, 0, data.length);
}
public static short geometricMean(short[] a, int offset, int length) {
if (a != null && a.length > 0) {
double logsum = 0.0;
for(int i = offset; i < length + offset; ++i) {
if (a[i] == 0) {
return 0;
}
double absa = (double)abs(a[i]);
logsum += log(absa);
}
return (short)((int)exp(logsum / (double)length));
} else {
return -1;
}
}
public static short geometricMean(short[] data) {
return geometricMean((short[])data, 0, data.length);
}
private static double[] invertAndSquare(double[] values) {
double[] ret = new double[values.length];
for(int i = 0; i < values.length; ++i) {
double val = values[i];
if (val == 0.0) {
ret[i] = Double.POSITIVE_INFINITY;
} else if (Double.isNaN(val)) {
ret[i] = Double.NaN;
} else if (Double.isInfinite(val)) {
ret[i] = 0.0;
} else {
ret[i] = 1.0 / MathBase.pow(val, 2.0);
}
}
return ret;
}
public static boolean isInside(double xp, double yp, int np, double[] x, double[] y) {
int inter = 0;
for(int i = 0; i < np; ++i) {
double xn;
double yn;
if (i < np - 1) {
xn = x[i + 1];
yn = y[i + 1];
} else {
xn = x[0];
yn = y[0];
}
if (y[i] != yn && (!(yp <= y[i]) || !(yp <= yn)) && (!(y[i] < yp) || !(yn < yp))) {
double xint = x[i] + (yp - y[i]) * (xn - x[i]) / (yn - y[i]);
if (xp < xint) {
++inter;
}
}
}
return inter % 2 > 0;
}
public static boolean isInside(float xp, float yp, int np, float[] x, float[] y) {
int inter = 0;
for(int i = 0; i < np; ++i) {
float xn;
float yn;
if (i < np - 1) {
xn = x[i + 1];
yn = y[i + 1];
} else {
xn = x[0];
yn = y[0];
}
if (y[i] != yn && (!(yp <= y[i]) || !(yp <= yn)) && (!(y[i] < yp) || !(yn < yp))) {
double xint = (double)(x[i] + (yp - y[i]) * (xn - x[i]) / (yn - y[i]));
if ((double)xp < xint) {
++inter;
}
}
}
return inter % 2 > 0;
}
public static boolean isInside(int xp, int yp, int np, int[] x, int[] y) {
int inter = 0;
for(int i = 0; i < np; ++i) {
int xn;
int yn;
if (i < np - 1) {
xn = x[i + 1];
yn = y[i + 1];
} else {
xn = x[0];
yn = y[0];
}
if (y[i] != yn && (yp > y[i] || yp > yn) && (y[i] >= yp || yn >= yp)) {
double xint = (double)x[i] + (double)((yp - y[i]) * (xn - x[i])) / (double)(yn - y[i]);
if ((double)xp < xint) {
++inter;
}
}
}
return inter % 2 > 0;
}
public static double landau(double x, double mpv, double sigma, boolean norm) {
if (sigma <= 0.0) {
return 0.0;
} else {
double[] p1 = new double[]{0.4259894875, -0.124976255, 0.039842437, -0.006298287635, 0.001511162253};
double[] q1 = new double[]{1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
double[] p2 = new double[]{0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 1.283617211E-4};
double[] q2 = new double[]{1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
double[] p3 = new double[]{0.1788544503, 0.09359161662, 0.006325387654, 6.611667319E-5, -2.031049101E-6};
double[] q3 = new double[]{1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
double[] p4 = new double[]{0.9874054407, 118.6723273, 849.279436, -743.7792444, 427.0262186};
double[] q4 = new double[]{1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511};
double[] p5 = new double[]{1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.9491};
double[] q5 = new double[]{1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357};
double[] p6 = new double[]{1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
double[] q6 = new double[]{1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939};
double[] a1 = new double[]{0.04166666667, -0.01996527778, 0.02709538966};
double[] a2 = new double[]{-1.84556867, -4.284640743};
double v = (x - mpv) / sigma;
double u;
double den;
if (v < -5.5) {
u = exp(v + 1.0);
if (u < 1.0E-10) {
return 0.0;
}
double ue = exp(-1.0 / u);
double us = sqrt(u);
den = 0.3989422803 * (ue / us) * (1.0 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
} else if (v < -1.0) {
u = exp(-v - 1.0);
den = exp(-u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
} else if (v < 1.0) {
den = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
} else if (v < 5.0) {
den = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
} else if (v < 12.0) {
u = 1.0 / v;
den = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
} else if (v < 50.0) {
u = 1.0 / v;
den = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
} else if (v < 300.0) {
u = 1.0 / v;
den = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
} else {
u = 1.0 / (v - v * log(v) / (v + 1.0));
den = u * u * (1.0 + (a2[0] + a2[1] * u) * u);
}
return !norm ? den : den / sigma;
}
}
public static double landauI(double x) {
double[] p1 = new double[]{0.2514091491, -0.06250580444, 0.0145838123, -0.002108817737, 7.41124729E-4};
double[] q1 = new double[]{1.0, -0.005571175625, 0.06225310236, -0.003137378427, 0.001931496439};
double[] p2 = new double[]{0.2868328584, 0.3564363231, 0.1523518695, 0.02251304883};
double[] q2 = new double[]{1.0, 0.6191136137, 0.1720721448, 0.02278594771};
double[] p3 = new double[]{0.2868329066, 0.3003828436, 0.09950951941, 0.008733827185};
double[] q3 = new double[]{1.0, 0.4237190502, 0.1095631512, 0.008693851567};
double[] p4 = new double[]{1.00035163, 4.503592498, 10.8588388, 7.536052269};
double[] q4 = new double[]{1.0, 5.539969678, 19.33581111, 27.21321508};
double[] p5 = new double[]{1.000006517, 49.09414111, 85.05544753, 153.2153455};
double[] q5 = new double[]{1.0, 50.09928881, 139.9819104, 420.0002909};
double[] p6 = new double[]{1.000000983, 132.9868456, 916.2149244, -960.5054274};
double[] q6 = new double[]{1.0, 133.9887843, 1055.990413, 553.2224619};
double[] a1 = new double[]{0.0, -0.4583333333, 0.6675347222, -1.641741416};
double[] a2 = new double[]{0.0, 1.0, -0.4227843351, -2.043403138};
double u;
double lan;
if (x < -5.5) {
u = exp(x + 1.0);
lan = 0.3989422803 * exp(-1.0 / u) * sqrt(u) * (1.0 + (a1[1] + (a1[2] + a1[3] * u) * u) * u);
} else if (x < -1.0) {
u = exp(-x - 1.0);
lan = exp(-u) / sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * x) * x) * x) * x) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * x) * x) * x) * x);
} else if (x < 1.0) {
lan = (p2[0] + (p2[1] + (p2[2] + p2[3] * x) * x) * x) / (q2[0] + (q2[1] + (q2[2] + q2[3] * x) * x) * x);
} else if (x < 4.0) {
lan = (p3[0] + (p3[1] + (p3[2] + p3[3] * x) * x) * x) / (q3[0] + (q3[1] + (q3[2] + q3[3] * x) * x) * x);
} else if (x < 12.0) {
u = 1.0 / x;
lan = (p4[0] + (p4[1] + (p4[2] + p4[3] * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + q4[3] * u) * u) * u);
} else if (x < 50.0) {
u = 1.0 / x;
lan = (p5[0] + (p5[1] + (p5[2] + p5[3] * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + q5[3] * u) * u) * u);
} else if (x < 300.0) {
u = 1.0 / x;
lan = (p6[0] + (p6[1] + (p6[2] + p6[3] * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + q6[3] * u) * u) * u);
} else {
u = 1.0 / (x - x * log(x) / (x + 1.0));
lan = 1.0 - (a2[1] + (a2[2] + a2[3] * u) * u) * u;
}
return lan;
}
public static double laplaceDist(double x, double alpha, double beta) {
double temp = exp(-abs((x - alpha) / beta));
temp /= 2.0 * beta;
return temp;
}
public static double laplaceDistI(double x, double alpha, double beta) {
double temp;
if (x <= alpha) {
temp = 0.5 * exp(-abs((x - alpha) / beta));
} else {
temp = 1.0 - 0.5 * exp(-abs((x - alpha) / beta));
}
return temp;
}
public static double lnGamma(double z) {
if (z <= 0.0) {
return 0.0;
} else {
double[] c = new double[]{2.5066282746310007, 76.18009172947146, -86.50532032941678, 24.01409824083091, -1.231739572450155, 0.001208650973866179, -5.395239384953E-6};
double y = z;
double tmp = z + 5.5;
tmp = (z + 0.5) * log(tmp) - tmp;
double ser = 1.000000000190015;
for(int i = 1; i < 7; ++i) {
++y;
ser += c[i] / y;
}
return tmp + log(c[0] * ser / z);
}
}
public static long locationMaximum(double[] a, int length) {
if (a != null && a.length > 0) {
double xmax = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmax < a[i]) {
xmax = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMinimum(double[] a, int length) {
if (a != null && a.length > 0) {
double xmin = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmin > a[i]) {
xmin = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMaximum(float[] a, int length) {
if (a != null && a.length > 0) {
float xmax = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmax < a[i]) {
xmax = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMinimum(float[] a, int length) {
if (a != null && a.length > 0) {
float xmin = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmin > a[i]) {
xmin = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMaximum(int[] a, int length) {
if (a != null && a.length > 0) {
int xmax = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmax < a[i]) {
xmax = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMinimum(int[] a, int length) {
if (a != null && a.length > 0) {
int xmin = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmin > a[i]) {
xmin = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMaximum(long[] a, int length) {
if (a != null && a.length > 0) {
long xmax = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmax < a[i]) {
xmax = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMinimum(long[] a, int length) {
if (a != null && a.length > 0) {
long xmin = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmin > a[i]) {
xmin = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMaximum(short[] a, int length) {
if (a != null && a.length > 0) {
short xmax = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmax < a[i]) {
xmax = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static long locationMinimum(short[] a, int length) {
if (a != null && a.length > 0) {
short xmin = a[0];
long location = 0L;
for(int i = 1; i < length; ++i) {
if (xmin > a[i]) {
xmin = a[i];
location = (long)i;
}
}
return location;
} else {
return -1L;
}
}
public static double logNormal(double x, double sigma, double theta, double m) {
if (!(x < theta) && !(sigma <= 0.0) && !(m <= 0.0)) {
double templog2 = log((x - theta) / m) * log((x - theta) / m);
double temp1 = exp(-templog2 / (2.0 * sigma * sigma));
double temp2 = (x - theta) * sigma * sqrt(6.283185307179586);
return temp1 / temp2;
} else {
System.err.println("Lognormal(): illegal parameter values");
return 0.0;
}
}
public static double maximum(double[] data) {
return maximum(data, data.length);
}
public static double maximum(double[] data, int length) {
double val = -1.7976931348623157E308;
for(int i = 0; i < length; ++i) {
val = max(val, data[i]);
}
return val;
}
public static double mean(double[] data) {
return mean(data, data.length);
}
public static double mean(double[] data, int length) {
double norm = 1.0 / (double)length;
double val = 0.0;
for(int i = 0; i < length; ++i) {
val += norm * data[i];
}
return val;
}
public static double median(double[] data) {
return median(data, data.length);
}
public static double median(double[] data, int length) {
double[] temp = sort(data, length, false);
return length % 2 == 0 ? 0.5 * (temp[length / 2] + temp[length / 2 + 1]) : temp[length / 2];
}
public static double minimum(double[] data) {
return minimum(data, data.length);
}
public static double minimum(double[] data, int length) {
double val = Double.MAX_VALUE;
for(int i = 0; i < length; ++i) {
val = min(val, data[i]);
}
return val;
}
public static float maximum(float[] data) {
return maximum(data, data.length);
}
public static float maximum(float[] data, int length) {
float val = -3.4028235E38F;
for(int i = 0; i < length; ++i) {
val = max(val, data[i]);
}
return val;
}
public static float mean(float[] data) {
return mean(data, data.length);
}
public static float mean(float[] data, int length) {
float norm = 1.0F / (float)length;
float val = 0.0F;
for(int i = 0; i < length; ++i) {
val += norm * data[i];
}
return val;
}
public static float median(float[] data) {
return median(data, data.length);
}
public static float median(float[] data, int length) {
float[] temp = sort(data, length, false);
return length % 2 == 0 ? (float)(0.5 * (double)(temp[length / 2] + temp[length / 2 + 1])) : temp[length / 2];
}
public static float minimum(float[] data) {
return minimum(data, data.length);
}
public static float minimum(float[] data, int length) {
float val = Float.MAX_VALUE;
for(int i = 0; i < length; ++i) {
val = min(val, data[i]);
}
return val;
}
public static int maximum(int[] data) {
return maximum(data, data.length);
}
public static int maximum(int[] data, int length) {
int val = -2147483647;
for(int i = 0; i < length; ++i) {
val = max(val, data[i]);
}
return val;
}
public static int mean(int[] data) {
return mean(data, data.length);
}
public static int mean(int[] data, int length) {
double norm = 1.0 / (double)length;
double val = 0.0;
for(int i = 0; i < length; ++i) {
val += norm * (double)data[i];
}
return (int)val;
}
public static int median(int[] data) {
return median(data, data.length);
}
public static int median(int[] data, int length) {
int[] temp = sort(data, length, false);
return length % 2 == 0 ? (int)(0.5 * (double)(temp[length / 2] + temp[length / 2 + 1])) : temp[length / 2];
}
public static int minimum(int[] data) {
return minimum(data, data.length);
}
public static int minimum(int[] data, int length) {
int val = Integer.MAX_VALUE;
for(int i = 0; i < length; ++i) {
val = min(val, data[i]);
}
return val;
}
public static long maximum(long[] data) {
return maximum(data, data.length);
}
public static long maximum(long[] data, int length) {
long val = -9223372036854775807L;
for(int i = 0; i < length; ++i) {
val = max(val, data[i]);
}
return val;
}
public static long mean(long[] data) {
return mean(data, data.length);
}
public static long mean(long[] data, int length) {
double norm = 1.0 / (double)length;
double val = 0.0;
for(int i = 0; i < length; ++i) {
val += norm * (double)data[i];
}
return (long)val;
}
public static long median(long[] data) {
return median(data, data.length);
}
public static long median(long[] data, int length) {
long[] temp = sort(data, length, false);
return length % 2 == 0 ? (long)(0.5 * (double)(temp[length / 2] + temp[length / 2 + 1])) : temp[length / 2];
}
public static long minimum(long[] data) {
return minimum(data, data.length);
}
public static long minimum(long[] data, int length) {
long val = Long.MAX_VALUE;
for(int i = 0; i < length; ++i) {
val = min(val, data[i]);
}
return val;
}
public static short maximum(short[] data) {
return maximum(data, data.length);
}
public static short maximum(short[] data, int length) {
short val = -32767;
for(int i = 0; i < length; ++i) {
val = max(val, data[i]);
}
return val;
}
public static short mean(short[] data) {
return mean(data, data.length);
}
public static short mean(short[] data, int length) {
double norm = 1.0 / (double)length;
double val = 0.0;
for(int i = 0; i < length; ++i) {
val += norm * (double)data[i];
}
return (short)((int)val);
}
public static short median(short[] data) {
return median(data, data.length);
}
public static short median(short[] data, int length) {
short[] temp = sort(data, length, false);
return length % 2 == 0 ? (short)((int)(0.5 * (double)(temp[length / 2] + temp[length / 2 + 1]))) : temp[length / 2];
}
public static short minimum(short[] data) {
return minimum(data, data.length);
}
public static short minimum(short[] data, int length) {
short val = 32767;
for(int i = 0; i < length; ++i) {
val = min(val, data[i]);
}
return val;
}
public static double[] normal2Plane(double[] p1, double[] p2, double[] p3, double[] normal) {
double[] v1 = new double[3];
double[] v2 = new double[3];
v1[0] = p2[0] - p1[0];
v1[1] = p2[1] - p1[1];
v1[2] = p2[2] - p1[2];
v2[0] = p3[0] - p1[0];
v2[1] = p3[1] - p1[1];
v2[2] = p3[2] - p1[2];
normCross(v1, v2, normal);
return normal;
}
public static float[] normal2Plane(float[] p1, float[] p2, float[] p3, float[] normal) {
float[] v1 = new float[3];
float[] v2 = new float[3];
v1[0] = p2[0] - p1[0];
v1[1] = p2[1] - p1[1];
v1[2] = p2[2] - p1[2];
v2[0] = p3[0] - p1[0];
v2[1] = p3[1] - p1[1];
v2[2] = p3[2] - p1[2];
normCross(v1, v2, normal);
return normal;
}
public static double normalize(double[] v) {
double av0 = abs(v[0]);
double av1 = abs(v[1]);
double av2 = abs(v[2]);
double amax;
double foo;
double bar;
if (av0 >= av1 && av0 >= av2) {
amax = av0;
foo = av1;
bar = av2;
} else if (av1 >= av0 && av1 >= av2) {
amax = av1;
foo = av0;
bar = av2;
} else {
amax = av2;
foo = av0;
bar = av1;
}
if (amax == 0.0) {
return 0.0;
} else {
double foofrac = foo / amax;
double barfrac = bar / amax;
double d = amax * sqrt(1.0 + foofrac * foofrac + barfrac * barfrac);
v[0] /= d;
v[1] /= d;
v[2] /= d;
return d;
}
}
public static float normalize(float[] v) {
float d = (float)sqrt((double)(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]));
if (d != 0.0F) {
v[0] /= d;
v[1] /= d;
v[2] /= d;
}
return d;
}
public static double normCross(double[] v1, double[] v2, double[] out) {
return normalize(cross(v1, v2, out));
}
public static float normCross(float[] v1, float[] v2, float[] out) {
return normalize(cross(v1, v2, out));
}
public static double normQuantile(double p) {
if (!(p <= 0.0) && !(p >= 1.0)) {
double a0 = 3.3871328727963665;
double a1 = 133.14166789178438;
double a2 = 1971.5909503065513;
double a3 = 13731.69376550946;
double a4 = 45921.95393154987;
double a5 = 67265.7709270087;
double a6 = 33430.57558358813;
double a7 = 2509.0809287301227;
double b1 = 42.31333070160091;
double b2 = 687.1870074920579;
double b3 = 5394.196021424751;
double b4 = 21213.794301586597;
double b5 = 39307.89580009271;
double b6 = 28729.085735721943;
double b7 = 5226.495278852854;
double c0 = 1.4234371107496835;
double c1 = 4.630337846156546;
double c2 = 5.769497221460691;
double c3 = 3.6478483247632045;
double c4 = 1.2704582524523684;
double c5 = 0.2417807251774506;
double c6 = 0.022723844989269184;
double c7 = 7.745450142783414E-4;
double d1 = 2.053191626637759;
double d2 = 1.6763848301838038;
double d3 = 0.6897673349851;
double d4 = 0.14810397642748008;
double d5 = 0.015198666563616457;
double d6 = 5.475938084995345E-4;
double d7 = 1.0507500716444169E-9;
double e0 = 6.657904643501103;
double e1 = 5.463784911164114;
double e2 = 1.7848265399172913;
double e3 = 0.29656057182850487;
double e4 = 0.026532189526576124;
double e5 = 0.0012426609473880784;
double e6 = 2.7115555687434876E-5;
double e7 = 2.0103343992922881E-7;
double f1 = 0.599832206555888;
double f2 = 0.1369298809227358;
double f3 = 0.014875361290850615;
double f4 = 7.868691311456133E-4;
double f5 = 1.8463183175100548E-5;
double f6 = 1.421511758316446E-7;
double f7 = 2.0442631033899397E-15;
double split1 = 0.425;
double split2 = 5.0;
double konst1 = 0.180625;
double konst2 = 1.6;
double q = p - 0.5;
double r;
double quantile;
if (abs(q) < 0.425) {
r = 0.180625 - q * q;
quantile = q * (((((((2509.0809287301227 * r + 33430.57558358813) * r + 67265.7709270087) * r + 45921.95393154987) * r + 13731.69376550946) * r + 1971.5909503065513) * r + 133.14166789178438) * r + 3.3871328727963665) / (((((((5226.495278852854 * r + 28729.085735721943) * r + 39307.89580009271) * r + 21213.794301586597) * r + 5394.196021424751) * r + 687.1870074920579) * r + 42.31333070160091) * r + 1.0);
} else {
if (q < 0.0) {
r = p;
} else {
r = 1.0 - p;
}
if (r <= 0.0) {
quantile = 0.0;
} else {
r = sqrt(-log(r));
if (r <= 5.0) {
--r;
quantile = (((((((7.745450142783414E-4 * r + 0.022723844989269184) * r + 0.2417807251774506) * r + 1.2704582524523684) * r + 3.6478483247632045) * r + 5.769497221460691) * r + 4.630337846156546) * r + 1.4234371107496835) / (((((((1.0507500716444169E-9 * r + 5.475938084995345E-4) * r + 0.015198666563616457) * r + 0.14810397642748008) * r + 0.6897673349851) * r + 1.6763848301838038) * r + 2.053191626637759) * r + 1.0);
} else {
r -= 5.0;
quantile = (((((((2.0103343992922881E-7 * r + 2.7115555687434876E-5) * r + 0.0012426609473880784) * r + 0.026532189526576124) * r + 0.29656057182850487) * r + 1.7848265399172913) * r + 5.463784911164114) * r + 6.657904643501103) / (((((((2.0442631033899397E-15 * r + 1.421511758316446E-7) * r + 1.8463183175100548E-5) * r + 7.868691311456133E-4) * r + 0.014875361290850615) * r + 0.1369298809227358) * r + 0.599832206555888) * r + 1.0);
}
if (q < 0.0) {
quantile = -quantile;
}
}
}
return quantile;
} else {
System.err.println("NormQuantile(): probability outside (0, 1)");
return 0.0;
}
}
public static double peakToPeak(double[] data) {
return peakToPeak(data, data.length);
}
public static double peakToPeak(double[] data, int length) {
return abs(maximum(data, length) - minimum(data, length));
}
public static float peakToPeak(float[] data) {
return peakToPeak(data, data.length);
}
public static float peakToPeak(float[] data, int length) {
return abs(maximum(data, length) - minimum(data, length));
}
public static int peakToPeak(int[] data) {
return peakToPeak(data, data.length);
}
public static int peakToPeak(int[] data, int length) {
return abs(maximum(data, length) - minimum(data, length));
}
public static long peakToPeak(long[] data) {
return peakToPeak(data, data.length);
}
public static long peakToPeak(long[] data, int length) {
return abs(maximum(data, length) - minimum(data, length));
}
public static short peakToPeak(short[] data) {
return peakToPeak(data, data.length);
}
public static short peakToPeak(short[] data, int length) {
return (short)abs(maximum(data, length) - minimum(data, length));
}
public static boolean permute(int n, int[] a) {
int i1 = -1;
int i;
for(i = n - 2; i > -1; --i) {
if (a[i] < a[i + 1]) {
i1 = i;
break;
}
}
if (i1 == -1) {
return false;
} else {
int itmp;
for(i = n - 1; i > i1; --i) {
if (a[i] > a[i1]) {
itmp = a[i1];
a[i1] = a[i];
a[i] = itmp;
break;
}
}
for(i = 0; i < (n - i1 - 1) / 2; ++i) {
itmp = a[i1 + i + 1];
a[i1 + i + 1] = a[n - i - 1];
a[n - i - 1] = itmp;
}
return true;
}
}
public static double poisson(double x, double par) {
if (x < 0.0) {
return 0.0;
} else if (x == 0.0) {
return 1.0 / exp(par);
} else {
double lnpoisson = x * log(par) - par - lnGamma(x + 1.0);
return exp(lnpoisson);
}
}
public static double poissonI(double x, double par) {
if (x < 0.0) {
return 0.0;
} else if (x < 1.0) {
return exp(-par);
} else {
int ix = (int)x;
double kMaxInt = 2000000.0;
double gam;
if (x < 2000000.0) {
gam = pow(par, (double)ix) / gamma((double)(ix + 1));
} else {
gam = pow(par, x) / gamma(x + 1.0);
}
return gam / exp(par);
}
}
public static double prob(double chi2, int ndf) {
if (ndf <= 0) {
return 0.0;
} else if (chi2 <= 0.0) {
return chi2 < 0.0 ? 0.0 : 1.0;
} else if (ndf == 1) {
return 1.0 - erf(sqrt(chi2) / sqrt(2.0));
} else {
double q = sqrt(2.0 * chi2) - sqrt((double)(2 * ndf - 1));
return ndf > 30 && q > 5.0 ? 0.5 * (1.0 - erf(q / sqrt(2.0))) : 1.0 - gamma(0.5 * (double)ndf, 0.5 * chi2);
}
}
public static double rms(double[] data) {
return rms(data, data.length);
}
public static double rms(double[] data, int length) {
if (length <= 0) {
return -1.0;
} else {
double norm = 1.0 / (double)length;
double val1 = 0.0;
double val2 = 0.0;
for(int i = 0; i < length; ++i) {
val1 += data[i];
val2 += data[i] * data[i];
}
val1 *= norm;
val2 *= norm;
return sqrt(abs(val2 - val1 * val1));
}
}
public static float rms(float[] data) {
return rms(data, data.length);
}
public static float rms(float[] data, int length) {
if (length <= 0) {
return -1.0F;
} else {
double norm = 1.0 / (double)length;
double val1 = 0.0;
double val2 = 0.0;
for(int i = 0; i < length; ++i) {
val1 += (double)data[i];
val2 += (double)(data[i] * data[i]);
}
val1 *= norm;
val2 *= norm;
return (float)sqrt(abs(val2 - val1 * val1));
}
}
public static int rms(int[] data) {
return rms(data, data.length);
}
public static int rms(int[] data, int length) {
if (length <= 0) {
return -1;
} else {
double norm = 1.0 / (double)length;
double val1 = 0.0;
double val2 = 0.0;
for(int i = 0; i < length; ++i) {
val1 += (double)data[i];
val2 += (double)(data[i] * data[i]);
}
val1 *= norm;
val2 *= norm;
return (int)sqrt(abs(val2 - val1 * val1));
}
}
public static long rms(long[] data) {
return rms(data, data.length);
}
public static long rms(long[] data, int length) {
if (length <= 0) {
return -1L;
} else {
double norm = 1.0 / (double)length;
double val1 = 0.0;
double val2 = 0.0;
for(int i = 0; i < length; ++i) {
val1 += (double)data[i];
val2 += (double)(data[i] * data[i]);
}
val1 *= norm;
val2 *= norm;
return (long)sqrt(abs(val2 - val1 * val1));
}
}
public static short rms(short[] data) {
return rms(data, data.length);
}
public static short rms(short[] data, int length) {
if (length <= 0) {
return -1;
} else {
double norm = 1.0 / (double)length;
double val1 = 0.0;
double val2 = 0.0;
for(int i = 0; i < length; ++i) {
val1 += (double)data[i];
val2 += (double)(data[i] * data[i]);
}
val1 *= norm;
val2 *= norm;
return (short)((int)sqrt(abs(val2 - val1 * val1)));
}
}
public static double sinc(double x, boolean norm) {
if (x == 0.0) {
return 1.0;
} else {
double val = norm ? java.lang.Math.PI : 1.0;
return MathBase.sin(val * x) / (val * x);
}
}
public static double[] sort(double[] a, int length, boolean down) {
if (a != null && a.length > 0) {
double[] index = Arrays.copyOf(a, length);
Arrays.sort(index);
if (down) {
int nlast = length - 1;
for(int i = 0; i < length / 2; ++i) {
double temp = index[i];
index[i] = index[nlast - i];
index[nlast - i] = temp;
}
}
return index;
} else {
return new double[0];
}
}
public static float[] sort(float[] a, int length, boolean down) {
if (a != null && a.length > 0) {
float[] index = Arrays.copyOf(a, length);
Arrays.sort(index);
if (down) {
int nlast = length - 1;
for(int i = 0; i < length / 2; ++i) {
float temp = index[i];
index[i] = index[nlast - i];
index[nlast - i] = temp;
}
}
return index;
} else {
return new float[0];
}
}
public static int[] sort(int[] a, int length, boolean down) {
if (a != null && a.length > 0) {
int[] index = Arrays.copyOf(a, length);
Arrays.sort(index);
if (down) {
int nlast = length - 1;
for(int i = 0; i < length / 2; ++i) {
int temp = index[i];
index[i] = index[nlast - i];
index[nlast - i] = temp;
}
}
return index;
} else {
return new int[0];
}
}
public static long[] sort(long[] a, int length, boolean down) {
if (a != null && a.length > 0) {
long[] index = Arrays.copyOf(a, length);
Arrays.sort(index);
if (down) {
int nlast = length - 1;
for(int i = 0; i < length / 2; ++i) {
long temp = index[i];
index[i] = index[nlast - i];
index[nlast - i] = temp;
}
}
return index;
} else {
return new long[0];
}
}
public static short[] sort(short[] a, int length, boolean down) {
if (a != null && a.length > 0) {
short[] index = Arrays.copyOf(a, length);
Arrays.sort(index);
if (down) {
int nlast = length - 1;
for(int i = 0; i < length / 2; ++i) {
short temp = index[i];
index[i] = index[nlast - i];
index[nlast - i] = temp;
}
}
return index;
} else {
return new short[0];
}
}
public static double struveH0(double x) {
boolean n1 = true;
boolean n2 = true;
double[] c1 = new double[]{1.00215845609912, -1.6396929268130915, 1.502369396182928, -0.7248511530212187, 0.18955327371093136, -0.03067052022988, 0.00337561447375194, -2.6965014312602E-4, 1.637461692612E-5, -7.8244408508E-7, 3.021593188E-8, -9.6326645E-10, 2.579337E-11, -5.8854E-13, 1.158E-14, -2.0E-16};
double[] c2 = new double[]{0.9928372757642394, -0.00696891281138625, 1.8205103787037E-4, -1.063258252844E-5, 9.8198294287E-7, -1.2250645445E-7, 1.894083312E-8, -3.44358226E-9, 7.1119102E-10, -1.6288744E-10, 4.065681E-11, -1.091505E-11, 3.12005E-12, -9.4202E-13, 2.9848E-13, -9.872E-14, 3.394E-14, -1.208E-14, 4.44E-15, -1.68E-15, 6.5E-16, -2.6E-16, 1.1E-16, -4.0E-17, 2.0E-17, -1.0E-17};
double c0 = 0.6366197723675814;
double v = abs(x);
v = abs(x);
int i;
double alfa;
double h;
double b0;
double b1;
double b2;
if (v < 8.0) {
double y = v / 8.0;
h = 2.0 * y * y - 1.0;
alfa = h + h;
b0 = 0.0;
b1 = 0.0;
b2 = 0.0;
for(i = 15; i >= 0; --i) {
b0 = c1[i] + alfa * b1 - b2;
b2 = b1;
b1 = b0;
}
h = y * (b0 - h * b2);
} else {
double r = 1.0 / v;
h = 128.0 * r * r - 1.0;
alfa = h + h;
b0 = 0.0;
b1 = 0.0;
b2 = 0.0;
for(i = 25; i >= 0; --i) {
b0 = c2[i] + alfa * b1 - b2;
b2 = b1;
b1 = b0;
}
h = besselY0(v) + r * 0.6366197723675814 * (b0 - h * b2);
}
if (x < 0.0) {
h = -h;
}
return h;
}
public static double struveH1(double x) {
boolean n3 = true;
boolean n4 = true;
double[] c3 = new double[]{0.5578891446481605, -0.11188325726569816, -0.16337958125200938, 0.322569320724059, -0.14581632367244243, 0.03292677399374035, -0.00460372142093573, 4.434706163314E-4, -3.142099529341E-5, 1.7123719938E-6, -7.416987005E-8, 2.61837671E-9, -7.685839E-11, 1.9067E-12, -4.052E-14, 7.5E-16, -1.0E-17};
double[] c4 = new double[]{1.0075764729386565, 0.00750316051248257, -7.043933264519E-5, 2.66205393382E-6, -1.8841157753E-7, 1.949014958E-8, -2.6126199E-9, 4.236269E-10, -7.955156E-11, 1.679973E-11, -3.9072E-12, 9.8543E-13, -2.6636E-13, 7.645E-14, -2.313E-14, 7.33E-15, -2.42E-15, 8.3E-16, -3.0E-16, 1.1E-16, -4.0E-17, 2.0E-17, -1.0E-17};
double c0 = 0.6366197723675814;
double cc = 0.2122065907891938;
double v = abs(x);
double h;
if (v == 0.0) {
h = 0.0;
} else {
int i;
if (v <= 0.3) {
double y = v * v;
double r = 1.0;
h = 1.0;
int i1 = (int)(-8.0 / log10(v));
for(i = 1; i <= i1; ++i) {
h = -h * y / (double)((2 * i + 1) * (2 * i + 3));
r += h;
}
h = 0.2122065907891938 * y * r;
} else {
double alfa;
double b0;
double b1;
double b2;
if (v < 8.0) {
h = v * v / 32.0 - 1.0;
alfa = h + h;
b0 = 0.0;
b1 = 0.0;
b2 = 0.0;
for(i = 16; i >= 0; --i) {
b0 = c3[i] + alfa * b1 - b2;
b2 = b1;
b1 = b0;
}
h = b0 - h * b2;
} else {
h = 128.0 / (v * v) - 1.0;
alfa = h + h;
b0 = 0.0;
b1 = 0.0;
b2 = 0.0;
for(i = 22; i >= 0; --i) {
b0 = c4[i] + alfa * b1 - b2;
b2 = b1;
b1 = b0;
}
h = besselY1(v) + 0.6366197723675814 * (b0 - h * b2);
}
}
}
return h;
}
public static double struveL0(double x) {
double pi = java.lang.Math.PI;
double s = 1.0;
double r = 1.0;
int i;
double sl0;
if (x <= 20.0) {
double a0 = 2.0 * x / java.lang.Math.PI;
for(i = 1; i <= 60; ++i) {
r *= x / (double)(2 * i + 1) * (x / (double)(2 * i + 1));
s += r;
if (abs(r / s) < 1.0E-12) {
break;
}
}
sl0 = a0 * s;
} else {
int km = (int)(5.0 * (x + 1.0));
if (x >= 50.0) {
km = 25;
}
for(i = 1; i <= km; ++i) {
r *= (double)((2 * i - 1) * (2 * i - 1)) / x / x;
s += r;
if (abs(r / s) < 1.0E-12) {
break;
}
}
double a1 = exp(x) / sqrt(6.283185307179586 * x);
r = 1.0;
double bi0 = 1.0;
for(i = 1; i <= 16; ++i) {
r = 0.125 * r * (2.0 * (double)i - 1.0) * (2.0 * (double)i - 1.0) / ((double)i * x);
bi0 += r;
if (abs(r / bi0) < 1.0E-12) {
break;
}
}
bi0 *= a1;
sl0 = -2.0 / (java.lang.Math.PI * x) * s + bi0;
}
return sl0;
}
public static double struveL1(double x) {
double pi = java.lang.Math.PI;
double r = 1.0;
double sl1;
double s;
int i;
if (x <= 20.0) {
s = 0.0;
for(i = 1; i <= 60; ++i) {
r *= x * x / (4.0 * (double)i * (double)i - 1.0);
s += r;
if (abs(r) < abs(s) * 1.0E-12) {
break;
}
}
sl1 = 0.6366197723675814 * s;
} else {
s = 1.0;
int km = (int)(0.5 * x);
if (x > 50.0) {
km = 25;
}
for(i = 1; i <= km; ++i) {
r *= (double)((2 * i + 3) * (2 * i + 1)) / x / x;
s += r;
if (abs(r / s) < 1.0E-12) {
break;
}
}
sl1 = 0.6366197723675814 * (-1.0 + 1.0 / (x * x) + 3.0 * s / (x * x * x * x));
double a1 = exp(x) / sqrt(6.283185307179586 * x);
r = 1.0;
double bi1 = 1.0;
for(i = 1; i <= 16; ++i) {
r = -0.125 * r * (4.0 - (2.0 * (double)i - 1.0) * (2.0 * (double)i - 1.0)) / ((double)i * x);
bi1 += r;
if (abs(r / bi1) < 1.0E-12) {
break;
}
}
sl1 += a1 * bi1;
}
return sl1;
}
public static double student(double T, double ndf) {
if (ndf < 1.0) {
return 0.0;
} else {
double rh = 0.5 * ndf;
double rh1 = rh + 0.5;
double denom = sqrt(ndf * java.lang.Math.PI) * gamma(rh) * pow(1.0 + T * T / ndf, rh1);
return gamma(rh1) / denom;
}
}
public static double studentI(double T, double ndf) {
double si = T > 0.0 ? 1.0 - 0.5 * betaIncomplete(ndf / (ndf + T * T), ndf * 0.5, 0.5) : 0.5 * betaIncomplete(ndf / (ndf + T * T), ndf * 0.5, 0.5);
return si;
}
public static double studentQuantile(double p, double ndf, boolean lower_tail) {
if (!(ndf < 1.0) && !(p >= 1.0) && !(p <= 0.0)) {
boolean neg;
double q;
if ((!lower_tail || !(p > 0.5)) && (lower_tail || !(p < 0.5))) {
neg = true;
q = 2.0 * (lower_tail ? p : 1.0 - p);
} else {
neg = false;
q = 2.0 * (lower_tail ? 1.0 - p : p);
}
double quantile;
if (ndf - 1.0 < 1.0E-8) {
double temp = 1.5707963267948966 * q;
quantile = cos(temp) / sin(temp);
} else if (ndf - 2.0 < 1.0E-8) {
quantile = sqrt(2.0 / (q * (2.0 - q)) - 2.0);
} else {
double a = 1.0 / (ndf - 0.5);
double b = 48.0 / (a * a);
double c = ((20700.0 * a / b - 98.0) * a - 16.0) * a + 96.36;
double d = ((94.5 / (b + c) - 3.0) / b + 1.0) * sqrt(a * 1.5707963267948966) * ndf;
double x = q * d;
double y = pow(x, 2.0 / ndf);
if (y > 0.05 + a) {
x = normQuantile(q * 0.5);
y = x * x;
if (ndf < 5.0) {
c += 0.3 * (ndf - 4.5) * (x + 0.6);
}
c += (((0.05 * d * x - 5.0) * x - 7.0) * x - 2.0) * x + b;
y = (((((0.4 * y + 6.3) * y + 36.0) * y + 94.5) / c - y - 3.0) / b + 1.0) * x;
y = a * y * y;
if (y > 0.002) {
y = exp(y) - 1.0;
} else {
y += 0.5 * y * y;
}
} else {
y = ((1.0 / (((ndf + 6.0) / (ndf * y) - 0.089 * d - 0.822) * (ndf + 2.0) * 3.0) + 0.5 / (ndf + 4.0)) * y - 1.0) * (ndf + 1.0) / (ndf + 2.0) + 1.0 / y;
}
quantile = sqrt(ndf * y);
}
if (neg) {
quantile = -quantile;
}
return quantile;
} else {
System.err.println("StudentQuantile() - illegal parameter values");
return 0.0;
}
}
public static double variance(double[] aa, double[] ww) {
int n = aa.length;
if (n != ww.length) {
throw new IllegalArgumentException("length of variable array, " + n + " and length of weight array, " + ww.length + " are different");
} else {
double nn = effectiveSampleNumber(ww);
double nterm = nn / (nn - 1.0);
double sumx = 0.0;
double sumw = 0.0;
double mean = 0.0;
double[] weight = invertAndSquare(ww);
int i;
for(i = 0; i < n; ++i) {
sumx += aa[i] * weight[i];
sumw += weight[i];
}
mean = sumx / sumw;
sumx = 0.0;
for(i = 0; i < n; ++i) {
sumx += weight[i] * MathBase.sqr(aa[i] - mean);
}
return sumx * nterm / sumw;
}
}
public static double vavilov(double x, double kappa, double beta2) {
double[] ac = new double[14];
double[] hc = new double[9];
int[] itype = new int[1];
int[] npt = new int[1];
vavilovSet(kappa, beta2, false, (double[])null, ac, hc, itype, npt);
double v = vavilovDenEval(x, ac, hc, itype[0]);
return v;
}
protected static double vavilovDenEval(double rlam, double[] AC, double[] HC, int itype) {
if (!(rlam < AC[0]) && !(rlam > AC[8])) {
double v = 0.0;
double[] h = new double[10];
double x;
if (itype == 1) {
double fn = 1.0;
x = (rlam + HC[0]) * HC[1];
h[1] = x;
h[2] = x * x - 1.0;
int k;
for(k = 2; k <= 8; ++k) {
++fn;
h[k + 1] = x * h[k] - fn * h[k - 1];
}
double s = 1.0 + HC[7] * h[9];
for(k = 2; k <= 6; ++k) {
s += HC[k] * h[k + 1];
}
v = HC[8] * exp(-0.5 * x * x) * max(s, 0.0);
} else if (itype == 2) {
x = rlam * rlam;
v = AC[1] * exp(-AC[2] * (rlam + AC[5] * x) - AC[3] * exp(-AC[4] * (rlam + AC[6] * x)));
} else if (itype == 3) {
if (rlam < AC[7]) {
x = rlam * rlam;
v = AC[1] * exp(-AC[2] * (rlam + AC[5] * x) - AC[3] * exp(-AC[4] * (rlam + AC[6] * x)));
} else {
x = 1.0 / rlam;
v = (AC[11] * x + AC[12]) * x;
}
} else if (itype == 4) {
}
return v;
} else {
return 0.0;
}
}
public static double vavilovI(double x, double kappa, double beta2) {
double[] ac = new double[14];
double[] hc = new double[9];
double[] wcm = new double[200];
int[] itype = new int[0];
int[] npt = new int[0];
vavilovSet(kappa, beta2, true, wcm, ac, hc, itype, npt);
double v;
if (x < ac[0]) {
v = 0.0;
} else if (x >= ac[8]) {
v = 1.0;
} else {
double xx = x - ac[0];
int k = (int)(xx * ac[10]);
v = min(wcm[k] + (xx - (double)k * ac[9]) * (wcm[k + 1] - wcm[k]) * ac[10], 1.0);
}
return v;
}
public static void vavilovSet(double rkappa, double beta2, boolean mode, double[] WCM, double[] AC, double[] HC, int[] itype, int[] npt) {
if (!(rkappa < 0.01) && !(rkappa > 12.0)) {
double BKMNX1 = 0.02;
double BKMNY1 = 0.05;
double BKMNX2 = 0.12;
double BKMNY2 = 0.05;
double BKMNX3 = 0.22;
double BKMNY3 = 0.05;
double BKMXX1 = 0.1;
double BKMXY1 = 1.0;
double BKMXX2 = 0.2;
double BKMXY2 = 1.0;
double BKMXX3 = 0.3;
double BKMXY3 = 1.0;
double FBKX1 = 25.0;
double FBKX2 = 24.999999999999996;
double FBKX3 = 25.000000000000004;
double FBKY1 = 2.1052631578947367;
double FBKY2 = 2.1052631578947367;
double FBKY3 = 2.1052631578947367;
double[] FNINV = new double[]{0.0, 1.0, 0.5, 0.33333333, 0.25, 0.2};
double[] EDGEC = new double[]{0.0, 0.0, 0.16666667, 0.041666667, 0.0083333333, 0.013888889, 0.0069444444, 7.7160493E-4};
double[] U1 = new double[]{0.0, 0.25850868, 0.032477982, -0.0059020496, 0.0, 0.024880692, 0.0047404356, -7.444513E-4, 0.0073225731, 0.0, 0.0011668284, 0.0, -0.0015727318, -0.0011210142};
double[] U2 = new double[]{0.0, 0.43142611, 0.040797543, -0.0091490215, 0.0, 0.042127077, 0.0073167928, -0.0014026047, 0.016195241, 0.0024714789, 0.0020751278, 0.0, -0.0025141668, -0.0014064022};
double[] U3 = new double[]{0.0, 0.25225955, 0.064820468, -0.023615759, 0.0, 0.023834176, 0.0021624675, -0.0026865597, -0.0054891384, 0.0039800522, 0.0048447456, -0.0089439554, -0.0062756944, -0.0024655436};
double[] U4 = new double[]{0.0, 1.2593231, -0.20374501, 0.095055662, -0.020771531, -0.04686518, -0.0077222986, 0.0032241039, 0.008988292, -0.0067167236, -0.013049241, 0.018786468, 0.014484097};
double[] U5 = new double[]{0.0, -0.024864376, -0.0010368495, 0.0014330117, 2.005273E-4, 0.0018751903, 0.0012668869, 4.8736023E-4, 0.0034850854, 0.0, -3.6597173E-4, 0.0019372124, 7.0761825E-4, 4.6898375E-4};
double[] U6 = new double[]{0.0, 0.035855696, -0.027542114, 0.012631023, -0.0030188807, -8.4479939E-4, 0.0, 4.5675843E-4, -0.0069836141, 0.0039876546, -0.0036055679, 0.0, 0.0015298434, 0.0019247256};
double[] U7 = new double[]{0.0, 10.234691, -3.5619655, 0.69387764, -0.14047599, -1.995239, -0.45679694, 0.0, 0.50505298};
double[] U8 = new double[]{0.0, 21.487518, -11.825253, 4.3133087, -1.4500543, -3.4343169, -1.1063164, -0.21000819, 1.7891643, -0.89601916, 0.39120793, 0.73410606, 0.0, -0.32454506};
double[] V1 = new double[]{0.0, 0.27827257, -0.0014227603, 0.0024848327, 0.0, 0.045091424, 0.0080559636, -0.0038974523, 0.0, -0.0030634124, 7.5633702E-4, 0.0054730726, 0.0019792507};
double[] V2 = new double[]{0.0, 0.41421789, -0.030061649, 0.0052249697, 0.0, 0.12693873, 0.022999801, -0.0086792801, 0.031875584, -0.0061757928, 0.0, 0.019716857, 0.0032596742};
double[] V3 = new double[]{0.0, 0.20191056, -0.046831422, 0.0096777473, -0.0017995317, 0.053921588, 0.003506874, -0.012621494, -0.0054996531, -0.0090029985, 0.0034958743, 0.018513506, 0.0068332334, -0.0012940502};
double[] V4 = new double[]{0.0, 1.3206081, 0.10036618, -0.022015201, 0.0061667091, -0.14986093, -0.012720568, 0.024972042, -0.0097751962, 0.026087455, -0.011399062, -0.048282515, -0.0098552378};
double[] V5 = new double[]{0.0, 0.016435243, 0.0360514, 0.002303652, -6.1666343E-4, -0.010775802, 0.0051476061, 0.0056856517, -0.013438433, 0.0, 0.0, -0.0025421507, 0.0020169108, -0.0015144931};
double[] V6 = new double[]{0.0, 0.033432405, 0.0060583916, -0.0023381379, 8.3846081E-4, -0.013346861, -0.0017402116, 0.0021052496, 0.0015528195, 0.002190067, -0.0013202847, -0.0045124157, -0.0015629454, 2.2499176E-4};
double[] V7 = new double[]{0.0, 5.4529572, -0.90906096, 0.086122438, 0.0, -1.2218009, -0.3232412, -0.027373591, 0.12173464, 0.0, 0.0, 0.040917471};
double[] V8 = new double[]{0.0, 9.3841352, -1.6276904, 0.16571423, 0.0, -1.8160479, -0.50919193, -0.051384654, 0.21413992, 0.0, 0.0, 0.066596366};
double[] W1 = new double[]{0.0, 0.29712951, 0.0097572934, 0.0, -0.0015291686, 0.035707399, 0.0096221631, -0.0018402821, -0.0049821585, 0.0018831112, 0.0043541673, 0.0020301312, -0.0018723311, -7.3403108E-4};
double[] W2 = new double[]{0.0, 0.40882635, 0.014474912, 0.0025023704, -0.0037707379, 0.18719727, 0.056954987, 0.0, 0.023020158, 0.0050574313, 0.009455014, 0.019300232};
double[] W3 = new double[]{0.0, 0.16861629, 0.0, 0.0036317285, -0.0043657818, 0.030144338, 0.013891826, -0.0058030495, -0.0038717547, 0.0085359607, 0.014507659, 0.0082387775, -0.010116105, -0.005513567};
double[] W4 = new double[]{0.0, 1.3493891, -0.0026863185, -0.003521604, 0.024434909, -0.083447911, -0.04806136, 0.0076473951, 0.02449443, -0.0162092, -0.037768479, -0.047890063, 0.017778596, 0.013179324};
double[] W5 = new double[]{0.0, 0.10264945, 0.032738857, 0.0, 0.0043608779, -0.043097757, -0.0022647176, 0.009453129, -0.012442571, -0.0032283517, -0.0075640352, -0.0088293329, 0.0052537299, 0.0013340546};
double[] W6 = new double[]{0.0, 0.029568177, -0.001630006, -2.1119745E-4, 0.0023599053, -0.0048515387, -0.0040797531, 4.0403265E-4, 0.0018200105, -0.0014346306, -0.0039165276, -0.0037432073, 0.001995038, 0.0012222675};
double[] W8 = new double[]{0.0, 6.6184645, -0.73866379, 0.044693973, 0.0, -1.4540925, -0.39529833, -0.044293243, 0.088741049};
itype[0] = 0;
double[] DRK = new double[6];
double[] DSIGM = new double[6];
double[] ALFA = new double[8];
double x;
double y;
double p2;
double fl;
if (rkappa >= 0.29) {
itype[0] = 1;
npt[0] = 100;
fl = 1.0 / sqrt(rkappa);
AC[0] = (-0.032227 * beta2 - 0.074275) * rkappa + (0.24533 * beta2 + 0.070152) * fl + (-0.5561 * beta2 - 3.1579);
AC[8] = (-0.013483 * beta2 - 0.048801) * rkappa + (-1.6921 * beta2 + 8.3656) * fl + (-0.73275 * beta2 - 3.5226);
DRK[1] = fl * fl;
DSIGM[1] = sqrt(rkappa / (1.0 - 0.5 * beta2));
int j;
for(j = 1; j <= 4; ++j) {
DRK[j + 1] = DRK[1] * DRK[j];
DSIGM[j + 1] = DSIGM[1] * DSIGM[j];
ALFA[j + 1] = (FNINV[j] - beta2 * FNINV[j + 1]) * DRK[j];
}
HC[0] = log(rkappa) + beta2 + 0.42278434;
HC[1] = DSIGM[1];
HC[2] = ALFA[3] * DSIGM[3];
HC[3] = (3.0 * ALFA[2] * ALFA[2] + ALFA[4]) * DSIGM[4] - 3.0;
HC[4] = (10.0 * ALFA[2] * ALFA[3] + ALFA[5]) * DSIGM[5] - 10.0 * HC[2];
HC[5] = HC[2] * HC[2];
HC[6] = HC[2] * HC[3];
HC[7] = HC[2] * HC[5];
for(j = 2; j <= 7; ++j) {
HC[j] *= EDGEC[j];
}
HC[8] = 0.39894228 * HC[1];
} else {
double xx;
double yy;
double x2;
double x3;
double y2;
double y3;
double xy;
double p3;
double q2;
double q3;
double pq;
if (rkappa >= 0.22) {
itype[0] = 2;
npt[0] = 150;
x = 1.0 + (rkappa - 0.3) * 25.000000000000004;
y = 1.0 + (sqrt(beta2) - 1.0) * 2.1052631578947367;
xx = 2.0 * x;
yy = 2.0 * y;
x2 = xx * x - 1.0;
x3 = xx * x2 - x;
y2 = yy * y - 1.0;
y3 = yy * y2 - y;
xy = x * y;
p2 = x2 * y;
p3 = x3 * y;
q2 = y2 * x;
q3 = y3 * x;
pq = x2 * y2;
AC[1] = W1[1] + W1[2] * x + W1[4] * x3 + W1[5] * y + W1[6] * y2 + W1[7] * y3 + W1[8] * xy + W1[9] * p2 + W1[10] * p3 + W1[11] * q2 + W1[12] * q3 + W1[13] * pq;
AC[2] = W2[1] + W2[2] * x + W2[3] * x2 + W2[4] * x3 + W2[5] * y + W2[6] * y2 + W2[8] * xy + W2[9] * p2 + W2[10] * p3 + W2[11] * q2;
AC[3] = W3[1] + W3[3] * x2 + W3[4] * x3 + W3[5] * y + W3[6] * y2 + W3[7] * y3 + W3[8] * xy + W3[9] * p2 + W3[10] * p3 + W3[11] * q2 + W3[12] * q3 + W3[13] * pq;
AC[4] = W4[1] + W4[2] * x + W4[3] * x2 + W4[4] * x3 + W4[5] * y + W4[6] * y2 + W4[7] * y3 + W4[8] * xy + W4[9] * p2 + W4[10] * p3 + W4[11] * q2 + W4[12] * q3 + W4[13] * pq;
AC[5] = W5[1] + W5[2] * x + W5[4] * x3 + W5[5] * y + W5[6] * y2 + W5[7] * y3 + W5[8] * xy + W5[9] * p2 + W5[10] * p3 + W5[11] * q2 + W5[12] * q3 + W5[13] * pq;
AC[6] = W6[1] + W6[2] * x + W6[3] * x2 + W6[4] * x3 + W6[5] * y + W6[6] * y2 + W6[7] * y3 + W6[8] * xy + W6[9] * p2 + W6[10] * p3 + W6[11] * q2 + W6[12] * q3 + W6[13] * pq;
AC[8] = W8[1] + W8[2] * x + W8[3] * x2 + W8[5] * y + W8[6] * y2 + W8[7] * y3 + W8[8] * xy;
AC[0] = -3.05;
} else if (rkappa >= 0.12) {
itype[0] = 3;
npt[0] = 200;
x = 1.0 + (rkappa - 0.2) * 24.999999999999996;
y = 1.0 + (sqrt(beta2) - 1.0) * 2.1052631578947367;
xx = 2.0 * x;
yy = 2.0 * y;
x2 = xx * x - 1.0;
x3 = xx * x2 - x;
y2 = yy * y - 1.0;
y3 = yy * y2 - y;
xy = x * y;
p2 = x2 * y;
p3 = x3 * y;
q2 = y2 * x;
q3 = y3 * x;
pq = x2 * y2;
AC[1] = V1[1] + V1[2] * x + V1[3] * x2 + V1[5] * y + V1[6] * y2 + V1[7] * y3 + V1[9] * p2 + V1[10] * p3 + V1[11] * q2 + V1[12] * q3;
AC[2] = V2[1] + V2[2] * x + V2[3] * x2 + V2[5] * y + V2[6] * y2 + V2[7] * y3 + V2[8] * xy + V2[9] * p2 + V2[11] * q2 + V2[12] * q3;
AC[3] = V3[1] + V3[2] * x + V3[3] * x2 + V3[4] * x3 + V3[5] * y + V3[6] * y2 + V3[7] * y3 + V3[8] * xy + V3[9] * p2 + V3[10] * p3 + V3[11] * q2 + V3[12] * q3 + V3[13] * pq;
AC[4] = V4[1] + V4[2] * x + V4[3] * x2 + V4[4] * x3 + V4[5] * y + V4[6] * y2 + V4[7] * y3 + V4[8] * xy + V4[9] * p2 + V4[10] * p3 + V4[11] * q2 + V4[12] * q3;
AC[5] = V5[1] + V5[2] * x + V5[3] * x2 + V5[4] * x3 + V5[5] * y + V5[6] * y2 + V5[7] * y3 + V5[8] * xy + V5[11] * q2 + V5[12] * q3 + V5[13] * pq;
AC[6] = V6[1] + V6[2] * x + V6[3] * x2 + V6[4] * x3 + V6[5] * y + V6[6] * y2 + V6[7] * y3 + V6[8] * xy + V6[9] * p2 + V6[10] * p3 + V6[11] * q2 + V6[12] * q3 + V6[13] * pq;
AC[7] = V7[1] + V7[2] * x + V7[3] * x2 + V7[5] * y + V7[6] * y2 + V7[7] * y3 + V7[8] * xy + V7[11] * q2;
AC[8] = V8[1] + V8[2] * x + V8[3] * x2 + V8[5] * y + V8[6] * y2 + V8[7] * y3 + V8[8] * xy + V8[11] * q2;
AC[0] = -3.04;
} else {
itype[0] = 4;
if (rkappa >= 0.02) {
itype[0] = 3;
}
npt[0] = 200;
x = 1.0 + (rkappa - 0.1) * 25.0;
y = 1.0 + (sqrt(beta2) - 1.0) * 2.1052631578947367;
xx = 2.0 * x;
yy = 2.0 * y;
x2 = xx * x - 1.0;
x3 = xx * x2 - x;
y2 = yy * y - 1.0;
y3 = yy * y2 - y;
xy = x * y;
p2 = x2 * y;
p3 = x3 * y;
q2 = y2 * x;
q3 = y3 * x;
pq = x2 * y2;
if (itype[0] == 3) {
AC[1] = U1[1] + U1[2] * x + U1[3] * x2 + U1[5] * y + U1[6] * y2 + U1[7] * y3 + U1[8] * xy + U1[10] * p3 + U1[12] * q3 + U1[13] * pq;
AC[2] = U2[1] + U2[2] * x + U2[3] * x2 + U2[5] * y + U2[6] * y2 + U2[7] * y3 + U2[8] * xy + U2[9] * p2 + U2[10] * p3 + U2[12] * q3 + U2[13] * pq;
AC[3] = U3[1] + U3[2] * x + U3[3] * x2 + U3[5] * y + U3[6] * y2 + U3[7] * y3 + U3[8] * xy + U3[9] * p2 + U3[10] * p3 + U3[11] * q2 + U3[12] * q3 + U3[13] * pq;
AC[4] = U4[1] + U4[2] * x + U4[3] * x2 + U4[4] * x3 + U4[5] * y + U4[6] * y2 + U4[7] * y3 + U4[8] * xy + U4[9] * p2 + U4[10] * p3 + U4[11] * q2 + U4[12] * q3;
AC[5] = U5[1] + U5[2] * x + U5[3] * x2 + U5[4] * x3 + U5[5] * y + U5[6] * y2 + U5[7] * y3 + U5[8] * xy + U5[10] * p3 + U5[11] * q2 + U5[12] * q3 + U5[13] * pq;
AC[6] = U6[1] + U6[2] * x + U6[3] * x2 + U6[4] * x3 + U6[5] * y + U6[7] * y3 + U6[8] * xy + U6[9] * p2 + U6[10] * p3 + U6[12] * q3 + U6[13] * pq;
AC[7] = U7[1] + U7[2] * x + U7[3] * x2 + U7[4] * x3 + U7[5] * y + U7[6] * y2 + U7[8] * xy;
}
AC[8] = U8[1] + U8[2] * x + U8[3] * x2 + U8[4] * x3 + U8[5] * y + U8[6] * y2 + U8[7] * y3 + U8[8] * xy + U8[9] * p2 + U8[10] * p3 + U8[11] * q2 + U8[13] * pq;
AC[0] = -3.03;
}
}
AC[9] = (AC[8] - AC[0]) / (double)npt[0];
AC[10] = 1.0 / AC[9];
if (itype[0] == 3) {
x = (AC[7] - AC[8]) / (AC[7] * AC[8]);
y = 1.0 / log(AC[8] / AC[7]);
p2 = AC[7] * AC[7];
AC[11] = p2 * (AC[1] * exp(-AC[2] * (AC[7] + AC[5] * p2) - AC[3] * exp(-AC[4] * (AC[7] + AC[6] * p2))) - 0.045 * y / AC[7]) / (1.0 + x * y * AC[7]);
AC[12] = (0.045 + x * AC[11]) * y;
}
if (itype[0] == 4) {
AC[13] = 0.995 / landauI(AC[8]);
}
if (mode) {
x = AC[0];
if (WCM != null && WCM.length != 0) {
WCM[0] = 0.0;
fl = vavilovDenEval(x, AC, HC, itype[0]);
int k;
for(k = 1; k <= npt[0]; ++k) {
x += AC[9];
double fu = vavilovDenEval(x, AC, HC, itype[0]);
WCM[k] = WCM[k - 1] + fl + fu;
fl = fu;
}
x = 0.5 * AC[9];
for(k = 1; k <= npt[0]; ++k) {
WCM[k] *= x;
}
} else {
throw new IllegalArgumentException("WCM must not be null or zero-length");
}
}
} else {
System.err.println("Vavilov distribution - illegal value of kappa");
}
}
}
|
Beta Was this translation helpful? Give feedback.
0 replies
-
package io.fair_acc.math;
import static io.fair_acc.dataset.DataSet.DIM_X;
import static io.fair_acc.dataset.DataSet.DIM_Y;
import static io.fair_acc.dataset.DataSet.DIM_Z;
import static io.fair_acc.dataset.Histogram.Boundary.LOWER;
import static io.fair_acc.dataset.Histogram.Boundary.UPPER;
import static io.fair_acc.math.DataSetMath.ErrType.EXP;
import static io.fair_acc.math.DataSetMath.ErrType.EYN;
import static io.fair_acc.math.DataSetMath.ErrType.EYP;
import static io.fair_acc.math.MathBase.max;
import static io.fair_acc.math.MathBase.min;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.Objects;
import org.jetbrains.annotations.NotNull;
import org.jtransforms.fft.DoubleFFT_1D;
import io.fair_acc.dataset.*;
import io.fair_acc.dataset.spi.DoubleDataSet;
import io.fair_acc.dataset.spi.DoubleErrorDataSet;
import io.fair_acc.dataset.spi.Histogram;
import io.fair_acc.dataset.spi.utils.DoublePointError;
import io.fair_acc.dataset.utils.NoDuplicatesList;
import io.fair_acc.math.spectra.Apodization;
import io.fair_acc.math.spectra.SpectrumTools;
/**
* Some math operation on DataSet, DataSetError and Histogram
*
* @author rstein
*/
public final class DataSetMath { // NOPMD - nomen est omen
private static final char INTEGRAL_SYMBOL = 0x222B;
private static final char DIFFERENTIAL_SYMBOL = 0x2202;
private static final char MULTIPLICATION_SYMBOL = 0x00B7;
private static final String DIFFERENTIAL = DIFFERENTIAL_SYMBOL + "/" + DIFFERENTIAL_SYMBOL + "x";
public static Formatter<Number> DEFAULT_FORMATTER = new DefaultNumberFormatter(); // NOSONAR NOPMD -- explicitly not getter/setter
/**
*
*/
private DataSetMath() {
// private function, never called
}
public static DoubleErrorDataSet applyMathOperation(final DoubleErrorDataSet ret, final MathOp op, final double x1, final double y1, final double y2, final double eyn1, final double eyp1, final double eyn2, final double eyp2) { // NOPMD NOSONAR
// switch through math operations
switch (op) {
case ADD:
return ret.add(x1, y1 + y2, MathBase.hypot(eyn1, eyn2), MathBase.hypot(eyp1, eyp2));
case SUBTRACT:
return ret.add(x1, y1 - y2, MathBase.hypot(eyn1, eyn2), MathBase.hypot(eyp1, eyp2));
case MULTIPLY:
return ret.add(x1, y1 * y2, MathBase.hypot(y2 * eyn1, y1 * eyn2), MathBase.hypot(y2 * eyp1, y1 * eyp2));
case DIVIDE:
final double newY = y1 / y2;
final double nEYN = MathBase.hypot(eyn1 / y2, newY * eyn2 / y2);
final double nEYP = MathBase.hypot(eyp1 / y2, newY * eyp2 / y2);
return ret.add(x1, newY, nEYN, nEYP);
case SQR:
return ret.add(x1, MathBase.sqr(y1 + y2), 2 * MathBase.abs(y1 + y2) * MathBase.hypot(eyn1, eyn2), 2 * MathBase.abs(y1 + y2) * MathBase.hypot(eyp1, eyp2));
case SQRT:
return ret.add(x1, MathBase.sqrt(y1 + y2), MathBase.sqrt(MathBase.abs(y1 + y2)) * MathBase.hypot(eyn1, eyn2), MathBase.sqrt(MathBase.abs(y1 + y2)) * MathBase.hypot(eyp1, eyp2));
case LOG10:
double norm = 1.0 / MathBase.log(10);
final double nEYNLog = y1 + y2 > 0 ? norm / MathBase.abs(y1 + y2) * MathBase.hypot(eyn1, eyn2) : Double.NaN;
final double nEYPLog = y1 + y2 > 0 ? norm / MathBase.abs(y1 + y2) * MathBase.hypot(eyp1, eyp2) : Double.NaN;
return ret.add(x1, 10 * MathBase.log10(y1 + y2), nEYNLog, nEYPLog);
case DB:
final double normDb = 20.0 / MathBase.log(10);
final double nEYNDb = y1 + y2 > 0 ? normDb / MathBase.abs(y1 + y2) * MathBase.hypot(eyn1, eyn2) : Double.NaN;
final double nEYPDb = y1 + y2 > 0 ? normDb / MathBase.abs(y1 + y2) * MathBase.hypot(eyp1, eyp2) : Double.NaN;
return ret.add(x1, 20 * MathBase.log10(y1 + y2), nEYNDb, nEYPDb);
case IDENTITY:
default:
return ret.add(x1, y1 + y2, eyn1, eyp1);
}
}
// convenience short-hand notation for getting error variables (if defined for dataset)
private static double error(final DataSet dataSet, final ErrType eType, final int index, final double x,
final boolean interpolate) {
if (!(dataSet instanceof DataSetError)) {
// data set does not have any error definition
return 0.0;
}
final DataSetError ds = (DataSetError) dataSet;
if (interpolate) {
switch (eType) {
case EXN:
return ds.getErrorNegative(DIM_X, x);
case EXP:
return ds.getErrorPositive(DIM_X, x);
case EYN:
return ds.getErrorNegative(DIM_Y, x);
case EYP:
return ds.getErrorPositive(DIM_Y, x);
}
} else {
switch (eType) {
case EXN:
return ds.getErrorNegative(DIM_X, index);
case EXP:
return ds.getErrorPositive(DIM_X, index);
case EYN:
return ds.getErrorNegative(DIM_Y, index);
case EYP:
return ds.getErrorPositive(DIM_Y, index);
}
}
return 0;
}
@SafeVarargs
public static DataSet addFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.ADD, format);
}
@SafeVarargs
public static DataSet addFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.ADD, format);
}
@SafeVarargs
public static DataSet addGaussianNoise(final DataSet function, final double sigma, @NotNull final Formatter<Number>... format) {
final int nLength = function.getDataCount();
final var dataSetName = getFormatter(format).format("{0}+noise({1})", function.getName(), sigma);
final var ret = new DoubleErrorDataSet(dataSetName, nLength);
for (var i = 0; i < nLength; i++) {
final double x = function.get(DIM_X, i);
final double y = function.get(DIM_Y, i) + TRandom.Gaus(0, sigma);
ret.add(x, y, sigma, sigma);
}
return ret;
}
@SafeVarargs
public static DataSet averageDataSetsFIR(@NotNull final List<DataSet> dataSets, final int nUpdates, @NotNull final Formatter<Number>... format) {
if (dataSets.isEmpty()) {
return new DoubleErrorDataSet(getFormatter(format).format("LP({0}, FIR)", "<empty>"));
}
final String functionName = getFormatter(format).format("LP({0}, FIR)", dataSets.get(0).getName());
if (dataSets.size() <= 1) {
final var newFunction = dataSets.get(0);
if (newFunction instanceof DataSetError) {
return new DoubleErrorDataSet(functionName, newFunction.getValues(DIM_X), newFunction.getValues(DIM_Y),
errors(newFunction, EYN), errors(newFunction, EYP), newFunction.getDataCount(), true);
}
final int ncount = newFunction.getDataCount();
return new DoubleErrorDataSet(functionName, newFunction.getValues(DIM_X), newFunction.getValues(DIM_Y),
new double[ncount], new double[ncount], ncount, true);
}
final int nAvg = min(nUpdates, dataSets.size());
final var newFunction = dataSets.get(dataSets.size() - 1);
final var retFunction = new DoubleErrorDataSet(functionName, newFunction.getDataCount() + 2);
for (var i = 0; i < newFunction.getDataCount(); i++) {
final double newX = newFunction.get(DIM_X, i);
var mean = 0.0;
var variance = 0.0;
var eyn = 0.0;
var eyp = 0.0;
var count = 0;
for (int j = max(0, dataSets.size() - nAvg); j < dataSets.size(); j++) {
final var oldFunction = dataSets.get(j);
final double oldX = oldFunction.get(DIM_X, i);
final double oldY = oldX == newX ? oldFunction.get(DIM_Y, i) : oldFunction.getValue(DIM_X, newX);
mean += oldY;
variance += oldY * oldY;
// whether we need to interpolate
final boolean inter = oldX != newX;
eyn += error(oldFunction, EYN, i, newX, inter);
eyp += error(oldFunction, EYP, i, newX, inter);
count++;
}
if (count == 0) {
// cannot compute average
retFunction.add(newX, Double.NaN, Double.NaN, Double.NaN);
} else {
mean /= count;
eyn /= count;
eyp /= count;
variance /= count;
final double mean2 = mean * mean;
final double diff = MathBase.abs(variance - mean2);
eyn = MathBase.sqrt(eyn * eyn + diff);
eyp = MathBase.sqrt(eyp * eyp + diff);
retFunction.add(newX, mean, eyn, eyp);
}
}
return retFunction;
}
@SafeVarargs
public static DataSet averageDataSetsIIR(final DataSet prevAverage, final DataSet prevAverage2, final DataSet newDataSet, final int nUpdates, @NotNull final Formatter<Number>... format) {
final String functionName = getFormatter(format).format("LP({0}, IIR)", newDataSet.getName());
if (prevAverage == null || prevAverage2 == null || prevAverage.getDataCount() == 0 || prevAverage2.getDataCount() == 0) {
final double[] yValues = newDataSet.getValues(DIM_Y);
final double[] eyn = errors(newDataSet, EYN);
final double[] eyp = errors(newDataSet, EYP);
if (prevAverage2 instanceof DoubleErrorDataSet) {
((DoubleErrorDataSet) prevAverage2).set(newDataSet.getValues(DIM_X), ArrayMath.sqr(yValues), ArrayMath.sqr(eyn), ArrayMath.sqr(eyp));
} else if (prevAverage2 instanceof DoubleDataSet) {
((DoubleDataSet) prevAverage2).set(newDataSet.getValues(DIM_X), ArrayMath.sqr(yValues));
}
return new DoubleErrorDataSet(functionName, newDataSet.getValues(DIM_X), yValues, eyn, eyp,
newDataSet.getDataCount(), true);
}
final int dataCount1 = prevAverage.getDataCount();
final int dataCount2 = prevAverage2.getDataCount();
final DoubleErrorDataSet retFunction = dataCount1 == 0
? new DoubleErrorDataSet(functionName, newDataSet.getValues(DIM_X), newDataSet.getValues(DIM_Y),
errors(newDataSet, EYN), errors(newDataSet, EYP), newDataSet.getDataCount(), true)
: new DoubleErrorDataSet(prevAverage.getName(), prevAverage.getValues(DIM_X), prevAverage.getValues(DIM_Y),
errors(prevAverage, EYN), errors(prevAverage, EYP), newDataSet.getDataCount(), true);
final double alpha = 1.0 / (1.0 + nUpdates);
final boolean avg2Empty = dataCount2 == 0;
for (var i = 0; i < dataCount1; i++) {
final double oldX = prevAverage.get(DIM_X, i);
final double oldY = prevAverage.get(DIM_Y, i);
final double oldY2 = avg2Empty ? oldY * oldY : prevAverage2.get(DIM_Y, i);
final double newX = newDataSet.get(DIM_X, i);
// whether we need to interpolate
final boolean inter = oldX != newX;
final double y = inter ? newDataSet.getValue(DIM_Y, oldX) : newDataSet.get(DIM_Y, i);
final double newVal = (1 - alpha) * oldY + alpha * y;
final double newVal2 = (1 - alpha) * oldY2 + alpha * (y * y);
final double eyn = error(newDataSet, EYN, i, newX, inter);
final double eyp = error(newDataSet, EYP, i, newX, inter);
if (prevAverage2 instanceof DoubleErrorDataSet) {
if (avg2Empty) {
((DoubleErrorDataSet) prevAverage2).add(newX, newVal2, eyn, eyp);
} else {
((DoubleErrorDataSet) prevAverage2).set(i, newX, newVal2, eyn, eyp);
}
}
final double newEYN = MathBase.sqrt(MathBase.abs(newVal2 - MathBase.pow(newVal, 2)) + eyn * eyn);
final double newEYP = MathBase.sqrt(MathBase.abs(newVal2 - MathBase.pow(newVal, 2)) + eyp * eyp);
retFunction.set(i, oldX, newVal, newEYN, newEYP);
}
return retFunction;
}
@SafeVarargs
public static DataSet dbFunction(final DataSet function, @NotNull final Formatter<Number>... format) {
return mathFunction(function, 0.0, MathOp.DB, format);
}
@SafeVarargs
public static DataSet dbFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.DB, format);
}
@SafeVarargs
public static DataSet derivativeFunction(final DataSet function, @NotNull final Formatter<Number>... format) {
return derivativeFunction(function, +1.0, format);
}
@SafeVarargs
public static DataSet derivativeFunction(final DataSet function, final double sign, @NotNull final Formatter<Number>... format) {
final String signAdd = sign == 1.0 ? "" : Double.toString(sign) + MULTIPLICATION_SYMBOL;
final int ncount = function.getDataCount();
final String functionName = getFormatter(format).format("{0}{1}({2})", signAdd, DIFFERENTIAL, function.getName());
final var retFunction = new DoubleErrorDataSet(functionName, ncount);
if (ncount <= 3) {
return retFunction;
}
// TODO: check error estimate for derivative ...
// // derivative for first point
// final double y0 = function.get(DIM_Y, 0);
// final double y1 = function.get(DIM_Y, 0);
// final double yen0 = error(function, EYN, 0);
// final double yep0 = error(function, EYP, 0);
// final double yen1 = error(function, EYN, 1);
// final double yep1 = error(function, EYP, 1);
// final double deltaY0 = y1 - y0;
// final double deltaYEN = MathBase.sqrt(MathBase.pow(yen0, 2) + MathBase.pow(yen1,
// 2));
// final double deltaYEP = MathBase.sqrt(MathBase.pow(yep0, 2) + MathBase.pow(yep1,
// 2));
// final double deltaX0 = function.get(DIM_X, 1) - function.get(DIM_X, 0);
// retFunction.add(function.get(DIM_X, 0), sign * (deltaY0 / deltaX0),
// deltaYEN, deltaYEP);
for (var i = 0; i < 2; i++) {
final double x0 = function.get(DIM_X, i);
retFunction.add(x0, 0, 0, 0);
}
for (var i = 2; i < ncount - 2; i++) {
final double x0 = function.get(DIM_X, i);
final double stepL = x0 - function.get(DIM_X, i - 1);
final double stepR = function.get(DIM_X, i + 1) - x0;
final double valL = function.get(DIM_Y, i - 1);
final double valC = function.get(DIM_Y, i);
final double valR = function.get(DIM_Y, i + 1);
final double yenL = error(function, EYN, i - 1);
final double yenC = error(function, EYN, i);
final double yenR = error(function, EYN, i + 1);
final double yen = MathBase.sqrt(MathBase.sqr(yenL) + MathBase.sqr(yenC) + MathBase.sqr(yenR))
/ 4;
final double yepL = error(function, EYP, i - 1);
final double yepC = error(function, EYP, i);
final double yepR = error(function, EYP, i + 1);
final double yep = MathBase.sqrt(MathBase.sqr(yepL) + MathBase.sqr(yepC) + MathBase.sqr(yepR))
/ 4;
// simple derivative computation
final double derivative = 0.5 * ((valC - valL) / stepL + (valR - valC) / stepR);
retFunction.add(x0, sign * derivative, yen, yep);
}
for (int i = ncount - 2; i < ncount; i++) {
final double x0 = function.get(DIM_X, i);
retFunction.add(x0, 0, 0, 0);
}
// // derivative for last point
// final int last = ncount - 1;
// final double deltaYN = function.get(DIM_Y, last) - function.get(DIM_Y, last - 1);
// final double deltaXN = function.get(DIM_X, last) - function.get(DIM_X, last - 1);
// retFunction.add(function.get(DIM_X, last), sign * (deltaYN / deltaXN), 0,
// 0);
return retFunction;
}
@SafeVarargs
public static DataSet divideFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.DIVIDE, format);
}
@SafeVarargs
public static DataSet divideFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.DIVIDE, format);
}
/**
* convenience short-hand notation for getting error variables (if defined for dataset)
*
* @param dataSet the source data set
* @param eType the error type
* @param x the data set x-value for which the error should be interpolated
* @return the given interpolated error
*/
public static double error(final DataSet dataSet, final ErrType eType, final double x) {
return error(dataSet, eType, -1, x, true);
}
/**
* convenience short-hand notation for getting error variables (if defined for dataset)
*
* @param dataSet the source data set
* @param eType the error type
* @param index the data set index
* @return the given error
*/
public static double error(final DataSet dataSet, final ErrType eType, final int index) {
return error(dataSet, eType, index, 0.0, false);
}
/**
* convenience short-hand notation for getting error variables (if defined for dataset)
*
* @param dataSet the source data set
* @param eType the error type
* @return the given error array (cropped to data set length if necessary)
*/
public static double[] errors(final DataSet dataSet, final ErrType eType) {
final int nDim = dataSet.getDataCount();
if (!(dataSet instanceof DataSetError)) {
// data set does not have any error definition
return new double[nDim];
}
final DataSetError ds = (DataSetError) dataSet;
switch (eType) {
case EXN:
return cropToLength(ds.getErrorsNegative(DIM_X), nDim);
case EXP:
return cropToLength(ds.getErrorsPositive(DIM_X), nDim);
case EYN:
return cropToLength(ds.getErrorsNegative(DIM_Y), nDim);
case EYP:
default:
return cropToLength(ds.getErrorsPositive(DIM_Y), nDim);
}
}
@SafeVarargs
public static DataSet filterFunction(final DataSet function, final double width, final Filter filterType, @NotNull final Formatter<Number>... format) {
final int n = function.getDataCount();
final String dataSetName = getFormatter(format).format("{0}({1},{2})", filterType.getTag(), function.getName(), width);
final var filteredFunction = new DoubleErrorDataSet(dataSetName, n);
for (var dim = 0; dim < filteredFunction.getDimension(); dim++) {
final var refAxisDescription = function.getAxisDescription(dim);
filteredFunction.getAxisDescription(dim).set(refAxisDescription.getName(), refAxisDescription.getUnit());
}
final var subArrayY = new double[n];
final var subArrayYn = new double[n];
final var subArrayYp = new double[n];
final double[] xValues = function.getValues(DIM_X);
final double[] yValues = function.getValues(DIM_Y);
final double[] yen = errors(function, EYN);
final double[] yep = errors(function, EYN);
for (var i = 0; i < n; i++) {
final double time0 = xValues[i];
var count = 0;
for (var j = 0; j < n; j++) {
final double time = xValues[j];
if (MathBase.abs(time0 - time) <= width) {
subArrayY[count] = yValues[j];
subArrayYn[count] = yen[j];
subArrayYp[count] = yep[j];
count++;
}
}
final double norm = count > 0 ? 1.0 / MathBase.sqrt(count) : 0.0;
switch (filterType) {
case MEDIAN:
filteredFunction.add(time0, Math.median(subArrayY, count), Math.median(subArrayYn, count),
Math.median(subArrayYp, count));
break;
case MIN:
filteredFunction.add(time0, Math.minimum(subArrayY, count), Math.minimum(subArrayYn, count),
Math.minimum(subArrayYp, count));
break;
case MAX:
filteredFunction.add(time0, Math.maximum(subArrayY, count), Math.maximum(subArrayYn, count),
Math.maximum(subArrayYp, count));
break;
case P2P:
filteredFunction.add(time0, Math.peakToPeak(subArrayY, count), Math.peakToPeak(subArrayYn, count),
Math.peakToPeak(subArrayYp, count));
break;
case RMS:
filteredFunction.add(time0, Math.rms(subArrayY, count), Math.rms(subArrayYn, count),
Math.rms(subArrayYp, count));
break;
case GEOMMEAN:
filteredFunction.add(time0, Math.geometricMean(subArrayY, 0, count),
Math.geometricMean(subArrayYn, 0, count), Math.geometricMean(subArrayYp, 0, count));
break;
case MEAN:
default:
filteredFunction.add(time0, Math.mean(subArrayY, count), Math.mean(subArrayYn, count) * norm,
Math.mean(subArrayYp, count) * norm);
break;
}
}
return filteredFunction;
}
@SafeVarargs
public static DataSet geometricMeanFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.GEOMMEAN, format);
}
@SafeVarargs
public static DataSet getSubRange(final DataSet function, final double xMin, final double xMax, @NotNull final Formatter<Number>... format) {
final int nLength = function.getDataCount();
final String dataSetName = getFormatter(format).format("subRange({0}, {1})", xMin, xMax);
final var ret = new DoubleErrorDataSet(dataSetName, nLength);
for (var i = 0; i < nLength; i++) {
final double x = function.get(DIM_X, i);
final double y = function.get(DIM_Y, i);
final double ex = error(function, EXP, i);
final double ey = error(function, EYP, i);
if (x >= xMin && x <= xMax) {
ret.add(x, y, ex, ey);
}
}
return ret;
}
@SafeVarargs
public static DataSet iirLowPassFilterFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
final int n = function.getDataCount();
final String dataSetName = getFormatter(format).format("iir{0}({1},{2})", Filter.MEAN.getTag(), function.getName(), width);
final var filteredFunction = new DoubleErrorDataSet(dataSetName, n);
if (n <= 1) {
if (!(function instanceof GridDataSet)) {
filteredFunction.set(function);
return filteredFunction;
}
filteredFunction.set(function.getValues(DIM_X), function.getValues(DIM_Y), errors(function, EYN),
errors(function, EYP));
for (var index = 0; index < function.getDataCount(); index++) {
final String label = function.getDataLabel(index);
if (label != null && !label.isEmpty()) {
filteredFunction.addDataLabel(index, label);
}
}
for (var index = 0; index < function.getDataCount(); index++) {
final String style = function.getStyle(index);
if (style != null && !style.isEmpty()) {
filteredFunction.addDataStyle(index, style);
}
}
filteredFunction.setStyle(function.getStyle());
return filteredFunction;
}
final double[] xValues = function.getValues(DIM_X);
final double[] yValues = function.getValues(DIM_Y);
final double[] yen = errors(function, EYN);
final double[] yep = errors(function, EYN);
final var yUp = new double[n];
final var yDown = new double[n];
final var ye1 = new double[n];
final var ye2 = new double[n];
final double smoothing = 0.5 * width;
// IIR smoothing algorithm:
// smoothed += elapsedTime * ( newValue - smoothed ) / smoothing
double smoothed = yValues[0];
double smoothed2 = smoothed * smoothed;
// for (int i = 1; i < n; i++) {
// final double x0 = xValues[i - 1];
// final double x1 = xValues[i];
// final double y = yValues[i];
// smoothed += (x1 - x0) * (y - smoothed) / smoothing;
// smoothed2 += (x1 - x0) * (y * y - smoothed2) / smoothing;
// final double newEYN = MathBase.sqrt(MathBase.abs(smoothed2 - smoothed *
// smoothed) + yen[i] * yen[i]);
// final double newEYP = MathBase.sqrt(MathBase.abs(smoothed2 - smoothed *
// smoothed) + yep[i] * yep[i]);
//
// filteredFunction.add(x1 - smoothing, smoothed, newEYN, newEYP);
// }
// calculate forward/backward to compensate for the IIR group-delay
for (var i = 1; i < n; i++) {
final double x0 = xValues[i - 1];
final double x1 = xValues[i];
final double y = yValues[i];
smoothed += (x1 - x0) * (y - smoothed) / smoothing;
smoothed2 += (x1 - x0) * (y * y - smoothed2) / smoothing;
yUp[i] = smoothed;
ye1[i] = smoothed2;
}
smoothed = yValues[n - 1];
smoothed2 = smoothed * smoothed;
for (var i = n - 2; i >= 0; i--) {
final double x0 = xValues[i];
final double x1 = xValues[i + 1];
final double y = yValues[i];
smoothed += (x1 - x0) * (y - smoothed) / smoothing;
smoothed2 += (x1 - x0) * (y * y - smoothed2) / smoothing;
yDown[i] = smoothed;
ye2[i] = smoothed2;
}
filteredFunction.add(xValues[0], yValues[0], yen[0], yep[0]);
for (var i = 1; i < n; i++) {
final double x1 = xValues[i];
final double y = 0.5 * (yUp[i] + yDown[i]);
final double mean2 = y * y;
final double y2 = 0.5 * MathBase.pow(ye1[i] + ye2[i], 1);
final double avgError2 = MathBase.abs(y2 - mean2);
final double newEYN = MathBase.sqrt(avgError2 + yen[i] * yen[i]);
final double newEYP = MathBase.sqrt(avgError2 + yep[i] * yep[i]);
filteredFunction.add(x1, y, newEYN, newEYP);
}
return filteredFunction;
}
public static DoublePointError integral(final DataSet function) {
final DataSet integratedFunction = integrateFunction(function);
final var lastPoint = integratedFunction.getDataCount() - 1;
if (lastPoint <= 0) {
return new DoublePointError(0.0, 0.0, 0.0, 0.0);
}
final double x = integratedFunction.get(DIM_X, lastPoint);
final double y = integratedFunction.get(DIM_Y, lastPoint);
final double yen = error(integratedFunction, EYN, lastPoint);
final double yep = error(integratedFunction, EYP, lastPoint);
final double ye = 0.5 * (yen + yep);
return new DoublePointError(x, y, 0.0, ye);
}
public static DoublePointError integral(final DataSet function, final double xMin, final double xMax) {
final DataSet integratedFunction = integrateFunction(function, xMin, xMax);
final var lastPoint = integratedFunction.getDataCount() - 1;
final double yen = error(integratedFunction, EYN, lastPoint);
final double yep = error(integratedFunction, EYP, lastPoint);
final double ye = 0.5 * (yen + yep);
return new DoublePointError(integratedFunction.get(DIM_X, lastPoint), integratedFunction.get(DIM_Y, lastPoint), 0.0, ye);
}
public static double integralSimple(final DataSet function) {
var integral1 = 0.0;
var integral2 = 0.0;
final int nCount = function.getDataCount();
if (nCount <= 1) {
return 0.0;
}
for (var i = 1; i < nCount; i++) {
final double step = function.get(DIM_X, i) - function.get(DIM_X, i - 1);
final double val1 = function.get(DIM_Y, i - 1);
final double val2 = function.get(DIM_Y, i);
integral1 += step * val1;
integral2 += step * val2;
}
return 0.5 * (integral1 + integral2);
}
@SafeVarargs
public static DataSet integrateFunction(final DataSet function, @NotNull final Formatter<Number>... format) {
return integrateFunction(function, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, format);
}
@SafeVarargs
public static DataSet integrateFunction(final DataSet function, final double xMin, final double xMax, @NotNull final Formatter<Number>... format) {
final int nLength = function.getDataCount();
final var pattern = "{0}({1})dyn|_'{'{2}'}'^'{'{3}'}'";
final String dataSetName = getFormatter(format).format(pattern, INTEGRAL_SYMBOL, function.getName(), xMin, xMax);
if (nLength <= 0) {
if (!(function instanceof GridDataSet) || function.getDimension() > 2) {
final int ncount = function.getDataCount();
final var emptyVector = new double[ncount];
return new DoubleErrorDataSet(dataSetName, function.getValues(DIM_X), emptyVector, emptyVector,
emptyVector, ncount, true);
}
throw new IllegalStateException("not yet implemented for non 2D dataSets");
}
if (!function.getAxisDescription(DIM_X).isDefined()) {
function.recomputeLimits(DIM_X);
}
double xMinLocal = function.getAxisDescription(DIM_X).getMin();
double xMaxLocal = function.getAxisDescription(DIM_X).getMax();
var sign = 1.0;
if (Double.isFinite(xMin) && Double.isFinite(xMax)) {
xMinLocal = min(xMin, xMax);
xMaxLocal = max(xMin, xMax);
if (xMin > xMax) {
sign = -1;
}
} else if (Double.isFinite(xMin)) {
xMinLocal = xMin;
} else if (Double.isFinite(xMax)) {
xMaxLocal = xMax;
}
final var retFunction = new DoubleErrorDataSet(dataSetName, nLength);
if (nLength <= 1) {
retFunction.add(function.get(DIM_X, 0), 0, 0, 0);
return retFunction;
}
var integral = 0.0;
var integralEN = 0.0;
var integralEP = 0.0;
if (Double.isFinite(xMin) && xMin <= function.get(DIM_X, 0)) {
// interpolate before range where discrete function is defined
final double val1 = function.getValue(DIM_Y, xMin);
final double x1 = function.get(DIM_X, 0);
final double val2 = function.get(DIM_Y, 0);
final double step = x1 - xMin;
integral += sign * 0.5 * step * (val1 + val2);
final double en1 = error(function, EYN, 0);
final double ep1 = error(function, EYP, 0);
// assuming uncorrelated errors between bins
integralEN = MathBase.hypot(integralEN, step * en1);
integralEP = MathBase.hypot(integralEP, step * ep1);
retFunction.add(xMin, integral, 0, 0);
}
retFunction.add(function.get(DIM_X, 0), integral, integralEN, integralEP);
for (var i = 1; i < nLength; i++) {
final double x0 = function.get(DIM_X, i - 1);
final double x1 = function.get(DIM_X, i);
double step = x1 - x0;
final double y0 = function.get(DIM_Y, i - 1);
final double en1 = error(function, EYN, i - 1);
final double ep1 = error(function, EYP, i - 1);
final double y1 = function.get(DIM_Y, i);
final double en2 = error(function, EYN, i);
final double ep2 = error(function, EYP, i);
// simple triangulation integration
if (x0 >= xMinLocal && x1 <= xMaxLocal) {
integral += sign * 0.5 * step * (y0 + y1);
// assuming uncorrelated errors between bins
integralEN = MathBase.hypot(integralEN, 0.5 * step * (en1 + en2));
integralEP = MathBase.hypot(integralEP, 0.5 * step * (ep1 + ep2));
} else if (x1 < xMinLocal && x0 < xMinLocal) { // NOSONAR NOPMD
// see below
} else if (x0 < xMinLocal && x1 > xMinLocal) {
retFunction.add(xMin, integral, integralEN, integralEP);
step = x1 - xMinLocal;
integral += sign * 0.5 * step * (function.getValue(DIM_Y, xMinLocal) + y1);
// assuming uncorrelated errors between bins
integralEN = MathBase.hypot(integralEN, 0.5 * step * (en1 + en2));
integralEP = MathBase.hypot(integralEP, 0.5 * step * (ep1 + ep2));
} else if (x0 < xMaxLocal && x1 > xMaxLocal) {
step = xMaxLocal - x0;
final double yAtMax = function.getValue(DIM_Y, xMaxLocal);
integral += sign * 0.5 * step * (y0 + yAtMax);
// assuming uncorrelated errors between bins
integralEN = MathBase.hypot(integralEN, 0.5 * step * (en1 + en2));
integralEP = MathBase.hypot(integralEP, 0.5 * step * (ep1 + ep2));
retFunction.add(xMaxLocal, integral, integralEN, integralEP);
}
retFunction.add(x1, integral, integralEN, integralEP);
}
if (Double.isFinite(xMax) && xMax > function.get(DIM_X, nLength - 1)) {
// interpolate after range where discrete function is defined
final double x0 = function.get(DIM_X, nLength - 1);
final double val1 = function.get(DIM_Y, nLength - 1);
final double val2 = function.getValue(DIM_Y, xMax);
final double step = xMax - x0;
final double en1 = error(function, EYN, nLength - 1);
final double ep1 = error(function, EYP, nLength - 1);
// assuming uncorrelated errors between bins
integralEN = MathBase.hypot(integralEN, step * en1);
integralEP = MathBase.hypot(integralEP, step * ep1);
integral += 0.5 * step * (val1 + val2);
retFunction.add(xMax, integral, integralEN, integralEP);
}
return retFunction;
}
@SafeVarargs
public static DataSet integrateFromCentre(@NotNull DataSet function, final double centre, final double width, final boolean normalise, @NotNull final Formatter<Number>... format) {
final double xMinLocal = function.getAxisDescription(DIM_X).getMin();
final double xMaxLocal = function.getAxisDescription(DIM_X).getMax();
final double centreLocal = Double.isFinite(centre) ? centre : SimpleDataSetEstimators.computeCentreOfMass(function, 0, function.getDataCount());
final var pattern = "{0}({1},c={2})dyn|_'{'c-{3}'}'^'{'c+{3}'}'";
final int nLength = function.getDataCount();
final String dataSetName = getFormatter(format).format(pattern, INTEGRAL_SYMBOL, function.getName(), centreLocal, width);
if (nLength < 2) { // need at least two
if (!(function instanceof GridDataSet) || function.getDimension() > 2) {
final int ncount = function.getDataCount();
final var emptyVector = new double[ncount];
return new DoubleErrorDataSet(dataSetName, function.getValues(DIM_X), emptyVector, emptyVector, emptyVector, ncount, true);
}
throw new IllegalStateException("not yet implemented for non 2D dataSets");
}
final var retFunction = new DoubleErrorDataSet(dataSetName, nLength / 2);
if (width <= 0.0) {
return retFunction;
}
if (centreLocal <= xMinLocal || centreLocal >= xMaxLocal) {
throw new IllegalArgumentException(String.format("centre %f is outside DataSetRange [%f,%f]", centreLocal, xMinLocal, xMaxLocal));
}
final double scaleLocal = normalise ? 1.0 / integral(function, max(centreLocal - width, xMinLocal), min(centreLocal + width, xMaxLocal)).getY() : 1.0;
final var centreLocalIndex = function.getIndex(DIM_X, centreLocal);
final var maxIntIndex = min(centreLocalIndex, function.getDataCount() - 1 - centreLocalIndex);
var integral = 0.0;
for (var i = 0; i < maxIntIndex; i++) {
// x-values
final double xL0 = function.get(DIM_X, centreLocalIndex - i - 1);
final double xL1 = function.get(DIM_X, centreLocalIndex - i);
final double stepL = xL1 - xL0;
final double xR0 = function.get(DIM_X, centreLocalIndex + i);
final double xR1 = function.get(DIM_X, centreLocalIndex + i + 1);
final double stepR = xR1 - xR0;
// y-values
final double yL0 = function.get(DIM_Y, centreLocalIndex - i - 1);
final double yL1 = function.get(DIM_Y, centreLocalIndex - i);
final double yR0 = function.get(DIM_Y, centreLocalIndex + i);
final double yR1 = function.get(DIM_Y, centreLocalIndex + i + 1);
// simple triangulation integration
integral += scaleLocal * 0.5 * (stepL * (yL0 + yL1) + stepR * (yR0 + yR1));
final double distanceFromCentre = 0.5 * (xR1 - xL0);
retFunction.add(distanceFromCentre, integral, 0, 0);
}
return retFunction;
}
public static double integralWidth(@NotNull DataSet function, final double centre, final double maxWidth, final double threshold) {
return SimpleDataSetEstimators.getZeroCrossing(integrateFromCentre(function, centre, maxWidth, true), threshold);
}
@SafeVarargs
public static DataSet inversedbFunction(final DataSet function, @NotNull final Formatter<Number>... format) {
return mathFunction(function, 1.0, MathOp.INV_DB, format);
}
@SafeVarargs
public static DataSet log10Function(final DataSet function, @NotNull final Formatter<Number>... format) {
return mathFunction(function, 0.0, MathOp.LOG10, format);
}
@SafeVarargs
public static DataSet log10Function(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.LOG10, format);
}
@SafeVarargs
public static DataSet lowPassFilterFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.MEAN, format);
}
@SafeVarargs
public static DataSet magnitudeSpectrum(final DataSet function, @NotNull final Formatter<Number>... format) {
return magnitudeSpectrum(function, Apodization.Hann, false, false, format);
}
@SafeVarargs
public static DataSet magnitudeSpectrum(final DataSet function, final Apodization apodization, final boolean dbScale, final boolean normalisedFrequency, @NotNull final Formatter<Number>... format) {
final String functionName = getFormatter(format).format("Mag{0}({1})", dbScale ? "[dB]" : "", function.getName());
final int n = function.getDataCount();
if (n == 0) {
return new DoubleErrorDataSet(functionName, 0);
}
final var fastFourierTrafo = new DoubleFFT_1D(n);
// N.B. since realForward computes the FFT in-place -> generate a copy
final var fftSpectra = new double[n];
for (var i = 0; i < n; i++) {
final double window = apodization.getIndex(i, n);
fftSpectra[i] = function.get(DIM_Y, i) * window;
}
fastFourierTrafo.realForward(fftSpectra);
final var mag = dbScale ? SpectrumTools.computeMagnitudeSpectrum_dB(fftSpectra, true) : SpectrumTools.computeMagnitudeSpectrum(fftSpectra, true);
final var dt = function.get(DIM_X, function.getDataCount() - 1) - function.get(DIM_X, 0);
final var fsampling = normalisedFrequency || dt <= 0 ? 0.5 / mag.length : 1.0 / dt;
final var ret = new DoubleErrorDataSet(functionName, mag.length);
for (var i = 0; i < mag.length; i++) {
// TODO: consider magnitude error estimate
ret.add(i * fsampling, mag[i], 0, 0);
}
return ret;
}
@SafeVarargs
public static DataSet magnitudeSpectrumComplex(final DataSet function, @NotNull final Formatter<Number>... format) {
return magnitudeSpectrumComplex(function, Apodization.Hann, false, false, format);
}
@SafeVarargs
public static DataSet magnitudeSpectrumComplex(final DataSet function, final Apodization apodization,
final boolean dbScale, final boolean normalisedFrequency, @NotNull final Formatter<Number>... format) {
final int n = function.getDataCount();
final var fastFourierTrafo = new DoubleFFT_1D(n);
// N.B. since realForward computes the FFT in-place -> generate a copy
final var fftSpectra = new double[2 * n];
for (var i = 0; i < n; i++) {
final double window = apodization.getIndex(i, n);
fftSpectra[2 * i] = function.get(DIM_Y, i) * window;
fftSpectra[2 * i + 1] = function.get(DIM_Z, i) * window;
}
fastFourierTrafo.complexForward(fftSpectra);
final var mag = dbScale ? SpectrumTools.computeMagnitudeSpectrum_dB(fftSpectra, true)
: SpectrumTools.computeMagnitudeSpectrum(fftSpectra, true);
final var dt = function.get(DIM_X, function.getDataCount() - 1) - function.get(DIM_X, 0);
final var fsampling = normalisedFrequency || dt <= 0 ? 0.5 / mag.length : 1.0 / dt;
final var functionName = getFormatter(format).format("Mag{0}({1})", dbScale ? "[dB]" : "", function.getName());
final var ret = new DoubleErrorDataSet(functionName, mag.length);
for (var i = 0; i < mag.length; i++) {
// TODO: consider magnitude error estimate
if (i < mag.length / 2) {
ret.add((i - mag.length / 2.0) * fsampling, mag[i + mag.length / 2], 0, 0);
} else {
ret.add((i - mag.length / 2.0) * fsampling, mag[i - mag.length / 2], 0, 0);
}
}
return ret;
}
public static DataSet magnitudeSpectrumDecibel(final DataSet function) {
return magnitudeSpectrum(function, Apodization.Hann, true, false);
}
public static boolean sameHorizontalBase(final DataSet function1, final DataSet function2) {
if (function1.getDataCount() != function2.getDataCount()) {
return false;
}
for (var i = 0; i < function1.getDataCount(); i++) {
final double X1 = function1.get(DIM_X, i);
final double X2 = function2.get(DIM_X, i);
if (X1 != X2) {
return false;
}
}
return true;
}
@SafeVarargs
public static DataSet mathFunction(final DataSet function1, final DataSet function2, final MathOp op, @NotNull final Formatter<Number>... format) {
final String newDataSetName = getFormatter(format).format("{0}{1}{2}", function1.getName(), op.getTag(), function2.getName());
final var ret = new DoubleErrorDataSet(newDataSetName, function1.getDataCount());
ret.getAxisDescription(DIM_X).set(function1.getAxisDescription(DIM_X));
ret.getAxisDescription(DIM_Y).set(function1.getAxisDescription(DIM_Y).getName(), function1.getAxisDescription(DIM_Y).getUnit());
final boolean needsInterpolation = !sameHorizontalBase(function1, function2);
if (needsInterpolation) {
final List<Double> xValues = getCommonBase(function1, function2);
for (double x : xValues) {
final double Y1 = function1.getValue(DIM_Y, x);
final double Y2 = function2.getValue(DIM_Y, x);
final double eyn1 = error(function1, EYN, 0, x, needsInterpolation);
final double eyp1 = error(function1, EYP, 0, x, needsInterpolation);
final double eyn2 = error(function2, EYN, 0, x, needsInterpolation);
final double eyp2 = error(function2, EYP, 0, x, needsInterpolation);
applyMathOperation(ret, op, x, Y1, Y2, eyn1, eyp1, eyn2, eyp2);
}
return ret;
}
for (var i = 0; i < function1.getDataCount(); i++) {
final double X1 = function1.get(DIM_X, i);
// not needed : final double X2 = function1.get(DIM_X, i) ...
final double Y1 = function1.get(DIM_Y, i);
final double Y2 = function2.get(DIM_Y, i);
final double eyn1 = error(function1, EYN, i);
final double eyp1 = error(function1, EYP, i);
final double eyn2 = error(function2, EYN, i, X1, needsInterpolation);
final double eyp2 = error(function2, EYP, i, X1, needsInterpolation);
applyMathOperation(ret, op, X1, Y1, Y2, eyn1, eyp1, eyn2, eyp2);
}
return ret;
}
public static List<Double> getCommonBase(final DataSet... functions) {
final List<Double> xValues = new NoDuplicatesList<>();
for (DataSet function : functions) {
for (var i = 0; i < function.getDataCount(); i++) {
if (function instanceof Histogram) {
xValues.add(((Histogram) function).getBinLimits(DIM_X, LOWER, i));
xValues.add(((Histogram) function).getBinCenter(DIM_X, i));
xValues.add(((Histogram) function).getBinLimits(DIM_X, UPPER, i));
} else {
xValues.add(function.get(DIM_X, i));
}
}
}
Collections.sort(xValues);
return xValues;
}
@SafeVarargs
public static DataSet mathFunction(final DataSet function, final double value, final MathOp op, @NotNull final Formatter<Number>... format) {
final String functionName = getFormatter(format).format("{0}({1})", op.getTag(), function.getName());
final double[] y = function.getValues(DIM_Y);
final double[] eyn = Arrays.copyOf(errors(function, EYN), y.length);
final double[] eyp = Arrays.copyOf(errors(function, EYP), y.length);
final int ncount = function.getDataCount();
switch (op) {
case ADD:
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.add(y, value), eyn, eyp,
ncount, true);
case SUBTRACT:
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.subtract(y, value), eyn, eyp,
ncount, true);
case MULTIPLY:
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.multiply(y, value),
ArrayMath.multiply(eyn, value), ArrayMath.multiply(eyp, value), ncount, true);
case DIVIDE:
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.divide(y, value),
ArrayMath.divide(eyn, value), ArrayMath.divide(eyp, value), ncount, true);
case SQR:
for (var i = 0; i < eyn.length; i++) {
eyn[i] = 2 * MathBase.abs(y[i] + value) * eyn[i];
eyp[i] = 2 * MathBase.abs(y[i] + value) * eyp[i];
}
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), value == 0.0 ? ArrayMath.sqr(y) : ArrayMath.sqr(ArrayMath.add(y, value)), eyn, eyp, ncount,
true);
case SQRT:
for (var i = 0; i < eyn.length; i++) {
eyn[i] = MathBase.sqrt(MathBase.abs(y[i] + value)) * eyn[i];
eyp[i] = MathBase.sqrt(MathBase.abs(y[i] + value)) * eyp[i];
}
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), value == 0.0 ? ArrayMath.sqrt(y) : ArrayMath.sqrt(ArrayMath.add(y, value)), eyn, eyp, ncount,
true);
case LOG10:
for (var i = 0; i < eyn.length; i++) {
eyn[i] = 0.0; // 0.0 as a work-around
eyp[i] = 0.0;
}
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.tenLog10(y), eyn, eyp,
ncount, true);
case DB:
for (var i = 0; i < eyn.length; i++) {
eyn[i] = 0.0; // 0.0 as a work-around
eyp[i] = 0.0;
}
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.decibel(y), eyn, eyp, ncount,
true);
case INV_DB:
for (var i = 0; i < eyn.length; i++) {
eyn[i] = 0.0; // 0.0 as a work-around
eyp[i] = 0.0;
}
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), ArrayMath.inverseDecibel(y), eyn, eyp,
ncount, true);
case IDENTITY:
default:
// return copy if nothing else matches
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), function.getValues(DIM_Y),
errors(function, EYN), errors(function, EYP), ncount, true);
}
}
@SafeVarargs
public static DataSet maxFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.MAX, format);
}
@SafeVarargs
public static DataSet medianFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.MEDIAN, format);
}
@SafeVarargs
public static DataSet minFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.MIN, format);
}
@SafeVarargs
public static DataSet multiplyFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.MULTIPLY, format);
}
@SafeVarargs
public static DataSet multiplyFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.MULTIPLY, format);
}
@SafeVarargs
public static DataSet normalisedFunction(final DataSet function, @NotNull final Formatter<Number>... format) {
return normalisedFunction(function, 1.0, format);
}
@SafeVarargs
public static DataSet normalisedFunction(final DataSet function, final double requiredIntegral, @NotNull final Formatter<Number>... format) {
final DoublePointError complexInt = integral(function);
final double integral = complexInt.getY() / requiredIntegral;
final int ncount = function.getDataCount();
// final double integralErr = complexInt.getErrorY() / requiredIntegral;
// TODO: add error propagation to normalised function error estimate
final String functionName = getFormatter(format).format("{0}", function.getName());
if (integral == 0) {
return new DoubleErrorDataSet(functionName, function.getValues(DIM_X), new double[ncount],
new double[ncount], new double[ncount], ncount, true);
}
final var xValues = function.getValues(DIM_X);
final var yValues = ArrayMath.divide(function.getValues(DIM_Y), integral);
final var eyp = ArrayMath.divide(errors(function, EYN), integral);
final var eyn = ArrayMath.divide(errors(function, EYP), integral);
return new DoubleErrorDataSet(functionName, xValues, yValues, eyp, eyn, ncount, true);
}
@SafeVarargs
public static DataSet normalisedMagnitudeSpectrumDecibel(final DataSet function, @NotNull final Formatter<Number>... format) {
return magnitudeSpectrum(function, Apodization.Hann, true, true, format);
}
@SafeVarargs
public static DataSet peakToPeakFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.P2P, format);
}
@SafeVarargs
public static DataSet rmsFilteredFunction(final DataSet function, final double width, @NotNull final Formatter<Number>... format) {
return filterFunction(function, width, Filter.RMS, format);
}
public static EditableDataSet setFunction(final EditableDataSet function, final double value, final double xMin, final double xMax) {
final int nLength = function.getDataCount();
double xMinLocal = function.get(DIM_X, 0);
double xMaxLocal = function.get(DIM_X, nLength - 1);
if (Double.isFinite(xMin) && Double.isFinite(xMax)) {
xMinLocal = min(xMin, xMax);
xMaxLocal = max(xMin, xMax);
} else if (Double.isFinite(xMin)) {
xMinLocal = xMin;
} else if (Double.isFinite(xMax)) {
xMaxLocal = xMax;
}
for (var i = 0; i < nLength; i++) {
final double x = function.get(DIM_X, i);
if (x >= xMinLocal && x <= xMaxLocal) {
function.set(i, x, value);
}
}
return function;
}
@SafeVarargs
public static DataSet sqrFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.SQR, format);
}
@SafeVarargs
public static DataSet sqrFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.SQR, format);
}
@SafeVarargs
public static DataSet sqrtFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.SQRT, format);
}
@SafeVarargs
public static DataSet sqrtFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.SQRT, format);
}
@SafeVarargs
public static DataSet subtractFunction(final DataSet function1, final DataSet function2, @NotNull final Formatter<Number>... format) {
return mathFunction(function1, function2, MathOp.SUBTRACT, format);
}
@SafeVarargs
public static DataSet subtractFunction(final DataSet function, final double value, @NotNull final Formatter<Number>... format) {
return mathFunction(function, value, MathOp.SUBTRACT, format);
}
@SafeVarargs
private static Formatter<Number> getFormatter(@NotNull final Formatter<Number>... format) {
return Objects.requireNonNull(format, "user-supplied format").length > 0 ? format[0] : DEFAULT_FORMATTER;
}
private static double[] cropToLength(final double[] in, final int length) {
// small helper routine to crop data array in case it's to long
if (in.length == length) {
return in;
}
return Arrays.copyOf(in, length);
}
public enum ErrType {
EXN,
EXP,
EYN,
EYP
}
public enum Filter {
MEAN("LowPass"),
MEDIAN("Median"),
MIN("Min"),
MAX("Max"),
P2P("PeakToPeak"),
RMS("RMS"),
GEOMMEAN("GeometricMean");
private final String tag;
Filter(final String tag) {
this.tag = tag;
}
public String getTag() {
return tag;
}
}
public enum MathOp {
ADD("+"),
SUBTRACT("-"),
MULTIPLY("*"),
DIVIDE("*"),
SQR("SQR"),
SQRT("SQRT"),
LOG10("Log10"),
DB("dB"),
INV_DB("dB^{-1}"),
IDENTITY("Copy");
private final String tag;
MathOp(final String tag) {
this.tag = tag;
}
public String getTag() {
return tag;
}
}
} |
Beta Was this translation helpful? Give feedback.
0 replies
-
|
For #641, seems it suffices to have <dependency>
<groupId>io.fair-acc</groupId>
<artifactId>dataset</artifactId>
<version>11.3.0</version>
</dependency>
<dependency>
<groupId>io.fair-acc</groupId>
<artifactId>math</artifactId>
<version>11.3.0</version>
</dependency>
<dependency>
<groupId>io.fair-acc</groupId>
<artifactId>bench</artifactId>
<version>11.3.0</version>
</dependency>while still debugging. |
Beta Was this translation helpful? Give feedback.
0 replies
Answer selected by
yezhengli-Mr9
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
For #641, seems it suffices to have
while still debugging.