You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by ra...@apache.org on 2018/06/27 14:51:30 UTC
[02/51] [partial] mahout git commit: MAHOUT-2042 and MAHOUT-2045
Delete directories which were moved/no longer in use
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/random/Multinomial.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/random/Multinomial.java b/math/src/main/java/org/apache/mahout/math/random/Multinomial.java
deleted file mode 100644
index d79c32c..0000000
--- a/math/src/main/java/org/apache/mahout/math/random/Multinomial.java
+++ /dev/null
@@ -1,202 +0,0 @@
-/*
- * 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.mahout.math.random;
-
-import java.util.Iterator;
-import java.util.List;
-import java.util.Map;
-import java.util.Random;
-
-import com.google.common.base.Preconditions;
-import com.google.common.collect.AbstractIterator;
-import com.google.common.collect.Iterables;
-import com.google.common.collect.Lists;
-import com.google.common.collect.Maps;
-import com.google.common.collect.Multiset;
-import org.apache.mahout.common.RandomUtils;
-import org.apache.mahout.math.list.DoubleArrayList;
-
-/**
- * Multinomial sampler that allows updates to element probabilities. The basic idea is that sampling is
- * done by using a simple balanced tree. Probabilities are kept in the tree so that we can navigate to
- * any leaf in log N time. Updates are simple because we can just propagate them upwards.
- * <p/>
- * In order to facilitate access by value, we maintain an additional map from value to tree node.
- */
-public final class Multinomial<T> implements Sampler<T>, Iterable<T> {
- // these lists use heap ordering. Thus, the root is at location 1, first level children at 2 and 3, second level
- // at 4, 5 and 6, 7.
- private final DoubleArrayList weight = new DoubleArrayList();
- private final List<T> values = Lists.newArrayList();
- private final Map<T, Integer> items = Maps.newHashMap();
- private Random rand = RandomUtils.getRandom();
-
- public Multinomial() {
- weight.add(0);
- values.add(null);
- }
-
- public Multinomial(Multiset<T> counts) {
- this();
- Preconditions.checkArgument(!counts.isEmpty(), "Need some data to build sampler");
- rand = RandomUtils.getRandom();
- for (T t : counts.elementSet()) {
- add(t, counts.count(t));
- }
- }
-
- public Multinomial(Iterable<WeightedThing<T>> things) {
- this();
- for (WeightedThing<T> thing : things) {
- add(thing.getValue(), thing.getWeight());
- }
- }
-
- public void add(T value, double w) {
- Preconditions.checkNotNull(value);
- Preconditions.checkArgument(!items.containsKey(value));
-
- int n = this.weight.size();
- if (n == 1) {
- weight.add(w);
- values.add(value);
- items.put(value, 1);
- } else {
- // parent comes down
- weight.add(weight.get(n / 2));
- values.add(values.get(n / 2));
- items.put(values.get(n / 2), n);
- n++;
-
- // new item goes in
- items.put(value, n);
- this.weight.add(w);
- values.add(value);
-
- // parents get incremented all the way to the root
- while (n > 1) {
- n /= 2;
- this.weight.set(n, this.weight.get(n) + w);
- }
- }
- }
-
- public double getWeight(T value) {
- if (items.containsKey(value)) {
- return weight.get(items.get(value));
- } else {
- return 0;
- }
- }
-
- public double getProbability(T value) {
- if (items.containsKey(value)) {
- return weight.get(items.get(value)) / weight.get(1);
- } else {
- return 0;
- }
- }
-
- public double getWeight() {
- if (weight.size() > 1) {
- return weight.get(1);
- } else {
- return 0;
- }
- }
-
- public void delete(T value) {
- set(value, 0);
- }
-
- public void set(T value, double newP) {
- Preconditions.checkArgument(items.containsKey(value));
- int n = items.get(value);
- if (newP <= 0) {
- // this makes the iterator not see such an element even though we leave a phantom in the tree
- // Leaving the phantom behind simplifies tree maintenance and testing, but isn't really necessary.
- items.remove(value);
- }
- double oldP = weight.get(n);
- while (n > 0) {
- weight.set(n, weight.get(n) - oldP + newP);
- n /= 2;
- }
- }
-
- @Override
- public T sample() {
- Preconditions.checkArgument(!weight.isEmpty());
- return sample(rand.nextDouble());
- }
-
- public T sample(double u) {
- u *= weight.get(1);
-
- int n = 1;
- while (2 * n < weight.size()) {
- // children are at 2n and 2n+1
- double left = weight.get(2 * n);
- if (u <= left) {
- n = 2 * n;
- } else {
- u -= left;
- n = 2 * n + 1;
- }
- }
- return values.get(n);
- }
-
- /**
- * Exposed for testing only. Returns a list of the leaf weights. These are in an
- * order such that probing just before and after the cumulative sum of these weights
- * will touch every element of the tree twice and thus will make it possible to test
- * every possible left/right decision in navigating the tree.
- */
- List<Double> getWeights() {
- List<Double> r = Lists.newArrayList();
- int i = Integer.highestOneBit(weight.size());
- while (i < weight.size()) {
- r.add(weight.get(i));
- i++;
- }
- i /= 2;
- while (i < Integer.highestOneBit(weight.size())) {
- r.add(weight.get(i));
- i++;
- }
- return r;
- }
-
- @Override
- public Iterator<T> iterator() {
- return new AbstractIterator<T>() {
- Iterator<T> valuesIterator = Iterables.skip(values, 1).iterator();
- @Override
- protected T computeNext() {
- while (valuesIterator.hasNext()) {
- T next = valuesIterator.next();
- if (items.containsKey(next)) {
- return next;
- }
- }
- return endOfData();
- }
- };
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/random/Normal.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/random/Normal.java b/math/src/main/java/org/apache/mahout/math/random/Normal.java
deleted file mode 100644
index c162f26..0000000
--- a/math/src/main/java/org/apache/mahout/math/random/Normal.java
+++ /dev/null
@@ -1,40 +0,0 @@
-/*
- * 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.mahout.math.random;
-
-import org.apache.mahout.common.RandomUtils;
-
-import java.util.Random;
-
-public final class Normal extends AbstractSamplerFunction {
- private final Random rand = RandomUtils.getRandom();
- private double mean = 0;
- private double sd = 1;
-
- public Normal() {}
-
- public Normal(double mean, double sd) {
- this.mean = mean;
- this.sd = sd;
- }
-
- @Override
- public Double sample() {
- return rand.nextGaussian() * sd + mean;
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/random/PoissonSampler.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/random/PoissonSampler.java b/math/src/main/java/org/apache/mahout/math/random/PoissonSampler.java
deleted file mode 100644
index e4e49f8..0000000
--- a/math/src/main/java/org/apache/mahout/math/random/PoissonSampler.java
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * 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.mahout.math.random;
-
-import com.google.common.collect.Lists;
-import org.apache.commons.math3.distribution.PoissonDistribution;
-import org.apache.mahout.common.RandomUtils;
-import org.apache.mahout.common.RandomWrapper;
-
-import java.util.List;
-
-/**
- * Samples from a Poisson distribution. Should probably not be used for lambda > 1000 or so.
- */
-public final class PoissonSampler extends AbstractSamplerFunction {
-
- private double limit;
- private Multinomial<Integer> partial;
- private final RandomWrapper gen;
- private final PoissonDistribution pd;
-
- public PoissonSampler(double lambda) {
- limit = 1;
- gen = RandomUtils.getRandom();
- pd = new PoissonDistribution(gen.getRandomGenerator(),
- lambda,
- PoissonDistribution.DEFAULT_EPSILON,
- PoissonDistribution.DEFAULT_MAX_ITERATIONS);
- }
-
- @Override
- public Double sample() {
- return sample(gen.nextDouble());
- }
-
- double sample(double u) {
- if (u < limit) {
- List<WeightedThing<Integer>> steps = Lists.newArrayList();
- limit = 1;
- int i = 0;
- while (u / 20 < limit) {
- double pdf = pd.probability(i);
- limit -= pdf;
- steps.add(new WeightedThing<>(i, pdf));
- i++;
- }
- steps.add(new WeightedThing<>(steps.size(), limit));
- partial = new Multinomial<>(steps);
- }
- return partial.sample(u);
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/random/Sampler.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/random/Sampler.java b/math/src/main/java/org/apache/mahout/math/random/Sampler.java
deleted file mode 100644
index 51460fa..0000000
--- a/math/src/main/java/org/apache/mahout/math/random/Sampler.java
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * 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.mahout.math.random;
-
-/**
- * Samples from a generic type.
- */
-public interface Sampler<T> {
- T sample();
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/random/WeightedThing.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/random/WeightedThing.java b/math/src/main/java/org/apache/mahout/math/random/WeightedThing.java
deleted file mode 100644
index 20f6df3..0000000
--- a/math/src/main/java/org/apache/mahout/math/random/WeightedThing.java
+++ /dev/null
@@ -1,71 +0,0 @@
-/*
- * 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.mahout.math.random;
-
-import com.google.common.base.Preconditions;
-import org.apache.mahout.common.RandomUtils;
-
-/**
- * Handy for creating multinomial distributions of things.
- */
-public final class WeightedThing<T> implements Comparable<WeightedThing<T>> {
- private double weight;
- private final T value;
-
- public WeightedThing(T thing, double weight) {
- this.value = Preconditions.checkNotNull(thing);
- this.weight = weight;
- }
-
- public WeightedThing(double weight) {
- this.value = null;
- this.weight = weight;
- }
-
- public T getValue() {
- return value;
- }
-
- public double getWeight() {
- return weight;
- }
-
- public void setWeight(double weight) {
- this.weight = weight;
- }
-
- @Override
- public int compareTo(WeightedThing<T> other) {
- return Double.compare(this.weight, other.weight);
- }
-
- @Override
- public boolean equals(Object o) {
- if (o instanceof WeightedThing) {
- @SuppressWarnings("unchecked")
- WeightedThing<T> other = (WeightedThing<T>) o;
- return weight == other.weight && value.equals(other.value);
- }
- return false;
- }
-
- @Override
- public int hashCode() {
- return 31 * RandomUtils.hashDouble(weight) + value.hashCode();
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/set/AbstractSet.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/set/AbstractSet.java b/math/src/main/java/org/apache/mahout/math/set/AbstractSet.java
deleted file mode 100644
index 7691420..0000000
--- a/math/src/main/java/org/apache/mahout/math/set/AbstractSet.java
+++ /dev/null
@@ -1,188 +0,0 @@
-/**
- * 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.
- */
-/*
-Copyright 1999 CERN - European Organization for Nuclear Research.
-Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
-is hereby granted without fee, provided that the above copyright notice appear in all copies and
-that both that copyright notice and this permission notice appear in supporting documentation.
-CERN makes no representations about the suitability of this software for any purpose.
-It is provided "as is" without expressed or implied warranty.
-*/
-package org.apache.mahout.math.set;
-
-import org.apache.mahout.math.PersistentObject;
-import org.apache.mahout.math.map.PrimeFinder;
-
-public abstract class AbstractSet extends PersistentObject {
- //public static boolean debug = false; // debug only
-
- /** The number of distinct associations in the map; its "size()". */
- protected int distinct;
-
- /**
- * The table capacity c=table.length always satisfies the invariant <tt>c * minLoadFactor <= s <= c *
- * maxLoadFactor</tt>, where s=size() is the number of associations currently contained. The term "c * minLoadFactor"
- * is called the "lowWaterMark", "c * maxLoadFactor" is called the "highWaterMark". In other words, the table capacity
- * (and proportionally the memory used by this class) oscillates within these constraints. The terms are precomputed
- * and cached to avoid recalculating them each time put(..) or removeKey(...) is called.
- */
- protected int lowWaterMark;
- protected int highWaterMark;
-
- /** The minimum load factor for the hashtable. */
- protected double minLoadFactor;
-
- /** The maximum load factor for the hashtable. */
- protected double maxLoadFactor;
-
- // these are public access for unit tests.
- public static final int DEFAULT_CAPACITY = 277;
- public static final double DEFAULT_MIN_LOAD_FACTOR = 0.2;
- public static final double DEFAULT_MAX_LOAD_FACTOR = 0.5;
-
- /**
- * Chooses a new prime table capacity optimized for growing that (approximately) satisfies the invariant <tt>c *
- * minLoadFactor <= size <= c * maxLoadFactor</tt> and has at least one FREE slot for the given size.
- */
- protected int chooseGrowCapacity(int size, double minLoad, double maxLoad) {
- return nextPrime(Math.max(size + 1, (int) ((4 * size / (3 * minLoad + maxLoad)))));
- }
-
- /**
- * Returns new high water mark threshold based on current capacity and maxLoadFactor.
- *
- * @return int the new threshold.
- */
- protected int chooseHighWaterMark(int capacity, double maxLoad) {
- return Math.min(capacity - 2, (int) (capacity * maxLoad)); //makes sure there is always at least one FREE slot
- }
-
- /**
- * Returns new low water mark threshold based on current capacity and minLoadFactor.
- *
- * @return int the new threshold.
- */
- protected int chooseLowWaterMark(int capacity, double minLoad) {
- return (int) (capacity * minLoad);
- }
-
- /**
- * Chooses a new prime table capacity neither favoring shrinking nor growing, that (approximately) satisfies the
- * invariant <tt>c * minLoadFactor <= size <= c * maxLoadFactor</tt> and has at least one FREE slot for the given
- * size.
- */
- protected int chooseMeanCapacity(int size, double minLoad, double maxLoad) {
- return nextPrime(Math.max(size + 1, (int) ((2 * size / (minLoad + maxLoad)))));
- }
-
- /**
- * Chooses a new prime table capacity optimized for shrinking that (approximately) satisfies the invariant <tt>c *
- * minLoadFactor <= size <= c * maxLoadFactor</tt> and has at least one FREE slot for the given size.
- */
- protected int chooseShrinkCapacity(int size, double minLoad, double maxLoad) {
- return nextPrime(Math.max(size + 1, (int) ((4 * size / (minLoad + 3 * maxLoad)))));
- }
-
- /** Removes all (key,value) associations from the receiver. */
- public abstract void clear();
-
- /**
- * Ensures that the receiver can hold at least the specified number of elements without needing to allocate new
- * internal memory. If necessary, allocates new internal memory and increases the capacity of the receiver. <p> This
- * method never need be called; it is for performance tuning only. Calling this method before <tt>put()</tt>ing a
- * large number of associations boosts performance, because the receiver will grow only once instead of potentially
- * many times. <p> <b>This default implementation does nothing.</b> Override this method if necessary.
- *
- * @param minCapacity the desired minimum capacity.
- */
- public void ensureCapacity(int minCapacity) {
- }
-
- /**
- * Returns <tt>true</tt> if the receiver contains no (key,value) associations.
- *
- * @return <tt>true</tt> if the receiver contains no (key,value) associations.
- */
- public boolean isEmpty() {
- return distinct == 0;
- }
-
- /**
- * Returns a prime number which is <code>>= desiredCapacity</code> and very close to <code>desiredCapacity</code>
- * (within 11% if <code>desiredCapacity >= 1000</code>).
- *
- * @param desiredCapacity the capacity desired by the user.
- * @return the capacity which should be used for a hashtable.
- */
- protected int nextPrime(int desiredCapacity) {
- return PrimeFinder.nextPrime(desiredCapacity);
- }
-
- /**
- * Initializes the receiver. You will almost certainly need to override this method in subclasses to initialize the
- * hash table.
- *
- * @param initialCapacity the initial capacity of the receiver.
- * @param minLoadFactor the minLoadFactor of the receiver.
- * @param maxLoadFactor the maxLoadFactor of the receiver.
- * @throws IllegalArgumentException if <tt>initialCapacity < 0 || (minLoadFactor < 0.0 || minLoadFactor >= 1.0) ||
- * (maxLoadFactor <= 0.0 || maxLoadFactor >= 1.0) || (minLoadFactor >=
- * maxLoadFactor)</tt>.
- */
- protected void setUp(int initialCapacity, double minLoadFactor, double maxLoadFactor) {
- if (initialCapacity < 0) {
- throw new IllegalArgumentException("Initial Capacity must not be less than zero: " + initialCapacity);
- }
- if (minLoadFactor < 0.0 || minLoadFactor >= 1.0) {
- throw new IllegalArgumentException("Illegal minLoadFactor: " + minLoadFactor);
- }
- if (maxLoadFactor <= 0.0 || maxLoadFactor >= 1.0) {
- throw new IllegalArgumentException("Illegal maxLoadFactor: " + maxLoadFactor);
- }
- if (minLoadFactor >= maxLoadFactor) {
- throw new IllegalArgumentException(
- "Illegal minLoadFactor: " + minLoadFactor + " and maxLoadFactor: " + maxLoadFactor);
- }
- }
-
- /**
- * Returns the number of (key,value) associations currently contained.
- *
- * @return the number of (key,value) associations currently contained.
- */
- public int size() {
- return distinct;
- }
-
- /**
- * Trims the capacity of the receiver to be the receiver's current size. Releases any superfluous internal memory. An
- * application can use this operation to minimize the storage of the receiver. <p> This default implementation does
- * nothing. Override this method if necessary.
- */
- public void trimToSize() {
- }
-
- protected static boolean equalsMindTheNull(Object a, Object b) {
- if (a == null && b == null) {
- return true;
- }
- if (a == null || b == null) {
- return false;
- }
- return a.equals(b);
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/set/HashUtils.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/set/HashUtils.java b/math/src/main/java/org/apache/mahout/math/set/HashUtils.java
deleted file mode 100644
index f5dfeb0..0000000
--- a/math/src/main/java/org/apache/mahout/math/set/HashUtils.java
+++ /dev/null
@@ -1,56 +0,0 @@
-/*
- * 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.mahout.math.set;
-
-/**
- * Computes hashes of primitive values. Providing these as statics allows the templated code
- * to compute hashes of sets.
- */
-public final class HashUtils {
-
- private HashUtils() {
- }
-
- public static int hash(byte x) {
- return x;
- }
-
- public static int hash(short x) {
- return x;
- }
-
- public static int hash(char x) {
- return x;
- }
-
- public static int hash(int x) {
- return x;
- }
-
- public static int hash(float x) {
- return Float.floatToIntBits(x) >>> 3 + Float.floatToIntBits((float) (Math.PI * x));
- }
-
- public static int hash(double x) {
- return hash(17 * Double.doubleToLongBits(x));
- }
-
- public static int hash(long x) {
- return (int) ((x * 11) >>> 32 ^ x);
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/set/OpenHashSet.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/set/OpenHashSet.java b/math/src/main/java/org/apache/mahout/math/set/OpenHashSet.java
deleted file mode 100644
index 285b5a5..0000000
--- a/math/src/main/java/org/apache/mahout/math/set/OpenHashSet.java
+++ /dev/null
@@ -1,548 +0,0 @@
-/**
- * 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.mahout.math.set;
-
-import java.nio.ByteBuffer;
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.Collection;
-import java.util.Iterator;
-import java.util.List;
-import java.util.Set;
-
-import org.apache.mahout.math.MurmurHash;
-import org.apache.mahout.math.function.ObjectProcedure;
-import org.apache.mahout.math.map.PrimeFinder;
-
-/**
- * Open hashing alternative to java.util.HashSet.
- **/
-public class OpenHashSet<T> extends AbstractSet implements Set<T> {
- protected static final byte FREE = 0;
- protected static final byte FULL = 1;
- protected static final byte REMOVED = 2;
- protected static final char NO_KEY_VALUE = 0;
-
- /** The hash table keys. */
- private Object[] table;
-
- /** The state of each hash table entry (FREE, FULL, REMOVED). */
- private byte[] state;
-
- /** The number of table entries in state==FREE. */
- private int freeEntries;
-
-
- /** Constructs an empty map with default capacity and default load factors. */
- public OpenHashSet() {
- this(DEFAULT_CAPACITY);
- }
-
- /**
- * Constructs an empty map with the specified initial capacity and default load factors.
- *
- * @param initialCapacity the initial capacity of the map.
- * @throws IllegalArgumentException if the initial capacity is less than zero.
- */
- public OpenHashSet(int initialCapacity) {
- this(initialCapacity, DEFAULT_MIN_LOAD_FACTOR, DEFAULT_MAX_LOAD_FACTOR);
- }
-
- /**
- * Constructs an empty map with the specified initial capacity and the specified minimum and maximum load factor.
- *
- * @param initialCapacity the initial capacity.
- * @param minLoadFactor the minimum load factor.
- * @param maxLoadFactor the maximum load factor.
- * @throws IllegalArgumentException if <tt>initialCapacity < 0 || (minLoadFactor < 0.0 || minLoadFactor >= 1.0) ||
- * (maxLoadFactor <= 0.0 || maxLoadFactor >= 1.0) || (minLoadFactor >=
- * maxLoadFactor)</tt>.
- */
- public OpenHashSet(int initialCapacity, double minLoadFactor, double maxLoadFactor) {
- setUp(initialCapacity, minLoadFactor, maxLoadFactor);
- }
-
- /** Removes all values associations from the receiver. Implicitly calls <tt>trimToSize()</tt>. */
- @Override
- public void clear() {
- Arrays.fill(this.state, 0, state.length - 1, FREE);
- distinct = 0;
- freeEntries = table.length; // delta
- trimToSize();
- }
-
- /**
- * Returns a deep copy of the receiver.
- *
- * @return a deep copy of the receiver.
- */
- @SuppressWarnings("unchecked")
- @Override
- public Object clone() {
- OpenHashSet<T> copy = (OpenHashSet<T>) super.clone();
- copy.table = copy.table.clone();
- copy.state = copy.state.clone();
- return copy;
- }
-
- /**
- * Returns <tt>true</tt> if the receiver contains the specified key.
- *
- * @return <tt>true</tt> if the receiver contains the specified key.
- */
- @Override
- @SuppressWarnings("unchecked")
- public boolean contains(Object key) {
- return indexOfKey((T)key) >= 0;
- }
-
- /**
- * Ensures that the receiver can hold at least the specified number of associations without needing to allocate new
- * internal memory. If necessary, allocates new internal memory and increases the capacity of the receiver. <p> This
- * method never need be called; it is for performance tuning only. Calling this method before <tt>add()</tt>ing a
- * large number of associations boosts performance, because the receiver will grow only once instead of potentially
- * many times and hash collisions get less probable.
- *
- * @param minCapacity the desired minimum capacity.
- */
- @Override
- public void ensureCapacity(int minCapacity) {
- if (table.length < minCapacity) {
- int newCapacity = nextPrime(minCapacity);
- rehash(newCapacity);
- }
- }
-
- /**
- * Applies a procedure to each key of the receiver, if any. Note: Iterates over the keys in no particular order.
- * Subclasses can define a particular order, for example, "sorted by key". All methods which <i>can</i> be expressed
- * in terms of this method (most methods can) <i>must guarantee</i> to use the <i>same</i> order defined by this
- * method, even if it is no particular order. This is necessary so that, for example, methods <tt>keys</tt> and
- * <tt>values</tt> will yield association pairs, not two uncorrelated lists.
- *
- * @param procedure the procedure to be applied. Stops iteration if the procedure returns <tt>false</tt>, otherwise
- * continues.
- * @return <tt>false</tt> if the procedure stopped before all keys where iterated over, <tt>true</tt> otherwise.
- */
- @SuppressWarnings("unchecked")
- public boolean forEachKey(ObjectProcedure<T> procedure) {
- for (int i = table.length; i-- > 0;) {
- if (state[i] == FULL) {
- if (!procedure.apply((T)table[i])) {
- return false;
- }
- }
- }
- return true;
- }
-
- /**
- * @param key the key to be added to the receiver.
- * @return the index where the key would need to be inserted, if it is not already contained. Returns -index-1 if the
- * key is already contained at slot index. Therefore, if the returned index < 0, then it is already contained
- * at slot -index-1. If the returned index >= 0, then it is NOT already contained and should be inserted at
- * slot index.
- */
- protected int indexOfInsertion(T key) {
- Object[] tab = table;
- byte[] stat = state;
- int length = tab.length;
-
- int hash = key.hashCode() & 0x7FFFFFFF;
- int i = hash % length;
- int decrement = hash % (length - 2); // double hashing, see http://www.eece.unm.edu/faculty/heileman/hash/node4.html
- //int decrement = (hash / length) % length;
- if (decrement == 0) {
- decrement = 1;
- }
-
- // stop if we find a removed or free slot, or if we find the key itself
- // do NOT skip over removed slots (yes, open addressing is like that...)
- while (stat[i] == FULL && tab[i] != key) {
- i -= decrement;
- //hashCollisions++;
- if (i < 0) {
- i += length;
- }
- }
-
- if (stat[i] == REMOVED) {
- // stop if we find a free slot, or if we find the key itself.
- // do skip over removed slots (yes, open addressing is like that...)
- // assertion: there is at least one FREE slot.
- int j = i;
- while (stat[i] != FREE && (stat[i] == REMOVED || tab[i] != key)) {
- i -= decrement;
- //hashCollisions++;
- if (i < 0) {
- i += length;
- }
- }
- if (stat[i] == FREE) {
- i = j;
- }
- }
-
-
- if (stat[i] == FULL) {
- // key already contained at slot i.
- // return a negative number identifying the slot.
- return -i - 1;
- }
- // not already contained, should be inserted at slot i.
- // return a number >= 0 identifying the slot.
- return i;
- }
-
- /**
- * @param key the key to be searched in the receiver.
- * @return the index where the key is contained in the receiver, returns -1 if the key was not found.
- */
- protected int indexOfKey(T key) {
- Object[] tab = table;
- byte[] stat = state;
- int length = tab.length;
-
- int hash = key.hashCode() & 0x7FFFFFFF;
- int i = hash % length;
- int decrement = hash % (length - 2); // double hashing, see http://www.eece.unm.edu/faculty/heileman/hash/node4.html
- //int decrement = (hash / length) % length;
- if (decrement == 0) {
- decrement = 1;
- }
-
- // stop if we find a free slot, or if we find the key itself.
- // do skip over removed slots (yes, open addressing is like that...)
- while (stat[i] != FREE && (stat[i] == REMOVED || (!key.equals(tab[i])))) {
- i -= decrement;
- //hashCollisions++;
- if (i < 0) {
- i += length;
- }
- }
-
- if (stat[i] == FREE) {
- return -1;
- } // not found
- return i; //found, return index where key is contained
- }
-
- /**
- * Fills all keys contained in the receiver into the specified list. Fills the list, starting at index 0. After this
- * call returns the specified list has a new size that equals <tt>this.size()</tt>.
- * This method can be used
- * to iterate over the keys of the receiver.
- *
- * @param list the list to be filled, can have any size.
- */
- @SuppressWarnings("unchecked")
- public void keys(List<T> list) {
- list.clear();
-
-
- Object [] tab = table;
- byte[] stat = state;
-
- for (int i = tab.length; i-- > 0;) {
- if (stat[i] == FULL) {
- list.add((T)tab[i]);
- }
- }
- }
-
- @SuppressWarnings("unchecked")
- @Override
- public boolean add(Object key) {
- int i = indexOfInsertion((T)key);
- if (i < 0) { //already contained
- return false;
- }
-
- if (this.distinct > this.highWaterMark) {
- int newCapacity = chooseGrowCapacity(this.distinct + 1, this.minLoadFactor, this.maxLoadFactor);
- rehash(newCapacity);
- return add(key);
- }
-
- this.table[i] = key;
- if (this.state[i] == FREE) {
- this.freeEntries--;
- }
- this.state[i] = FULL;
- this.distinct++;
-
- if (this.freeEntries < 1) { //delta
- int newCapacity = chooseGrowCapacity(this.distinct + 1, this.minLoadFactor, this.maxLoadFactor);
- rehash(newCapacity);
- return add(key);
- }
-
- return true;
- }
-
- /**
- * Rehashes the contents of the receiver into a new table with a smaller or larger capacity. This method is called
- * automatically when the number of keys in the receiver exceeds the high water mark or falls below the low water
- * mark.
- */
- @SuppressWarnings("unchecked")
- protected void rehash(int newCapacity) {
- int oldCapacity = table.length;
- //if (oldCapacity == newCapacity) return;
-
- Object[] oldTable = table;
- byte[] oldState = state;
-
- Object[] newTable = new Object[newCapacity];
- byte[] newState = new byte[newCapacity];
-
- this.lowWaterMark = chooseLowWaterMark(newCapacity, this.minLoadFactor);
- this.highWaterMark = chooseHighWaterMark(newCapacity, this.maxLoadFactor);
-
- this.table = newTable;
- this.state = newState;
- this.freeEntries = newCapacity - this.distinct; // delta
-
- for (int i = oldCapacity; i-- > 0;) {
- if (oldState[i] == FULL) {
- Object element = oldTable[i];
- int index = indexOfInsertion((T)element);
- newTable[index] = element;
- newState[index] = FULL;
- }
- }
- }
-
- /**
- * Removes the given key with its associated element from the receiver, if present.
- *
- * @param key the key to be removed from the receiver.
- * @return <tt>true</tt> if the receiver contained the specified key, <tt>false</tt> otherwise.
- */
- @SuppressWarnings("unchecked")
- @Override
- public boolean remove(Object key) {
- int i = indexOfKey((T)key);
- if (i < 0) {
- return false;
- } // key not contained
-
- this.state[i] = REMOVED;
- this.distinct--;
-
- if (this.distinct < this.lowWaterMark) {
- int newCapacity = chooseShrinkCapacity(this.distinct, this.minLoadFactor, this.maxLoadFactor);
- rehash(newCapacity);
- }
-
- return true;
- }
-
- /**
- * Initializes the receiver.
- *
- * @param initialCapacity the initial capacity of the receiver.
- * @param minLoadFactor the minLoadFactor of the receiver.
- * @param maxLoadFactor the maxLoadFactor of the receiver.
- * @throws IllegalArgumentException if <tt>initialCapacity < 0 || (minLoadFactor < 0.0 || minLoadFactor >= 1.0) ||
- * (maxLoadFactor <= 0.0 || maxLoadFactor >= 1.0) || (minLoadFactor >=
- * maxLoadFactor)</tt>.
- */
- @Override
- protected final void setUp(int initialCapacity, double minLoadFactor, double maxLoadFactor) {
- int capacity = initialCapacity;
- super.setUp(capacity, minLoadFactor, maxLoadFactor);
- capacity = nextPrime(capacity);
- if (capacity == 0) {
- capacity = 1;
- } // open addressing needs at least one FREE slot at any time.
-
- this.table = new Object[capacity];
- this.state = new byte[capacity];
-
- // memory will be exhausted long before this pathological case happens, anyway.
- this.minLoadFactor = minLoadFactor;
- if (capacity == PrimeFinder.LARGEST_PRIME) {
- this.maxLoadFactor = 1.0;
- } else {
- this.maxLoadFactor = maxLoadFactor;
- }
-
- this.distinct = 0;
- this.freeEntries = capacity; // delta
-
- // lowWaterMark will be established upon first expansion.
- // establishing it now (upon instance construction) would immediately make the table shrink upon first put(...).
- // After all the idea of an "initialCapacity" implies violating lowWaterMarks when an object is young.
- // See ensureCapacity(...)
- this.lowWaterMark = 0;
- this.highWaterMark = chooseHighWaterMark(capacity, this.maxLoadFactor);
- }
-
- /**
- * Trims the capacity of the receiver to be the receiver's current size. Releases any superfluous internal memory. An
- * application can use this operation to minimize the storage of the receiver.
- */
- @Override
- public void trimToSize() {
- // * 1.2 because open addressing's performance exponentially degrades beyond that point
- // so that even rehashing the table can take very long
- int newCapacity = nextPrime((int) (1 + 1.2 * size()));
- if (table.length > newCapacity) {
- rehash(newCapacity);
- }
- }
-
- /**
- * Access for unit tests.
- * @param capacity
- * @param minLoadFactor
- * @param maxLoadFactor
- */
- void getInternalFactors(int[] capacity,
- double[] minLoadFactor,
- double[] maxLoadFactor) {
- capacity[0] = table.length;
- minLoadFactor[0] = this.minLoadFactor;
- maxLoadFactor[0] = this.maxLoadFactor;
- }
-
- @Override
- public boolean isEmpty() {
- return size() == 0;
- }
-
- /**
- * OpenHashSet instances are only equal to other OpenHashSet instances, not to
- * any other collection. Hypothetically, we should check for and permit
- * equals on other Sets.
- */
- @Override
- @SuppressWarnings("unchecked")
- public boolean equals(Object obj) {
- if (obj == this) {
- return true;
- }
-
- if (!(obj instanceof OpenHashSet)) {
- return false;
- }
- final OpenHashSet<T> other = (OpenHashSet<T>) obj;
- if (other.size() != size()) {
- return false;
- }
-
- return forEachKey(new ObjectProcedure<T>() {
- @Override
- public boolean apply(T key) {
- return other.contains(key);
- }
- });
- }
-
- @Override
- public int hashCode() {
- ByteBuffer buf = ByteBuffer.allocate(size());
- for (int i = 0; i < table.length; i++) {
- Object v = table[i];
- if (state[i] == FULL) {
- buf.putInt(v.hashCode());
- }
- }
- return MurmurHash.hash(buf, this.getClass().getName().hashCode());
- }
-
- /**
- * Implement the standard Java Collections iterator. Note that 'remove' is silently
- * ineffectual here. This method is provided for convenience, only.
- */
- @Override
- public Iterator<T> iterator() {
- List<T> keyList = new ArrayList<>();
- keys(keyList);
- return keyList.iterator();
- }
-
- @Override
- public Object[] toArray() {
- List<T> keyList = new ArrayList<>();
- keys(keyList);
- return keyList.toArray();
- }
-
- @Override
- public boolean addAll(Collection<? extends T> c) {
- boolean anyAdded = false;
- for (T o : c) {
- boolean added = add(o);
- anyAdded |= added;
- }
- return anyAdded;
- }
-
- @Override
- public boolean containsAll(Collection<?> c) {
- for (Object o : c) {
- if (!contains(o)) {
- return false;
- }
- }
- return true;
- }
-
- @Override
- public boolean removeAll(Collection<?> c) {
- boolean anyRemoved = false;
- for (Object o : c) {
- boolean removed = remove(o);
- anyRemoved |= removed;
- }
- return anyRemoved;
- }
-
- @Override
- public boolean retainAll(Collection<?> c) {
- final Collection<?> finalCollection = c;
- final boolean[] modified = new boolean[1];
- modified[0] = false;
- forEachKey(new ObjectProcedure<T>() {
- @Override
- public boolean apply(T element) {
- if (!finalCollection.contains(element)) {
- remove(element);
- modified[0] = true;
- }
- return true;
- }
- });
- return modified[0];
- }
-
- @Override
- public <T1> T1[] toArray(T1[] a) {
- return keys().toArray(a);
- }
-
- public List<T> keys() {
- List<T> keys = new ArrayList<>();
- keys(keys);
- return keys;
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/solver/ConjugateGradientSolver.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/solver/ConjugateGradientSolver.java b/math/src/main/java/org/apache/mahout/math/solver/ConjugateGradientSolver.java
deleted file mode 100644
index 02bde9b..0000000
--- a/math/src/main/java/org/apache/mahout/math/solver/ConjugateGradientSolver.java
+++ /dev/null
@@ -1,213 +0,0 @@
-/**
- * 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.mahout.math.solver;
-
-import org.apache.mahout.math.CardinalityException;
-import org.apache.mahout.math.DenseVector;
-import org.apache.mahout.math.Vector;
-import org.apache.mahout.math.VectorIterable;
-import org.apache.mahout.math.function.Functions;
-import org.apache.mahout.math.function.PlusMult;
-import org.slf4j.Logger;
-import org.slf4j.LoggerFactory;
-
-/**
- * <p>Implementation of a conjugate gradient iterative solver for linear systems. Implements both
- * standard conjugate gradient and pre-conditioned conjugate gradient.
- *
- * <p>Conjugate gradient requires the matrix A in the linear system Ax = b to be symmetric and positive
- * definite. For convenience, this implementation could be extended relatively easily to handle the
- * case where the input matrix to be be non-symmetric, in which case the system A'Ax = b would be solved.
- * Because this requires only one pass through the matrix A, it is faster than explicitly computing A'A,
- * then passing the results to the solver.
- *
- * <p>For inputs that may be ill conditioned (often the case for highly sparse input), this solver
- * also accepts a parameter, lambda, which adds a scaled identity to the matrix A, solving the system
- * (A + lambda*I)x = b. This obviously changes the solution, but it will guarantee solvability. The
- * ridge regression approach to linear regression is a common use of this feature.
- *
- * <p>If only an approximate solution is required, the maximum number of iterations or the error threshold
- * may be specified to end the algorithm early at the expense of accuracy. When the matrix A is ill conditioned,
- * it may sometimes be necessary to increase the maximum number of iterations above the default of A.numCols()
- * due to numerical issues.
- *
- * <p>By default the solver will run a.numCols() iterations or until the residual falls below 1E-9.
- *
- * <p>For more information on the conjugate gradient algorithm, see Golub & van Loan, "Matrix Computations",
- * sections 10.2 and 10.3 or the <a href="http://en.wikipedia.org/wiki/Conjugate_gradient">conjugate gradient
- * wikipedia article</a>.
- */
-
-public class ConjugateGradientSolver {
-
- public static final double DEFAULT_MAX_ERROR = 1.0e-9;
-
- private static final Logger log = LoggerFactory.getLogger(ConjugateGradientSolver.class);
- private static final PlusMult PLUS_MULT = new PlusMult(1.0);
-
- private int iterations;
- private double residualNormSquared;
-
- public ConjugateGradientSolver() {
- this.iterations = 0;
- this.residualNormSquared = Double.NaN;
- }
-
- /**
- * Solves the system Ax = b with default termination criteria. A must be symmetric, square, and positive definite.
- * Only the squareness of a is checked, since testing for symmetry and positive definiteness are too expensive. If
- * an invalid matrix is specified, then the algorithm may not yield a valid result.
- *
- * @param a The linear operator A.
- * @param b The vector b.
- * @return The result x of solving the system.
- * @throws IllegalArgumentException if a is not square or if the size of b is not equal to the number of columns of a.
- *
- */
- public Vector solve(VectorIterable a, Vector b) {
- return solve(a, b, null, b.size() + 2, DEFAULT_MAX_ERROR);
- }
-
- /**
- * Solves the system Ax = b with default termination criteria using the specified preconditioner. A must be
- * symmetric, square, and positive definite. Only the squareness of a is checked, since testing for symmetry
- * and positive definiteness are too expensive. If an invalid matrix is specified, then the algorithm may not
- * yield a valid result.
- *
- * @param a The linear operator A.
- * @param b The vector b.
- * @param precond A preconditioner to use on A during the solution process.
- * @return The result x of solving the system.
- * @throws IllegalArgumentException if a is not square or if the size of b is not equal to the number of columns of a.
- *
- */
- public Vector solve(VectorIterable a, Vector b, Preconditioner precond) {
- return solve(a, b, precond, b.size() + 2, DEFAULT_MAX_ERROR);
- }
-
-
- /**
- * Solves the system Ax = b, where A is a linear operator and b is a vector. Uses the specified preconditioner
- * to improve numeric stability and possibly speed convergence. This version of solve() allows control over the
- * termination and iteration parameters.
- *
- * @param a The matrix A.
- * @param b The vector b.
- * @param preconditioner The preconditioner to apply.
- * @param maxIterations The maximum number of iterations to run.
- * @param maxError The maximum amount of residual error to tolerate. The algorithm will run until the residual falls
- * below this value or until maxIterations are completed.
- * @return The result x of solving the system.
- * @throws IllegalArgumentException if the matrix is not square, if the size of b is not equal to the number of
- * columns of A, if maxError is less than zero, or if maxIterations is not positive.
- */
-
- public Vector solve(VectorIterable a,
- Vector b,
- Preconditioner preconditioner,
- int maxIterations,
- double maxError) {
-
- if (a.numRows() != a.numCols()) {
- throw new IllegalArgumentException("Matrix must be square, symmetric and positive definite.");
- }
-
- if (a.numCols() != b.size()) {
- throw new CardinalityException(a.numCols(), b.size());
- }
-
- if (maxIterations <= 0) {
- throw new IllegalArgumentException("Max iterations must be positive.");
- }
-
- if (maxError < 0.0) {
- throw new IllegalArgumentException("Max error must be non-negative.");
- }
-
- Vector x = new DenseVector(b.size());
-
- iterations = 0;
- Vector residual = b.minus(a.times(x));
- residualNormSquared = residual.dot(residual);
-
- log.info("Conjugate gradient initial residual norm = {}", Math.sqrt(residualNormSquared));
- double previousConditionedNormSqr = 0.0;
- Vector updateDirection = null;
- while (Math.sqrt(residualNormSquared) > maxError && iterations < maxIterations) {
- Vector conditionedResidual;
- double conditionedNormSqr;
- if (preconditioner == null) {
- conditionedResidual = residual;
- conditionedNormSqr = residualNormSquared;
- } else {
- conditionedResidual = preconditioner.precondition(residual);
- conditionedNormSqr = residual.dot(conditionedResidual);
- }
-
- ++iterations;
-
- if (iterations == 1) {
- updateDirection = new DenseVector(conditionedResidual);
- } else {
- double beta = conditionedNormSqr / previousConditionedNormSqr;
-
- // updateDirection = residual + beta * updateDirection
- updateDirection.assign(Functions.MULT, beta);
- updateDirection.assign(conditionedResidual, Functions.PLUS);
- }
-
- Vector aTimesUpdate = a.times(updateDirection);
-
- double alpha = conditionedNormSqr / updateDirection.dot(aTimesUpdate);
-
- // x = x + alpha * updateDirection
- PLUS_MULT.setMultiplicator(alpha);
- x.assign(updateDirection, PLUS_MULT);
-
- // residual = residual - alpha * A * updateDirection
- PLUS_MULT.setMultiplicator(-alpha);
- residual.assign(aTimesUpdate, PLUS_MULT);
-
- previousConditionedNormSqr = conditionedNormSqr;
- residualNormSquared = residual.dot(residual);
-
- log.info("Conjugate gradient iteration {} residual norm = {}", iterations, Math.sqrt(residualNormSquared));
- }
- return x;
- }
-
- /**
- * Returns the number of iterations run once the solver is complete.
- *
- * @return The number of iterations run.
- */
- public int getIterations() {
- return iterations;
- }
-
- /**
- * Returns the norm of the residual at the completion of the solver. Usually this should be close to zero except in
- * the case of a non positive definite matrix A, which results in an unsolvable system, or for ill conditioned A, in
- * which case more iterations than the default may be needed.
- *
- * @return The norm of the residual in the solution.
- */
- public double getResidualNorm() {
- return Math.sqrt(residualNormSquared);
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/solver/EigenDecomposition.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/solver/EigenDecomposition.java b/math/src/main/java/org/apache/mahout/math/solver/EigenDecomposition.java
deleted file mode 100644
index 871ba44..0000000
--- a/math/src/main/java/org/apache/mahout/math/solver/EigenDecomposition.java
+++ /dev/null
@@ -1,892 +0,0 @@
-/*
- * 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.
- */
-
-/**
- * Adapted from the public domain Jama code.
- */
-
-package org.apache.mahout.math.solver;
-
-import org.apache.mahout.math.DenseMatrix;
-import org.apache.mahout.math.DenseVector;
-import org.apache.mahout.math.Matrix;
-import org.apache.mahout.math.Vector;
-import org.apache.mahout.math.function.Functions;
-
-/**
- * Eigenvalues and eigenvectors of a real matrix.
- * <p/>
- * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector
- * matrix V is orthogonal. I.e. A = V.times(D.times(V.transpose())) and V.times(V.transpose())
- * equals the identity matrix.
- * <p/>
- * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues
- * in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu,
- * lambda]. The columns of V represent the eigenvectors in the sense that A*V = V*D, i.e.
- * A.times(V) equals V.times(D). The matrix V may be badly conditioned, or even singular, so the
- * validity of the equation A = V*D*inverse(V) depends upon V.cond().
- */
-public class EigenDecomposition {
-
- /** Row and column dimension (square matrix). */
- private final int n;
- /** Arrays for internal storage of eigenvalues. */
- private final Vector d;
- private final Vector e;
- /** Array for internal storage of eigenvectors. */
- private final Matrix v;
-
- public EigenDecomposition(Matrix x) {
- this(x, isSymmetric(x));
- }
-
- public EigenDecomposition(Matrix x, boolean isSymmetric) {
- n = x.columnSize();
- d = new DenseVector(n);
- e = new DenseVector(n);
- v = new DenseMatrix(n, n);
-
- if (isSymmetric) {
- v.assign(x);
-
- // Tridiagonalize.
- tred2();
-
- // Diagonalize.
- tql2();
-
- } else {
- // Reduce to Hessenberg form.
- // Reduce Hessenberg to real Schur form.
- hqr2(orthes(x));
- }
- }
-
- /**
- * Return the eigenvector matrix
- *
- * @return V
- */
- public Matrix getV() {
- return v.like().assign(v);
- }
-
- /**
- * Return the real parts of the eigenvalues
- */
- public Vector getRealEigenvalues() {
- return d;
- }
-
- /**
- * Return the imaginary parts of the eigenvalues
- */
- public Vector getImagEigenvalues() {
- return e;
- }
-
- /**
- * Return the block diagonal eigenvalue matrix
- *
- * @return D
- */
- public Matrix getD() {
- Matrix x = new DenseMatrix(n, n);
- x.assign(0);
- x.viewDiagonal().assign(d);
- for (int i = 0; i < n; i++) {
- double v = e.getQuick(i);
- if (v > 0) {
- x.setQuick(i, i + 1, v);
- } else if (v < 0) {
- x.setQuick(i, i - 1, v);
- }
- }
- return x;
- }
-
- // Symmetric Householder reduction to tridiagonal form.
- private void tred2() {
- // This is derived from the Algol procedures tred2 by
- // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
- // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutine in EISPACK.
-
- d.assign(v.viewColumn(n - 1));
-
- // Householder reduction to tridiagonal form.
-
- for (int i = n - 1; i > 0; i--) {
-
- // Scale to avoid under/overflow.
-
- double scale = d.viewPart(0, i).norm(1);
- double h = 0.0;
-
-
- if (scale == 0.0) {
- e.setQuick(i, d.getQuick(i - 1));
- for (int j = 0; j < i; j++) {
- d.setQuick(j, v.getQuick(i - 1, j));
- v.setQuick(i, j, 0.0);
- v.setQuick(j, i, 0.0);
- }
- } else {
-
- // Generate Householder vector.
-
- for (int k = 0; k < i; k++) {
- d.setQuick(k, d.getQuick(k) / scale);
- h += d.getQuick(k) * d.getQuick(k);
- }
- double f = d.getQuick(i - 1);
- double g = Math.sqrt(h);
- if (f > 0) {
- g = -g;
- }
- e.setQuick(i, scale * g);
- h -= f * g;
- d.setQuick(i - 1, f - g);
- for (int j = 0; j < i; j++) {
- e.setQuick(j, 0.0);
- }
-
- // Apply similarity transformation to remaining columns.
-
- for (int j = 0; j < i; j++) {
- f = d.getQuick(j);
- v.setQuick(j, i, f);
- g = e.getQuick(j) + v.getQuick(j, j) * f;
- for (int k = j + 1; k <= i - 1; k++) {
- g += v.getQuick(k, j) * d.getQuick(k);
- e.setQuick(k, e.getQuick(k) + v.getQuick(k, j) * f);
- }
- e.setQuick(j, g);
- }
- f = 0.0;
- for (int j = 0; j < i; j++) {
- e.setQuick(j, e.getQuick(j) / h);
- f += e.getQuick(j) * d.getQuick(j);
- }
- double hh = f / (h + h);
- for (int j = 0; j < i; j++) {
- e.setQuick(j, e.getQuick(j) - hh * d.getQuick(j));
- }
- for (int j = 0; j < i; j++) {
- f = d.getQuick(j);
- g = e.getQuick(j);
- for (int k = j; k <= i - 1; k++) {
- v.setQuick(k, j, v.getQuick(k, j) - (f * e.getQuick(k) + g * d.getQuick(k)));
- }
- d.setQuick(j, v.getQuick(i - 1, j));
- v.setQuick(i, j, 0.0);
- }
- }
- d.setQuick(i, h);
- }
-
- // Accumulate transformations.
-
- for (int i = 0; i < n - 1; i++) {
- v.setQuick(n - 1, i, v.getQuick(i, i));
- v.setQuick(i, i, 1.0);
- double h = d.getQuick(i + 1);
- if (h != 0.0) {
- for (int k = 0; k <= i; k++) {
- d.setQuick(k, v.getQuick(k, i + 1) / h);
- }
- for (int j = 0; j <= i; j++) {
- double g = 0.0;
- for (int k = 0; k <= i; k++) {
- g += v.getQuick(k, i + 1) * v.getQuick(k, j);
- }
- for (int k = 0; k <= i; k++) {
- v.setQuick(k, j, v.getQuick(k, j) - g * d.getQuick(k));
- }
- }
- }
- for (int k = 0; k <= i; k++) {
- v.setQuick(k, i + 1, 0.0);
- }
- }
- d.assign(v.viewRow(n - 1));
- v.viewRow(n - 1).assign(0);
- v.setQuick(n - 1, n - 1, 1.0);
- e.setQuick(0, 0.0);
- }
-
- // Symmetric tridiagonal QL algorithm.
- private void tql2() {
-
- // This is derived from the Algol procedures tql2, by
- // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
- // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutine in EISPACK.
-
- e.viewPart(0, n - 1).assign(e.viewPart(1, n - 1));
- e.setQuick(n - 1, 0.0);
-
- double f = 0.0;
- double tst1 = 0.0;
- double eps = Math.pow(2.0, -52.0);
- for (int l = 0; l < n; l++) {
-
- // Find small subdiagonal element
-
- tst1 = Math.max(tst1, Math.abs(d.getQuick(l)) + Math.abs(e.getQuick(l)));
- int m = l;
- while (m < n) {
- if (Math.abs(e.getQuick(m)) <= eps * tst1) {
- break;
- }
- m++;
- }
-
- // If m == l, d.getQuick(l) is an eigenvalue,
- // otherwise, iterate.
-
- if (m > l) {
- do {
- // Compute implicit shift
-
- double g = d.getQuick(l);
- double p = (d.getQuick(l + 1) - g) / (2.0 * e.getQuick(l));
- double r = Math.hypot(p, 1.0);
- if (p < 0) {
- r = -r;
- }
- d.setQuick(l, e.getQuick(l) / (p + r));
- d.setQuick(l + 1, e.getQuick(l) * (p + r));
- double dl1 = d.getQuick(l + 1);
- double h = g - d.getQuick(l);
- for (int i = l + 2; i < n; i++) {
- d.setQuick(i, d.getQuick(i) - h);
- }
- f += h;
-
- // Implicit QL transformation.
-
- p = d.getQuick(m);
- double c = 1.0;
- double c2 = c;
- double c3 = c;
- double el1 = e.getQuick(l + 1);
- double s = 0.0;
- double s2 = 0.0;
- for (int i = m - 1; i >= l; i--) {
- c3 = c2;
- c2 = c;
- s2 = s;
- g = c * e.getQuick(i);
- h = c * p;
- r = Math.hypot(p, e.getQuick(i));
- e.setQuick(i + 1, s * r);
- s = e.getQuick(i) / r;
- c = p / r;
- p = c * d.getQuick(i) - s * g;
- d.setQuick(i + 1, h + s * (c * g + s * d.getQuick(i)));
-
- // Accumulate transformation.
-
- for (int k = 0; k < n; k++) {
- h = v.getQuick(k, i + 1);
- v.setQuick(k, i + 1, s * v.getQuick(k, i) + c * h);
- v.setQuick(k, i, c * v.getQuick(k, i) - s * h);
- }
- }
- p = -s * s2 * c3 * el1 * e.getQuick(l) / dl1;
- e.setQuick(l, s * p);
- d.setQuick(l, c * p);
-
- // Check for convergence.
-
- } while (Math.abs(e.getQuick(l)) > eps * tst1);
- }
- d.setQuick(l, d.getQuick(l) + f);
- e.setQuick(l, 0.0);
- }
-
- // Sort eigenvalues and corresponding vectors.
-
- for (int i = 0; i < n - 1; i++) {
- int k = i;
- double p = d.getQuick(i);
- for (int j = i + 1; j < n; j++) {
- if (d.getQuick(j) > p) {
- k = j;
- p = d.getQuick(j);
- }
- }
- if (k != i) {
- d.setQuick(k, d.getQuick(i));
- d.setQuick(i, p);
- for (int j = 0; j < n; j++) {
- p = v.getQuick(j, i);
- v.setQuick(j, i, v.getQuick(j, k));
- v.setQuick(j, k, p);
- }
- }
- }
- }
-
- // Nonsymmetric reduction to Hessenberg form.
- private Matrix orthes(Matrix x) {
- // Working storage for nonsymmetric algorithm.
- Vector ort = new DenseVector(n);
- Matrix hessenBerg = new DenseMatrix(n, n).assign(x);
-
- // This is derived from the Algol procedures orthes and ortran,
- // by Martin and Wilkinson, Handbook for Auto. Comp.,
- // Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutines in EISPACK.
-
- int low = 0;
- int high = n - 1;
-
- for (int m = low + 1; m <= high - 1; m++) {
-
- // Scale column.
-
- Vector hColumn = hessenBerg.viewColumn(m - 1).viewPart(m, high - m + 1);
- double scale = hColumn.norm(1);
-
- if (scale != 0.0) {
- // Compute Householder transformation.
-
- ort.viewPart(m, high - m + 1).assign(hColumn, Functions.plusMult(1 / scale));
- double h = ort.viewPart(m, high - m + 1).getLengthSquared();
-
- double g = Math.sqrt(h);
- if (ort.getQuick(m) > 0) {
- g = -g;
- }
- h -= ort.getQuick(m) * g;
- ort.setQuick(m, ort.getQuick(m) - g);
-
- // Apply Householder similarity transformation
- // H = (I-u*u'/h)*H*(I-u*u')/h)
-
- Vector ortPiece = ort.viewPart(m, high - m + 1);
- for (int j = m; j < n; j++) {
- double f = ortPiece.dot(hessenBerg.viewColumn(j).viewPart(m, high - m + 1)) / h;
- hessenBerg.viewColumn(j).viewPart(m, high - m + 1).assign(ortPiece, Functions.plusMult(-f));
- }
-
- for (int i = 0; i <= high; i++) {
- double f = ortPiece.dot(hessenBerg.viewRow(i).viewPart(m, high - m + 1)) / h;
- hessenBerg.viewRow(i).viewPart(m, high - m + 1).assign(ortPiece, Functions.plusMult(-f));
- }
- ort.setQuick(m, scale * ort.getQuick(m));
- hessenBerg.setQuick(m, m - 1, scale * g);
- }
- }
-
- // Accumulate transformations (Algol's ortran).
-
- v.assign(0);
- v.viewDiagonal().assign(1);
-
- for (int m = high - 1; m >= low + 1; m--) {
- if (hessenBerg.getQuick(m, m - 1) != 0.0) {
- ort.viewPart(m + 1, high - m).assign(hessenBerg.viewColumn(m - 1).viewPart(m + 1, high - m));
- for (int j = m; j <= high; j++) {
- double g = ort.viewPart(m, high - m + 1).dot(v.viewColumn(j).viewPart(m, high - m + 1));
- // Double division avoids possible underflow
- g = g / ort.getQuick(m) / hessenBerg.getQuick(m, m - 1);
- v.viewColumn(j).viewPart(m, high - m + 1).assign(ort.viewPart(m, high - m + 1), Functions.plusMult(g));
- }
- }
- }
- return hessenBerg;
- }
-
-
- // Complex scalar division.
- private double cdivr;
- private double cdivi;
-
- private void cdiv(double xr, double xi, double yr, double yi) {
- double r;
- double d;
- if (Math.abs(yr) > Math.abs(yi)) {
- r = yi / yr;
- d = yr + r * yi;
- cdivr = (xr + r * xi) / d;
- cdivi = (xi - r * xr) / d;
- } else {
- r = yr / yi;
- d = yi + r * yr;
- cdivr = (r * xr + xi) / d;
- cdivi = (r * xi - xr) / d;
- }
- }
-
-
- // Nonsymmetric reduction from Hessenberg to real Schur form.
-
- private void hqr2(Matrix h) {
-
- // This is derived from the Algol procedure hqr2,
- // by Martin and Wilkinson, Handbook for Auto. Comp.,
- // Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutine in EISPACK.
-
- // Initialize
-
- int nn = this.n;
- int n = nn - 1;
- int low = 0;
- int high = nn - 1;
- double eps = Math.pow(2.0, -52.0);
- double exshift = 0.0;
- double p = 0;
- double q = 0;
- double r = 0;
- double s = 0;
- double z = 0;
- double w;
- double x;
- double y;
-
- // Store roots isolated by balanc and compute matrix norm
-
- double norm = h.aggregate(Functions.PLUS, Functions.ABS);
-
- // Outer loop over eigenvalue index
-
- int iter = 0;
- while (n >= low) {
-
- // Look for single small sub-diagonal element
-
- int l = n;
- while (l > low) {
- s = Math.abs(h.getQuick(l - 1, l - 1)) + Math.abs(h.getQuick(l, l));
- if (s == 0.0) {
- s = norm;
- }
- if (Math.abs(h.getQuick(l, l - 1)) < eps * s) {
- break;
- }
- l--;
- }
-
- // Check for convergence
-
- if (l == n) {
- // One root found
- h.setQuick(n, n, h.getQuick(n, n) + exshift);
- d.setQuick(n, h.getQuick(n, n));
- e.setQuick(n, 0.0);
- n--;
- iter = 0;
-
-
- } else if (l == n - 1) {
- // Two roots found
- w = h.getQuick(n, n - 1) * h.getQuick(n - 1, n);
- p = (h.getQuick(n - 1, n - 1) - h.getQuick(n, n)) / 2.0;
- q = p * p + w;
- z = Math.sqrt(Math.abs(q));
- h.setQuick(n, n, h.getQuick(n, n) + exshift);
- h.setQuick(n - 1, n - 1, h.getQuick(n - 1, n - 1) + exshift);
- x = h.getQuick(n, n);
-
- // Real pair
- if (q >= 0) {
- if (p >= 0) {
- z = p + z;
- } else {
- z = p - z;
- }
- d.setQuick(n - 1, x + z);
- d.setQuick(n, d.getQuick(n - 1));
- if (z != 0.0) {
- d.setQuick(n, x - w / z);
- }
- e.setQuick(n - 1, 0.0);
- e.setQuick(n, 0.0);
- x = h.getQuick(n, n - 1);
- s = Math.abs(x) + Math.abs(z);
- p = x / s;
- q = z / s;
- r = Math.sqrt(p * p + q * q);
- p /= r;
- q /= r;
-
- // Row modification
-
- for (int j = n - 1; j < nn; j++) {
- z = h.getQuick(n - 1, j);
- h.setQuick(n - 1, j, q * z + p * h.getQuick(n, j));
- h.setQuick(n, j, q * h.getQuick(n, j) - p * z);
- }
-
- // Column modification
-
- for (int i = 0; i <= n; i++) {
- z = h.getQuick(i, n - 1);
- h.setQuick(i, n - 1, q * z + p * h.getQuick(i, n));
- h.setQuick(i, n, q * h.getQuick(i, n) - p * z);
- }
-
- // Accumulate transformations
-
- for (int i = low; i <= high; i++) {
- z = v.getQuick(i, n - 1);
- v.setQuick(i, n - 1, q * z + p * v.getQuick(i, n));
- v.setQuick(i, n, q * v.getQuick(i, n) - p * z);
- }
-
- // Complex pair
-
- } else {
- d.setQuick(n - 1, x + p);
- d.setQuick(n, x + p);
- e.setQuick(n - 1, z);
- e.setQuick(n, -z);
- }
- n -= 2;
- iter = 0;
-
- // No convergence yet
-
- } else {
-
- // Form shift
-
- x = h.getQuick(n, n);
- y = 0.0;
- w = 0.0;
- if (l < n) {
- y = h.getQuick(n - 1, n - 1);
- w = h.getQuick(n, n - 1) * h.getQuick(n - 1, n);
- }
-
- // Wilkinson's original ad hoc shift
-
- if (iter == 10) {
- exshift += x;
- for (int i = low; i <= n; i++) {
- h.setQuick(i, i, x);
- }
- s = Math.abs(h.getQuick(n, n - 1)) + Math.abs(h.getQuick(n - 1, n - 2));
- x = y = 0.75 * s;
- w = -0.4375 * s * s;
- }
-
- // MATLAB's new ad hoc shift
-
- if (iter == 30) {
- s = (y - x) / 2.0;
- s = s * s + w;
- if (s > 0) {
- s = Math.sqrt(s);
- if (y < x) {
- s = -s;
- }
- s = x - w / ((y - x) / 2.0 + s);
- for (int i = low; i <= n; i++) {
- h.setQuick(i, i, h.getQuick(i, i) - s);
- }
- exshift += s;
- x = y = w = 0.964;
- }
- }
-
- iter++; // (Could check iteration count here.)
-
- // Look for two consecutive small sub-diagonal elements
-
- int m = n - 2;
- while (m >= l) {
- z = h.getQuick(m, m);
- r = x - z;
- s = y - z;
- p = (r * s - w) / h.getQuick(m + 1, m) + h.getQuick(m, m + 1);
- q = h.getQuick(m + 1, m + 1) - z - r - s;
- r = h.getQuick(m + 2, m + 1);
- s = Math.abs(p) + Math.abs(q) + Math.abs(r);
- p /= s;
- q /= s;
- r /= s;
- if (m == l) {
- break;
- }
- double hmag = Math.abs(h.getQuick(m - 1, m - 1)) + Math.abs(h.getQuick(m + 1, m + 1));
- double threshold = eps * Math.abs(p) * (Math.abs(z) + hmag);
- if (Math.abs(h.getQuick(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < threshold) {
- break;
- }
- m--;
- }
-
- for (int i = m + 2; i <= n; i++) {
- h.setQuick(i, i - 2, 0.0);
- if (i > m + 2) {
- h.setQuick(i, i - 3, 0.0);
- }
- }
-
- // Double QR step involving rows l:n and columns m:n
-
- for (int k = m; k <= n - 1; k++) {
- boolean notlast = k != n - 1;
- if (k != m) {
- p = h.getQuick(k, k - 1);
- q = h.getQuick(k + 1, k - 1);
- r = notlast ? h.getQuick(k + 2, k - 1) : 0.0;
- x = Math.abs(p) + Math.abs(q) + Math.abs(r);
- if (x != 0.0) {
- p /= x;
- q /= x;
- r /= x;
- }
- }
- if (x == 0.0) {
- break;
- }
- s = Math.sqrt(p * p + q * q + r * r);
- if (p < 0) {
- s = -s;
- }
- if (s != 0) {
- if (k != m) {
- h.setQuick(k, k - 1, -s * x);
- } else if (l != m) {
- h.setQuick(k, k - 1, -h.getQuick(k, k - 1));
- }
- p += s;
- x = p / s;
- y = q / s;
- z = r / s;
- q /= p;
- r /= p;
-
- // Row modification
-
- for (int j = k; j < nn; j++) {
- p = h.getQuick(k, j) + q * h.getQuick(k + 1, j);
- if (notlast) {
- p += r * h.getQuick(k + 2, j);
- h.setQuick(k + 2, j, h.getQuick(k + 2, j) - p * z);
- }
- h.setQuick(k, j, h.getQuick(k, j) - p * x);
- h.setQuick(k + 1, j, h.getQuick(k + 1, j) - p * y);
- }
-
- // Column modification
-
- for (int i = 0; i <= Math.min(n, k + 3); i++) {
- p = x * h.getQuick(i, k) + y * h.getQuick(i, k + 1);
- if (notlast) {
- p += z * h.getQuick(i, k + 2);
- h.setQuick(i, k + 2, h.getQuick(i, k + 2) - p * r);
- }
- h.setQuick(i, k, h.getQuick(i, k) - p);
- h.setQuick(i, k + 1, h.getQuick(i, k + 1) - p * q);
- }
-
- // Accumulate transformations
-
- for (int i = low; i <= high; i++) {
- p = x * v.getQuick(i, k) + y * v.getQuick(i, k + 1);
- if (notlast) {
- p += z * v.getQuick(i, k + 2);
- v.setQuick(i, k + 2, v.getQuick(i, k + 2) - p * r);
- }
- v.setQuick(i, k, v.getQuick(i, k) - p);
- v.setQuick(i, k + 1, v.getQuick(i, k + 1) - p * q);
- }
- } // (s != 0)
- } // k loop
- } // check convergence
- } // while (n >= low)
-
- // Backsubstitute to find vectors of upper triangular form
-
- if (norm == 0.0) {
- return;
- }
-
- for (n = nn - 1; n >= 0; n--) {
- p = d.getQuick(n);
- q = e.getQuick(n);
-
- // Real vector
-
- double t;
- if (q == 0) {
- int l = n;
- h.setQuick(n, n, 1.0);
- for (int i = n - 1; i >= 0; i--) {
- w = h.getQuick(i, i) - p;
- r = 0.0;
- for (int j = l; j <= n; j++) {
- r += h.getQuick(i, j) * h.getQuick(j, n);
- }
- if (e.getQuick(i) < 0.0) {
- z = w;
- s = r;
- } else {
- l = i;
- if (e.getQuick(i) == 0.0) {
- if (w == 0.0) {
- h.setQuick(i, n, -r / (eps * norm));
- } else {
- h.setQuick(i, n, -r / w);
- }
-
- // Solve real equations
-
- } else {
- x = h.getQuick(i, i + 1);
- y = h.getQuick(i + 1, i);
- q = (d.getQuick(i) - p) * (d.getQuick(i) - p) + e.getQuick(i) * e.getQuick(i);
- t = (x * s - z * r) / q;
- h.setQuick(i, n, t);
- if (Math.abs(x) > Math.abs(z)) {
- h.setQuick(i + 1, n, (-r - w * t) / x);
- } else {
- h.setQuick(i + 1, n, (-s - y * t) / z);
- }
- }
-
- // Overflow control
-
- t = Math.abs(h.getQuick(i, n));
- if (eps * t * t > 1) {
- for (int j = i; j <= n; j++) {
- h.setQuick(j, n, h.getQuick(j, n) / t);
- }
- }
- }
- }
-
- // Complex vector
-
- } else if (q < 0) {
- int l = n - 1;
-
- // Last vector component imaginary so matrix is triangular
-
- if (Math.abs(h.getQuick(n, n - 1)) > Math.abs(h.getQuick(n - 1, n))) {
- h.setQuick(n - 1, n - 1, q / h.getQuick(n, n - 1));
- h.setQuick(n - 1, n, -(h.getQuick(n, n) - p) / h.getQuick(n, n - 1));
- } else {
- cdiv(0.0, -h.getQuick(n - 1, n), h.getQuick(n - 1, n - 1) - p, q);
- h.setQuick(n - 1, n - 1, cdivr);
- h.setQuick(n - 1, n, cdivi);
- }
- h.setQuick(n, n - 1, 0.0);
- h.setQuick(n, n, 1.0);
- for (int i = n - 2; i >= 0; i--) {
- double ra = 0.0;
- double sa = 0.0;
- for (int j = l; j <= n; j++) {
- ra += h.getQuick(i, j) * h.getQuick(j, n - 1);
- sa += h.getQuick(i, j) * h.getQuick(j, n);
- }
- w = h.getQuick(i, i) - p;
-
- if (e.getQuick(i) < 0.0) {
- z = w;
- r = ra;
- s = sa;
- } else {
- l = i;
- if (e.getQuick(i) == 0) {
- cdiv(-ra, -sa, w, q);
- h.setQuick(i, n - 1, cdivr);
- h.setQuick(i, n, cdivi);
- } else {
-
- // Solve complex equations
-
- x = h.getQuick(i, i + 1);
- y = h.getQuick(i + 1, i);
- double vr = (d.getQuick(i) - p) * (d.getQuick(i) - p) + e.getQuick(i) * e.getQuick(i) - q * q;
- double vi = (d.getQuick(i) - p) * 2.0 * q;
- if (vr == 0.0 && vi == 0.0) {
- double hmag = Math.abs(x) + Math.abs(y);
- vr = eps * norm * (Math.abs(w) + Math.abs(q) + hmag + Math.abs(z));
- }
- cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
- h.setQuick(i, n - 1, cdivr);
- h.setQuick(i, n, cdivi);
- if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
- h.setQuick(i + 1, n - 1, (-ra - w * h.getQuick(i, n - 1) + q * h.getQuick(i, n)) / x);
- h.setQuick(i + 1, n, (-sa - w * h.getQuick(i, n) - q * h.getQuick(i, n - 1)) / x);
- } else {
- cdiv(-r - y * h.getQuick(i, n - 1), -s - y * h.getQuick(i, n), z, q);
- h.setQuick(i + 1, n - 1, cdivr);
- h.setQuick(i + 1, n, cdivi);
- }
- }
-
- // Overflow control
-
- t = Math.max(Math.abs(h.getQuick(i, n - 1)), Math.abs(h.getQuick(i, n)));
- if (eps * t * t > 1) {
- for (int j = i; j <= n; j++) {
- h.setQuick(j, n - 1, h.getQuick(j, n - 1) / t);
- h.setQuick(j, n, h.getQuick(j, n) / t);
- }
- }
- }
- }
- }
- }
-
- // Vectors of isolated roots
-
- for (int i = 0; i < nn; i++) {
- if (i < low || i > high) {
- for (int j = i; j < nn; j++) {
- v.setQuick(i, j, h.getQuick(i, j));
- }
- }
- }
-
- // Back transformation to get eigenvectors of original matrix
-
- for (int j = nn - 1; j >= low; j--) {
- for (int i = low; i <= high; i++) {
- z = 0.0;
- for (int k = low; k <= Math.min(j, high); k++) {
- z += v.getQuick(i, k) * h.getQuick(k, j);
- }
- v.setQuick(i, j, z);
- }
- }
- }
-
- private static boolean isSymmetric(Matrix a) {
- /*
- Symmetry flag.
- */
- int n = a.columnSize();
-
- boolean isSymmetric = true;
- for (int j = 0; (j < n) && isSymmetric; j++) {
- for (int i = 0; (i < n) && isSymmetric; i++) {
- isSymmetric = a.getQuick(i, j) == a.getQuick(j, i);
- }
- }
- return isSymmetric;
- }
-}
http://git-wip-us.apache.org/repos/asf/mahout/blob/e0573de3/math/src/main/java/org/apache/mahout/math/solver/JacobiConditioner.java
----------------------------------------------------------------------
diff --git a/math/src/main/java/org/apache/mahout/math/solver/JacobiConditioner.java b/math/src/main/java/org/apache/mahout/math/solver/JacobiConditioner.java
deleted file mode 100644
index 7524564..0000000
--- a/math/src/main/java/org/apache/mahout/math/solver/JacobiConditioner.java
+++ /dev/null
@@ -1,47 +0,0 @@
-/**
- * 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.mahout.math.solver;
-
-import org.apache.mahout.math.DenseVector;
-import org.apache.mahout.math.Matrix;
-import org.apache.mahout.math.Vector;
-
-/**
- * Implements the Jacobi preconditioner for a matrix A. This is defined as inv(diag(A)).
- */
-public final class JacobiConditioner implements Preconditioner {
-
- private final DenseVector inverseDiagonal;
-
- public JacobiConditioner(Matrix a) {
- if (a.numCols() != a.numRows()) {
- throw new IllegalArgumentException("Matrix must be square.");
- }
-
- inverseDiagonal = new DenseVector(a.numCols());
- for (int i = 0; i < a.numCols(); ++i) {
- inverseDiagonal.setQuick(i, 1.0 / a.getQuick(i, i));
- }
- }
-
- @Override
- public Vector precondition(Vector v) {
- return v.times(inverseDiagonal);
- }
-
-}