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

Copyright (C) 2020, 2021

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      bump_function - infinitely differentiable with compact support

Syntax        #include "libscl.h"
              REAL phi(REAL x, REAL x0, REAL r)
              REAL d1phi(REAL x, REAL x0, REAL r)
              REAL d2phi(REAL x, REAL x0, REAL r)
              REAL phi(const realmat x, const realmat& x0, const realmat& r)
              realmat d1phi(const realmat x,const realmat& x0,const realmat& r)
              realmat d2phi(const realmat x,const realmat& x0,const realmat& r)
              void bfquad(REAL x0, REAL r, INTEGER n, realmat& a, realmat& w)
              void bfquad(const realmat& x0, const realmat& r, INTEGER n, 
                vector(realmat)& a, realmat& w)

Prototype in  libscl.h

Description   x0 is the center and r is the radius of the support of the
              bump function. phi evaluates the bump function, d1phi returns 
              its first derivative, and d2phi its second. The realmat 
              inputs x, x0, r must be column vectors of the same order. 
              This requirement is not checked. 
              The REAL function bfquad returns a quadrature rule of length 
              n with abcissae a and weights w.  It is a trapezoid rule with 
              unequal spacing. The multivariate version is a Riemann sum and 
              n is the number of points along the diameter of the support 
              of the bump function before trimming to the support. The 
              number of quadrature points is w.size(). See 
              https://people.sc.fsu.edu/~jburkardt/c_src/stroud/stroud.html
              for potentially better rules.

Remarks       Infinitely differentiable functions with compact support
              are used to define distributions, weak derivatives, and 
              solutions to differential equations. This by facilitating
              integration by parts: integral_S (D^a u(x)) phi(x) dx 
              = (-1)^|a| integral_S u(x) (D^a phi(x)) dx where S, includes 
              the support of phi and a is a multi-index.
              Note that for the multivariate quadrature rule, a indexes
              from 0 to one less than npoints=a.size=w.size whereas w 
              indexes from 1 to npoints. See the example.

Example       heat equation: (d/dt)u(p,t) = alpha*(d^2/p^2)u(p,t), alpha = 1
              variables    : p is position on a bar, t is time, x = (p,t)
              domain       : 0 < p < pi, 0 < t < 3
              boundary cond: u(0,t) = u(pi,t) = T, u(p,0) = sin(p), T = 50
              solution     : u(p,t) = T + exp(-alpha*t)*sin(p)

              realmat x0(2,1,1.0);
              realmat radius(2,1,0.75);
              INTEGER n = 150;
              vector<realmat> a;
              realmat w;
              bfquad(x0,radius,n,a,w);
              INTEGER npoints = w.size();
              REAL lhs = 0.0; REAL rhs = 0.0;
              for (INTEGER j=1; j<=npoints; ++j) {
                realmat x = a[j-1];
                REAL u = 50.0 + exp(-x[2])*sin(x[1]);
                realmat H = d2phi(x,x0,radius);
                rhs += u*H(1,1)*w[j];
                realmat D = d1phi(x,x0,radius);
                lhs -= u*D[2]*w[j];
              }
              cout << "lhs = "<< lhs << "\nrhs = " << rhs << '\n';

              answer:
              lhs = -0.0812326
              rhs = -0.0812346

Return value  The bump function or its first and second derivatives as a
              REAL or realmat as indicated by the syntax above.

Functions     Library: pow
called        libscl: (none)

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

#include "libscl.h"
using namespace scl;
using namespace std;

namespace {
  REAL z2(REAL x, REAL x0, REAL r)  
  {
    return pow((x - x0)/r, 2);
  }
  
  REAL d1z2(REAL x, REAL x0, REAL r)  
  {
    return 2.0*(x - x0)/pow(r,2);
  }
    
  REAL d2z2(REAL x, REAL x0, REAL r)  
  {
    return 2.0/pow(r,2);
  }

