<commit_before>/*
*   Copyright (c) 2010-2016, MIT Probabilistic Computing Project
*
*   Lead Developers: <NAME> and <NAME>
*   Authors: <NAME>, <NAME>, <NAME>, <NAME>
*   Research Leads: <NAME>, <NAME>
*
*   Licensed under the Apache License, Version 2.0 (the "License");
*   you may not use this file except in compliance with the License.
*   You may obtain a copy of the License at
*
*       http://www.apache.org/licenses/LICENSE-2.0
*
*   Unless required by applicable law or agreed to in writing, software
*   distributed under the License is distributed on an "AS IS" BASIS,
*   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*   See the License for the specific language governing permissions and
*   limitations under the License.
*/
#include "numerics.h"

#include <boost/math/special_functions/bessel.hpp>

using namespace std;

namespace numerics {

// return sign of val (signum function)
template <typename T> int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}


double estimate_vonmises_kappa(const vector<double>& X) {
    // Newton's method solution for ML estimate of kappa

    double N = (double) X.size();
    double sum_sin_x = 0;
    double sum_cos_x = 0;
    vector<double>::const_iterator it = X.begin();
    for (; it != X.end(); it++) {
        double x = *it;
        sum_sin_x += sin(x);
        sum_cos_x += cos(x);
    }
    double R2 = (sum_sin_x/N)*(sum_sin_x/N) + (sum_cos_x/N)*(sum_cos_x/N);
    double R = sqrt(R2) ;

    double kappa = R*(2.0-R2)/(1.0-R2);

    double Ap, kappa_1;
    Ap = boost::math::cyl_bessel_i(1, kappa)/boost::math::cyl_bessel_i(0, kappa);
    kappa_1 = kappa - (Ap-R)/(1-Ap*Ap-Ap/kappa);
    Ap = boost::math::cyl_bessel_i(1, kappa_1)/boost::math::cyl_bessel_i(0, kappa_1);
    kappa = kappa_1 - (Ap-R)/(1-Ap*Ap-Ap/kappa_1);
    Ap = boost::math::cyl_bessel_i(1, kappa)/boost::math::cyl_bessel_i(0, kappa);
    kappa_1 = kappa - (Ap-R)/(1-Ap*Ap-Ap/kappa);
    Ap = boost::math::cyl_bessel_i(1, kappa_1)/boost::math::cyl_bessel_i(0, kappa_1);
    kappa = kappa_1 - (Ap-R)/(1-Ap*Ap-Ap/kappa_1);

    return (kappa > 0) ? kappa : 1.0/X.size();
}

double log_bessel_0(double x){
    // calclate bessel function. Overflow will be a problem in which case we
    // use an approximation
    double i0;
    try{
      i0 = log(boost::math::cyl_bessel_i(0, x));  
    }catch(std::exception const&  ex){
      return x - .5*log(2*M_PI*x);
    }
    return i0;
  }

double calc_crp_alpha_hyperprior(double alpha) {
    double logp = 0;
    // invert the effect of log gridding
    // logp += +log(alpha);
    return logp;
}

double calc_continuous_hyperprior(double r, double nu, double s) {
    double logp = 0;
    // invert the effect of log gridding
    // MAYBE THIS IS THE WRONG PLACE FOR THIS
    // logp += log(r) + log(nu) + log(s);
    return logp;
}

double logaddexp(const vector<double>& logs) {
    double maximum = *std::max_element(logs.begin(), logs.end());
    double result = 0;
    vector<double>::const_iterator it;
    for (it = logs.begin(); it != logs.end(); it++) {
        result += exp(*it - maximum);
    }
    return log(result) + maximum;
}
<commit_msg>Local, boost-free definition of modified Bessel functions.<commit_after>/*
*   Copyright (c) 2010-2016, MIT Probabilistic Computing Project
*
*   Lead Developers: <NAME> and <NAME>
*   Authors: <NAME>, <NAME>, <NAME>, <NAME>
*   Research Leads: <NAME>, <NAME>
*
*   Licensed under the Apache License, Version 2.0 (the "License");
*   you may not use this file except in compliance with the License.
*   You may obtain a copy of the License at
*
*       http://www.apache.org/licenses/LICENSE-2.0
*
*   Unless required by applicable law or agreed to in writing, software
*   distributed under the License is distributed on an "AS IS" BASIS,
*   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*   See the License for the specific language governing permissions and
*   limitations under the License.
*/
#include "numerics.h"

