# include <cmath>
# include <cstdlib>

# include "sigmoid_derivative.hpp"

//****************************************************************************80

int nchoosek ( int n, int k )

//****************************************************************************80
//
//  Purpose:
//
//    nchoosek() computes the binomial coefficient C(N,K).
//
//  Discussion:
//
//    The value is calculated in such a way as to avoid overflow and
//    roundoff.  The calculation is done in int arithmetic.
//
//    The formula used is:
//
//      C(N,K) = N! / ( K! * (N-K)! )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    09 December 2013
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    ML Wolfson, HV Wright,
//    Algorithm 160:
//    Combinatorial of M Things Taken N at a Time,
//    Communications of the ACM,
//    Volume 6, Number 4, April 1963, page 161.
//
//  Input:
//
//    int N, K, are the values of N and K.
//
//  Output:
//
//    int nchoosek: the number of combinations of N
//    things taken K at a time.
//
{
  int i;
  int mn;
  int mx;
  int value;

  if ( k < n - k )
  {
    mn = k;
    mx = n - k;
  }
  else
  {
    mn = n - k;
    mx = k;
  }

  if ( mn < 0 )
  {
    value = 0;
  }
  else if ( mn == 0 )
  {
    value = 1;
  }
  else
  {
    value = mx + 1;

    for ( i = 2; i <= mn; i++ )
    {
      value = ( value * ( mx + i ) ) / i;
    }
  }

  return value;
}
//****************************************************************************80

double *sigmoid_derivative_coef ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    sigmoid_derivative_coef(): sigmoid derivative expansion coefficients.
//
//  Discussion:
//
//    s(x) = 1 / ( 1 + exp ( - x ) )
//
//    s(x)^(n) = sum ( 1 <= j <= n + 1 ) coef(j) * s(x)^j 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Joe McKenna,
//    Derivatives of the sigmoid function,
//    https://joepatmckenna.github.io/calculus/derivative/sigmoid!20function/linear!20albegra/2018/01/20/sigmoid-derivs/
//
//  Input:
//
//    int n: the order of the derivative.
//    0 <= n.
//
//  Output:
//
//    double coef[n+2]: the expansion coefficients.
//    coef[0] is always 0, and coef[1] is always 1.
//
{
  double cnk;
  double *coef;
  int j;
  int k;
  int mop;
  
  coef = new double[n+2];

  for ( k = 0; k < n + 2; k++ )
  {
    coef[k] = 0.0;
  }

  for ( k = 0; k <= n; k++ )
  {
    cnk = 0.0;
    mop = -1;
    for ( j = 0; j <= k; j++ )
    {
      mop = - mop;
      cnk = cnk + mop * pow ( j + 1, n ) * nchoosek ( k, j );
    }
    coef[k+1] = cnk;
  }

  return coef;
}
//****************************************************************************80

double sigmoid_derivative_value ( int n, double x )

//****************************************************************************80
//
//  Purpose:
//
//    sigmoid_derivative_value(): evaluate sigmoid derivative at x.
//
//  Discussion:
//
//    s(x) = 1 / ( 1 + exp ( - x ) )
//
//    s(x)^(n) = sum ( 1 <= j <= n + 1 ) coef(j) * s(x)^j 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 May 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Joe McKenna,
//    Derivatives of the sigmoid function,
//    https://joepatmckenna.github.io/calculus/derivative/sigmoid!20function/linear!20albegra/2018/01/20/sigmoid-derivs/
//
//  Input:
//
//    int n: the order of the derivative.
//
//    double x: the evaluation point.
//
//  Output:
//
//    double sigmoid_derivative_value: the value of s^(n)(x).
//
{
  double cnk;
  int j;
  int k;
  int mop;
  double sig;
  double sigk;
  double value;

  sig = 1.0 / ( 1.0 + exp ( - x ) );
  value = 0.0;
  sigk = 1.0;

  for ( k = 0; k <= n; k++ )
  {
    sigk = sigk * sig;
    cnk = 0.0;
    mop = -1;
    for ( j = 0; j <= k; j++ )
    {
      mop = - mop;
      cnk = cnk + mop * pow ( j + 1, n ) * nchoosek ( k, j );
    }
    value = value + cnk * sigk;
  }

  return value;
}