"""
.. index:: multidimensional scaling (mds)

.. index::
   single: projection; multidimensional scaling (mds)

**********************************
Multidimensional scaling (``mds``)
**********************************

The functionality to perform multidimensional scaling
(http://en.wikipedia.org/wiki/Multidimensional_scaling).

The main class to perform multidimensional scaling is
:class:`Orange.projection.mds.MDS`

.. autoclass:: Orange.projection.mds.MDS
   :members:
   :exclude-members: Torgerson, get_distance, get_stress, calc_stress, run

   .. automethod:: calc_stress(stress_func=SgnRelStress)
   .. automethod:: run(iter, stress_func=SgnRelStress, eps=1e-3, progress_callback=None)

Stress functions
================

Stress functions that can be used for MDS have to be implemented as functions
or callable classes:

    .. method:: \ __call__(correct, current, weight=1.0)
       
       Compute the stress using the correct and the current distance value (the
       :obj:`Orange.projection.mds.MDS.distances` and
       :obj:`Orange.projection.mds.MDS.projected_distances` elements).
       
       :param correct: correct (actual) distance between elements, represented by
           the two points.
       :type correct: float
       
       :param current: current distance between the points in the MDS space.
       :type current: float

This module provides the following stress functions:

   * :obj:`SgnRelStress`
   * :obj:`KruskalStress`
   * :obj:`SammonStress`
   * :obj:`SgnSammonStress`

Examples
========

MDS Scatterplot
---------------

The following script computes the Euclidean distance between the data
instances and runs MDS. Final coordinates are plotted with matplotlib
(not included with orange, http://matplotlib.sourceforge.net/).

Example (:download:`mds-scatterplot.py <code/mds-scatterplot.py>`)

.. literalinclude:: code/mds-scatterplot.py
    :lines: 7-

The script produces a file *mds-scatterplot.py.png*. Color denotes
the class. Iris is a relatively simple data set with respect to
classification; to no surprise we see that MDS finds such instance
placement in 2D where instances of different classes are well separated.
Note that MDS has no knowledge of points' classes.

.. image:: files/mds-scatterplot.png


A more advanced example
-----------------------

The following script performs 10 steps of Smacof optimization before computing
the stress. This is suitable if you have a large dataset and want to save some
time.

Example (:download:`mds-advanced.py <code/mds-advanced.py>`)

.. literalinclude:: code/mds-advanced.py
    :lines: 7-

A few representative lines of the output are::

    <-0.633911848068, 0.112218663096> [5.1, 3.5, 1.4, 0.2, 'Iris-setosa']
    <-0.624193906784, -0.111143872142> [4.9, 3.0, 1.4, 0.2, 'Iris-setosa']
    ...
    <0.265250980854, 0.237793982029> [7.0, 3.2, 4.7, 1.4, 'Iris-versicolor']
    <0.208580598235, 0.116296850145> [6.4, 3.2, 4.5, 1.5, 'Iris-versicolor']
    ...
    <0.635814905167, 0.238721415401> [6.3, 3.3, 6.0, 2.5, 'Iris-virginica']
    <0.356859534979, -0.175976261497> [5.8, 2.7, 5.1, 1.9, 'Iris-virginica']
    ...


"""

import numpy
from numpy.linalg import svd

import Orange.core
from Orange import orangeom as orangemds
from Orange.utils import deprecated_keywords
from Orange.utils import deprecated_members

KruskalStress = orangemds.KruskalStress()
SammonStress = orangemds.SammonStress()
SgnSammonStress = orangemds.SgnSammonStress()
SgnRelStress = orangemds.SgnRelStress()

PointList = Orange.core.FloatListList
FloatListList = Orange.core.FloatListList

def _mycompare((a,aa),(b,bb)):
    if a == b:
        return 0
    if a < b:
        return -1
    else:
        return 1
            
class PivotMDS(object):
    def __init__(self, distances=None, pivots=50, dim=2, **kwargs):
        self.dst = numpy.array([m for m in distances])
        self.n = len(self.dst)

        if type(pivots) == type(1):
            self.k = pivots
            self.pivots = numpy.random.permutation(len(self.dst))[:pivots]
            #self.pivots.sort()
        elif type(pivots) == type([]):
            self.pivots = pivots
            #self.pivots.sort()
            self.k = len(self.pivots)
        else:
            raise AttributeError('pivots')
        
    def optimize(self):