  REAL
  z2(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    INTEGER d = x.size();
    REAL sumsqr = 0.0;
    for (INTEGER i=1; i<=d; ++i) {
      REAL e = (x[i] - x0[i])/r[i];
      sumsqr += pow(e,2);
    }
    return sumsqr;
  }
  
  scl::realmat 
  d1z2(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    INTEGER d = x.size();
    scl::realmat del(d,1);
    for (INTEGER i=1; i<=d; ++i) {
      del[i] = 2.0*(x[i] - x0[i])/pow(r[i],2);
    }
    return del;
  }
  
  scl::realmat 
  d2z2(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    INTEGER d = x.size();
    scl::realmat H(d,d,0.0);
    for (INTEGER i=1; i<=d; ++i) {
      H(i,i) = 2.0/pow(r[i],2);
    }
    return H;
  }

  REAL piecewise_linear(REAL x)
  {
    REAL a1 = -1.0;
    REAL b1 = -1.0/3.0;
    REAL c1 = -1.0;
    REAL d1 = -0.4;

    REAL a2 = b1;
    REAL b2 = 1.0/3.0;
    REAL c2 = d1;
    REAL d2 = 0.4;

    REAL a3 = b2;
    REAL b3 = 1.0;
    REAL c3 = d2;
    REAL d3 = 1.0;

    if ( (a1 <= x) && (x < b1) ) {
      return c1 + ((d1-c1)/(b1-a1))*(x-a1);
    }
    else if ( (a2 <= x) && (x < b2) ) {
      return c2 + ((d2-c2)/(b2-a2))*(x-a2);
    }
    else {
      return c3 + ((d3-c3)/(b3-a3))*(x-a3);
    }
  }

  void grd
  (INTEGER d, const realmat& inc, const realmat& r, const realmat& x0, 
    realmat& x, vector<realmat>& vgrid)
  {
    INTEGER dim = inc.nrow();

    realmat lo = x0 - r;
    realmat hi = x0 + r;
  
    if (d == dim) {  // external call must have d = dim
  
      //vgrid.clear();
  
      x.resize(dim,1);
  
      INTEGER j = 1;
  
      for (REAL t = lo[j]; t<=hi[j]; t += inc[j]) {
        x[j] = t;
        if (d == 1) {
          REAL sumsqr = 0.0;
          for (INTEGER i=1; i<=dim; ++i) {
            REAL e = (x[i] - x0[i])/r[i];
            sumsqr += pow(e,2);
          }
          if (sumsqr <= 1.0) {
            vgrid.push_back(x);
          }
        }
        else {
          grd(d-1,inc,r,x0,x,vgrid);
        }
      }
    } 
    else {
      
      INTEGER j = dim - d + 1;
  
      for (REAL t = lo[j]; t<=hi[j]; t += inc[j]) {
        x[j] = t;
        if (d == 1) {
          REAL sumsqr = 0.0;
          for (INTEGER i=1; i<=dim; ++i) {
            REAL e = (x[i] - x0[i])/r[i];
            sumsqr += pow(e,2);
          }
          if (sumsqr <= 1.0) {
            vgrid.push_back(x);
          }    
        }
        else {
          grd(d-1,inc,r,x0,x,vgrid);
        }
      }
    }
  }
}

namespace scl {
  REAL phi(REAL x, REAL x0, REAL r) 
  {
    if (z2(x,x0,r) >= 1.0) return 0.0;
    return exp( 1.0/(z2(x,x0,r)-1.0) );
  }
  
  REAL d1phi(REAL x, REAL x0, REAL r)
  {
    if (z2(x,x0,r) >= 1.0) return 0.0;
    return ( -1.0/pow(z2(x,x0,r)-1.0,2) ) * phi(x, x0, r) * d1z2(x,x0,r);
  }
  
