/*-----------------------------------------------------------------------------

Copyright (C) 2020

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-------------------------------------------------------------------------------

Function      ball_rand - uniform on a ball psuedo random number generator.
              ball_volume - volume of a d-ball with specified radius.
              log_ball_volume - log volume of a d-ball with specified radius.

Syntax        #include "libscl.h"
              realmat ball_rand(INT_32BIT* ix, INTEGER d, REAL radius)
              realmat ball_rand(INT_32BIT& ix, INTEGER d, REAL radius)
              REAL ball_volume(INTEGER d, REAL radius)
              REAL log_ball_volume(INTEGER d, REAL radius)

Prototype in  libscl.h

Description   ix is an INT_32BIT seed on entry, d is dimension. Typical
              usage is:
              static INT_32BIT ix; ix=983874; INTEGER d=3; REAL r =3.0;
              realmat x;
              for (i=0; i<n; i++) {x=ball_rand(&ix,d,r); do_something};
              REAL volume = ball_volume(d,r);

Remark        If Y is drawn from the uncorrelated multivariate normal 
              distribution, then S = Y/||Y|| has the uniform distribution 
              on the unit d-sphere.  Multiplying S by pow(U,1/d), where U has 
              the uniform distribution on the unit interval (0,1), creates 
              the uniform distribution in the unit d-dimensional ball.
              https://blogs.sas.com/content/iml/2016/04/06/
              generate-points-uniformly-in-ball.html

References    Harman, Radoslav, and Vladimir Lacko (2010), "On decompositional 
              algorithms of uniform sampling from n-spheres and n-balls," 
              Journal of Multivariate Analysis 101, 2297-2304.
              https://en.wikipedia.org/wiki/
              Unit_sphere#Unit_spheres_and_balls_in_Euclidean_space

Return value  A realmat distributed uniformly on the d-ball of specified radius.
              A REAL giving the volume of the d-ball of specified radius.

Functions     Library: pow, log
called        libscl:  unsk, ran, error

------------------------------------------------------------------------------*/

#include "libscl.h"

namespace scl {

  REAL ball_volume(INTEGER d, REAL radius)
  {
    if (d < 0 || radius < 0.0) error("Error, ball_volume, bad input");
  
    const REAL pi = M_PI;
    // const REAL pi = 4.0*atan(1.0);
  
    REAL factor = 1.0;
    if (d > 0) factor = std::pow(radius,d);
  
    REAL V[d+1];
    if (d == 0) return 1.0;
    if (d == 1) return 2.0*factor;
    if (d == 2) return pi*factor;
    V[0] = 1.0;
    V[1] = 2.0;
    for (INTEGER n=2; n<=d; ++n) V[n] = (2.0*pi/REAL(n))*V[n-2];
    return factor*V[d];
  }

  REAL log_ball_volume(INTEGER d, REAL radius)
  {
    if (d < 0 || radius < 0.0) error("Error, log_ball_volume, bad input");

    const REAL pi = M_PI;
    const REAL ln2 = M_LN2;
    const REAL lpi = std::log(pi);

    // const REAL pi = 4.0*atan(1.0);
    // const REAL ln2 = std::log(2.0);

    REAL log_factor = 0.0;
    if (d > 0) log_factor = REAL(d)*std::log(radius);

    REAL lV[d+1];
    if (d == 0) return 0.0;
    if (d == 1) return ln2 + log_factor;
    if (d == 2) return lpi + log_factor;
    lV[0] = 0.0;
    lV[1] = ln2;
    for (INTEGER n=2; n<=d; ++n) {
      lV[n] = ln2 + lpi - std::log(REAL(n)) + lV[n-2];
    }
    lV[d] += log_factor;
    return lV[d];
  }
  
  realmat ball_rand(INT_32BIT& ix, INTEGER d, REAL radius)
  {
    INT_32BIT iy = ix;
    realmat rv = ball_rand(&iy,d,radius);
    ix = iy;
    return rv;
  }
  
  realmat ball_rand(INT_32BIT* ix, INTEGER d, REAL radius)
  {
    realmat y(d,1);
    REAL sum = 0.0;
    for (INTEGER i=1; i<=d; ++i) {
      y[i] = unsk(ix);
      sum += pow(y[i],2);
    }
    y = y/sqrt(sum);

    REAL u = ran(ix);

    REAL v = pow(u,1.0/REAL(d));

    y = (radius*v)*y;
    
    return y;
  }

}
