---

CompareHistogram.cpp

#include <LEDA/d_array.h>
#include <LEDA/graph_alg.h>
#include "SmartArray.h"
#include "Emd.h"
#include "CompareHistograms.h"

double CompareHistograms(const SmartArray<double> hist_values1, const SmartArray<double> weights1,
                                const SmartArray<double> hist_values2, const SmartArray<double> weights2){

 // variables to fill later.
 graph G[2];
 node root;
 node_array <double> node_weight[2];
 edge_array <double> edge_weight[2];
 double emd_return;
 Matrix embedded[2];
 int ncols[2];
 int cdim[2];
 node_array <int> column_map[2]; // a mapping of nodes into matrix columns.
 ColumnVector **gvec[2]; // column vectors corresponding to matrix columns.
 double *n_weight[2]; // normalized node weights (to sum up to 1).
 ColumnVector mu[2];
 double sigma[2];

         
 // read in two histograms.
 for( i = 0; i < 2; i++ ) {
         if (i==0)
                 embedded[i] = Get_Graph(G[i], node_weight[i], column_map[i], hist_values1, weights1);
         if (i==1)
                 embedded[i] = Get_Graph(G[i], node_weight[i], column_map[i], hist_values2, weights2);
         
         ncols[i] = embedded[i].Ncols();
         cdim[i] = embedded[i].Nrows();
            
         // copy weights into a column-weight array.
         double orig_weight[ncols[i]];
         node v;
         forall_nodes( v, G[i] ) {
                 column_map[i][v] -= 1;
                 orig_weight[column_map[i][v]] = node_weight[i][v];
         }

         // normalize column-weights.
         n_weight[i] = (double *) calloc( ncols[i], sizeof(double) );
         normalize_weights( orig_weight, ncols[i], n_weight[i] );

         // copy the matrix into an array of column vectors -- for "emd.c".
         gvec[i] = (ColumnVector **) calloc( ncols[i], sizeof(ColumnVector *) );
         for( int j=0; j<ncols[i]; j++ )
                 gvec[i][j] = new ColumnVector( embedded[i].Column(j+1) );

         // normalize column-vectors.
         normalize_column_vectors(
                         gvec[i], ncols[i], n_weight[i], mu[i], sigma[i] );
 }

        // int min_dim = DMIN( cdim[0], cdim[1] );
        int min_dim = 3;

        // set up structures for emd.
        signature_t s[2] = { {ncols[0], gvec[0], n_weight[0]},
                {ncols[1], gvec[1], n_weight[1]}};
        int flow_size;  
        flow_t flow[ncols[0]+ncols[1]-1];

        // set up structures for least squares.
        d_array <int, ColumnVector *> X;
        d_array <int, ColumnVector *> Y;
        d_array <int, double> pair_weight;

        // set up emd/lset iteration.
        double old_emd_value = 999999.9;
        double curr_emd_value = old_emd_value - 1;
        while( curr_emd_value < old_emd_value ) {
                // call emd.
                old_emd_value = curr_emd_value;
                curr_emd_value = emd( s, s+1, vector_dist, flow, &flow_size );
                emd_return = curr_emd_value;

                // set up LSET.
                X.clear(); Y.clear(); pair_weight.clear();
                for( int k=0; k<flow_size; k++ ) {
                        X[k] = gvec[0][flow[k].from];
                        Y[k] = gvec[1][flow[k].to];
                        pair_weight[k] = flow[k].amount;
                }

                // call least squares.
                double error = least_squares_estimation( X, Y, pair_weight )/flow_size;
        }
                
 return emd_return;     
}


/**************************************************************************************/