  REAL d2phi(REAL x, REAL x0, REAL r)
  {
    if (z2(x,x0,r) >= 1.0) return 0.0;
    REAL sum = 0.0;
    sum += 2.0/pow(z2(x,x0,r)-1.0,3) * pow(d1z2(x,x0,r),2);
    sum += pow(1.0/pow(z2(x,x0,r)-1.0,2),2) * pow(d1z2(x,x0,r),2);
    sum += ( -1.0/pow(z2(x,x0,r)-1.0,2) ) * d2z2(x,x0,r);
    sum *= phi(x,x0,r);
    return sum;
  }
  
  REAL 
  phi(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    if (z2(x,x0,r) >= 1.0) return 0.0;
    return exp( 1.0/(z2(x,x0,r)-1.0) );
  }
  
  scl::realmat
  d1phi(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    INTEGER d = x.size();
    if (z2(x,x0,r) >= 1.0) {
      return scl::realmat(d,d,0.0);
    }
    return ( -1.0/pow(z2(x,x0,r)-1.0,2) ) * phi(x, x0, r) * d1z2(x,x0,r);
  }
  
  scl::realmat
  d2phi(const scl::realmat x, const scl::realmat& x0, const scl::realmat& r)
  {
    INTEGER d = x.size();
    if (z2(x,x0,r) >= 1.0) {
      return scl::realmat(d,d,0.0);
    }
    scl::realmat sum(d,d,0.0);
    sum += 2.0/pow(z2(x,x0,r)-1.0,3) * (d1z2(x,x0,r)*T(d1z2(x,x0,r)));
    sum += pow(1.0/pow(z2(x,x0,r)-1.0,2),2) * (d1z2(x,x0,r)*T(d1z2(x,x0,r)));
    sum += ( -1.0/pow(z2(x,x0,r)-1.0,2) ) * d2z2(x,x0,r);
    REAL ph = phi(x,x0,r);
    for (INTEGER j=1; j<=d; ++j) {
      for (INTEGER i=1; i<=d; ++i) {
        sum(i,j) *= ph;
      }
    }
    return sum;
  }
  
  void bfquad(REAL x0, REAL r, INTEGER n, realmat& a, realmat& w)
  {
    if (n < 24) error("Error, bump_function, bfquad, n must be 24 or larger");
  
    a.resize(n,1);
    w.resize(n,1);
  
    for (INTEGER i=1; i<=n; ++i) {
      REAL x = -1.0 + 2.0*REAL(i - 1)/REAL(n - 1);
      a[i] = piecewise_linear(x);
    }
  
    w[n] = 0.0;
    for (INTEGER i=2; i<=n; ++i) {
      REAL delta = a[i]-a[i-1];
      w[i-1] = 0.5*delta;
    }
    for (INTEGER i=2; i<=n; ++i) {
      REAL delta = a[i]-a[i-1];
      w[i] += 0.5*delta;
    }
  
    for (INTEGER i=1; i<=n; ++i) {
      a[i] *= r;
      a[i] += x0;
      w[i] *= r;
    }
    return;
  }
  
  extern void
  bfquad(const realmat& x0, const realmat& r, INTEGER n, 
  vector<realmat>& a, realmat& w)
  {
    INTEGER d = x0.nrow();
  
    if (x0.ncol() != 1) error("Error, bfquad, x0 not a column vector");
    if (r.nrow() != d || r.ncol() !=1) error("Error, grid, bad r dimension");
  
    for (INTEGER i=1; i<=d; ++i) {
      if (r[i] <= 0.0) error("Error, bfquad, elements of r not positive");
    }
  
    realmat inc = 2.0*r/REAL(n);
  
    a.clear();
  
    realmat x;
  
    grd(d, inc, r, x0, x, a);
  
    INTEGER ngrid = a.size();
  
    REAL wi = 1.0;
    for (INTEGER i=1; i<=d; ++i) wi *= inc[i];
  
    w.resize(ngrid,1);
    for (INTEGER i=1; i<=ngrid; ++i) w[i] = wi;
  }
}