using namespace std;

namespace numerics {

// return sign of val (signum function)
template <typename T> int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}

#define arraycount(A) (sizeof(A)/sizeof(*A))

double polyeval(const double a[], size_t n, double x) {
    size_t i = n;
    double y = a[--i];

    while (0 < i--) {
        y *= x;
        y += a[i];
    }

    return y;
}

// Modified Bessel functions of the first and second kinds, derived
// from
//
//      <NAME> and <NAME>, `Stable rational minimax
//      approximations to the modified Bessel functions I_0(x) and
//      I_1(x)', AECL-4928, Chalk River Nuclear Laboratories, Chalk
//      River, Ontario, October 1974.
//
//      http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/06/178/6178667.pdf

double i_0(double x) {
    static const double P1[] = {
        -2.2335582639474375249e+15,
        -5.5050369673018427753e+14,
        -3.2940087627407749166e+13,
        -8.4925101247114157499e+11,
        -1.1912746104985237192e+10,
        -1.0313066708737980747e+08,
        -5.9545626019847898221e+05,
        -2.4125195876041896775e+03,
        -7.0935347449210549190e+00,
        -1.5453977791786851041e-02,
        -2.5172644670688975051e-05,
        -3.0517226450451067446e-08,
        -2.6843448573468483278e-11,
        -1.5982226675653184646e-14,
        -5.2487866627945699800e-18,
    };
    static const double Q1[] = {
        -2.2335582639474375245e+15,
        +7.8858692566751002988e+12,
        -1.2207067397808979846e+10,
        +1.0377081058062166144e+07,
        -4.8527560179962773045e+03,
        +1.0,
    };
    static const double P2[] = {
        -2.2210262233306573296e-04,
        +1.3067392038106924055e-02,
        -4.4700805721174453923e-01,
        +5.5674518371240761397e+00,
        -2.3517945679239481621e+01,
        +3.1611322818701131207e+01,
        -9.6090021968656180000e+00,
    };
    static const double Q2[] = {
        -5.5194330231005480228e-04,
        +3.2547697594819615062e-02,
        -1.1151759188741312645e+00,
        +1.3982595353892851542e+01,
        -6.0228002066743340583e+01,
        +8.5539563258012929600e+01,
        -3.1446690275135491500e+01,
        +1.0,
    };
    double y, q, v;

    if (x < 0)
        x = -x;

    if (x == 0)
        return 1;

    if (x <= 15) {
        y = x*x;
        q = polyeval(P1, arraycount(P1), y)/polyeval(Q1, arraycount(Q1), y);
        v = q;
    } else {
        y = 1/x - static_cast<double>(1)/15;
        q = polyeval(P2, arraycount(P2), y)/polyeval(Q2, arraycount(Q2), y);
        v = (exp(x)/sqrt(x)) * q;
    }

    return v;
}

double i_1(double x) {
    static const double P1[] = {
        -1.4577180278143463643e+15,
        -1.7732037840791591320e+14,
        -6.9876779648010090070e+12,
        -1.3357437682275493024e+11,
        -1.4828267606612366099e+09,
        -1.0588550724769347106e+07,
        -5.1894091982308017540e+04,
        -1.8225946631657315931e+02,
        -4.7207090827310162436e-01,
        -9.1746443287817501309e-04,
        -1.3466829827635152875e-06,
        -1.4831904935994647675e-09,
        -1.1928788903603238754e-12,
        -6.5245515583151902910e-16,
        -1.9705291802535139930e-19,
    };
    static const double Q1[] = {
        -2.9154360556286927285e+15,
        +9.7887501377547640438e+12,
        -1.4386907088588283434e+10,
        +1.1594225856856884006e+07,
        -5.1326864679904189920e+03,
        +1.0,
    };
    static const double P2[] = {
        +1.4582087408985668208e-05,
        -8.9359825138577646443e-04,
        +2.9204895411257790122e-02,
        -3.4198728018058047439e-01,
        +1.3960118277609544334e+00,
        -1.9746376087200685843e+00,
        +8.5591872901933459000e-01,
        -6.0437159056137599999e-02,
    };
    static const double Q2[] = {
        +3.7510433111922824643e-05,
        -2.2835624489492512649e-03,
        +7.4212010813186530069e-02,
        -8.5017476463217924408e-01,
        +3.2593714889036996297e+00,
        -3.8806586721556593450e+00,
        +1.0,
    };
    double xabs, y, q, v;

    if (x == 0)
        return 0;

    xabs = fabs(x);
    if (xabs <= 15) {
        y = x*x;
        q = polyeval(P1, arraycount(P1), y)/polyeval(Q1, arraycount(Q1), y);
        v = xabs*q;
    } else {
        y = 1/xabs - static_cast<double>(1)/15;
        q = polyeval(P2, arraycount(P2), y)/polyeval(Q2, arraycount(Q2), y);
        v = (exp(xabs)/sqrt(xabs)) * q;
    }

    if (x < 0)
        v *= -v;
    return v;
}