#        # Classical MDS (Torgerson)
#        J = identity(self.n) - (1/float(self.n))
#        B = -1/2. * dot(dot(J, self.dst**2), J)
#        w,v = linalg.eig(B)
#        tmp = zip([float(val) for val in w], range(self.n))
#        tmp.sort()
#        w1, w2 = tmp[-1][0], tmp[-2][0]
#        v1, v2 = v[:, tmp[-1][1]], v[:, tmp[-2][1]]
#        return v1 * sqrt(w1), v2 * sqrt(w2) 
        
        # Pivot MDS
        d = self.dst[[self.pivots]].T
        C = d**2
        # double-center d
        cavg = numpy.sum(d, axis=0)/(self.k+0.0)      # column sum
        ravg = numpy.sum(d, axis=1)/(self.n+0.0)    # row sum
        tavg = numpy.sum(cavg)/(self.n+0.0)   # total sum
        # TODO: optimize
        for i in xrange(self.n):
            for j in xrange(self.k):
                C[i,j] += -ravg[i] - cavg[j]
        
        C = -0.5 * (C + tavg)
        w,v = numpy.linalg.eig(numpy.dot(C.T, C))
        tmp = zip([float(val) for val in w], range(self.n))
        tmp.sort()
        w1, w2 = tmp[-1][0], tmp[-2][0]
        v1, v2 = v[:, tmp[-1][1]], v[:, tmp[-2][1]]
        x = numpy.dot(C, v1)
        y = numpy.dot(C, v2)
        return x, y


@deprecated_members(
    {"projectedDistances": "projected_distances",
     "originalDistances": "original_distances",
     "avgStress": "avg_stress",
     "progressCallback": "progress_callback",
     "getStress": "calc_stress",
     "get_stress": "calc_stress",
     "calcStress": "calc_stress",
     "getDistance": "calc_distance",
     "get_distance": "calc_distance",
     "calcDistance": "calc_distance",
     "Torgerson": "torgerson",
     "SMACOFstep": "smacof_step",
     "LSMT": "lsmt"})
class MDS(object):
    """
    Main class for performing multidimensional scaling.
    
    :param distances: original dissimilarity - a distance matrix to operate on.
    :type distances: :class:`Orange.misc.SymMatrix`
    
    :param dim: dimension of the projected space.
    :type dim: int
    
    :param points: an initial configuration of points (optional)
    :type points: :class:`Orange.core.FloatListList`
    
    An instance of MDS object has the following attributes and functions:
    
    .. attribute:: points
       
       Holds the current configuration of projected points in an
       :class:`Orange.core.FloatListList` object.
       
    .. attribute:: distances
    
       An :class:`Orange.misc.SymMatrix` containing the distances that we
       want to achieve (lsmt changes these).
       
    .. attribute:: projected_distances

       An :class:`Orange.misc.SymMatrix` containing the distances between
       projected points.
       
    .. attribute:: original_distances

       An :class:`Orange.misc.SymMatrix` containing the original distances
       between points.
       
    .. attribute:: stress
       
       An :class:`Orange.misc.SymMatrix` holding the stress.
    
    .. attribute:: dim

       An integer holding the dimension of the projected space.
       
    .. attribute:: n

       An integer holding the number of elements (points).
       
    .. attribute:: avg_stress

       A float holding the average stress in the :obj:`stress` matrix.
       
    .. attribute:: progress_callback

       A function that gets called after each optimization step in the
       :func:`run` method.
    
    """
    
    def __init__(self, distances=None, dim=2, **kwargs):
        self.mds=orangemds.MDS(distances, dim, **kwargs)
        self.original_distances=Orange.misc.SymMatrix([m for m in self.distances])

    def __getattr__(self, name):
        if name in ["points", "projected_distances", "distances" ,"stress",
                    "progress_callback", "n", "dim", "avg_stress"]:
            #print "rec:",name            
            return self.__dict__["mds"].__dict__[name]
        else:
            raise AttributeError(name)

    def __setattr__(self, name, value):
        #print "setattr"
        if name=="points":
            for i in range(len(value)):
                for j in range(len(value[i])):
                    self.mds.points[i][j]=value[i][j]
            return
            
        if name in ["projected_distances", "distances" ,"stress",
                    "progress_callback"]:
            self.mds.__setattr__(name, value)
        else:
            self.__dict__[name]=value
            
    def __nonzero__(self):
        return True
            
    def smacof_step(self):
        """
        Perform a single iteration of a Smacof algorithm that optimizes
        :obj:`stress` and updates the :obj:`points`.
        """
        self.mds.SMACOFstep()

    def calc_distance(self):
        """
        Compute the distances between points and update the
        :obj:`projected_distances` matrix.
        
        """
        self.mds.get_distance()
        

    @deprecated_keywords({"stressFunc": "stress_func"})
    def calc_stress(self, stress_func=SgnRelStress):
        """
        Compute the stress between the current :obj:`projected_distances` and
        :obj:`distances` matrix using *stress_func* and update the
        :obj:`stress` matrix and :obj:`avgStress` accordingly.
        
        """
        self.mds.getStress(stress_func)

    @deprecated_keywords({"stressFunc": "stress_func"})
    def optimize(self, iter, stress_func=SgnRelStress, eps=1e-3,
                 progress_callback=None):
        self.mds.progress_callback=progress_callback
        self.mds.optimize(iter, stress_func, eps)

    @deprecated_keywords({"stressFunc": "stress_func"})
    def run(self, iter, stress_func=SgnRelStress, eps=1e-3,
            progress_callback=None):
        """
        Perform optimization until stopping conditions are met.
        Stopping conditions are:
           
           * optimization runs for *iter* iterations of smacof_step function, or
           * stress improvement (old stress minus new stress) is smaller than
             eps * old stress.
        
        :param iter: maximum number of optimization iterations.
        :type iter: int
        
        :param stress_func: stress function.
        """
        self.optimize(iter, stress_func, eps, progress_callback)

    def torgerson(self):
        """
        Run the Torgerson algorithm that computes an initial analytical
        solution of the problem.
        
        """
        # Torgerson's initial approximation
        O = numpy.array([m for m in self.distances])
        
