You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@ignite.apache.org by sb...@apache.org on 2017/04/18 05:25:52 UTC
[37/54] [abbrv] ignite git commit: IGNITE-5000 Rename Ignite Math
module to Ignite ML module
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/Tracer.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/Tracer.java b/modules/ml/src/main/java/org/apache/ignite/math/Tracer.java
new file mode 100644
index 0000000..89d4669
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/Tracer.java
@@ -0,0 +1,456 @@
+/*
+ * 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.ignite.math;
+
+import java.awt.Color;
+import java.awt.Desktop;
+import java.io.BufferedReader;
+import java.io.BufferedWriter;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.io.InputStreamReader;
+import java.nio.file.Files;
+import java.nio.file.Paths;
+import java.nio.file.StandardOpenOption;
+import java.util.Locale;
+import java.util.function.Function;
+import java.util.stream.Collectors;
+import org.apache.ignite.IgniteLogger;
+import org.apache.ignite.lang.IgniteUuid;
+
+/**
+ * Utility methods to support output of {@link Vector} and {@link Matrix} instances to plain text or HTML.
+ */
+public class Tracer {
+ /**
+ * Double to color mapper.
+ */
+ public interface ColorMapper extends Function<Double, Color> {
+ }
+
+ /** Continuous red-to-blue color mapping. */
+ static private ColorMapper defaultColorMapper(double min, double max) {
+ double range = max - min;
+
+ return new ColorMapper() {
+ /** {@inheritDoc} */
+ @Override public Color apply(Double d) {
+ int r = (int)Math.round(255 * d);
+ int g = 0;
+ int b = (int)Math.round(255 * (1 - d));
+
+ return new Color(r, g, b);
+ }
+ };
+ }
+
+ /**
+ * Default vector color mapper implementation that map given double value
+ * to continuous red-blue (R_B) specter.
+ *
+ * @param vec Vector to map.
+ * @return {@link ColorMapper} for the given vector.
+ */
+ static private ColorMapper mkVectorColorMapper(Vector vec) {
+ return defaultColorMapper(vec.minValue(), vec.maxValue());
+ }
+
+ /** Default matrix color mapper implementation that map given double value
+ * to continuous red-blue (R_B) specter.
+ * @param mtx Matrix to be mapped.
+ * @return Color mapper for given matrix.
+ */
+ static private ColorMapper mkMatrixColorMapper(Matrix mtx) {
+ return defaultColorMapper(mtx.minValue(), mtx.maxValue());
+ }
+
+ /**
+ * @param vec Vector to show.
+ * @param log {@link IgniteLogger} instance for output.
+ * @param fmt Format string for vector elements.
+ */
+ public static void showAscii(Vector vec, IgniteLogger log, String fmt) {
+ String cls = vec.getClass().getSimpleName();
+
+ log.info(String.format("%s(%d) [%s]", cls, vec.size(), mkString(vec, fmt)));
+ }
+
+ /**
+ * @param vec Vector to show as plain text.
+ * @param log {@link IgniteLogger} instance for output.
+ */
+ public static void showAscii(Vector vec, IgniteLogger log) {
+ showAscii(vec, log, "%4f");
+ }
+
+ /**
+ * @param vec Vector to show as plain text.
+ * @param fmt Format string for vector elements.
+ */
+ public static void showAscii(Vector vec, String fmt) {
+ String cls = vec.getClass().getSimpleName();
+
+ System.out.println(String.format("%s(%d) [%s]", cls, vec.size(), mkString(vec, fmt)));
+ }
+
+ /**
+ * @param mtx Matrix to show as plain text.
+ */
+ public static void showAscii(Matrix mtx) {
+ showAscii(mtx, "%4f");
+ }
+
+ /**
+ * @param mtx Matrix to show.
+ * @param row Matrix row to output.
+ * @param fmt Format string for matrix elements in the row.
+ * @return String representation of given matrix row according to given format.
+ */
+ static private String rowStr(Matrix mtx, int row, String fmt) {
+ StringBuilder buf = new StringBuilder();
+
+ boolean first = true;
+
+ int cols = mtx.columnSize();
+
+ for (int col = 0; col < cols; col++) {
+ String s = String.format(fmt, mtx.get(row, col));
+
+ if (!first)
+ buf.append(", ");
+
+ buf.append(s);
+
+ first = false;
+ }
+
+ return buf.toString();
+ }
+
+ /**
+ * @param mtx {@link Matrix} object to show as a plain text.
+ * @param fmt Format string for matrix rows.
+ */
+ public static void showAscii(Matrix mtx, String fmt) {
+ String cls = mtx.getClass().getSimpleName();
+
+ int rows = mtx.rowSize();
+ int cols = mtx.columnSize();
+
+ System.out.println(String.format("%s(%dx%d)", cls, rows, cols));
+
+ for (int row = 0; row < rows; row++)
+ System.out.println(rowStr(mtx, row, fmt));
+ }
+
+ /**
+ * @param mtx {@link Matrix} object to show as a plain text.
+ * @param log {@link IgniteLogger} instance to output the logged matrix.
+ * @param fmt Format string for matrix rows.
+ */
+ public static void showAscii(Matrix mtx, IgniteLogger log, String fmt) {
+ String cls = mtx.getClass().getSimpleName();
+
+ int rows = mtx.rowSize();
+ int cols = mtx.columnSize();
+
+ log.info(String.format("%s(%dx%d)", cls, rows, cols));
+
+ for (int row = 0; row < rows; row++)
+ log.info(rowStr(mtx, row, fmt));
+ }
+
+ /**
+ * @param vec {@link Vector} object to show as a plain text.
+ */
+ public static void showAscii(Vector vec) {
+ showAscii(vec, "%4f");
+ }
+
+ /**
+ * Saves given vector as CSV file.
+ *
+ * @param vec Vector to save.
+ * @param fmt Format to use.
+ * @param filePath Path of the file to save to.
+ */
+ public static void saveAsCsv(Vector vec, String fmt, String filePath) throws IOException {
+ String s = mkString(vec, fmt);
+
+ Files.write(Paths.get(filePath), s.getBytes(), StandardOpenOption.CREATE, StandardOpenOption.WRITE);
+ }
+
+ /**
+ * Saves given matrix as CSV file.
+ *
+ * @param mtx Matrix to save.
+ * @param fmt Format to use.
+ * @param filePath Path of the file to save to.
+ */
+ public static void saveAsCsv(Matrix mtx, String fmt, String filePath) throws IOException {
+ String s = mkString(mtx, fmt);
+
+ Files.write(Paths.get(filePath), s.getBytes(), StandardOpenOption.CREATE, StandardOpenOption.WRITE);
+ }
+
+ /**
+ * Shows given matrix in the browser with D3-based visualization.
+ *
+ * @param mtx Matrix to show.
+ * @throws IOException Thrown in case of any errors.
+ */
+ public static void showHtml(Matrix mtx) throws IOException {
+ showHtml(mtx, mkMatrixColorMapper(mtx));
+ }
+
+ /**
+ * Shows given matrix in the browser with D3-based visualization.
+ *
+ * @param mtx Matrix to show.
+ * @param cm Optional color mapper. If not provided - red-to-blue (R_B) mapper will be used.
+ * @throws IOException Thrown in case of any errors.
+ */
+ public static void showHtml(Matrix mtx, ColorMapper cm) throws IOException {
+ // Read it every time so that we can change it at runtime.
+ String tmpl = fileToString("d3-matrix-template.html");
+
+ String cls = mtx.getClass().getSimpleName();
+
+ double min = mtx.minValue();
+ double max = mtx.maxValue();
+
+ openHtmlFile(tmpl.
+ replaceAll("/\\*@NAME@\\*/.*\n", "var name = \"" + cls + "\";\n").
+ replaceAll("/\\*@MIN@\\*/.*\n", "var min = " + dataColorJson(min, cm.apply(min)) + ";\n").
+ replaceAll("/\\*@MAX@\\*/.*\n", "var max = " + dataColorJson(max, cm.apply(max)) + ";\n").
+ replaceAll("/\\*@DATA@\\*/.*\n", "var data = " + mkJsArrayString(mtx, cm) + ";\n")
+ );
+ }
+
+ /**
+ * Shows given vector in the browser with D3-based visualization.
+ *
+ * @param vec Vector to show.
+ * @throws IOException Thrown in case of any errors.
+ */
+ public static void showHtml(Vector vec) throws IOException {
+ showHtml(vec, mkVectorColorMapper(vec));
+ }
+
+ /**
+ * @param d Value of {@link Matrix} or {@link Vector} element.
+ * @param clr {@link Color} to paint.
+ * @return JSON representation for given value and color.
+ */
+ static private String dataColorJson(double d, Color clr) {
+ return "{" +
+ "d: " + String.format("%4f", d) +
+ ", r: " + clr.getRed() +
+ ", g: " + clr.getGreen() +
+ ", b: " + clr.getBlue() +
+ "}";
+ }
+
+ /**
+ * Shows given vector in the browser with D3-based visualization.
+ *
+ * @param vec Vector to show.
+ * @param cm Optional color mapper. If not provided - red-to-blue (R_B) mapper will be used.
+ * @throws IOException Thrown in case of any errors.
+ */
+ public static void showHtml(Vector vec, ColorMapper cm) throws IOException {
+ // Read it every time so that we can change it at runtime.
+ String tmpl = fileToString("d3-vector-template.html");
+
+ String cls = vec.getClass().getSimpleName();
+
+ double min = vec.minValue();
+ double max = vec.maxValue();
+
+ openHtmlFile(tmpl.
+ replaceAll("/\\*@NAME@\\*/.*\n", "var name = \"" + cls + "\";\n").
+ replaceAll("/\\*@MIN@\\*/.*\n", "var min = " + dataColorJson(min, cm.apply(min)) + ";\n").
+ replaceAll("/\\*@MAX@\\*/.*\n", "var max = " + dataColorJson(max, cm.apply(max)) + ";\n").
+ replaceAll("/\\*@DATA@\\*/.*\n", "var data = " + mkJsArrayString(vec, cm) + ";\n")
+ );
+ }
+
+ /**
+ * Reads file content into the string.
+ *
+ * @param fileName Name of the file (on classpath) to read.
+ * @return Content of the file.
+ * @throws IOException If an I/O error of some sort has occurred.
+ */
+ private static String fileToString(String fileName) throws IOException {
+ assert Tracer.class.getResourceAsStream(fileName) != null : "Can't get resource: " + fileName;
+
+ InputStreamReader is = new InputStreamReader(Tracer.class.getResourceAsStream(fileName));
+
+ String str = new BufferedReader(is).lines().collect(Collectors.joining("\n"));
+
+ is.close();
+
+ return str;
+ }
+
+ /**
+ * Opens file in the browser with given HTML content.
+ *
+ * @param html HTML content.
+ * @throws IOException Thrown in case of any errors.
+ */
+ static private void openHtmlFile(String html) throws IOException {
+ File temp = File.createTempFile(IgniteUuid.randomUuid().toString(), ".html");
+
+ BufferedWriter bw = new BufferedWriter(new FileWriter(temp));
+
+ bw.write(html);
+
+ bw.close();
+
+ Desktop.getDesktop().browse(temp.toURI());
+ }
+
+ /**
+ * Gets string presentation of this vector.
+ *
+ * @param vec Vector to string-ify.
+ * @param fmt {@link String#format(Locale, String, Object...)} format.
+ */
+ private static String mkString(Vector vec, String fmt) {
+ boolean first = true;
+
+ StringBuilder buf = new StringBuilder();
+
+ for (Vector.Element x : vec.all()) {
+ String s = String.format(Locale.US, fmt, x.get());
+
+ if (!first) {
+ buf.append(", ");
+ buf.append(s);
+ }
+ else {
+ buf.append(s);
+ first = false;
+ }
+ }
+
+ return buf.toString();
+ }
+
+ /**
+ * Gets JavaScript array presentation of this vector.
+ *
+ * @param vec Vector to JavaScript-ify.
+ * @param cm Color mapper to user.
+ */
+ private static String mkJsArrayString(Vector vec, ColorMapper cm) {
+ boolean first = true;
+
+ StringBuilder buf = new StringBuilder();
+
+ for (Vector.Element x : vec.all()) {
+ double d = x.get();
+
+ String s = dataColorJson(d, cm.apply(d));
+
+ if (!first)
+ buf.append(", ");
+
+ buf.append(s);
+
+ first = false;
+ }
+
+ return '[' + buf.toString() + ']';
+ }
+
+ /**
+ * Gets JavaScript array presentation of this vector.
+ *
+ * @param mtx Matrix to JavaScript-ify.
+ * @param cm Color mapper to user.
+ */
+ private static String mkJsArrayString(Matrix mtx, ColorMapper cm) {
+ boolean first = true;
+
+ StringBuilder buf = new StringBuilder();
+
+ int rows = mtx.rowSize();
+ int cols = mtx.columnSize();
+
+ for (int row = 0; row < rows; row++) {
+ StringBuilder rowBuf = new StringBuilder();
+
+ boolean rowFirst = true;
+
+ for (int col = 0; col < cols; col++) {
+ double d = mtx.get(row, col);
+
+ String s = dataColorJson(d, cm.apply(d));
+
+ if (!rowFirst)
+ rowBuf.append(", ");
+
+ rowBuf.append(s);
+
+ rowFirst = false;
+ }
+
+ if (!first)
+ buf.append(", ");
+
+ buf.append('[').append(rowBuf.toString()).append(']');
+
+ first = false;
+ }
+
+ return '[' + buf.toString() + ']';
+ }
+
+ /**
+ * @param mtx Matrix to log.
+ * @param fmt Output format.
+ * @return Formatted representation of a matrix.
+ */
+ private static String mkString(Matrix mtx, String fmt) {
+ StringBuilder buf = new StringBuilder();
+
+ int rows = mtx.rowSize();
+ int cols = mtx.columnSize();
+
+ for (int row = 0; row < rows; row++) {
+ for (int col = 0; col < cols; col++) {
+ String s = String.format(Locale.US, fmt, mtx.get(row, col));
+
+ if (col != 0)
+ buf.append(", ");
+
+ buf.append(s);
+
+ if (col == cols - 1 && row != rows - 1)
+ buf.append(",\n");
+
+ }
+ }
+
+ return buf.toString();
+ }
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/ValueMapper.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/ValueMapper.java b/modules/ml/src/main/java/org/apache/ignite/math/ValueMapper.java
new file mode 100644
index 0000000..9459bd1
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/ValueMapper.java
@@ -0,0 +1,27 @@
+// @java.file.header
+
+/* _________ _____ __________________ _____
+ * __ ____/___________(_)______ /__ ____/______ ____(_)_______
+ * _ / __ __ ___/__ / _ __ / _ / __ _ __ `/__ / __ __ \
+ * / /_/ / _ / _ / / /_/ / / /_/ / / /_/ / _ / _ / / /
+ * \____/ /_/ /_/ \_,__/ \____/ \__,_/ /_/ /_/ /_/
+ */
+
+package org.apache.ignite.math;
+
+import java.io.Serializable;
+
+/**
+ * Utility mapper that can be used to map arbitrary values types to and from double.
+ */
+public interface ValueMapper<V> extends Serializable {
+ /**
+ * @param v
+ */
+ public V fromDouble(double v);
+
+ /**
+ * @param v
+ */
+ public double toDouble(V v);
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/Vector.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/Vector.java b/modules/ml/src/main/java/org/apache/ignite/math/Vector.java
new file mode 100644
index 0000000..ac2a6c7
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/Vector.java
@@ -0,0 +1,498 @@
+/*
+ * 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.ignite.math;
+
+import java.io.Externalizable;
+import java.util.Spliterator;
+import java.util.function.IntToDoubleFunction;
+import org.apache.ignite.lang.IgniteUuid;
+import org.apache.ignite.math.exceptions.CardinalityException;
+import org.apache.ignite.math.exceptions.IndexException;
+import org.apache.ignite.math.exceptions.UnsupportedOperationException;
+import org.apache.ignite.math.functions.IgniteBiFunction;
+import org.apache.ignite.math.functions.IgniteDoubleFunction;
+
+/**
+ * A vector interface.
+ *
+ * Based on its flavor it can have vastly different implementations tailored for
+ * for different types of data (e.g. dense vs. sparse), different sizes of data or different operation
+ * optimizations.
+ *
+ * Note also that not all operations can be supported by all underlying implementations. If an operation is not
+ * supported a {@link UnsupportedOperationException} is thrown. This exception can also be thrown in partial cases
+ * where an operation is unsupported only in special cases, e.g. where a given operation cannot be deterministically
+ * completed in polynomial time.
+ *
+ * Based on ideas from <a href="http://mahout.apache.org/">Apache Mahout</a>.
+ */
+public interface Vector extends MetaAttributes, Externalizable, StorageOpsMetrics, Destroyable {
+ /**
+ * Holder for vector's element.
+ */
+ interface Element {
+ /**
+ * Gets element's value.
+ *
+ * @return The value of this vector element.
+ */
+ double get();
+
+ /**
+ * Gets element's index in the vector.
+ *
+ * @return The index of this vector element.
+ */
+ int index();
+
+ /**
+ * Sets element's value.
+ *
+ * @param val Value to set.
+ */
+ void set(double val);
+ }
+
+ /**
+ * Gets cardinality of this vector (maximum number of the elements).
+ *
+ * @return This vector's cardinality.
+ */
+ public int size();
+
+ /**
+ * Creates new copy of this vector.
+ *
+ * @return New copy vector.
+ */
+ public Vector copy();
+
+ /**
+ * Gets iterator over all elements in this vector.
+ *
+ * NOTE: implementation can choose to reuse {@link Element} instance so you need to copy it
+ * if you want to retain it outside of iteration.
+ *
+ * @return Iterator.
+ */
+ public Iterable<Element> all();
+
+ /**
+ * Iterates ove all non-zero elements in this vector.
+ *
+ * NOTE: implementation can choose to reuse {@link Element} instance so you need to copy it
+ * if you want to retain it outside of iteration.
+ *
+ * @return Iterator.
+ */
+ public Iterable<Element> nonZeroes();
+
+ /**
+ * Gets spliterator for all values in this vector.
+ *
+ * @return Spliterator for all values.
+ */
+ public Spliterator<Double> allSpliterator();
+
+ /**
+ * Gets spliterator for all non-zero values in this vector.
+ *
+ * @return Spliterator for all non-zero values.
+ */
+ public Spliterator<Double> nonZeroSpliterator();
+
+ /**
+ * Sorts this vector in ascending order.
+ */
+ public Vector sort();
+
+ /**
+ * Gets element at the given index.
+ *
+ * NOTE: implementation can choose to reuse {@link Element} instance so you need to copy it
+ * if you want to retain it outside of iteration.
+ *
+ * @param idx Element's index.
+ * @return Vector's element at the given index.
+ * @throws IndexException Throw if index is out of bounds.
+ */
+ public Element getElement(int idx);
+
+ /**
+ * Assigns given value to all elements of this vector.
+ *
+ * @param val Value to assign.
+ * @return This vector.
+ */
+ public Vector assign(double val);
+
+ /**
+ * Assigns values from given array to this vector.
+ *
+ * @param vals Values to assign.
+ * @return This vector.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector assign(double[] vals);
+
+ /**
+ * Copies values from the argument vector to this one.
+ *
+ * @param vec Argument vector.
+ * @return This vector.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector assign(Vector vec);
+
+ /**
+ * Assigns each vector element to the value generated by given function.
+ *
+ * @param fun Function that takes the index and returns value.
+ * @return This vector.
+ */
+ public Vector assign(IntToDoubleFunction fun);
+
+ /**
+ * Maps all values in this vector through a given function.
+ *
+ * @param fun Mapping function.
+ * @return This vector.
+ */
+ public Vector map(IgniteDoubleFunction<Double> fun);
+
+ /**
+ * Maps all values in this vector through a given function.
+ *
+ * For this vector <code>A</code>, argument vector <code>B</code> and the
+ * function <code>F</code> this method maps every element <code>x</code> as:
+ * <code>A(x) = F(A(x), B(x))</code>
+ *
+ * @param vec Argument vector.
+ * @param fun Mapping function.
+ * @return This function.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector map(Vector vec, IgniteBiFunction<Double, Double, Double> fun);
+
+ /**
+ * Maps all elements of this vector by applying given function to each element with a constant
+ * second parameter <code>y</code>.
+ *
+ * @param fun Mapping function.
+ * @param y Second parameter for mapping function.
+ * @return This vector.
+ */
+ public Vector map(IgniteBiFunction<Double, Double, Double> fun, double y);
+
+ /**
+ * Creates new vector containing values from this vector divided by the argument.
+ *
+ * @param x Division argument.
+ * @return New vector.
+ */
+ public Vector divide(double x);
+
+ /**
+ * Gets dot product of two vectors.
+ *
+ * @param vec Argument vector.
+ * @return Dot product of two vectors.
+ */
+ public double dot(Vector vec);
+
+ /**
+ * Gets the value at specified index.
+ *
+ * @param idx Vector index.
+ * @return Vector value.
+ * @throws IndexException Throw if index is out of bounds.
+ */
+ public double get(int idx);
+
+ /**
+ * Gets the value at specified index without checking for index boundaries.
+ *
+ * @param idx Vector index.
+ * @return Vector value.
+ */
+ public double getX(int idx);
+
+ /**
+ * Creates new empty vector of the same underlying class but of different cardinality.
+ *
+ * @param crd Cardinality for new vector.
+ * @return New vector.
+ */
+ public Vector like(int crd);
+
+ /**
+ * Creates new matrix of compatible flavor with given size.
+ *
+ * @param rows Number of rows.
+ * @param cols Number of columns.
+ * @return New matrix.
+ */
+ public Matrix likeMatrix(int rows, int cols);
+
+ /**
+ * Converts this vector into [N x 1] or [1 x N] matrix where N is this vector cardinality.
+ *
+ * @param rowLike {@code true} for rowLike [N x 1], or {@code false} for column [1 x N] matrix.
+ * @return Newly created matrix.
+ */
+ public Matrix toMatrix(boolean rowLike);
+
+ /**
+ * Converts this vector into [N+1 x 1] or [1 x N+1] matrix where N is this vector cardinality.
+ * (0,0) element of this matrix will be {@code zeroVal} parameter.
+ *
+ * @param rowLike {@code true} for rowLike [N+1 x 1], or {@code false} for column [1 x N+1] matrix.
+ * @return Newly created matrix.
+ */
+ public Matrix toMatrixPlusOne(boolean rowLike, double zeroVal);
+
+ /**
+ * Creates new vector containing element by element difference between this vector and the argument one.
+ *
+ * @param vec Argument vector.
+ * @return New vector.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector minus(Vector vec);
+
+ /**
+ * Creates new vector containing the normalized (L_2 norm) values of this vector.
+ *
+ * @return New vector.
+ */
+ public Vector normalize();
+
+ /**
+ * Creates new vector containing the normalized (L_power norm) values of this vector.
+ * See http://en.wikipedia.org/wiki/Lp_space for details.
+ *
+ * @param power The power to use. Must be >= 0. May also be {@link Double#POSITIVE_INFINITY}.
+ * @return New vector {@code x} such that {@code norm(x, power) == 1}
+ */
+ public Vector normalize(double power);
+
+ /**
+ * Creates new vector containing the {@code log(1 + entry) / L_2 norm} values of this vector.
+ *
+ * @return New vector.
+ */
+ public Vector logNormalize();
+
+ /**
+ * Creates new vector with a normalized value calculated as {@code log_power(1 + entry) / L_power norm}.
+ *
+ * @param power The power to use. Must be > 1. Cannot be {@link Double#POSITIVE_INFINITY}.
+ * @return New vector
+ */
+ public Vector logNormalize(double power);
+
+ /**
+ * Gets the k-norm of the vector. See http://en.wikipedia.org/wiki/Lp_space for more details.
+ *
+ * @param power The power to use.
+ * @see #normalize(double)
+ */
+ public double kNorm(double power);
+
+ /**
+ * Gets minimal value in this vector.
+ *
+ * @return Minimal value.
+ */
+ public double minValue();
+
+ /**
+ * Gets maximum value in this vector.
+ *
+ * @return Maximum c.
+ */
+ public double maxValue();
+
+ /**
+ * Gets minimal element in this vector.
+ *
+ * @return Minimal element.
+ */
+ public Element minElement();
+
+ /**
+ * Gets maximum element in this vector.
+ *
+ * @return Maximum element.
+ */
+ public Element maxElement();
+
+ /**
+ * Creates new vector containing sum of each element in this vector and argument.
+ *
+ * @param x Argument value.
+ * @return New vector.
+ */
+ public Vector plus(double x);
+
+ /**
+ * Creates new vector containing element by element sum from both vectors.
+ *
+ * @param vec Other argument vector to add.
+ * @return New vector.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector plus(Vector vec);
+
+ /**
+ * Sets value.
+ *
+ * @param idx Vector index to set value at.
+ * @param val Value to set.
+ * @return This vector.
+ * @throws IndexException Throw if index is out of bounds.
+ */
+ public Vector set(int idx, double val);
+
+ /**
+ * Sets value without checking for index boundaries.
+ *
+ * @param idx Vector index to set value at.
+ * @param val Value to set.
+ * @return This vector.
+ */
+ public Vector setX(int idx, double val);
+
+ /**
+ * Increments value at given index without checking for index boundaries.
+ *
+ * @param idx Vector index.
+ * @param val Increment value.
+ * @return This vector.
+ */
+ public Vector incrementX(int idx, double val);
+
+ /**
+ * Increments value at given index.
+ *
+ * @param idx Vector index.
+ * @param val Increment value.
+ * @return This vector.
+ * @throws IndexException Throw if index is out of bounds.
+ */
+ public Vector increment(int idx, double val);
+
+ /**
+ * Gets number of non-zero elements in this vector.
+ *
+ * @return Number of non-zero elements in this vector.
+ */
+ public int nonZeroElements();
+
+ /**
+ * Gets a new vector that contains product of each element and the argument.
+ *
+ * @param x Multiply argument.
+ * @return New vector.
+ */
+ public Vector times(double x);
+
+ /**
+ * Gets a new vector that is an element-wie product of this vector and the argument.
+ *
+ * @param vec Vector to multiply by.
+ * @return New vector.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public Vector times(Vector vec);
+
+ /**
+ * @param off Offset into parent vector.
+ * @param len Length of the view.
+ */
+ public Vector viewPart(int off, int len);
+
+ /**
+ * Gets vector storage model.
+ */
+ public VectorStorage getStorage();
+
+ /**
+ * Gets the sum of all elements in this vector.
+ *
+ * @return Vector's sum
+ */
+ public double sum();
+
+ /**
+ * Gets the cross product of this vector and the other vector.
+ *
+ * @param vec Second vector.
+ * @return New matrix as a cross product of two vectors.
+ */
+ public Matrix cross(Vector vec);
+
+ /**
+ * Folds this vector into a single value.
+ *
+ * @param foldFun Folding function that takes two parameters: accumulator and the current value.
+ * @param mapFun Mapping function that is called on each vector element before its passed to the accumulator (as its
+ * second parameter).
+ * @param <T> Type of the folded value.
+ * @param zeroVal Zero value for fold operation.
+ * @return Folded value of this vector.
+ */
+ public <T> T foldMap(IgniteBiFunction<T, Double, T> foldFun, IgniteDoubleFunction<Double> mapFun, T zeroVal);
+
+ /**
+ * Combines & maps two vector and folds them into a single value.
+ *
+ * @param vec Another vector to combine with.
+ * @param foldFun Folding function.
+ * @param combFun Combine function.
+ * @param <T> Type of the folded value.
+ * @param zeroVal Zero value for fold operation.
+ * @return Folded value of these vectors.
+ * @throws CardinalityException Thrown when cardinalities mismatch.
+ */
+ public <T> T foldMap(Vector vec, IgniteBiFunction<T, Double, T> foldFun, IgniteBiFunction<Double, Double, Double> combFun,
+ T zeroVal);
+
+ /**
+ * Gets the sum of squares of all elements in this vector.
+ *
+ * @return Length squared value.
+ */
+ public double getLengthSquared();
+
+ /**
+ * Get the square of the distance between this vector and the argument vector.
+ *
+ * @param vec Another vector.
+ * @return Distance squared.
+ * @throws CardinalityException Thrown if cardinalities mismatch.
+ */
+ public double getDistanceSquared(Vector vec);
+
+ /**
+ * Auto-generated globally unique vector ID.
+ *
+ * @return Vector GUID.
+ */
+ public IgniteUuid guid();
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/VectorKeyMapper.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/VectorKeyMapper.java b/modules/ml/src/main/java/org/apache/ignite/math/VectorKeyMapper.java
new file mode 100644
index 0000000..17d76f5
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/VectorKeyMapper.java
@@ -0,0 +1,29 @@
+/*
+ * 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.ignite.math;
+
+/**
+ * Maps {@link Vector} element index to cache key.
+ */
+public interface VectorKeyMapper<K> extends KeyMapper<K> {
+ /**
+ * @param i Vector element index.
+ * @return Cache key for given element index.
+ */
+ public K apply(int i);
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/VectorStorage.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/VectorStorage.java b/modules/ml/src/main/java/org/apache/ignite/math/VectorStorage.java
new file mode 100644
index 0000000..f410254
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/VectorStorage.java
@@ -0,0 +1,53 @@
+/*
+ * 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.ignite.math;
+
+import java.io.Externalizable;
+
+/**
+ * Data storage support for {@link Vector}.
+ */
+public interface VectorStorage extends Externalizable, StorageOpsMetrics, Destroyable {
+ /**
+ *
+ *
+ */
+ public int size();
+
+ /**
+ * @param i Vector element index.
+ * @return Value obtained for given element index.
+ */
+ public double get(int i);
+
+ /**
+ * @param i Vector element index.
+ * @param v Value to set at given index.
+ */
+ public void set(int i, double v);
+
+ /**
+ * Gets underlying array if {@link StorageOpsMetrics#isArrayBased()} returns {@code true}.
+ * Returns {@code null} if in other cases.
+ *
+ * @see StorageOpsMetrics#isArrayBased()
+ */
+ public default double[] data() {
+ return null;
+ }
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/decompositions/CholeskyDecomposition.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/CholeskyDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/CholeskyDecomposition.java
new file mode 100644
index 0000000..9554737
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/CholeskyDecomposition.java
@@ -0,0 +1,306 @@
+/*
+ * 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.ignite.math.decompositions;
+
+import org.apache.ignite.math.Matrix;
+import org.apache.ignite.math.Vector;
+import org.apache.ignite.math.exceptions.CardinalityException;
+import org.apache.ignite.math.exceptions.NonPositiveDefiniteMatrixException;
+import org.apache.ignite.math.exceptions.NonSymmetricMatrixException;
+
+/**
+ * Calculates the Cholesky decomposition of a matrix.
+ *
+ * This class inspired by class from Apache Common Math with similar name.
+ *
+ * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
+ * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
+ */
+public class CholeskyDecomposition extends DecompositionSupport {
+ /**
+ * Default threshold above which off-diagonal elements are considered too different
+ * and matrix not symmetric.
+ */
+ public static final double DFLT_REL_SYMMETRY_THRESHOLD = 1.0e-15;
+
+ /**
+ * Default threshold below which diagonal elements are considered null
+ * and matrix not positive definite.
+ */
+ public static final double DFLT_ABS_POSITIVITY_THRESHOLD = 1.0e-10;
+
+ /** Row-oriented storage for L<sup>T</sup> matrix data. */
+ private double[][] lTData;
+ /** Cached value of L. */
+ private Matrix cachedL;
+ /** Cached value of LT. */
+ private Matrix cachedLT;
+ /** Origin matrix */
+ private Matrix origin;
+
+ /**
+ * Calculates the Cholesky decomposition of the given matrix.
+ *
+ * Calling this constructor is equivalent to call {@link #CholeskyDecomposition(Matrix, double, double)} with the
+ * thresholds set to the default values {@link #DFLT_REL_SYMMETRY_THRESHOLD} and
+ * {@link #DFLT_ABS_POSITIVITY_THRESHOLD}.
+ *
+ * @param mtx the matrix to decompose.
+ * @throws CardinalityException if matrix is not square.
+ * @see #CholeskyDecomposition(Matrix, double, double)
+ * @see #DFLT_REL_SYMMETRY_THRESHOLD
+ * @see #DFLT_ABS_POSITIVITY_THRESHOLD
+ */
+ public CholeskyDecomposition(final Matrix mtx) {
+ this(mtx, DFLT_REL_SYMMETRY_THRESHOLD, DFLT_ABS_POSITIVITY_THRESHOLD);
+ }
+
+ /**
+ * Calculates the Cholesky decomposition of the given matrix.
+ *
+ * @param mtx the matrix to decompose.
+ * @param relSymmetryThreshold threshold above which off-diagonal elements are considered too different and matrix
+ * not symmetric
+ * @param absPositivityThreshold threshold below which diagonal elements are considered null and matrix not positive
+ * definite
+ * @see #CholeskyDecomposition(Matrix)
+ * @see #DFLT_REL_SYMMETRY_THRESHOLD
+ * @see #DFLT_ABS_POSITIVITY_THRESHOLD
+ */
+ public CholeskyDecomposition(final Matrix mtx, final double relSymmetryThreshold,
+ final double absPositivityThreshold) {
+ assert mtx != null;
+
+ if (mtx.columnSize() != mtx.rowSize())
+ throw new CardinalityException(mtx.rowSize(), mtx.columnSize());
+
+ origin = mtx;
+
+ final int order = mtx.rowSize();
+
+ lTData = toDoubleArr(mtx);
+ cachedL = null;
+ cachedLT = null;
+
+ // Check the matrix before transformation.
+ for (int i = 0; i < order; ++i) {
+ final double[] lI = lTData[i];
+
+ // Check off-diagonal elements (and reset them to 0).
+ for (int j = i + 1; j < order; ++j) {
+ final double[] lJ = lTData[j];
+
+ final double lIJ = lI[j];
+ final double lJI = lJ[i];
+
+ final double maxDelta = relSymmetryThreshold * Math.max(Math.abs(lIJ), Math.abs(lJI));
+
+ if (Math.abs(lIJ - lJI) > maxDelta)
+ throw new NonSymmetricMatrixException(i, j, relSymmetryThreshold);
+
+ lJ[i] = 0;
+ }
+ }
+
+ // Transform the matrix.
+ for (int i = 0; i < order; ++i) {
+ final double[] ltI = lTData[i];
+
+ // Check diagonal element.
+ if (ltI[i] <= absPositivityThreshold)
+ throw new NonPositiveDefiniteMatrixException(ltI[i], i, absPositivityThreshold);
+
+ ltI[i] = Math.sqrt(ltI[i]);
+ final double inverse = 1.0 / ltI[i];
+
+ for (int q = order - 1; q > i; --q) {
+ ltI[q] *= inverse;
+ final double[] ltQ = lTData[q];
+
+ for (int p = q; p < order; ++p)
+ ltQ[p] -= ltI[q] * ltI[p];
+ }
+ }
+ }
+
+ /** */
+ @Override public void destroy() {
+ if (cachedL != null)
+ cachedL.destroy();
+ if (cachedLT != null)
+ cachedLT.destroy();
+ }
+
+ /**
+ * Returns the matrix L of the decomposition.
+ * <p>L is an lower-triangular matrix</p>
+ *
+ * @return the L matrix
+ */
+ public Matrix getL() {
+ if (cachedL == null)
+ cachedL = getLT().transpose();
+
+ return cachedL;
+ }
+
+ /**
+ * Returns the transpose of the matrix L of the decomposition.
+ * <p>L<sup>T</sup> is an upper-triangular matrix</p>
+ *
+ * @return the transpose of the matrix L of the decomposition
+ */
+ public Matrix getLT() {
+
+ if (cachedLT == null) {
+ Matrix like = like(origin, origin.rowSize(), origin.columnSize());
+ like.assign(lTData);
+
+ cachedLT = like;
+ }
+
+ // return the cached matrix
+ return cachedLT;
+ }
+
+ /**
+ * Return the determinant of the matrix
+ *
+ * @return determinant of the matrix
+ */
+ public double getDeterminant() {
+ double determinant = 1.0;
+
+ for (int i = 0; i < lTData.length; ++i) {
+ double lTii = lTData[i][i];
+ determinant *= lTii * lTii;
+ }
+
+ return determinant;
+ }
+
+ /**
+ * Solve the linear equation A × X = B for matrices A.
+ *
+ * @param b right-hand side of the equation A × X = B
+ * @return a vector X that minimizes the two norm of A × X - B
+ * @throws CardinalityException if the vectors dimensions do not match
+ */
+ public Vector solve(final Vector b) {
+ final int m = lTData.length;
+
+ if (b.size() != m)
+ throw new CardinalityException(b.size(), m);
+
+ final double[] x = b.getStorage().data();
+
+ // Solve LY = b
+ for (int j = 0; j < m; j++) {
+ final double[] lJ = lTData[j];
+
+ x[j] /= lJ[j];
+
+ final double xJ = x[j];
+
+ for (int i = j + 1; i < m; i++)
+ x[i] -= xJ * lJ[i];
+ }
+
+ // Solve LTX = Y
+ for (int j = m - 1; j >= 0; j--) {
+ x[j] /= lTData[j][j];
+
+ final double xJ = x[j];
+
+ for (int i = 0; i < j; i++)
+ x[i] -= xJ * lTData[i][j];
+ }
+
+ return likeVector(origin, m).assign(x);
+ }
+
+ /**
+ * Solve the linear equation A × X = B for matrices A.
+ *
+ * @param b right-hand side of the equation A × X = B
+ * @return a matrix X that minimizes the two norm of A × X - B
+ * @throws CardinalityException if the matrices dimensions do not match
+ */
+ public Matrix solve(final Matrix b) {
+ final int m = lTData.length;
+
+ if (b.rowSize() != m)
+ throw new CardinalityException(b.rowSize(), m);
+
+ final int nColB = b.columnSize();
+ final double[][] x = b.getStorage().data();
+
+ // Solve LY = b
+ for (int j = 0; j < m; j++) {
+ final double[] lJ = lTData[j];
+ final double lJJ = lJ[j];
+ final double[] xJ = x[j];
+
+ for (int k = 0; k < nColB; ++k)
+ xJ[k] /= lJJ;
+
+ for (int i = j + 1; i < m; i++) {
+ final double[] xI = x[i];
+ final double lJI = lJ[i];
+
+ for (int k = 0; k < nColB; ++k)
+ xI[k] -= xJ[k] * lJI;
+ }
+ }
+
+ // Solve LTX = Y
+ for (int j = m - 1; j >= 0; j--) {
+ final double lJJ = lTData[j][j];
+ final double[] xJ = x[j];
+
+ for (int k = 0; k < nColB; ++k)
+ xJ[k] /= lJJ;
+
+ for (int i = 0; i < j; i++) {
+ final double[] xI = x[i];
+ final double lIJ = lTData[i][j];
+
+ for (int k = 0; k < nColB; ++k)
+ xI[k] -= xJ[k] * lIJ;
+ }
+ }
+
+ return like(origin, m, b.columnSize()).assign(x);
+ }
+
+ /** */
+ private double[][] toDoubleArr(Matrix mtx) {
+ if (mtx.isArrayBased())
+ return mtx.getStorage().data();
+
+ double[][] res = new double[mtx.rowSize()][];
+
+ for (int row = 0; row < mtx.rowSize(); row++) {
+ res[row] = new double[mtx.columnSize()];
+ for (int col = 0; col < mtx.columnSize(); col++)
+ res[row][col] = mtx.get(row, col);
+ }
+
+ return res;
+ }
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/decompositions/DecompositionSupport.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/DecompositionSupport.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/DecompositionSupport.java
new file mode 100644
index 0000000..2c76284
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/DecompositionSupport.java
@@ -0,0 +1,105 @@
+/*
+ * 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.ignite.math.decompositions;
+
+import org.apache.ignite.math.Destroyable;
+import org.apache.ignite.math.Matrix;
+import org.apache.ignite.math.Vector;
+import org.apache.ignite.math.impls.matrix.CacheMatrix;
+import org.apache.ignite.math.impls.matrix.DenseLocalOnHeapMatrix;
+import org.apache.ignite.math.impls.matrix.PivotedMatrixView;
+import org.apache.ignite.math.impls.matrix.RandomMatrix;
+import org.apache.ignite.math.impls.vector.DenseLocalOnHeapVector;
+
+/**
+ * Helper methods to support decomposition of matrix types having some functionality limited.
+ */
+public abstract class DecompositionSupport implements Destroyable {
+ /**
+ * Create the like matrix with read-only matrices support.
+ *
+ * @param matrix Matrix for like.
+ * @return Like matrix.
+ */
+ protected Matrix like(Matrix matrix) {
+ if (isCopyLikeSupport(matrix))
+ return new DenseLocalOnHeapMatrix(matrix.rowSize(), matrix.columnSize());
+ else
+ return matrix.like(matrix.rowSize(), matrix.columnSize());
+ }
+
+ /**
+ * Create the like matrix with specified size with read-only matrices support.
+ *
+ * @param matrix Matrix for like.
+ * @return Like matrix.
+ */
+ protected Matrix like(Matrix matrix, int rows, int cols) {
+ if (isCopyLikeSupport(matrix))
+ return new DenseLocalOnHeapMatrix(rows, cols);
+ else
+ return matrix.like(rows, cols);
+ }
+
+ /**
+ * Create the like vector with read-only matrices support.
+ *
+ * @param matrix Matrix for like.
+ * @param crd Cardinality of the vector.
+ * @return Like vector.
+ */
+ protected Vector likeVector(Matrix matrix, int crd) {
+ if (isCopyLikeSupport(matrix))
+ return new DenseLocalOnHeapVector(crd);
+ else
+ return matrix.likeVector(crd);
+ }
+
+ /**
+ * Create the like vector with read-only matrices support.
+ *
+ * @param matrix Matrix for like.
+ * @return Like vector.
+ */
+ protected Vector likeVector(Matrix matrix) {
+ return likeVector(matrix, matrix.rowSize());
+ }
+
+ /**
+ * Create the copy of matrix with read-only matrices support.
+ *
+ * @param matrix Matrix for copy.
+ * @return Copy.
+ */
+ protected Matrix copy(Matrix matrix) {
+ if (isCopyLikeSupport(matrix)) {
+ DenseLocalOnHeapMatrix cp = new DenseLocalOnHeapMatrix(matrix.rowSize(), matrix.columnSize());
+
+ cp.assign(matrix);
+
+ return cp;
+ }
+ else
+ return matrix.copy();
+ }
+
+ /** */
+ private boolean isCopyLikeSupport(Matrix matrix) {
+ return matrix instanceof RandomMatrix || matrix instanceof PivotedMatrixView || matrix instanceof CacheMatrix;
+ }
+}
http://git-wip-us.apache.org/repos/asf/ignite/blob/732dfea9/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java
----------------------------------------------------------------------
diff --git a/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java
new file mode 100644
index 0000000..66fe13c
--- /dev/null
+++ b/modules/ml/src/main/java/org/apache/ignite/math/decompositions/EigenDecomposition.java
@@ -0,0 +1,923 @@
+/*
+ * 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.ignite.math.decompositions;
+
+import org.apache.ignite.math.Matrix;
+import org.apache.ignite.math.Vector;
+import org.apache.ignite.math.functions.Functions;
+
+/**
+ * This class provides EigenDecomposition of given matrix. The class is based on
+ * class with similar name from <a href="http://mahout.apache.org/">Apache Mahout</a> library.
+ *
+ * @see <a href=http://mathworld.wolfram.com/EigenDecomposition.html>MathWorld</a>
+ */
+public class EigenDecomposition extends DecompositionSupport {
+ /** Row and column dimension (square matrix). */
+ private final int n;
+
+ /** Array for internal storage of eigen vectors. */
+ private final Matrix v;
+
+ /** Array for internal storage of eigenvalues. */
+ private final Vector d;
+ /** Array for internal storage of eigenvalues. */
+ private final Vector e;
+
+ /** */
+ public EigenDecomposition(Matrix matrix) {
+ this(matrix, isSymmetric(matrix));
+ }
+
+ /** */
+ public EigenDecomposition(Matrix matrix, boolean isSymmetric) {
+ n = matrix.columnSize();
+
+ d = likeVector(matrix);
+ e = likeVector(matrix);
+ v = like(matrix);
+
+ if (isSymmetric) {
+ v.assign(matrix);
+
+ // Tridiagonalize.
+ tred2();
+
+ // Diagonalize.
+ tql2();
+
+ }
+ else
+ // Reduce to Hessenberg form.
+ // Reduce Hessenberg to real Schur form.
+ hqr2(orthes(matrix));
+ }
+
+ /**
+ * Return the eigen vector matrix
+ *
+ * @return V
+ */
+ public Matrix getV() {
+ return like(v).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 res = like(v, d.size(), d.size());
+ res.assign(0);
+ res.viewDiagonal().assign(d);
+ for (int i = 0; i < n; i++) {
+ double v = e.getX(i);
+ if (v > 0)
+ res.setX(i, i + 1, v);
+ else if (v < 0)
+ res.setX(i, i - 1, v);
+ }
+ return res;
+ }
+
+ /**
+ * Destroys decomposition components and other internal components of decomposition.
+ */
+ @Override public void destroy() {
+ e.destroy();
+ v.destroy();
+ d.destroy();
+ }
+
+ /** */
+ 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).kNorm(1);
+ double h = 0.0;
+
+ if (scale == 0.0) {
+ e.setX(i, d.getX(i - 1));
+ for (int j = 0; j < i; j++) {
+ d.setX(j, v.getX(i - 1, j));
+ v.setX(i, j, 0.0);
+ v.setX(j, i, 0.0);
+ }
+ }
+ else {
+
+ // Generate Householder vector.
+
+ for (int k = 0; k < i; k++) {
+ d.setX(k, d.getX(k) / scale);
+ h += d.getX(k) * d.getX(k);
+ }
+
+ double f = d.getX(i - 1);
+ double g = Math.sqrt(h);
+
+ if (f > 0)
+ g = -g;
+
+ e.setX(i, scale * g);
+ h -= f * g;
+ d.setX(i - 1, f - g);
+
+ for (int j = 0; j < i; j++)
+ e.setX(j, 0.0);
+
+ // Apply similarity transformation to remaining columns.
+
+ for (int j = 0; j < i; j++) {
+ f = d.getX(j);
+ v.setX(j, i, f);
+ g = e.getX(j) + v.getX(j, j) * f;
+
+ for (int k = j + 1; k <= i - 1; k++) {
+ g += v.getX(k, j) * d.getX(k);
+ e.setX(k, e.getX(k) + v.getX(k, j) * f);
+ }
+
+ e.setX(j, g);
+ }
+
+ f = 0.0;
+
+ for (int j = 0; j < i; j++) {
+ e.setX(j, e.getX(j) / h);
+ f += e.getX(j) * d.getX(j);
+ }
+
+ double hh = f / (h + h);
+
+ for (int j = 0; j < i; j++)
+ e.setX(j, e.getX(j) - hh * d.getX(j));
+
+ for (int j = 0; j < i; j++) {
+ f = d.getX(j);
+ g = e.getX(j);
+
+ for (int k = j; k <= i - 1; k++)
+ v.setX(k, j, v.getX(k, j) - (f * e.getX(k) + g * d.getX(k)));
+
+ d.setX(j, v.getX(i - 1, j));
+ v.setX(i, j, 0.0);
+ }
+ }
+
+ d.setX(i, h);
+ }
+ }
+
+ /** */
+ private Matrix orthes(Matrix matrix) {
+ // Working storage for nonsymmetric algorithm.
+ Vector ort = likeVector(matrix);
+ Matrix hessenBerg = like(matrix).assign(matrix);
+
+ // 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 hCol = hessenBerg.viewColumn(m - 1).viewPart(m, high - m + 1);
+ double scale = hCol.kNorm(1);
+
+ if (scale != 0.0) {
+ // Compute Householder transformation.
+ ort.viewPart(m, high - m + 1).map(hCol, Functions.plusMult(1 / scale));
+ double h = ort.viewPart(m, high - m + 1).getLengthSquared();
+
+ double g = Math.sqrt(h);
+
+ if (ort.getX(m) > 0)
+ g = -g;
+
+ h -= ort.getX(m) * g;
+ ort.setX(m, ort.getX(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).map(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).map(ortPiece, Functions.plusMult(-f));
+ }
+
+ ort.setX(m, scale * ort.getX(m));
+ hessenBerg.setX(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.getX(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.getX(m) / hessenBerg.getX(m, m - 1);
+ v.viewColumn(j).viewPart(m, high - m + 1).map(ort.viewPart(m, high - m + 1), Functions.plusMult(g));
+ }
+ }
+ }
+
+ return hessenBerg;
+ }
+
+ /** 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.setX(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.getX(l)) + Math.abs(e.getX(l)));
+ int m = l;
+
+ while (m < n) {
+ if (Math.abs(e.getX(m)) <= eps * tst1)
+ break;
+
+ m++;
+ }
+
+ // If m == l, d.getX(l) is an eigenvalue,
+ // otherwise, iterate.
+
+ if (m > l) {
+ do {
+ // Compute implicit shift
+
+ double g = d.getX(l);
+ double p = (d.getX(l + 1) - g) / (2.0 * e.getX(l));
+ double r = Math.hypot(p, 1.0);
+
+ if (p < 0)
+ r = -r;
+
+ d.setX(l, e.getX(l) / (p + r));
+ d.setX(l + 1, e.getX(l) * (p + r));
+ double dl1 = d.getX(l + 1);
+ double h = g - d.getX(l);
+
+ for (int i = l + 2; i < n; i++)
+ d.setX(i, d.getX(i) - h);
+
+ f += h;
+
+ // Implicit QL transformation.
+
+ p = d.getX(m);
+ double c = 1.0;
+ double c2 = c;
+ double c3 = c;
+ double el1 = e.getX(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.getX(i);
+ h = c * p;
+ r = Math.hypot(p, e.getX(i));
+ e.setX(i + 1, s * r);
+ s = e.getX(i) / r;
+ c = p / r;
+ p = c * d.getX(i) - s * g;
+ d.setX(i + 1, h + s * (c * g + s * d.getX(i)));
+
+ // Accumulate transformation.
+
+ for (int k = 0; k < n; k++) {
+ h = v.getX(k, i + 1);
+ v.setX(k, i + 1, s * v.getX(k, i) + c * h);
+ v.setX(k, i, c * v.getX(k, i) - s * h);
+ }
+ }
+
+ p = -s * s2 * c3 * el1 * e.getX(l) / dl1;
+ e.setX(l, s * p);
+ d.setX(l, c * p);
+
+ // Check for convergence.
+
+ }
+ while (Math.abs(e.getX(l)) > eps * tst1);
+ }
+
+ d.setX(l, d.getX(l) + f);
+ e.setX(l, 0.0);
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+
+ for (int i = 0; i < n - 1; i++) {
+ int k = i;
+ double p = d.getX(i);
+
+ for (int j = i + 1; j < n; j++)
+ if (d.getX(j) > p) {
+ k = j;
+ p = d.getX(j);
+ }
+
+ if (k != i) {
+ d.setX(k, d.getX(i));
+ d.setX(i, p);
+
+ for (int j = 0; j < n; j++) {
+ p = v.getX(j, i);
+ v.setX(j, i, v.getX(j, k));
+ v.setX(j, k, p);
+ }
+ }
+ }
+ }
+
+ /** */
+ 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.foldMap(Functions.PLUS, Functions.ABS, 0.0);
+
+ // 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.getX(l - 1, l - 1)) + Math.abs(h.getX(l, l));
+
+ if (s == 0.0)
+ s = norm;
+
+ if (Math.abs(h.getX(l, l - 1)) < eps * s)
+ break;
+
+ l--;
+ }
+
+ // Check for convergence
+
+ if (l == n) {
+ // One root found
+ h.setX(n, n, h.getX(n, n) + exshift);
+ d.setX(n, h.getX(n, n));
+ e.setX(n, 0.0);
+ n--;
+ iter = 0;
+ }
+ else if (l == n - 1) {
+ // Two roots found
+ w = h.getX(n, n - 1) * h.getX(n - 1, n);
+ p = (h.getX(n - 1, n - 1) - h.getX(n, n)) / 2.0;
+ q = p * p + w;
+ z = Math.sqrt(Math.abs(q));
+ h.setX(n, n, h.getX(n, n) + exshift);
+ h.setX(n - 1, n - 1, h.getX(n - 1, n - 1) + exshift);
+ x = h.getX(n, n);
+
+ // Real pair
+ if (q >= 0) {
+ if (p >= 0)
+ z = p + z;
+ else
+ z = p - z;
+
+ d.setX(n - 1, x + z);
+ d.setX(n, d.getX(n - 1));
+
+ if (z != 0.0)
+ d.setX(n, x - w / z);
+
+ e.setX(n - 1, 0.0);
+ e.setX(n, 0.0);
+ x = h.getX(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.getX(n - 1, j);
+ h.setX(n - 1, j, q * z + p * h.getX(n, j));
+ h.setX(n, j, q * h.getX(n, j) - p * z);
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= n; i++) {
+ z = h.getX(i, n - 1);
+ h.setX(i, n - 1, q * z + p * h.getX(i, n));
+ h.setX(i, n, q * h.getX(i, n) - p * z);
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ z = v.getX(i, n - 1);
+ v.setX(i, n - 1, q * z + p * v.getX(i, n));
+ v.setX(i, n, q * v.getX(i, n) - p * z);
+ }
+
+ // Complex pair
+
+ }
+ else {
+ d.setX(n - 1, x + p);
+ d.setX(n, x + p);
+ e.setX(n - 1, z);
+ e.setX(n, -z);
+ }
+
+ n -= 2;
+ iter = 0;
+
+ // No convergence yet
+
+ }
+ else {
+ // Form shift
+ x = h.getX(n, n);
+ y = 0.0;
+ w = 0.0;
+
+ if (l < n) {
+ y = h.getX(n - 1, n - 1);
+ w = h.getX(n, n - 1) * h.getX(n - 1, n);
+ }
+
+ // Wilkinson's original ad hoc shift
+
+ if (iter == 10) {
+ exshift += x;
+
+ for (int i = low; i <= n; i++)
+ h.setX(i, i, x);
+
+ s = Math.abs(h.getX(n, n - 1)) + Math.abs(h.getX(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.setX(i, i, h.getX(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.getX(m, m);
+ r = x - z;
+ s = y - z;
+ p = (r * s - w) / h.getX(m + 1, m) + h.getX(m, m + 1);
+ q = h.getX(m + 1, m + 1) - z - r - s;
+ r = h.getX(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.getX(m - 1, m - 1)) + Math.abs(h.getX(m + 1, m + 1));
+ double threshold = eps * Math.abs(p) * (Math.abs(z) + hmag);
+
+ if (Math.abs(h.getX(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < threshold)
+ break;
+
+ m--;
+ }
+
+ for (int i = m + 2; i <= n; i++) {
+ h.setX(i, i - 2, 0.0);
+
+ if (i > m + 2)
+ h.setX(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.getX(k, k - 1);
+ q = h.getX(k + 1, k - 1);
+ r = notlast ? h.getX(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.setX(k, k - 1, -s * x);
+ else if (l != m)
+ h.setX(k, k - 1, -h.getX(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.getX(k, j) + q * h.getX(k + 1, j);
+
+ if (notlast) {
+ p += r * h.getX(k + 2, j);
+ h.setX(k + 2, j, h.getX(k + 2, j) - p * z);
+ }
+
+ h.setX(k, j, h.getX(k, j) - p * x);
+ h.setX(k + 1, j, h.getX(k + 1, j) - p * y);
+ }
+
+ // Column modification
+
+ for (int i = 0; i <= Math.min(n, k + 3); i++) {
+ p = x * h.getX(i, k) + y * h.getX(i, k + 1);
+
+ if (notlast) {
+ p += z * h.getX(i, k + 2);
+ h.setX(i, k + 2, h.getX(i, k + 2) - p * r);
+ }
+
+ h.setX(i, k, h.getX(i, k) - p);
+ h.setX(i, k + 1, h.getX(i, k + 1) - p * q);
+ }
+
+ // Accumulate transformations
+
+ for (int i = low; i <= high; i++) {
+ p = x * v.getX(i, k) + y * v.getX(i, k + 1);
+
+ if (notlast) {
+ p += z * v.getX(i, k + 2);
+ v.setX(i, k + 2, v.getX(i, k + 2) - p * r);
+ }
+
+ v.setX(i, k, v.getX(i, k) - p);
+ v.setX(i, k + 1, v.getX(i, k + 1) - p * q);
+ }
+ } // (s != 0)
+ } // k loop
+ } // check convergence
+ } // while (n >= low)
+
+ // Back substitute to find vectors of upper triangular form
+
+ if (norm == 0.0)
+ return;
+
+ for (n = nn - 1; n >= 0; n--) {
+ p = d.getX(n);
+ q = e.getX(n);
+
+ // Real vector
+
+ double t;
+
+ if (q == 0) {
+ int l = n;
+ h.setX(n, n, 1.0);
+
+ for (int i = n - 1; i >= 0; i--) {
+ w = h.getX(i, i) - p;
+ r = 0.0;
+
+ for (int j = l; j <= n; j++)
+ r += h.getX(i, j) * h.getX(j, n);
+
+ if (e.getX(i) < 0.0) {
+ z = w;
+ s = r;
+ }
+ else {
+ l = i;
+
+ if (e.getX(i) == 0.0) {
+ if (w == 0.0)
+ h.setX(i, n, -r / (eps * norm));
+ else
+ h.setX(i, n, -r / w);
+
+ // Solve real equations
+
+ }
+ else {
+ x = h.getX(i, i + 1);
+ y = h.getX(i + 1, i);
+ q = (d.getX(i) - p) * (d.getX(i) - p) + e.getX(i) * e.getX(i);
+ t = (x * s - z * r) / q;
+ h.setX(i, n, t);
+
+ if (Math.abs(x) > Math.abs(z))
+ h.setX(i + 1, n, (-r - w * t) / x);
+ else
+ h.setX(i + 1, n, (-s - y * t) / z);
+ }
+
+ // Overflow control
+
+ t = Math.abs(h.getX(i, n));
+
+ if (eps * t * t > 1) {
+ for (int j = i; j <= n; j++)
+ h.setX(j, n, h.getX(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.getX(n, n - 1)) > Math.abs(h.getX(n - 1, n))) {
+ h.setX(n - 1, n - 1, q / h.getX(n, n - 1));
+ h.setX(n - 1, n, -(h.getX(n, n) - p) / h.getX(n, n - 1));
+ }
+ else {
+ cdiv(0.0, -h.getX(n - 1, n), h.getX(n - 1, n - 1) - p, q);
+ h.setX(n - 1, n - 1, cdivr);
+ h.setX(n - 1, n, cdivi);
+ }
+
+ h.setX(n, n - 1, 0.0);
+ h.setX(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.getX(i, j) * h.getX(j, n - 1);
+ sa += h.getX(i, j) * h.getX(j, n);
+ }
+
+ w = h.getX(i, i) - p;
+
+ if (e.getX(i) < 0.0) {
+ z = w;
+ r = ra;
+ s = sa;
+ }
+ else {
+ l = i;
+
+ if (e.getX(i) == 0) {
+ cdiv(-ra, -sa, w, q);
+ h.setX(i, n - 1, cdivr);
+ h.setX(i, n, cdivi);
+ }
+ else {
+
+ // Solve complex equations
+
+ x = h.getX(i, i + 1);
+ y = h.getX(i + 1, i);
+
+ double vr = (d.getX(i) - p) * (d.getX(i) - p) + e.getX(i) * e.getX(i) - q * q;
+ double vi = (d.getX(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.setX(i, n - 1, cdivr);
+ h.setX(i, n, cdivi);
+
+ if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
+ h.setX(i + 1, n - 1, (-ra - w * h.getX(i, n - 1) + q * h.getX(i, n)) / x);
+ h.setX(i + 1, n, (-sa - w * h.getX(i, n) - q * h.getX(i, n - 1)) / x);
+ }
+ else {
+ cdiv(-r - y * h.getX(i, n - 1), -s - y * h.getX(i, n), z, q);
+
+ h.setX(i + 1, n - 1, cdivr);
+ h.setX(i + 1, n, cdivi);
+ }
+ }
+
+ // Overflow control
+
+ t = Math.max(Math.abs(h.getX(i, n - 1)), Math.abs(h.getX(i, n)));
+
+ if (eps * t * t > 1)
+ for (int j = i; j <= n; j++) {
+ h.setX(j, n - 1, h.getX(j, n - 1) / t);
+ h.setX(j, n, h.getX(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.setX(i, j, h.getX(i, j));
+ }
+
+ // Back transformation to get eigen vectors 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.getX(i, k) * h.getX(k, j);
+
+ v.setX(i, j, z);
+ }
+ }
+
+ /** */
+ private static boolean isSymmetric(Matrix matrix) {
+ int cols = matrix.columnSize();
+ int rows = matrix.rowSize();
+
+ if (cols != rows)
+ return false;
+
+ for (int i = 0; i < cols; i++)
+ for (int j = 0; j < rows; j++) {
+ if (matrix.getX(i, j) != matrix.get(j, i))
+ return false;
+ }
+
+ return true;
+ }
+
+ /** Complex scalar division - real part. */
+ private double cdivr;
+ /** Complex scalar division - imaginary part. */
+ 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;
+ }
+ }
+}