double estimate_vonmises_kappa(const vector<double>& X) {
    // Newton's method solution for ML estimate of kappa

    double N = (double) X.size();
    double sum_sin_x = 0;
    double sum_cos_x = 0;
    vector<double>::const_iterator it = X.begin();
    for (; it != X.end(); it++) {
        double x = *it;
        sum_sin_x += sin(x);
        sum_cos_x += cos(x);
    }
    double R2 = (sum_sin_x/N)*(sum_sin_x/N) + (sum_cos_x/N)*(sum_cos_x/N);
    double R = sqrt(R2) ;

    double kappa = R*(2.0-R2)/(1.0-R2);

    double Ap, kappa_1;
    Ap = i_1(kappa)/i_0(kappa);
    kappa_1 = kappa - (Ap-R)/(1-Ap*Ap-Ap/kappa);
    Ap = i_1(kappa_1)/i_0(kappa_1);
    kappa = kappa_1 - (Ap-R)/(1-Ap*Ap-Ap/kappa_1);
    Ap = i_1(kappa)/i_0(kappa);
    kappa_1 = kappa - (Ap-R)/(1-Ap*Ap-Ap/kappa);
    Ap = i_1(kappa_1)/i_0(kappa_1);
    kappa = kappa_1 - (Ap-R)/(1-Ap*Ap-Ap/kappa_1);

    return (kappa > 0) ? kappa : 1.0/X.size();
}

double log_bessel_0(double x) {
    double i0 = i_0(x);
    // As x grows without bound, i_0(x) ---> e^x / sqrt(2 pi x).  So if
    // naive calculation of i_0 overflows, approximate it with the
    // asymptotic expansion.  See NIST DLMF 10.40(i) (Modified Bessel
    // Functions, Asymptotic Expansions for Large Argument, Hankel's
    // Expansions), Eq. (10.40.1) <http://dlmf.nist.gov/10.40.E1>.  We
    // assume the summation to be negligible if x is large enough to
    // overflow.
    if (isinf(i0))
        return x - .5*log(2*M_PI*x);
    return log(i0);
}

double calc_crp_alpha_hyperprior(double alpha) {
    double logp = 0;
    // invert the effect of log gridding
    // logp += +log(alpha);
    return logp;
}

double calc_continuous_hyperprior(double r, double nu, double s) {
    double logp = 0;
    // invert the effect of log gridding
    // MAYBE THIS IS THE WRONG PLACE FOR THIS
    // logp += log(r) + log(nu) + log(s);
    return logp;
}

double logaddexp(const vector<double>& logs) {
    double maximum = *std::max_element(logs.begin(), logs.end());
    double result = 0;
    vector<double>::const_iterator it;
    for (it = logs.begin(); it != logs.end(); it++) {
        result += exp(*it - maximum);
    }
    return log(result) + maximum;
}

{'fraction_non_alphanumeric': 0.14223215116820256, 'fraction_numerical': 0.1589872475100065, 'mean_word_length': 2.5493888338288735, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5184316', 'n_tokens_mistral': 5145, 'n_tokens_neox': 4170, 'n_words': 913}