Matrix GetGraph(graph &G,node_array <double> &node_weight, node_array <int> &column_map,
                const SmartArray<double> hist_values, const SmartArray<double> weights){

 int number, i, node1, node2, root_number, mass_no=0, segment;
 node v, a, b; edge e;
 double mass, weight, x, y, z;
 char node_type;

 number = hist_values.GetSize();

 for (i=0; i<number; i++)      //Generating a graph with n nodes and 0 edges.
       G.new_node();


  double *Node_Mass = (double*)malloc(number * sizeof(double));  //Stores the Weight for each node
  if (Node_Mass == NULL){
          cerr << "Not enough memory !" << endl;
          exit(-1);
  }

  Matrix Vector_l2(3, G.number_of_nodes()); //There are only x,y,z coordinates

  for (i=0; i<number; i++){
          Node_Mass[i] = weights[i];
          Vector_l2(1,i+1) = hist_values[i];
          Vector_l2(2,i+1) = 0;
          Vector_l2(3,i+1) = 0;
  }

  node_array<double> node_mass;
  node_array<int> node_label;
  node_label.init(G);
  node_mass.init(G);
  node_weight.init(G);
  column_map.init(G);
  int in = 0;
  forall_nodes(v,G) {
    node_label[v] = in;
    node_mass[v] = Node_Mass[in];
    node_weight[v] = Node_Mass[in];
    column_map[v] = in + 1;
    in++;
  }

 return Vector_l2;
}

// distance defined on pairs of column vectors.
float vector_dist( ColumnVector **v1, ColumnVector **v2 ) {
        return (float) sqrt( (*(*v1) - *(*v2)).SumSquare() );
}

// normalizing weights to sum up to 1.
void normalize_weights( double weight[], int nv, double normed_weight[] ) {
        // compute sum of weights and normalize them.
        double weight_sum = 0;
        for( int i=0; i<nv; i++ )
                weight_sum += weight[i];
        for( int i=0; i<nv; i++ )
                normed_weight[i] = weight[i]/weight_sum;
}

//  "normalizing" a set of column vectors by finding their weighted average,
//  shifting column-vectors to make average 0, and dividing by the deviation.
void normalize_column_vectors( ColumnVector *mcol[], int nc,
                double normed_weight[], ColumnVector &mu, double &sigma ){
        int nr = (*mcol[0]).Nrows();

        // compute mu.
        mu = ColumnVector(nr);
        mu *= 0;
        for( int i=0; i<nc; i++ )
                mu += normed_weight[i]*(*mcol[i]);

        // subtract mu from everything.
        for( int i=0; i<nc; i++ )
                *mcol[i] -= mu;

        // compute sigma.
        sigma = 0;
        for( int i=0; i<nc; i++ )
                sigma += (double) normed_weight[i]*((*mcol[i]).SumSquare());
        sigma = sqrt(sigma);

        // normalize by sigma.
        for( int i=0; i<nc; i++ )
                (*mcol[i]) /= sigma;

#ifdef _DEBUG_
        cerr << "computed sigma: " << sigma << endl;
        cerr << "computed mu: " << mu.t();
#endif
}

//  Projecting X's onto Y's using linear least squares estimate.
double least_squares_estimation(
                d_array <int, ColumnVector *> &X, d_array <int, ColumnVector *> &Y,
                d_array <int, double> &normed_weight ) {

        int ncols = X.size();
        int nrows = X[0]->Nrows();

#ifdef _DEBUG_
        // compute initial error (debug).
        double init_error = 0;
        for( int i=0; i<ncols; i++ ) {
                init_error += (double) normed_weight[i]*sqrt( ((*X[i])-(*Y[i])).SumSquare() );
        }
        cerr << "initial error in LSET: " << init_error << endl;
#endif

        // compute Sigma_XY.
        Matrix Sigma_XY( nrows, nrows );
        Sigma_XY *= 0;
        for( int i=0; i<ncols; i++ ) {
                Sigma_XY += (normed_weight[i]*(*Y[i]))*(*X[i]).t();
        }

        // compute SVD.
        Matrix U(nrows, nrows), V(nrows, nrows);
        DiagonalMatrix D(nrows);
        SVD( Sigma_XY, D, U, V );

        // compute R.
        Matrix R = U * V.t();

        // update X's. Since they are represented by pointers,
        // need to update each pointer's value once - use a d_array to check it.
        d_array <long, ColumnVector *> ptr_map;
        for( int i=0; i<ncols; i++ ) {
                if( !ptr_map.defined( (long) X[i] ) ) {
                        *X[i] = R * (*X[i]);
                        ptr_map[(long) X[i]] = X[i];
                }
        }

        // return the error of the estimate.
        double error = 0;
        for( int i=0; i<ncols; i++ ) {
                error += (double) normed_weight[i]*sqrt( ((*X[i])-(*Y[i])).SumSquare() );
        }

#ifdef _DEBUG_
        cerr << "final error in LSET: " << error << endl;
#endif

        return error;
}

---