##        #B = matrixmultiply(O,O)
##        # bug!? B = O**2
##        B = dot(O,O)
##        # double-center B
##        cavg = sum(B, axis=0)/(self.n+0.0)      # column sum
##        ravg = sum(B, axis=1)/(self.n+0.0)    # row sum
##        tavg = sum(cavg)/(self.n+0.0)   # total sum
##        # B[row][column]
##        for i in xrange(self.n):
##            for j in xrange(self.n):
##                B[i,j] += -cavg[j]-ravg[i]
##        B = -0.5*(B+tavg)

        # B = double-center O**2 !!!
        J = numpy.identity(self.n) - (1/numpy.float(self.n))
        B = -0.5 * numpy.dot(numpy.dot(J, O**2), J)
        
        # SVD-solve B = ULU'
        #(U,L,V) = singular_value_decomposition(B)
        (U,L,V)=svd(B)
        # X = U(L^0.5)
        # # self.X = matrixmultiply(U,identity(self.n)*sqrt(L))
        # X is n-dimensional, we take the two dimensions with the largest singular values
        idx = numpy.argsort(L)[-self.dim:].tolist()
        idx.reverse()
        
        Lt = numpy.take(L,idx)   # take those singular values
        Ut = numpy.take(U,idx,axis=1) # take those columns that are enabled
        Dt = numpy.identity(self.dim)*numpy.sqrt(Lt)  # make a diagonal matrix, with squarooted values
        self.points = Orange.core.FloatListList(numpy.dot(Ut,Dt))
        self.freshD = 0
        
#        D = identity(self.n)*sqrt(L)  # make a diagonal matrix, with squarooted values
#        X = matrixmultiply(U,D)
#        self.X = take(X,idx,1)
    
    # Kruskal's monotone transformation
    def lsmt(self):
        """
        Execute Kruskal monotone transformation.
        
        """
        # optimize the distance transformation
        # build vector o
        effect = 0
        self.getDistance()
        o = []
        for i in xrange(1,self.n):
            for j in xrange(i):
                o.append((self.original_distances[i,j],(i,j)))
        o.sort(_mycompare)
        # find the ties in o, and construct the d vector sorting in order within ties
        d = []
        td = []
        uv = [] # numbers of consecutively tied o values
        (i,j) = o[0][1]
        distnorm = self.projected_distances[i,j]*self.projected_distances[i,j]
        td = [self.projected_distances[i,j]] # fetch distance
        for l in xrange(1,len(o)):
            # copy now sorted distances in an array
            # but sort distances within a tied o
            (i,j) = o[l][1]
            cd = self.projected_distances[i,j]
            distnorm += self.projected_distances[i,j]*self.projected_distances[i,j]
            if o[l][0] != o[l-1][0]:
                # differing value, flush
                sum = reduce(lambda x,y:x+y,td)+0.0
                d.append([sum,len(td),sum/len(td),td])
                td = []
            td.append(cd)
        sum = reduce(lambda x,y:x+y,td)+0.0
        d.append([sum,len(td),sum/len(td),td])
        ####
        # keep merging non-monotonous areas in d
        monotony = 0
        while not monotony and len(d) > 1:
            monotony = 1
            pi = 0 # index
            n = 1  # n-areas
            nd = []
            r = d[0] # current area
            for i in range(1,len(d)):
                tr = d[i]
                if r[2]>=tr[2]:
                    monotony = 0
                    effect = 1
                    r[0] += tr[0]
                    r[1] += tr[1]
                    r[2] = tr[0]/tr[1]
                    r[3] += tr[3]
                else:
                    nd.append(r)
                    r = tr
            nd.append(r)
            d = nd
        # normalizing multiplier
        sum = 0.0
        for i in d:
            sum += i[2]*i[2]*i[1]
        f = numpy.sqrt(distnorm/numpy.max(sum,1e-6))
        # transform O
        k = 0
        for i in d:
            for j in range(i[1]):
                (ii,jj) = o[k][1]
                self.distances[ii,jj] = f*i[2]
                k += 1
        assert(len(o) == k)
        self.freshD = 0
        return effect
