Jensen-Shannon divergence

//*****************************************************************************
// JSDivergence
// Created 20110703
// 
// Shawn Driscoll
// Salk Institute for Biological Studies
//
// Calculates the Jensen-Shannon divergence between the probability 
// distributions P and Q.  P and Q must span the same range of values
// and contain the same number of points.
//
// Jensen-Shannon Divergence is defined as
//
//   JSD(P||Q) = 1/2*D(P||M) + 1/2*D(Q||M) 
// 
// where M is 1/2*(P+Q) and D(P||M) is the Kullback-Leibler divergence.
//
// According to Lin (1991), the Jensen–Shannon divergence is bounded by 1.0.
//
//   0 <= JSD(P||Q) <= 1.0
// 
//*****************************************************************************
function JSDivergence(P,Q)
    //--------------------------------------------------------------------------
    // parameters
    //--------------------------------------------------------------------------
    wave P,Q;
 
    //--------------------------------------------------------------------------
    // variables
    //--------------------------------------------------------------------------
    variable n_jsd = 0;
    
    //--------------------------------------------------------------------------
    // main function
    //--------------------------------------------------------------------------
    
    ////
    // create mean distribution
    matrixop/free M = (P+Q)/2;
    
    ////
    // calculate squared jensen-shannon divergence
    n_jsd = 0.5*KLDivergence(P,M) + 0.5*KLDivergence(Q,M);
    
    ////
    // take square root since this is the actual metric
    n_jsd = sqrt(n_jsd);
    
    //
    // save to global variable
    variable/g v_JSDivergence = n_jsd;
    
    //
    // print result
    printf "v_JSDivergence= %g\r",n_jsd;
    
    // 
    // return result
    return n_jsd;
 
end