List of usage examples for org.apache.commons.math3.special Erf erfc
public static double erfc(double x)
From source file:etomica.potential.EwaldSummation.java
public double uReal() { int nAtoms = box.getLeafList().getAtomCount(); double uReal = 0.0; for (int i = 0; i < nAtoms; i++) {//H IAtom atomA = box.getLeafList().getAtom(i); double chargeA = atomAgentManager.getAgent(atomA).charge; if (chargeA == 0) continue; int aIndex = atomA.getParentGroup().getIndex(); IVectorMutable positionA = atomA.getPosition(); for (int j = i; j < nAtoms; j++) { IAtom atomB = box.getLeafList().getAtom(j); int bIndex = atomB.getParentGroup().getIndex(); if (aIndex == bIndex && nRealShells[0] == 0 && nRealShells[1] == 0 && nRealShells[2] == 0) continue;//Skip atom-pairs in the same molecule in the orig. cell. double chargeB = atomAgentManager.getAgent(atomB).charge; if (chargeB == 0) continue; IVectorMutable positionB = atomB.getPosition(); rAB.Ev1Mv2(positionA, positionB);// get vector rAB box.getBoundary().nearestImage(rAB);// minimum image boolean isSelf = i == j; for (int nx = -nRealShells[0]; nx <= nRealShells[0]; nx++) { Lxyz.setX(0, nx * boxSize[0]); for (int ny = -nRealShells[1]; ny <= nRealShells[1]; ny++) { Lxyz.setX(1, ny * boxSize[1]); for (int nz = -nRealShells[2]; nz <= nRealShells[2]; nz++) { boolean centerImage = nx * nx + ny * ny + nz * nz == 0; if (aIndex == bIndex && centerImage) continue;//Skip atom-pairs in the same molecule in the orig. cell & ignores self+centerImage too Lxyz.setX(2, nz * boxSize[2]); drTmp.Ev1Pv2(rAB, Lxyz); double r2 = drTmp.squared(); if (r2 > rCutSquared) continue; double drTmpM = Math.sqrt(r2); double tmepReal = chargeA * chargeB * Erf.erfc(alpha * drTmpM) / drTmpM;//Don't worry about 1/2 factor;j>i uReal += (isSelf ? 0.5 : 1.0) * tmepReal; }//from w ww.j a v a 2 s . c o m } } } // close for all sites in j-th molecule } // close for the outside loop return uReal; }
From source file:edu.duke.cs.osprey.voxq.QuadraticQFunction.java
private boolean erfcInvNumericsOK(double x, double refRange) { //some points near edges of bell curve make erfc inversion innaccurate. Detect this. //refRange used to determine tolerance x = Math.abs(x);//from w w w .java 2 s . co m double invResult = myErfcInv(Erf.erfc(x)); if (Double.isInfinite(invResult))//definitely going to be a problem return false; else if (Math.abs(x - invResult) > 0.01 * refRange) return false; return true; }
From source file:dsp.unige.figures.ChannelHelper.java
/** * Returns the bit error rate when NOT using any FEC code * // ww w . ja va 2s . co m * @param sta the station object * @param sat the satellite object * @param rate the information rate * @return the ber rate in (0,0.5) */ public static double getUncodedBER(Station sta, Satellite sat, double rate) { double SdBW, Eb, N0, EbN0; SdBW = getSdBW(sta, sat); Eb = 10 * Math.log10(Math.pow(10, SdBW / 10d) / (rate * 1000d)); N0 = getN0dBW(sta, sat); EbN0 = Eb - N0; double ber = 0.5 * Erf.erfc(Math.sqrt(EbN0)); if (!Double.isNaN(ber) && ber < 0.5) return ber; else return 0.5; }
From source file:edu.duke.cs.osprey.voxq.QuadraticQFunction.java
private double cumulDistrInv(double F) { //inverse of cumulDistr if (Math.abs(a) < 1e-14) {//linear case double ans = (Math.log(b * F + Math.exp(b * xLo + c)) - c) / b; if (ans < xLo - 1e-6 || ans > xHi + 1e-6)//DEBUG!!! System.out.println("Out of range QuadraticQFunction draw...");//DEBUG!!! //double Fcheck = cumulDistr(ans);//DEBUG!!! return ans; } else {//w w w. ja v a2s.c o m double C = 0.5 * Math.exp(c - 0.25 * b * b / a) * Math.sqrt(-Math.PI / a); double erfArg1 = cdErfArg(xLo); double erfArg2; if (erfArg1 > 1) {//these three options are mathematically equivalent, //but will be selected to best provide numerical stability double erfcVal = -F / C + Erf.erfc(erfArg1); erfArg2 = myErfcInv(erfcVal); } else if (erfArg1 < -1) { double oerfcVal = F / C + Erf.erfc(-erfArg1); erfArg2 = -myErfcInv(oerfcVal); } else { double erfVal = F / C + Erf.erf(erfArg1); erfArg2 = Erf.erfInv(erfVal); } double denom = 2 * Math.sqrt(-a); double ans = (-erfArg2 * denom - b) / (2 * a); //double Fcheck = cumulDistr(ans);//DEBUG!!! if (ans < xLo - 1 || ans > xHi + 1)//DEBUG!!! System.out.println("Out of range QuadraticQFunction draw...");//DEBUG!!! return ans; } }
From source file:com.intel.diceros.test.securerandom.DRNGTest.java
/** * Discrete Fourier Transform (Spectral) Test * <p>/* w w w .j av a 2 s . c om*/ * The focus of this test is the peak heights in the Discrete Fourier * Transform of the sequence. The purpose of this test is to detect * periodic features (i.e., repetitive patterns that are near each other) * in the tested sequence that would indicate a deviation from the * assumption of randomness. The intention is to detect whether the number * of peaks exceeding the 95% threshold is significantly different than 5%. */ private void testDiscreteFourierTransform(final int[] epsilon) { final int n = epsilon.length; double p_value, upperBound, N_l, N_o, d; double[] m = new double[n / 2 + 1], X = new double[n]; int i, count; for (i = 0; i < n; i++) X[i] = 2 * epsilon[i] - 1; double[] X1 = new double[n]; for (i = 0; i < X.length; i++) { X1[i] = X[i]; } FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD); Complex[] Xc = fft.transform(X, TransformType.FORWARD); m[0] = Math.sqrt(Xc[0].getReal() * Xc[0].getReal()); for (i = 0; i < n / 2; i++) m[i + 1] = Math.sqrt(Math.pow(Xc[2 * i].getReal(), 2) + Math.pow(Xc[2 * i + 1].getReal(), 2)); count = 0; upperBound = Math.sqrt(2.995732274 * n); for (i = 0; i < n / 2; i++) if (m[i] < upperBound) count++; N_l = (double) count; N_o = 0.95 * n / 2.0; d = (N_l - N_o) / Math.sqrt(n / 4.0 * 0.95 * 0.05); p_value = Erf.erfc(Math.abs(d) / Math.sqrt(2.0)); assertTrue("RNG test failed, test discrete fourier transform.", p_value >= 0.01); }
From source file:etomica.potential.EwaldSummation.java
public IVector[] gradient(IAtomList atoms) { int nAtoms = box.getLeafList().getAtomCount(); double coeff = 4.0 * Math.PI / volume; double kCutSquared = kCut * kCut; // criteria for spherical cutoff in fourier space if (gradient.length < nAtoms) { gradient = new IVectorMutable[nAtoms]; sinkrj = new double[nAtoms]; coskrj = new double[nAtoms]; for (int i = 0; i < nAtoms; i++) { gradient[i] = space.makeVector(); }/* ww w . jav a 2 s . co m*/ } else { for (int i = 0; i < nAtoms; i++) { gradient[i].E(0); //for Vecors and scalar sinkrj[i] = 0; coskrj[i] = 0; } } if (doRealSum) { //Real gradient //Cross Interaction for (int i = 0; i < nAtoms; i++) { IAtom atomA = box.getLeafList().getAtom(i); double chargeA = atomAgentManager.getAgent(atomA).charge; if (chargeA == 0) continue; int aIndex = atomA.getParentGroup().getIndex(); // molecule a IVectorMutable positionA = atomA.getPosition(); for (int j = i + 1; j < nAtoms; j++) { IAtom atomB = box.getLeafList().getAtom(j); int bIndex = atomB.getParentGroup().getIndex(); // molecule b if (nRealShells[0] == 0 && nRealShells[1] == 0 && nRealShells[2] == 0 && aIndex == bIndex) continue;//Skip same molecules! double chargeB = atomAgentManager.getAgent(atomB).charge; if (chargeB == 0) continue; IVectorMutable positionB = atomB.getPosition(); rAB.Ev1Mv2(positionA, positionB); //rAB == rA - rB box.getBoundary().nearestImage(rAB); for (int nx = -nRealShells[0]; nx <= nRealShells[0]; nx++) { Lxyz.setX(0, nx * boxSize[0]); for (int ny = -nRealShells[1]; ny <= nRealShells[1]; ny++) { Lxyz.setX(1, ny * boxSize[1]); for (int nz = -nRealShells[2]; nz <= nRealShells[2]; nz++) { if (aIndex == bIndex && nx * nx + ny * ny + nz * nz == 0) continue; Lxyz.setX(2, nz * boxSize[2]); drTmp.Ev1Pv2(rAB, Lxyz); double rAB2 = drTmp.squared(); if (rAB2 > rCutSquared) continue; double rABMagnitude = Math.sqrt(rAB2); double rAB3 = rABMagnitude * rAB2; double B = Erf.erfc(alpha * rABMagnitude) + 2.0 * alpha * rABMagnitude / sqrtPI * Math.exp(-alpha2 * rAB2); double realCoeff = -chargeA * chargeB * B / rAB3; // gradU = -F gradient[i].PEa1Tv1(realCoeff, drTmp); gradient[j].PEa1Tv1(-realCoeff, drTmp); } } } } } } //Fourier gradient Part for (int xAxis = -nKs[0]; xAxis < nKs[0] + 1; xAxis++) { kVector.setX(0, (xAxis * basis[0]));// assign value to the x-axis for (int yAxis = -nKs[1]; yAxis < nKs[1] + 1; yAxis++) { kVector.setX(1, (yAxis * basis[1]));// assign value to the y-axis for (int zAxis = -nKs[2]; zAxis < nKs[2] + 1; zAxis++) { if ((xAxis * xAxis + yAxis * yAxis + zAxis * zAxis) == 0) continue;// first check: k is a non-zero vector kVector.setX(2, (zAxis * basis[2]));// assign value to the z-axis, now the vector is specified double kSquared = kVector.squared(); if (kSquared > kCutSquared) continue;// k-vector should be within the sphere with kCutoff as the radius double sCoskr = 0.0, sSinkr = 0.0; for (int j = 0; j < nAtoms; j++) { // Loop over atoms (4*nBasis) IAtom atom = box.getLeafList().getAtom(j); IVectorMutable position = atom.getPosition(); sinkrj[j] = Math.sin(position.dot(kVector)); coskrj[j] = Math.cos(position.dot(kVector)); double chargej = atomAgentManager.getAgent(atom).charge; sCoskr += chargej * coskrj[j]; sSinkr += chargej * sinkrj[j]; } //End loop over j double coeffk = coeff / kSquared * Math.exp(-kSquared / 4.0 / alpha2); for (int i = 0; i < nAtoms; i++) { IAtom atom = box.getLeafList().getAtom(i); double chargei = atomAgentManager.getAgent(atom).charge; double coeffki = coeffk * chargei * (sinkrj[i] * sCoskr - coskrj[i] * sSinkr); gradient[i].PEa1Tv1(-coeffki, kVector); // gradU = -F } } //end of storing Sin and Cos } } //End loop over ks //Intra-Molecular gradient: for (int i = 0; i < numMolecules; i++) { IMolecule molecule = moleculeList.getMolecule(i); int numSites = molecule.getChildList().getAtomCount(); for (int siteA = 0; siteA < numSites; siteA++) { IAtom atomA = molecule.getChildList().getAtom(siteA); // index = 0, 1, 2, 3|||leafIndex=0...184 double chargeA = atomAgentManager.getAgent(atomA).charge; if (chargeA == 0) continue; IVectorMutable positionA = atomA.getPosition(); for (int siteB = siteA + 1; siteB < numSites; siteB++) { IAtom atomB = molecule.getChildList().getAtom(siteB); double chargeB = atomAgentManager.getAgent(atomB).charge; if (chargeB == 0) continue; IVectorMutable positionB = atomB.getPosition(); rAB.Ev1Mv2(positionA, positionB); box.getBoundary().nearestImage(rAB); double rAB2 = rAB.squared(); double rABMagnitude = Math.sqrt(rAB2); double B = 2 * alpha / sqrtPI * Math.exp(-alpha2 * rAB2) - Erf.erf(alpha * rABMagnitude) / rABMagnitude; double coeffAB = -chargeA * chargeB * B / rAB2; // gradU = -F gradient[atomA.getLeafIndex()].PEa1Tv1(coeffAB, rAB); gradient[atomB.getLeafIndex()].PEa1Tv1(-coeffAB, rAB); } } } return gradient; }
From source file:etomica.potential.EwaldSummation.java
public Tensor secondDerivative(IAtom atom0, IAtom atom1) { double kCutSquared = kCut * kCut; // criteria for spherical cutoff in fourier space IVectorMutable pos0 = atom0.getPosition(); IVectorMutable pos1 = atom1.getPosition(); rAB.Ev1Mv2(pos0, pos1);//from ww w. j av a 2 s.c o m box.getBoundary().nearestImage(rAB); double q0 = atomAgentManager.getAgent(atom0).charge; double q1 = atomAgentManager.getAgent(atom1).charge; secondDerivative.E(0); if (q0 * q1 == 0) return secondDerivative; double coeff = 4.0 * Math.PI / volume; for (int xAxis = -nKs[0]; xAxis < nKs[0] + 1; xAxis++) { kVector.setX(0, (xAxis * basis[0]));// assign value to the x-axis for (int yAxis = -nKs[1]; yAxis < nKs[1] + 1; yAxis++) { kVector.setX(1, (yAxis * basis[1]));// assign value to the y-axis for (int zAxis = -nKs[2]; zAxis < nKs[2] + 1; zAxis++) { if ((xAxis * xAxis + yAxis * yAxis + zAxis * zAxis) == 0) continue;// first check: k is a non-zero vector kVector.setX(2, (zAxis * basis[2]));// assign value to the z-axis, now the vector is specified double kSquared = kVector.squared(); if (kSquared > kCutSquared) continue;// k-vector should be within the sphere with kCutoff as the radius tempTensorkk.Ev1v2(kVector, kVector); tempTensorkk.TE(Math.exp(-kSquared / 4.0 / alpha2) * Math.cos(kVector.dot(rAB)) / kSquared); secondDerivative.PE(tempTensorkk); } } } secondDerivative.TE(coeff); boolean intraMolecular = atom0.getParentGroup() == atom1.getParentGroup(); for (int nx = -nRealShells[0]; nx <= nRealShells[0]; nx++) { Lxyz.setX(0, nx * boxSize[0]); for (int ny = -nRealShells[1]; ny <= nRealShells[1]; ny++) { Lxyz.setX(1, ny * boxSize[1]); for (int nz = -nRealShells[2]; nz <= nRealShells[2]; nz++) { Lxyz.setX(2, nz * boxSize[2]); drTmp.Ev1Pv2(rAB, Lxyz); double rAB2 = drTmp.squared(); if (rAB2 > rCutSquared && (!intraMolecular || (nx * nx + ny * ny + nz * nz > 0))) continue; double rABMagnitude = Math.sqrt(rAB2); double rAB3 = rABMagnitude * rAB2; double rAB4 = rAB2 * rAB2; double rAB5 = rABMagnitude * rAB4; double erfc = Erf.erfc(alpha * rABMagnitude); double exp_Alpha2r2 = Math.exp(-alpha2 * rAB2); double B0 = (6 * alpha / rAB4 + 4 * alpha3 / rAB2) / sqrtPI * exp_Alpha2r2; tempTensorkk.Ev1v2(drTmp, drTmp); if (intraMolecular && nx * nx + ny * ny + nz * nz == 0) { double A = (1.0 - erfc) / rAB3 - 2 * alpha / sqrtPI * exp_Alpha2r2 / rAB2; double B = B0 - 3 * (1 - erfc) / rAB5; tempTensorkk.TE(-B); tempTensorkk.PEa1Tt1(-A, identity);//ch } else { double A = erfc / rAB3 + 2 * alpha / sqrtPI * exp_Alpha2r2 / rAB2; double B = B0 + 3 * erfc / rAB5; tempTensorkk.TE(-B);//ch tempTensorkk.PEa1Tt1(A, identity); } secondDerivative.PE(tempTensorkk); } //nz } //ny } //nx secondDerivative.TE(q0 * q1); return secondDerivative; }
From source file:com.intel.diceros.test.securerandom.DRNGTest.java
/** * Maurer's "Universal Statistical" Test * <p>// w w w . ja v a 2 s. c o m * The focus of this test is the number of bits between matching patterns * (a measure that is related to the length of a compressed sequence). The * purpose of the test is to detect whether or not the sequence can be * significantly compressed without loss of information. A significantly * compressible sequence is considered to be non-random. */ private void testUniversal(final int[] epsilon) { int i, j, p, L, Q, K, n = epsilon.length; double arg, sqrt2, sigma, phi, sum, p_value, c; int decRep; long[] T; double[] expected_value = { 0, 0, 0, 0, 0, 0, 5.2177052, 6.1962507, 7.1836656, 8.1764248, 9.1723243, 10.170032, 11.168765, 12.168070, 13.167693, 14.167488, 15.167379 }; double[] variance = { 0, 0, 0, 0, 0, 0, 2.954, 3.125, 3.238, 3.311, 3.356, 3.384, 3.401, 3.410, 3.416, 3.419, 3.421 }; L = 5; if (n >= 387840) L = 6; if (n >= 904960) L = 7; if (n >= 2068480) L = 8; if (n >= 4654080) L = 9; if (n >= 10342400) L = 10; if (n >= 22753280) L = 11; if (n >= 49643520) L = 12; if (n >= 107560960) L = 13; if (n >= 231669760) L = 14; if (n >= 496435200) L = 15; if (n >= 1059061760) L = 16; Q = 10 * (int) Math.pow(2, L); K = (int) (Math.floor(n / L) - (double) Q); /* BLOCKS TO TEST */ p = (int) Math.pow(2, L); T = new long[p]; assertTrue("L is out of range.", L >= 6 && L <= 16); assertTrue("Q is less than " + (10 * Math.pow(2, L)), ((double) Q >= 10 * Math.pow(2, L))); c = 0.7 - 0.8 / (double) L + (4 + 32 / (double) L) * Math.pow(K, -3 / (double) L) / 15; sigma = c * Math.sqrt(variance[L] / (double) K); sqrt2 = Math.sqrt(2); sum = 0.0; for (i = 0; i < p; i++) T[i] = 0; for (i = 1; i <= Q; i++) { /* INITIALIZE TABLE */ decRep = 0; for (j = 0; j < L; j++) decRep += epsilon[(i - 1) * L + j] * (long) Math.pow(2, L - 1 - j); T[decRep] = i; } for (i = Q + 1; i <= Q + K; i++) { /* PROCESS BLOCKS */ decRep = 0; for (j = 0; j < L; j++) decRep += epsilon[(i - 1) * L + j] * (long) Math.pow(2, L - 1 - j); sum += Math.log(i - T[decRep]) / Math.log(2); T[decRep] = i; } phi = sum / (double) K; arg = Math.abs(phi - expected_value[L]) / (sqrt2 * sigma); p_value = Erf.erfc(arg); assertTrue("RNG test failed, test universal.", p_value >= 0.01); }
From source file:com.intel.diceros.test.securerandom.DRNGTest.java
/** * Random Excursions Variant Test/*w ww .ja v a2s .co m*/ * <p> * The focus of this test is the total number of times that a particular * state is visited (i.e., occurs) in a cumulative sum random walk. The * purpose of this test is to detect deviations from the expected number * of visits to various states in the random walk. This test is actually * a series of eighteen tests (and conclusions), one test and conclusion * for each of the states: -9, -8, , -1 and +1, +2, , +9. */ public void testRandomExcursionsVariant(final int[] epsilon) { final int n = epsilon.length; int[] stateX = { -9, -8, -7, -6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; int[] S_k = new int[n]; int J = 0, i; S_k[0] = 2 * epsilon[0] - 1; for (i = 1; i < n; i++) { S_k[i] = S_k[i - 1] + 2 * epsilon[i] - 1; if (0 == S_k[i]) { J++; } } if (0 != S_k[n - 1]) { J++; } // int constraint = (int) Math.max(0.005 * Math.pow(n, 0.5), 500); int p, x, count; double p_value; for (p = 0; p < 18; p++) { x = stateX[p]; count = 0; for (i = 0; i < n; i++) { if (S_k[i] == x) { count++; } } p_value = Erf.erfc(Math.abs(count - J) / (Math.sqrt(2.0 * J * (4.0 * Math.abs(x) - 2)))); assertTrue("RNG test failed, test random excursions variant.", p_value >= 0.01); } }
From source file:reflex.module.ReflexErf.java
public ReflexValue erfc(List<ReflexValue> params) { if (params.size() != 1) { throw new ReflexException(-1, "erfc needs one number parameter"); }//from w ww .j a v a2s .com if (!params.get(0).isNumber()) { throw new ReflexException(-1, "erfc needs one number parameter"); } double value = params.get(0).asDouble(); return new ReflexValue(Erf.erfc(value)); }