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