001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.ml.clustering.spectral;
031
032
033import java.util.Iterator;
034
035import org.openimaj.math.matrix.DiagonalMatrix;
036import org.openimaj.math.matrix.MatlibMatrixUtils;
037import org.openimaj.util.pair.DoubleObjectPair;
038
039import ch.akuhn.matrix.SparseMatrix;
040import ch.akuhn.matrix.Vector;
041import ch.akuhn.matrix.eigenvalues.Eigenvalues;
042
043/**
044 * Functions which turn a graph weight adjacency matrix into the Laplacian
045 * matrix.
046 * @author Sina Samangooei (ss@ecs.soton.ac.uk)
047 *
048 */
049public abstract class GraphLaplacian{
050
051
052        /**
053         * @param adj the adjanceny matrix should be square and symmetric
054         * @return the laplacian
055         */
056        public SparseMatrix laplacian(SparseMatrix adj){
057                DiagonalMatrix degree = new DiagonalMatrix(adj.rowCount());
058
059                int i = 0;
060                for (Vector row : adj.rows()) {
061//                      row.put(i, 0);
062                        degree.put(i, i, row.sum());
063                        i++;
064                }
065
066                return laplacian(adj,degree);
067        }
068
069        /**
070         * @param adj square and symmetric
071         * @param degree the sum of the adjacency for a node in the diagonals
072         * @return the laplacian
073         */
074        public abstract SparseMatrix laplacian(SparseMatrix adj, DiagonalMatrix degree);
075
076        /**
077         * @param evd
078         * @return provides an iterator over the (presumeably sorted)
079         */
080        public Iterator<DoubleObjectPair<Vector>> eigenIterator(Eigenvalues evd) {
081                return new FBEigenIterator(evd);
082        }
083        
084        /**
085         * The symmetric normalised Laplacian is defined as:
086         * L = D - W
087         * @author Sina Samangooei (ss@ecs.soton.ac.uk)
088         *
089         */
090        public static class Unnormalised extends GraphLaplacian{
091                @Override
092                public SparseMatrix laplacian(SparseMatrix adj, DiagonalMatrix degree) {
093                        SparseMatrix ret = MatlibMatrixUtils.plusInplace(
094                                DiagonalMatrix.ones(degree.rowCount()), 
095                                MatlibMatrixUtils.minusInplace(
096                                        degree,
097                                        adj
098                                )
099                        );
100                        return ret;
101                }
102        }
103
104        /**
105         * The inverted symmetric normalised Laplacian is defined as:
106         * L = D^-1/2 A D^-1/2
107         * @author Sina Samangooei (ss@ecs.soton.ac.uk)
108         *
109         */
110        public static class Normalised extends GraphLaplacian{
111
112                @Override
113                public SparseMatrix laplacian(SparseMatrix adj, DiagonalMatrix degree) {
114                        DiagonalMatrix invSqrtDegree = MatlibMatrixUtils.powInplace(degree,-1./2.);
115                        for (int i = 0; i < degree.rowCount(); i++) {
116                                if(Double.isNaN(degree.get(i, i)) || Double.isInfinite(degree.get(i, i)) ) 
117                                        invSqrtDegree.put(i, i, 0); 
118                        }
119                        SparseMatrix ret = MatlibMatrixUtils.times(
120                                MatlibMatrixUtils.times(
121                                        invSqrtDegree, adj
122                                ),
123                                invSqrtDegree
124                        );
125                        return ret;
126                }
127
128                @Override
129                public Iterator<DoubleObjectPair<Vector>> eigenIterator(Eigenvalues evd) {
130                        return new FBEigenIterator(evd);
131                }
132
133        }
134        
135        /**
136         * The inverted symmetric normalised Laplacian is defined as:
137         * L = D^-1 . W
138         * @author Sina Samangooei (ss@ecs.soton.ac.uk)
139         *
140         */
141        public static class Warped extends GraphLaplacian{
142
143                @Override
144                public SparseMatrix laplacian(SparseMatrix adj, DiagonalMatrix degree) {
145                        DiagonalMatrix invDegree = MatlibMatrixUtils.powInplace(degree,-1.);
146                        for (int i = 0; i < degree.rowCount(); i++) {
147                                if(Double.isNaN(degree.get(i, i)) || Double.isInfinite(degree.get(i, i)) ) 
148                                        invDegree.put(i, i, 0); 
149                        }
150                        SparseMatrix ret = MatlibMatrixUtils.times(
151                                MatlibMatrixUtils.times(
152                                        invDegree, adj
153                                ),
154                                invDegree
155                        );
156                        return ret;
157                }
158
159                @Override
160                public Iterator<DoubleObjectPair<Vector>> eigenIterator(Eigenvalues evd) {
161                        return new FBEigenIterator(evd);
162                }
163
164        }
165
166}