Return-Path: X-Original-To: archive-asf-public-internal@cust-asf2.ponee.io Delivered-To: archive-asf-public-internal@cust-asf2.ponee.io Received: from cust-asf.ponee.io (cust-asf.ponee.io [163.172.22.183]) by cust-asf2.ponee.io (Postfix) with ESMTP id D9836200C6B for ; Mon, 17 Apr 2017 18:01:48 +0200 (CEST) Received: by cust-asf.ponee.io (Postfix) id D82D7160B9C; Mon, 17 Apr 2017 16:01:48 +0000 (UTC) Delivered-To: archive-asf-public@cust-asf.ponee.io Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by cust-asf.ponee.io (Postfix) with SMTP id 5187A160BB7 for ; Mon, 17 Apr 2017 18:01:46 +0200 (CEST) Received: (qmail 64707 invoked by uid 500); 17 Apr 2017 16:01:45 -0000 Mailing-List: contact commits-help@ignite.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@ignite.apache.org Delivered-To: mailing list commits@ignite.apache.org Received: (qmail 64239 invoked by uid 99); 17 Apr 2017 16:01:44 -0000 Received: from git1-us-west.apache.org (HELO git1-us-west.apache.org) (140.211.11.23) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 17 Apr 2017 16:01:44 +0000 Received: by git1-us-west.apache.org (ASF Mail Server at git1-us-west.apache.org, from userid 33) id 7AC32E178A; Mon, 17 Apr 2017 16:01:43 +0000 (UTC) Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit From: av@apache.org To: commits@ignite.apache.org Date: Mon, 17 Apr 2017 16:01:53 -0000 Message-Id: In-Reply-To: <7bc1acfb896b40699a7b469ab9cc87c4@git.apache.org> References: <7bc1acfb896b40699a7b469ab9cc87c4@git.apache.org> X-Mailer: ASF-Git Admin Mailer Subject: [11/26] ignite git commit: IGNITE-5000 Rename Ignite Math module to Ignite ML module archived-at: Mon, 17 Apr 2017 16:01:49 -0000 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 { + } + + /** 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 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 Apache Mahout. + */ +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 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 nonZeroes(); + + /** + * Gets spliterator for all values in this vector. + * + * @return Spliterator for all values. + */ + public Spliterator allSpliterator(); + + /** + * Gets spliterator for all non-zero values in this vector. + * + * @return Spliterator for all non-zero values. + */ + public Spliterator 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 fun); + + /** + * Maps all values in this vector through a given function. + * + * For this vector A, argument vector B and the + * function F this method maps every element x as: + * A(x) = F(A(x), B(x)) + * + * @param vec Argument vector. + * @param fun Mapping function. + * @return This function. + * @throws CardinalityException Thrown if cardinalities mismatch. + */ + public Vector map(Vector vec, IgniteBiFunction fun); + + /** + * Maps all elements of this vector by applying given function to each element with a constant + * second parameter y. + * + * @param fun Mapping function. + * @param y Second parameter for mapping function. + * @return This vector. + */ + public Vector map(IgniteBiFunction 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 Type of the folded value. + * @param zeroVal Zero value for fold operation. + * @return Folded value of this vector. + */ + public T foldMap(IgniteBiFunction foldFun, IgniteDoubleFunction 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 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 foldMap(Vector vec, IgniteBiFunction foldFun, IgniteBiFunction 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 extends KeyMapper { + /** + * @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 MathWorld + * @see Wikipedia + */ +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 LT 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. + *

L is an lower-triangular matrix

+ * + * @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. + *

LT is an upper-triangular matrix

+ * + * @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 Apache Mahout library. + * + * @see MathWorld + */ +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; + } + } +}