Skip to content
Snippets Groups Projects
Commit 54304858 authored by Koen Wermer's avatar Koen Wermer
Browse files

added real world test files

parent d7c3d8af
No related branches found
No related tags found
No related merge requests found
module Settings where
testFile, postCondVoid, postCondRefType, postCondPrimType :: String
testFile = "arrays1"
testFile = "BaseSecantSolver"
-- The post condition may depend on the type of the method we are looking at
postCondVoid = "true"
......
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.solvers;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.MathInternalError;
/**
* Base class for all bracketing <em>Secant</em>-based methods for root-finding
* (approximating a zero of a univariate real function).
*
* <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and
* {@link IllinoisSolver <em>Illinois</em>} methods is based on the
* following article: M. Dowell and P. Jarratt,
* <em>A modified regula falsi method for computing the root of an
* equation</em>, BIT Numerical Mathematics, volume 11, number 2,
* pages 168-174, Springer, 1971.</p>
*
* <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is
* based on the following article: M. Dowell and P. Jarratt,
* <em>The "Pegasus" method for computing the root of an equation</em>,
* BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer,
* 1972.</p>
*
* <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a
* bracketing method, so it is not implemented here. It has a separate
* implementation.</p>
*
* @since 3.0
*/
public abstract class BaseSecantSolver
extends AbstractUnivariateSolver
implements BracketedUnivariateSolver<UnivariateFunction> {
/** Default absolute accuracy. */
protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/** The kinds of solutions that the algorithm may accept. */
private AllowedSolution allowed;
/** The <em>Secant</em>-based root-finding method to use. */
private final Method method;
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
* @param method <em>Secant</em>-based root-finding method to use.
*/
protected BaseSecantSolver(final double absoluteAccuracy, final Method method1) {
super(absoluteAccuracy);
this.allowed = AllowedSolution.ANY_SIDE;
this.method = method1;
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
* @param method <em>Secant</em>-based root-finding method to use.
*/
protected BaseSecantSolver(final double relativeAccuracy,
final double absoluteAccuracy1,
final Method method2) {
super(relativeAccuracy, absoluteAccuracy1);
this.allowed = AllowedSolution.ANY_SIDE;
this.method = method2;
}
/**
* Construct a solver.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
* @param functionValueAccuracy Maximum function value error.
* @param method <em>Secant</em>-based root-finding method to use
*/
protected BaseSecantSolver(final double relativeAccuracy2,
final double absoluteAccuracy2,
final double functionValueAccuracy,
final Method method3) {
super(relativeAccuracy2, absoluteAccuracy2, functionValueAccuracy);
this.allowed = AllowedSolution.ANY_SIDE;
this.method = method3;
}
protected final double doSolve()
throws ConvergenceException {
// Get initial solution
double x0 = getMin();
double x1 = getMax();
double f0 = computeObjectiveValue(x0);
double f1 = computeObjectiveValue(x1);
// If one of the bounds is the exact root, return it. Since these are
// not under-approximations or over-approximations, we can return them
// regardless of the allowed solutions.
if (f0 == 0.0) {
return x0;
}
if (f1 == 0.0) {
return x1;
}
// Verify bracketing of initial solution.
verifyBracketing(x0, x1);
// Get accuracies.
final double ftol = getFunctionValueAccuracy();
final double atol = getAbsoluteAccuracy();
final double rtol = getRelativeAccuracy();
// Keep track of inverted intervals, meaning that the left bound is
// larger than the right bound.
boolean inverted = false;
// Keep finding better approximations.
while (true) {
// Calculate the next approximation.
final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
final double fx = computeObjectiveValue(x);
// If the new approximation is the exact root, return it. Since
// this is not an under-approximation or an over-approximation,
// we can return it regardless of the allowed solutions.
if (fx == 0.0) {
return x;
}
// Update the bounds with the new approximation.
if (f1 * fx < 0) {
// The value of x1 has switched to the other bound, thus inverting
// the interval.
x0 = x1;
f0 = f1;
inverted = !inverted;
} else {
switch (this.method) {
case ILLINOIS:
f0 *= 0.5;
break;
case PEGASUS:
f0 *= f1 / (f1 + fx);
break;
case REGULA_FALSI:
// Detect early that algorithm is stuck, instead of waiting
// for the maximum number of iterations to be exceeded.
if (x == x1) {
throw new ConvergenceException();
}
break;
default:
// Should never happen.
throw new MathInternalError();
}
}
// Update from [x0, x1] to [x0, x].
x1 = x;
f1 = fx;
// If the function value of the last approximation is too small,
// given the function value accuracy, then we can't get closer to
// the root than we already are.
if (FastMath.abs(f1) <= ftol) {
switch (this.allowed) {
case ANY_SIDE:
return x1;
case LEFT_SIDE:
if (inverted) {
return x1;
}
break;
case RIGHT_SIDE:
if (!inverted) {
return x1;
}
break;
case BELOW_SIDE:
if (f1 <= 0) {
return x1;
}
break;
case ABOVE_SIDE:
if (f1 >= 0) {
return x1;
}
break;
default:
throw new MathInternalError();
}
}
// If the current interval is within the given accuracies, we
// are satisfied with the current approximation.
if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1),
atol)) {
switch (this.allowed) {
case ANY_SIDE:
return x1;
case LEFT_SIDE:
return inverted ? x1 : x0;
case RIGHT_SIDE:
return inverted ? x0 : x1;
case BELOW_SIDE:
return (f1 <= 0) ? x1 : x0;
case ABOVE_SIDE:
return (f1 >= 0) ? x1 : x0;
default:
throw new MathInternalError();
}
}
}
}
protected enum Method {
REGULA_FALSI,
ILLINOIS,
PEGASUS;
}
//
// UnivariateSolverUtils
//
public static double midpoint(double a, double b) {
return (a + b) * 0.5;
}
public static boolean isBracketing(UnivariateFunction function1,
final double lower1,
final double upper1)
throws NullArgumentException {
if (function1 == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
final double fLo = function1.value(lower1);
final double fHi = function1.value(upper1);
return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
}
public static boolean isSequence(final double start,
final double mid,
final double end) {
return (start < mid) && (mid < end);
}
public static void verifyInterval(final double lower2,
final double upper2)
throws NumberIsTooLargeException {
if (lower2 >= upper2) {
throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
lower2, upper2, false);
}
}
public static void verifySequence(final double lower3,
final double initial,
final double upper3)
throws NumberIsTooLargeException {
verifyInterval(lower3, initial);
verifyInterval(initial, upper3);
}
public static void verifyBracketing(UnivariateFunction function2,
final double lower4,
final double upper4)
throws NullArgumentException,
NoBracketingException {
if (function2 == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
verifyInterval(lower4, upper4);
if (!isBracketing(function2, lower4, upper4)) {
throw new NoBracketingException(lower4, upper4,
function2.value(lower4),
function2.value(upper4));
}
}
}
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
/** Class representing the gradient of a multivariate function.
* <p>
* The vectorial components of the function represent the derivatives
* with respect to each function parameters.
* </p>
* @since 3.1
*/
public class GradientFunction implements MultivariateVectorFunction {
/** Underlying real-valued function. */
private final MultivariateDifferentiableFunction f;
/** Simple constructor.
* @param f underlying real-valued function
*/
public GradientFunction(final MultivariateDifferentiableFunction f) {
this.f = f;
}
/** {@inheritDoc} */
public double[] value(double[] point) {
// set up parameters
final DerivativeStructure[] dsX = new DerivativeStructure[point.length];
for (int i = 0; i < point.length; ++i) {
dsX[i] = new DerivativeStructure(point.length, 1, i, point[i]);
}
// compute the derivatives
final DerivativeStructure dsY = f.value(dsX);
// extract the gradient
final double[] y = new double[point.length];
final int[] orders = new int[point.length];
for (int i = 0; i < point.length; ++i) {
orders[i] = 1;
y[i] = dsY.getPartialDerivative(orders);
orders[i] = 0;
}
return y;
}
}
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.util;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.Serializable;
import java.util.ConcurrentModificationException;
import java.util.NoSuchElementException;
/**
* Open addressed map from int to double.
* <p>This class provides a dedicated map from integers to doubles with a
* much smaller memory overhead than standard <code>java.util.Map</code>.</p>
* <p>This class is not synchronized. The specialized iterators returned by
* {@link #iterator()} are fail-fast: they throw a
* <code>ConcurrentModificationException</code> when they detect the map has been
* modified during iteration.</p>
* @since 2.0
*/
public class OpenIntToDoubleHashMap implements Serializable {
/** Status indicator for free table entries. */
protected static final byte FREE = 0;
/** Status indicator for full table entries. */
protected static final byte FULL = 1;
/** Status indicator for removed table entries. */
protected static final byte REMOVED = 2;
/** Serializable version identifier */
private static final long serialVersionUID = -3646337053166149105L;
/** Load factor for the map. */
private static final float LOAD_FACTOR = 0.5f;
/** Default starting size.
* <p>This must be a power of two for bit mask to work properly. </p>
*/
private static final int DEFAULT_EXPECTED_SIZE = 16;
/** Multiplier for size growth when map fills up.
* <p>This must be a power of two for bit mask to work properly. </p>
*/
private static final int RESIZE_MULTIPLIER = 2;
/** Number of bits to perturb the index when probing for collision resolution. */
private static final int PERTURB_SHIFT = 5;
/** Keys table. */
private int[] keys;
/** Values table. */
private double[] values;
/** States table. */
private byte[] states;
/** Return value for missing entries. */
private final double missingEntries;
/** Current size of the map. */
private int size;
/** Bit mask for hash values. */
private int mask;
/** Modifications count. */
private transient int count;
/**
* Build an empty map with default size and using NaN for missing entries.
*/
public OpenIntToDoubleHashMap() {
this(DEFAULT_EXPECTED_SIZE, Double.NaN);
}
/**
* Build an empty map with default size
* @param missingEntries value to return when a missing entry is fetched
*/
public OpenIntToDoubleHashMap(final double missingEntries) {
this(DEFAULT_EXPECTED_SIZE, missingEntries);
}
/**
* Build an empty map with specified size and using NaN for missing entries.
* @param expectedSize expected number of elements in the map
*/
public OpenIntToDoubleHashMap(final int expectedSize) {
this(expectedSize, Double.NaN);
}
/**
* Build an empty map with specified size.
* @param expectedSize expected number of elements in the map
* @param missingEntries value to return when a missing entry is fetched
*/
public OpenIntToDoubleHashMap(final int expectedSize,
final double missingEntries) {
final int capacity = computeCapacity(expectedSize);
keys = new int[capacity];
values = new double[capacity];
states = new byte[capacity];
this.missingEntries = missingEntries;
mask = capacity - 1;
}
/**
* Copy constructor.
* @param source map to copy
*/
public OpenIntToDoubleHashMap(final OpenIntToDoubleHashMap source) {
final int length = source.keys.length;
keys = new int[length];
System.arraycopy(source.keys, 0, keys, 0, length);
values = new double[length];
System.arraycopy(source.values, 0, values, 0, length);
states = new byte[length];
System.arraycopy(source.states, 0, states, 0, length);
missingEntries = source.missingEntries;
size = source.size;
mask = source.mask;
count = source.count;
}
/**
* Compute the capacity needed for a given size.
* @param expectedSize expected size of the map
* @return capacity to use for the specified size
*/
private static int computeCapacity(final int expectedSize) {
if (expectedSize == 0) {
return 1;
}
final int capacity = (int) FastMath.ceil(expectedSize / LOAD_FACTOR);
final int powerOfTwo = Integer.highestOneBit(capacity);
if (powerOfTwo == capacity) {
return capacity;
}
return nextPowerOfTwo(capacity);
}
/**
* Find the smallest power of two greater than the input value
* @param i input value
* @return smallest power of two greater than the input value
*/
private static int nextPowerOfTwo(final int i) {
return Integer.highestOneBit(i) << 1;
}
/**
* Get the stored value associated with the given key
* @param key key associated with the data
* @return data associated with the key
*/
public double get(final int key) {
final int hash = hashOf(key);
int index = hash & mask;
if (containsKey(key, index)) {
return values[index];
}
if (states[index] == FREE) {
return missingEntries;
}
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
return values[index];
}
}
return missingEntries;
}
/**
* Check if a value is associated with a key.
* @param key key to check
* @return true if a value is associated with key
*/
public boolean containsKey(final int key) {
final int hash = hashOf(key);
int index = hash & mask;
if (containsKey(key, index)) {
return true;
}
if (states[index] == FREE) {
return false;
}
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
return true;
}
}
return false;
}
/**
* Get an iterator over map elements.
* <p>The specialized iterators returned are fail-fast: they throw a
* <code>ConcurrentModificationException</code> when they detect the map
* has been modified during iteration.</p>
* @return iterator over the map elements
*/
public Iterator iterator() {
return new Iterator();
}
/**
* Perturb the hash for starting probing.
* @param hash initial hash
* @return perturbed hash
*/
private static int perturb(final int hash) {
return hash & 0x7fffffff;
}
/**
* Find the index at which a key should be inserted
* @param key key to lookup
* @return index at which key should be inserted
*/
private int findInsertionIndex(final int key) {
return findInsertionIndex(keys, states, key, mask);
}
/**
* Find the index at which a key should be inserted
* @param keys keys table
* @param states states table
* @param key key to lookup
* @param mask bit mask for hash values
* @return index at which key should be inserted
*/
private static int findInsertionIndex(final int[] keys, final byte[] states,
final int key, final int mask) {
final int hash = hashOf(key);
int index = hash & mask;
if (states[index] == FREE) {
return index;
} else if (states[index] == FULL && keys[index] == key) {
return changeIndexSign(index);
}
int perturb = perturb(hash);
int j = index;
if (states[index] == FULL) {
while (true) {
j = probe(perturb, j);
index = j & mask;
perturb >>= PERTURB_SHIFT;
if (states[index] != FULL || keys[index] == key) {
break;
}
}
}
if (states[index] == FREE) {
return index;
} else if (states[index] == FULL) {
// due to the loop exit condition,
// if (states[index] == FULL) then keys[index] == key
return changeIndexSign(index);
}
final int firstRemoved = index;
while (true) {
j = probe(perturb, j);
index = j & mask;
if (states[index] == FREE) {
return firstRemoved;
} else if (states[index] == FULL && keys[index] == key) {
return changeIndexSign(index);
}
perturb >>= PERTURB_SHIFT;
}
}
/**
* Compute next probe for collision resolution
* @param perturb perturbed hash
* @param j previous probe
* @return next probe
*/
private static int probe(final int perturb, final int j) {
return (j << 2) + j + perturb + 1;
}
/**
* Change the index sign
* @param index initial index
* @return changed index
*/
private static int changeIndexSign(final int index) {
return -index - 1;
}
/**
* Get the number of elements stored in the map.
* @return number of elements stored in the map
*/
public int size() {
return size;
}
/**
* Remove the value associated with a key.
* @param key key to which the value is associated
* @return removed value
*/
public double remove(final int key) {
final int hash = hashOf(key);
int index = hash & mask;
if (containsKey(key, index)) {
return doRemove(index);
}
if (states[index] == FREE) {
return missingEntries;
}
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
return doRemove(index);
}
}
return missingEntries;
}
/**
* Check if the tables contain an element associated with specified key
* at specified index.
* @param key key to check
* @param index index to check
* @return true if an element is associated with key at index
*/
private boolean containsKey(final int key, final int index) {
return (key != 0 || states[index] == FULL) && keys[index] == key;
}
/**
* Remove an element at specified index.
* @param index index of the element to remove
* @return removed value
*/
private double doRemove(int index) {
keys[index] = 0;
states[index] = REMOVED;
final double previous = values[index];
values[index] = missingEntries;
--size;
++count;
return previous;
}
/**
* Put a value associated with a key in the map.
* @param key key to which value is associated
* @param value value to put in the map
* @return previous value associated with the key
*/
public double put(final int key, final double value) {
int index = findInsertionIndex(key);
double previous = missingEntries;
boolean newMapping = true;
if (index < 0) {
index = changeIndexSign(index);
previous = values[index];
newMapping = false;
}
keys[index] = key;
states[index] = FULL;
values[index] = value;
if (newMapping) {
++size;
if (shouldGrowTable()) {
growTable();
}
++count;
}
return previous;
}
/**
* Grow the tables.
*/
private void growTable() {
final int oldLength = states.length;
final int[] oldKeys = keys;
final double[] oldValues = values;
final byte[] oldStates = states;
final int newLength = RESIZE_MULTIPLIER * oldLength;
final int[] newKeys = new int[newLength];
final double[] newValues = new double[newLength];
final byte[] newStates = new byte[newLength];
final int newMask = newLength - 1;
for (int i = 0; i < oldLength; ++i) {
if (oldStates[i] == FULL) {
final int key = oldKeys[i];
final int index = findInsertionIndex(newKeys, newStates, key, newMask);
newKeys[index] = key;
newValues[index] = oldValues[i];
newStates[index] = FULL;
}
}
mask = newMask;
keys = newKeys;
values = newValues;
states = newStates;
}
/**
* Check if tables should grow due to increased size.
* @return true if tables should grow
*/
private boolean shouldGrowTable() {
return size > (mask + 1) * LOAD_FACTOR;
}
/**
* Compute the hash value of a key
* @param key key to hash
* @return hash value of the key
*/
private static int hashOf(final int key) {
final int h = key ^ ((key >>> 20) ^ (key >>> 12));
return h ^ (h >>> 7) ^ (h >>> 4);
}
/** Iterator class for the map. */
public class Iterator {
/** Reference modification count. */
private final int referenceCount;
/** Index of current element. */
private int current;
/** Index of next element. */
private int next;
/**
* Simple constructor.
*/
private Iterator() {
// preserve the modification count of the map to detect concurrent modifications later
referenceCount = count;
// initialize current index
next = -1;
try {
advance();
} catch (NoSuchElementException nsee) { // NOPMD
// ignored
}
}
/**
* Check if there is a next element in the map.
* @return true if there is a next element
*/
public boolean hasNext() {
return next >= 0;
}
/**
* Get the key of current entry.
* @return key of current entry
* @exception ConcurrentModificationException if the map is modified during iteration
* @exception NoSuchElementException if there is no element left in the map
*/
public int key()
throws ConcurrentModificationException, NoSuchElementException {
if (referenceCount != count) {
throw new ConcurrentModificationException();
}
if (current < 0) {
throw new NoSuchElementException();
}
return keys[current];
}
/**
* Get the value of current entry.
* @return value of current entry
* @exception ConcurrentModificationException if the map is modified during iteration
* @exception NoSuchElementException if there is no element left in the map
*/
public double value()
throws ConcurrentModificationException, NoSuchElementException {
if (referenceCount != count) {
throw new ConcurrentModificationException();
}
if (current < 0) {
throw new NoSuchElementException();
}
return values[current];
}
/**
* Advance iterator one step further.
* @exception ConcurrentModificationException if the map is modified during iteration
* @exception NoSuchElementException if there is no element left in the map
*/
public void advance()
throws ConcurrentModificationException, NoSuchElementException {
if (referenceCount != count) {
throw new ConcurrentModificationException();
}
// advance on step
current = next;
// prepare next step
try {
while (states[++next] != FULL) { // NOPMD
// nothing to do
}
} catch (ArrayIndexOutOfBoundsException e) {
next = -2;
if (current < 0) {
throw new NoSuchElementException();
}
}
}
}
/**
* Read a serialized object.
* @param stream input stream
* @throws IOException if object cannot be read
* @throws ClassNotFoundException if the class corresponding
* to the serialized object cannot be found
*/
private void readObject(final ObjectInputStream stream)
throws IOException, ClassNotFoundException {
stream.defaultReadObject();
count = 0;
}
}
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment