/*******************************************************************************
* Copyright 2007, 2009 Stephen O'Rourke (stephen.orourke@sydney.edu.au)
*
* Licensed 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 tml.vectorspace.factorisation;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Random;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import tml.utils.DistanceLib;
import tml.utils.DistanceLib.DistanceMeasure;
import weka.core.Attribute;
import weka.core.FastVector;
import weka.core.Instance;
import weka.core.Instances;
import Jama.Matrix;
/**
* This class converts instances for plotting using Multidimensional Scaling. It
* use a Newton-Raphson algorithm to project instances into 2 dimensions.
*
* Details of this algorithm can be found in the paper:
* http://www.pavis.org/essay/multidimensional_scaling.html
*
* @author Stephen O'Rourke
*/
public class MultiDimensionalScalingNR {
public static final int X = 0; // attribute Y
public static final int Y = 1; // attribute X
private static final int p = 2; // number of dimensions
private double tolerence = 0.01;
private int maxIterations = 1000; // maximum iterations
private double error;
private Matrix d;
private Matrix d_hat;
private Instances initialX;
private DistanceMeasure lowDimensionalDistanceMeasure = DistanceMeasure.EUCLIDEAN;
private DistanceMeasure highDimensionlDistanceMeasure = DistanceMeasure.COSINE;
private final Log logger = LogFactory.getLog(getClass());
public Instances scale(Instances instances) {
// approximation error
error = 0.0;
double error_previous;
// number of points
int n = instances.numInstances();
// distance between points in the p-dimensional layout
d = new Matrix(n, n);
Matrix d_previous;
// dissimilarity between vectors
d_hat = new Matrix(n, n);
// points instances
FastVector attributes = new FastVector(p);
attributes.addElement(new Attribute("X"));
attributes.addElement(new Attribute("Y"));
Instances x = new Instances("MDS", attributes, instances.numInstances());
Instances x_previous;
// initialise points sequence
ArrayList<Integer> kseq = new ArrayList<Integer>();
for (int k = 0; k < n; k++) {
kseq.add(k);
}
// initialise x
if (initialX != null) {
x = new Instances(initialX);
} else {
Random rand = new Random();
for (int k = 0; k < n; k++) {
Instance x_inst = new Instance(p);
x_inst.setValue(X, rand.nextDouble() - rand.nextInt(1));
x_inst.setValue(Y, rand.nextDouble() - rand.nextInt(1));
x.add(x_inst);
}
}
// calculate d
for (int j = 0; j < n; j++) {
for (int i = 0; i < j; i++) {
double distance = this.distance(x.instance(i), x.instance(j));
d.set(i, j, distance);
d.set(j, i, distance);
double dissimilarity = this.dissimilarity(instances.instance(i), instances.instance(j));
d_hat.set(i, j, dissimilarity);
d_hat.set(j, i, dissimilarity);
if (d_hat.get(i, j) != 0) {
error += Math.pow(d.get(i, j) - d_hat.get(i, j), 2) / Math.pow(d_hat.get(i, j), 2);
}
}
}
// record previous results
error_previous = error;
d_previous = d.copy();
x_previous = new Instances(x);
// start of Newton-Raphson method
logger.info("Starting Newton-Raphson MDS.");
for (int iter = 0; iter < maxIterations; iter++) {
// randomise points sequence to ensure faster convergence
Collections.shuffle(kseq);
for (int k : kseq) {
Matrix gradient = new Matrix(p, 1);
Matrix hessian = new Matrix(p, p);
// calculate gradient vector
for (int a = 0; a < p; a++) {
double sum = 0;
for (int l = 0; l < n; l++) {
if (k != l) {
if (d.get(k, l) != 0 && d_hat.get(k, l) != 0) {
sum += ((d.get(k, l) - d_hat.get(k, l)) / (d.get(k, l) * Math.pow(d_hat.get(k, l), 2))) * (x.instance(k).value(a) - x.instance(l).value(a));
}
}
}
gradient.set(a, 0, 2 * sum);
}
// calculate hessian matrix
for (int a = 0; a < p; a++) {
for (int b = 0; b < p; b++) {
double sum = 0.0;
if (a != b) {
for (int l = 0; l < n; l++) {
if (k != l) {
if (d.get(k, l) != 0 && d_hat.get(k, l) != 0) {
sum += ((x.instance(k).value(a) - x.instance(l).value(a)) * (x.instance(k).value(b) - x.instance(l).value(b))) / (Math.pow(d.get(k, l), 3) * d_hat.get(k, l));
}
}
}
sum = 2 * sum;
} else {
for (int l = 0; l < n; l++) {
if (k != l) {
if (d_hat.get(k, l) != 0 && d.get(k, l) != 0) {
sum += (1.0 / Math.pow(d_hat.get(k, l), 2)) - (Math.pow(d.get(k, l), 2) - Math.pow((x.instance(k).value(a) - x.instance(l).value(a)), 2)) / (Math.pow(d.get(k, l), 3) * d_hat.get(k, l));
}
}
}
sum = 2 * sum;
}
hessian.set(a, b, sum);
}
}
// update x
Matrix x_k = new Matrix(x.instance(k).toDoubleArray(), p);
Matrix x_k_tilda = x_k.minus(hessian.inverse().times(gradient));
x.instance(k).setValue(X, x_k_tilda.get(X, 0));
x.instance(k).setValue(Y, x_k_tilda.get(Y, 0));
}
// calculate d and error
error = 0;
for (int j = 0; j < n; j++) {
for (int i = 0; i < j; i++) {
double distance = this.distance(x.instance(i), x.instance(j));
d.set(i, j, distance);
d.set(j, i, distance);
if (d_hat.get(i, j) != 0) {
error += Math.pow(d.get(i, j) - d_hat.get(i, j), 2) / Math.pow(d_hat.get(i, j), 2);
}
}
}
if (error < error_previous) {
logger.debug(iter + ".\t error " + error);
if (error_previous - error <= tolerence) {
break;
}
error_previous = error;
d_previous = d.copy();
x_previous = new Instances(x);
} else // invalidates last run
{
x = new Instances(x_previous);
d = d_previous.copy();
}
}
logger.info("Finished Newton-Raphson MDS.");
return x;
}
protected double distance(Instance inst1, Instance inst2) {
double distance = DistanceLib.distance(lowDimensionalDistanceMeasure, inst1, inst2);
return distance;
}
protected double dissimilarity(Instance inst1, Instance inst2) {
double distance = Math.sqrt(1 - DistanceLib.distance(highDimensionlDistanceMeasure, inst1, inst2));
return distance;
}
public double getTolerence() {
return tolerence;
}
public void setTolerence(double tolerence) {
this.tolerence = tolerence;
}
public int getMaxIterations() {
return maxIterations;
}
public void setMaxIterations(int maxIterations) {
this.maxIterations = maxIterations;
}
public DistanceMeasure getLowDimensionalDistanceMeasure() {
return lowDimensionalDistanceMeasure;
}
public void setLowDimensionalDistanceMeasure(DistanceMeasure lowDimensionalDistanceMeasure) {
this.lowDimensionalDistanceMeasure = lowDimensionalDistanceMeasure;
}
public DistanceMeasure getHighDimensionlDistanceMeasure() {
return highDimensionlDistanceMeasure;
}
public void setHighDimensionlDistanceMeasure(DistanceMeasure highDimensionlDistanceMeasure) {
this.highDimensionlDistanceMeasure = highDimensionlDistanceMeasure;
}
public Instances getInitialX() {
return initialX;
}
public void setInitialX(Instances initialX) {
this.initialX = initialX;
}
public Matrix d() {
return d;
}
public Matrix d_hat() {
return d_hat;
}
public double error() {
return error;
}
}