Add graph references
This commit is contained in:
Abdelrahman Said
@@ -0,0 +1,46 @@
|
||||
# Declare the files needed to compile our vendored plfit copy
|
||||
add_library(
|
||||
plfit_vendored
|
||||
OBJECT
|
||||
EXCLUDE_FROM_ALL
|
||||
gss.c
|
||||
hzeta.c
|
||||
kolmogorov.c
|
||||
lbfgs.c
|
||||
mt.c
|
||||
options.c
|
||||
plfit.c
|
||||
plfit_error.c
|
||||
rbinom.c
|
||||
sampling.c
|
||||
)
|
||||
|
||||
target_include_directories(
|
||||
plfit_vendored
|
||||
PRIVATE
|
||||
${PROJECT_SOURCE_DIR}/include
|
||||
${PROJECT_BINARY_DIR}/include
|
||||
PUBLIC
|
||||
${CMAKE_CURRENT_SOURCE_DIR}
|
||||
)
|
||||
|
||||
if (BUILD_SHARED_LIBS)
|
||||
set_property(TARGET plfit_vendored PROPERTY POSITION_INDEPENDENT_CODE ON)
|
||||
endif()
|
||||
|
||||
# Since these are included as object files, they should call the
|
||||
# function as is (without visibility specification)
|
||||
target_compile_definitions(plfit_vendored PRIVATE IGRAPH_STATIC)
|
||||
|
||||
use_all_warnings(plfit_vendored)
|
||||
|
||||
if (MSVC)
|
||||
target_compile_options(
|
||||
plfit_vendored PRIVATE
|
||||
/wd4100
|
||||
) # disable unreferenced parameter warning
|
||||
endif()
|
||||
|
||||
if(IGRAPH_OPENMP_SUPPORT)
|
||||
target_link_libraries(plfit_vendored PRIVATE OpenMP::OpenMP_C)
|
||||
endif()
|
||||
@@ -0,0 +1,133 @@
|
||||
/*
|
||||
* ANSI C implementation of vector operations.
|
||||
*
|
||||
* Copyright (c) 2007-2010 Naoaki Okazaki
|
||||
* All rights reserved.
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
* of this software and associated documentation files (the "Software"), to deal
|
||||
* in the Software without restriction, including without limitation the rights
|
||||
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
* copies of the Software, and to permit persons to whom the Software is
|
||||
* furnished to do so, subject to the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be included in
|
||||
* all copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
* THE SOFTWARE.
|
||||
*/
|
||||
|
||||
/* $Id$ */
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <memory.h>
|
||||
|
||||
#if LBFGS_FLOAT == 32 && LBFGS_IEEE_FLOAT
|
||||
#define fsigndiff(x, y) (((*(uint32_t*)(x)) ^ (*(uint32_t*)(y))) & 0x80000000U)
|
||||
#else
|
||||
#define fsigndiff(x, y) (*(x) * (*(y) / fabs(*(y))) < 0.)
|
||||
#endif/*LBFGS_IEEE_FLOAT*/
|
||||
|
||||
inline static void* vecalloc(size_t size)
|
||||
{
|
||||
void *memblock = malloc(size);
|
||||
if (memblock) {
|
||||
memset(memblock, 0, size);
|
||||
}
|
||||
return memblock;
|
||||
}
|
||||
|
||||
inline static void vecfree(void *memblock)
|
||||
{
|
||||
free(memblock);
|
||||
}
|
||||
|
||||
inline static void vecset(lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
x[i] = c;
|
||||
}
|
||||
}
|
||||
|
||||
inline static void veccpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
y[i] = x[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecncpy(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
y[i] = -x[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecadd(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const lbfgsfloatval_t c, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
y[i] += c * x[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecdiff(lbfgsfloatval_t *z, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
z[i] = x[i] - y[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecscale(lbfgsfloatval_t *y, const lbfgsfloatval_t c, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
y[i] *= c;
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecmul(lbfgsfloatval_t *y, const lbfgsfloatval_t *x, const int n)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0;i < n;++i) {
|
||||
y[i] *= x[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vecdot(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const lbfgsfloatval_t *y, const int n)
|
||||
{
|
||||
int i;
|
||||
*s = 0.;
|
||||
for (i = 0;i < n;++i) {
|
||||
*s += x[i] * y[i];
|
||||
}
|
||||
}
|
||||
|
||||
inline static void vec2norm(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)
|
||||
{
|
||||
vecdot(s, x, x, n);
|
||||
*s = (lbfgsfloatval_t)sqrt(*s);
|
||||
}
|
||||
|
||||
inline static void vec2norminv(lbfgsfloatval_t* s, const lbfgsfloatval_t *x, const int n)
|
||||
{
|
||||
vec2norm(s, x, n);
|
||||
*s = (lbfgsfloatval_t)(1.0 / *s);
|
||||
}
|
||||
@@ -0,0 +1,294 @@
|
||||
/*
|
||||
* SSE2 implementation of vector oprations (64bit double).
|
||||
*
|
||||
* Copyright (c) 2007-2010 Naoaki Okazaki
|
||||
* All rights reserved.
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
* of this software and associated documentation files (the "Software"), to deal
|
||||
* in the Software without restriction, including without limitation the rights
|
||||
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
* copies of the Software, and to permit persons to whom the Software is
|
||||
* furnished to do so, subject to the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be included in
|
||||
* all copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
* THE SOFTWARE.
|
||||
*/
|
||||
|
||||
/* $Id$ */
|
||||
|
||||
#include <stdlib.h>
|
||||
#ifndef __APPLE__
|
||||
#include <malloc.h>
|
||||
#endif
|
||||
#include <memory.h>
|
||||
|
||||
#if 1400 <= _MSC_VER
|
||||
#include <intrin.h>
|
||||
#endif/*1400 <= _MSC_VER*/
|
||||
|
||||
#if HAVE_EMMINTRIN_H
|
||||
#include <emmintrin.h>
|
||||
#endif/*HAVE_EMMINTRIN_H*/
|
||||
|
||||
inline static void* vecalloc(size_t size)
|
||||
{
|
||||
#if defined(_WIN32)
|
||||
void *memblock = _aligned_malloc(size, 16);
|
||||
#elif defined(__APPLE__) /* OS X always aligns on 16-byte boundaries */
|
||||
void *memblock = malloc(size);
|
||||
#else
|
||||
void *memblock = NULL, *p = NULL;
|
||||
if (posix_memalign(&p, 16, size) == 0) {
|
||||
memblock = p;
|
||||
}
|
||||
#endif
|
||||
if (memblock != NULL) {
|
||||
memset(memblock, 0, size);
|
||||
}
|
||||
return memblock;
|
||||
}
|
||||
|
||||
inline static void vecfree(void *memblock)
|
||||
{
|
||||
#ifdef _MSC_VER
|
||||
_aligned_free(memblock);
|
||||
#else
|
||||
free(memblock);
|
||||
#endif
|
||||
}
|
||||
|
||||
#define fsigndiff(x, y) \
|
||||
((_mm_movemask_pd(_mm_set_pd(*(x), *(y))) + 1) & 0x002)
|
||||
|
||||
#define vecset(x, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM0 = _mm_set1_pd(c); \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
_mm_store_pd((x)+i , XMM0); \
|
||||
_mm_store_pd((x)+i+2, XMM0); \
|
||||
_mm_store_pd((x)+i+4, XMM0); \
|
||||
_mm_store_pd((x)+i+6, XMM0); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define veccpy(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128d XMM0 = _mm_load_pd((x)+i ); \
|
||||
__m128d XMM1 = _mm_load_pd((x)+i+2); \
|
||||
__m128d XMM2 = _mm_load_pd((x)+i+4); \
|
||||
__m128d XMM3 = _mm_load_pd((x)+i+6); \
|
||||
_mm_store_pd((y)+i , XMM0); \
|
||||
_mm_store_pd((y)+i+2, XMM1); \
|
||||
_mm_store_pd((y)+i+4, XMM2); \
|
||||
_mm_store_pd((y)+i+6, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecncpy(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128d XMM0 = _mm_setzero_pd(); \
|
||||
__m128d XMM1 = _mm_setzero_pd(); \
|
||||
__m128d XMM2 = _mm_setzero_pd(); \
|
||||
__m128d XMM3 = _mm_setzero_pd(); \
|
||||
__m128d XMM4 = _mm_load_pd((x)+i ); \
|
||||
__m128d XMM5 = _mm_load_pd((x)+i+2); \
|
||||
__m128d XMM6 = _mm_load_pd((x)+i+4); \
|
||||
__m128d XMM7 = _mm_load_pd((x)+i+6); \
|
||||
XMM0 = _mm_sub_pd(XMM0, XMM4); \
|
||||
XMM1 = _mm_sub_pd(XMM1, XMM5); \
|
||||
XMM2 = _mm_sub_pd(XMM2, XMM6); \
|
||||
XMM3 = _mm_sub_pd(XMM3, XMM7); \
|
||||
_mm_store_pd((y)+i , XMM0); \
|
||||
_mm_store_pd((y)+i+2, XMM1); \
|
||||
_mm_store_pd((y)+i+4, XMM2); \
|
||||
_mm_store_pd((y)+i+6, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecadd(y, x, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM7 = _mm_set1_pd(c); \
|
||||
for (i = 0;i < (n);i += 4) { \
|
||||
__m128d XMM0 = _mm_load_pd((x)+i ); \
|
||||
__m128d XMM1 = _mm_load_pd((x)+i+2); \
|
||||
__m128d XMM2 = _mm_load_pd((y)+i ); \
|
||||
__m128d XMM3 = _mm_load_pd((y)+i+2); \
|
||||
XMM0 = _mm_mul_pd(XMM0, XMM7); \
|
||||
XMM1 = _mm_mul_pd(XMM1, XMM7); \
|
||||
XMM2 = _mm_add_pd(XMM2, XMM0); \
|
||||
XMM3 = _mm_add_pd(XMM3, XMM1); \
|
||||
_mm_store_pd((y)+i , XMM2); \
|
||||
_mm_store_pd((y)+i+2, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecdiff(z, x, y, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128d XMM0 = _mm_load_pd((x)+i ); \
|
||||
__m128d XMM1 = _mm_load_pd((x)+i+2); \
|
||||
__m128d XMM2 = _mm_load_pd((x)+i+4); \
|
||||
__m128d XMM3 = _mm_load_pd((x)+i+6); \
|
||||
__m128d XMM4 = _mm_load_pd((y)+i ); \
|
||||
__m128d XMM5 = _mm_load_pd((y)+i+2); \
|
||||
__m128d XMM6 = _mm_load_pd((y)+i+4); \
|
||||
__m128d XMM7 = _mm_load_pd((y)+i+6); \
|
||||
XMM0 = _mm_sub_pd(XMM0, XMM4); \
|
||||
XMM1 = _mm_sub_pd(XMM1, XMM5); \
|
||||
XMM2 = _mm_sub_pd(XMM2, XMM6); \
|
||||
XMM3 = _mm_sub_pd(XMM3, XMM7); \
|
||||
_mm_store_pd((z)+i , XMM0); \
|
||||
_mm_store_pd((z)+i+2, XMM1); \
|
||||
_mm_store_pd((z)+i+4, XMM2); \
|
||||
_mm_store_pd((z)+i+6, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecscale(y, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM7 = _mm_set1_pd(c); \
|
||||
for (i = 0;i < (n);i += 4) { \
|
||||
__m128d XMM0 = _mm_load_pd((y)+i ); \
|
||||
__m128d XMM1 = _mm_load_pd((y)+i+2); \
|
||||
XMM0 = _mm_mul_pd(XMM0, XMM7); \
|
||||
XMM1 = _mm_mul_pd(XMM1, XMM7); \
|
||||
_mm_store_pd((y)+i , XMM0); \
|
||||
_mm_store_pd((y)+i+2, XMM1); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecmul(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128d XMM0 = _mm_load_pd((x)+i ); \
|
||||
__m128d XMM1 = _mm_load_pd((x)+i+2); \
|
||||
__m128d XMM2 = _mm_load_pd((x)+i+4); \
|
||||
__m128d XMM3 = _mm_load_pd((x)+i+6); \
|
||||
__m128d XMM4 = _mm_load_pd((y)+i ); \
|
||||
__m128d XMM5 = _mm_load_pd((y)+i+2); \
|
||||
__m128d XMM6 = _mm_load_pd((y)+i+4); \
|
||||
__m128d XMM7 = _mm_load_pd((y)+i+6); \
|
||||
XMM4 = _mm_mul_pd(XMM4, XMM0); \
|
||||
XMM5 = _mm_mul_pd(XMM5, XMM1); \
|
||||
XMM6 = _mm_mul_pd(XMM6, XMM2); \
|
||||
XMM7 = _mm_mul_pd(XMM7, XMM3); \
|
||||
_mm_store_pd((y)+i , XMM4); \
|
||||
_mm_store_pd((y)+i+2, XMM5); \
|
||||
_mm_store_pd((y)+i+4, XMM6); \
|
||||
_mm_store_pd((y)+i+6, XMM7); \
|
||||
} \
|
||||
}
|
||||
|
||||
|
||||
|
||||
#if 3 <= __SSE__ || defined(__SSE3__)
|
||||
/*
|
||||
Horizontal add with haddps SSE3 instruction. The work register (rw)
|
||||
is unused.
|
||||
*/
|
||||
#define __horizontal_sum(r, rw) \
|
||||
r = _mm_hadd_ps(r, r); \
|
||||
r = _mm_hadd_ps(r, r);
|
||||
|
||||
#else
|
||||
/*
|
||||
Horizontal add with SSE instruction. The work register (rw) is used.
|
||||
*/
|
||||
#define __horizontal_sum(r, rw) \
|
||||
rw = r; \
|
||||
r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(1, 0, 3, 2)); \
|
||||
r = _mm_add_ps(r, rw); \
|
||||
rw = r; \
|
||||
r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(2, 3, 0, 1)); \
|
||||
r = _mm_add_ps(r, rw);
|
||||
|
||||
#endif
|
||||
|
||||
#define vecdot(s, x, y, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM0 = _mm_setzero_pd(); \
|
||||
__m128d XMM1 = _mm_setzero_pd(); \
|
||||
__m128d XMM2, XMM3, XMM4, XMM5; \
|
||||
for (i = 0;i < (n);i += 4) { \
|
||||
XMM2 = _mm_load_pd((x)+i ); \
|
||||
XMM3 = _mm_load_pd((x)+i+2); \
|
||||
XMM4 = _mm_load_pd((y)+i ); \
|
||||
XMM5 = _mm_load_pd((y)+i+2); \
|
||||
XMM2 = _mm_mul_pd(XMM2, XMM4); \
|
||||
XMM3 = _mm_mul_pd(XMM3, XMM5); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_pd(XMM1, XMM3); \
|
||||
} \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
_mm_store_sd((s), XMM0); \
|
||||
}
|
||||
|
||||
#define vec2norm(s, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM0 = _mm_setzero_pd(); \
|
||||
__m128d XMM1 = _mm_setzero_pd(); \
|
||||
__m128d XMM2, XMM3, XMM4, XMM5; \
|
||||
for (i = 0;i < (n);i += 4) { \
|
||||
XMM2 = _mm_load_pd((x)+i ); \
|
||||
XMM3 = _mm_load_pd((x)+i+2); \
|
||||
XMM4 = XMM2; \
|
||||
XMM5 = XMM3; \
|
||||
XMM2 = _mm_mul_pd(XMM2, XMM4); \
|
||||
XMM3 = _mm_mul_pd(XMM3, XMM5); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_pd(XMM1, XMM3); \
|
||||
} \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
XMM0 = _mm_sqrt_pd(XMM0); \
|
||||
_mm_store_sd((s), XMM0); \
|
||||
}
|
||||
|
||||
|
||||
#define vec2norminv(s, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128d XMM0 = _mm_setzero_pd(); \
|
||||
__m128d XMM1 = _mm_setzero_pd(); \
|
||||
__m128d XMM2, XMM3, XMM4, XMM5; \
|
||||
for (i = 0;i < (n);i += 4) { \
|
||||
XMM2 = _mm_load_pd((x)+i ); \
|
||||
XMM3 = _mm_load_pd((x)+i+2); \
|
||||
XMM4 = XMM2; \
|
||||
XMM5 = XMM3; \
|
||||
XMM2 = _mm_mul_pd(XMM2, XMM4); \
|
||||
XMM3 = _mm_mul_pd(XMM3, XMM5); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_pd(XMM1, XMM3); \
|
||||
} \
|
||||
XMM2 = _mm_set1_pd(1.0); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
XMM1 = _mm_shuffle_pd(XMM0, XMM0, _MM_SHUFFLE2(1, 1)); \
|
||||
XMM0 = _mm_add_pd(XMM0, XMM1); \
|
||||
XMM0 = _mm_sqrt_pd(XMM0); \
|
||||
XMM2 = _mm_div_pd(XMM2, XMM0); \
|
||||
_mm_store_sd((s), XMM2); \
|
||||
}
|
||||
@@ -0,0 +1,302 @@
|
||||
/*
|
||||
* SSE/SSE3 implementation of vector oprations (32bit float).
|
||||
*
|
||||
* Copyright (c) 2007-2010 Naoaki Okazaki
|
||||
* All rights reserved.
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
* of this software and associated documentation files (the "Software"), to deal
|
||||
* in the Software without restriction, including without limitation the rights
|
||||
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
* copies of the Software, and to permit persons to whom the Software is
|
||||
* furnished to do so, subject to the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be included in
|
||||
* all copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
* THE SOFTWARE.
|
||||
*/
|
||||
|
||||
/* $Id$ */
|
||||
|
||||
#include <stdlib.h>
|
||||
#ifndef __APPLE__
|
||||
#include <malloc.h>
|
||||
#endif
|
||||
#include <memory.h>
|
||||
|
||||
#if 1400 <= _MSC_VER
|
||||
#include <intrin.h>
|
||||
#endif/*_MSC_VER*/
|
||||
|
||||
#if HAVE_XMMINTRIN_H
|
||||
#include <xmmintrin.h>
|
||||
#endif/*HAVE_XMMINTRIN_H*/
|
||||
|
||||
#if LBFGS_FLOAT == 32 && LBFGS_IEEE_FLOAT
|
||||
#define fsigndiff(x, y) (((*(uint32_t*)(x)) ^ (*(uint32_t*)(y))) & 0x80000000U)
|
||||
#else
|
||||
#define fsigndiff(x, y) (*(x) * (*(y) / fabs(*(y))) < 0.)
|
||||
#endif/*LBFGS_IEEE_FLOAT*/
|
||||
|
||||
inline static void* vecalloc(size_t size)
|
||||
{
|
||||
#if defined(_MSC_VER)
|
||||
void *memblock = _aligned_malloc(size, 16);
|
||||
#elif defined(__APPLE__) /* OS X always aligns on 16-byte boundaries */
|
||||
void *memblock = malloc(size);
|
||||
#else
|
||||
void *memblock = NULL, *p = NULL;
|
||||
if (posix_memalign(&p, 16, size) == 0) {
|
||||
memblock = p;
|
||||
}
|
||||
#endif
|
||||
if (memblock != NULL) {
|
||||
memset(memblock, 0, size);
|
||||
}
|
||||
return memblock;
|
||||
}
|
||||
|
||||
inline static void vecfree(void *memblock)
|
||||
{
|
||||
#ifdef _MSC_VER
|
||||
_aligned_free(memblock);
|
||||
#else
|
||||
free(memblock);
|
||||
#endif
|
||||
}
|
||||
|
||||
#define vecset(x, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM0 = _mm_set_ps1(c); \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
_mm_store_ps((x)+i , XMM0); \
|
||||
_mm_store_ps((x)+i+ 4, XMM0); \
|
||||
_mm_store_ps((x)+i+ 8, XMM0); \
|
||||
_mm_store_ps((x)+i+12, XMM0); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define veccpy(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
__m128 XMM0 = _mm_load_ps((x)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((x)+i+ 4); \
|
||||
__m128 XMM2 = _mm_load_ps((x)+i+ 8); \
|
||||
__m128 XMM3 = _mm_load_ps((x)+i+12); \
|
||||
_mm_store_ps((y)+i , XMM0); \
|
||||
_mm_store_ps((y)+i+ 4, XMM1); \
|
||||
_mm_store_ps((y)+i+ 8, XMM2); \
|
||||
_mm_store_ps((y)+i+12, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecncpy(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
const uint32_t mask = 0x80000000; \
|
||||
__m128 XMM4 = _mm_load_ps1((float*)&mask); \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
__m128 XMM0 = _mm_load_ps((x)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((x)+i+ 4); \
|
||||
__m128 XMM2 = _mm_load_ps((x)+i+ 8); \
|
||||
__m128 XMM3 = _mm_load_ps((x)+i+12); \
|
||||
XMM0 = _mm_xor_ps(XMM0, XMM4); \
|
||||
XMM1 = _mm_xor_ps(XMM1, XMM4); \
|
||||
XMM2 = _mm_xor_ps(XMM2, XMM4); \
|
||||
XMM3 = _mm_xor_ps(XMM3, XMM4); \
|
||||
_mm_store_ps((y)+i , XMM0); \
|
||||
_mm_store_ps((y)+i+ 4, XMM1); \
|
||||
_mm_store_ps((y)+i+ 8, XMM2); \
|
||||
_mm_store_ps((y)+i+12, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecadd(y, x, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM7 = _mm_set_ps1(c); \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128 XMM0 = _mm_load_ps((x)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((x)+i+4); \
|
||||
__m128 XMM2 = _mm_load_ps((y)+i ); \
|
||||
__m128 XMM3 = _mm_load_ps((y)+i+4); \
|
||||
XMM0 = _mm_mul_ps(XMM0, XMM7); \
|
||||
XMM1 = _mm_mul_ps(XMM1, XMM7); \
|
||||
XMM2 = _mm_add_ps(XMM2, XMM0); \
|
||||
XMM3 = _mm_add_ps(XMM3, XMM1); \
|
||||
_mm_store_ps((y)+i , XMM2); \
|
||||
_mm_store_ps((y)+i+4, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecdiff(z, x, y, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
__m128 XMM0 = _mm_load_ps((x)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((x)+i+ 4); \
|
||||
__m128 XMM2 = _mm_load_ps((x)+i+ 8); \
|
||||
__m128 XMM3 = _mm_load_ps((x)+i+12); \
|
||||
__m128 XMM4 = _mm_load_ps((y)+i ); \
|
||||
__m128 XMM5 = _mm_load_ps((y)+i+ 4); \
|
||||
__m128 XMM6 = _mm_load_ps((y)+i+ 8); \
|
||||
__m128 XMM7 = _mm_load_ps((y)+i+12); \
|
||||
XMM0 = _mm_sub_ps(XMM0, XMM4); \
|
||||
XMM1 = _mm_sub_ps(XMM1, XMM5); \
|
||||
XMM2 = _mm_sub_ps(XMM2, XMM6); \
|
||||
XMM3 = _mm_sub_ps(XMM3, XMM7); \
|
||||
_mm_store_ps((z)+i , XMM0); \
|
||||
_mm_store_ps((z)+i+ 4, XMM1); \
|
||||
_mm_store_ps((z)+i+ 8, XMM2); \
|
||||
_mm_store_ps((z)+i+12, XMM3); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecscale(y, c, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM7 = _mm_set_ps1(c); \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
__m128 XMM0 = _mm_load_ps((y)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((y)+i+4); \
|
||||
XMM0 = _mm_mul_ps(XMM0, XMM7); \
|
||||
XMM1 = _mm_mul_ps(XMM1, XMM7); \
|
||||
_mm_store_ps((y)+i , XMM0); \
|
||||
_mm_store_ps((y)+i+4, XMM1); \
|
||||
} \
|
||||
}
|
||||
|
||||
#define vecmul(y, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
__m128 XMM0 = _mm_load_ps((x)+i ); \
|
||||
__m128 XMM1 = _mm_load_ps((x)+i+ 4); \
|
||||
__m128 XMM2 = _mm_load_ps((x)+i+ 8); \
|
||||
__m128 XMM3 = _mm_load_ps((x)+i+12); \
|
||||
__m128 XMM4 = _mm_load_ps((y)+i ); \
|
||||
__m128 XMM5 = _mm_load_ps((y)+i+ 4); \
|
||||
__m128 XMM6 = _mm_load_ps((y)+i+ 8); \
|
||||
__m128 XMM7 = _mm_load_ps((y)+i+12); \
|
||||
XMM4 = _mm_mul_ps(XMM4, XMM0); \
|
||||
XMM5 = _mm_mul_ps(XMM5, XMM1); \
|
||||
XMM6 = _mm_mul_ps(XMM6, XMM2); \
|
||||
XMM7 = _mm_mul_ps(XMM7, XMM3); \
|
||||
_mm_store_ps((y)+i , XMM4); \
|
||||
_mm_store_ps((y)+i+ 4, XMM5); \
|
||||
_mm_store_ps((y)+i+ 8, XMM6); \
|
||||
_mm_store_ps((y)+i+12, XMM7); \
|
||||
} \
|
||||
}
|
||||
|
||||
|
||||
|
||||
#if 3 <= __SSE__ || defined(__SSE3__)
|
||||
/*
|
||||
Horizontal add with haddps SSE3 instruction. The work register (rw)
|
||||
is unused.
|
||||
*/
|
||||
#define __horizontal_sum(r, rw) \
|
||||
r = _mm_hadd_ps(r, r); \
|
||||
r = _mm_hadd_ps(r, r);
|
||||
|
||||
#else
|
||||
/*
|
||||
Horizontal add with SSE instruction. The work register (rw) is used.
|
||||
*/
|
||||
#define __horizontal_sum(r, rw) \
|
||||
rw = r; \
|
||||
r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(1, 0, 3, 2)); \
|
||||
r = _mm_add_ps(r, rw); \
|
||||
rw = r; \
|
||||
r = _mm_shuffle_ps(r, rw, _MM_SHUFFLE(2, 3, 0, 1)); \
|
||||
r = _mm_add_ps(r, rw);
|
||||
|
||||
#endif
|
||||
|
||||
#define vecdot(s, x, y, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM0 = _mm_setzero_ps(); \
|
||||
__m128 XMM1 = _mm_setzero_ps(); \
|
||||
__m128 XMM2, XMM3, XMM4, XMM5; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
XMM2 = _mm_load_ps((x)+i ); \
|
||||
XMM3 = _mm_load_ps((x)+i+4); \
|
||||
XMM4 = _mm_load_ps((y)+i ); \
|
||||
XMM5 = _mm_load_ps((y)+i+4); \
|
||||
XMM2 = _mm_mul_ps(XMM2, XMM4); \
|
||||
XMM3 = _mm_mul_ps(XMM3, XMM5); \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_ps(XMM1, XMM3); \
|
||||
} \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM1); \
|
||||
__horizontal_sum(XMM0, XMM1); \
|
||||
_mm_store_ss((s), XMM0); \
|
||||
}
|
||||
|
||||
#define vec2norm(s, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM0 = _mm_setzero_ps(); \
|
||||
__m128 XMM1 = _mm_setzero_ps(); \
|
||||
__m128 XMM2, XMM3; \
|
||||
for (i = 0;i < (n);i += 8) { \
|
||||
XMM2 = _mm_load_ps((x)+i ); \
|
||||
XMM3 = _mm_load_ps((x)+i+4); \
|
||||
XMM2 = _mm_mul_ps(XMM2, XMM2); \
|
||||
XMM3 = _mm_mul_ps(XMM3, XMM3); \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_ps(XMM1, XMM3); \
|
||||
} \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM1); \
|
||||
__horizontal_sum(XMM0, XMM1); \
|
||||
XMM2 = XMM0; \
|
||||
XMM1 = _mm_rsqrt_ss(XMM0); \
|
||||
XMM3 = XMM1; \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM1); \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM3); \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM0); \
|
||||
XMM1 = _mm_mul_ss(XMM1, _mm_set_ss(-0.5f)); \
|
||||
XMM3 = _mm_mul_ss(XMM3, _mm_set_ss(1.5f)); \
|
||||
XMM3 = _mm_add_ss(XMM3, XMM1); \
|
||||
XMM3 = _mm_mul_ss(XMM3, XMM2); \
|
||||
_mm_store_ss((s), XMM3); \
|
||||
}
|
||||
|
||||
#define vec2norminv(s, x, n) \
|
||||
{ \
|
||||
int i; \
|
||||
__m128 XMM0 = _mm_setzero_ps(); \
|
||||
__m128 XMM1 = _mm_setzero_ps(); \
|
||||
__m128 XMM2, XMM3; \
|
||||
for (i = 0;i < (n);i += 16) { \
|
||||
XMM2 = _mm_load_ps((x)+i ); \
|
||||
XMM3 = _mm_load_ps((x)+i+4); \
|
||||
XMM2 = _mm_mul_ps(XMM2, XMM2); \
|
||||
XMM3 = _mm_mul_ps(XMM3, XMM3); \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM2); \
|
||||
XMM1 = _mm_add_ps(XMM1, XMM3); \
|
||||
} \
|
||||
XMM0 = _mm_add_ps(XMM0, XMM1); \
|
||||
__horizontal_sum(XMM0, XMM1); \
|
||||
XMM2 = XMM0; \
|
||||
XMM1 = _mm_rsqrt_ss(XMM0); \
|
||||
XMM3 = XMM1; \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM1); \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM3); \
|
||||
XMM1 = _mm_mul_ss(XMM1, XMM0); \
|
||||
XMM1 = _mm_mul_ss(XMM1, _mm_set_ss(-0.5f)); \
|
||||
XMM3 = _mm_mul_ss(XMM3, _mm_set_ss(1.5f)); \
|
||||
XMM3 = _mm_add_ss(XMM3, XMM1); \
|
||||
_mm_store_ss((s), XMM3); \
|
||||
}
|
||||
+152
@@ -0,0 +1,152 @@
|
||||
/* gss.c
|
||||
*
|
||||
* Copyright (C) 2012 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include "plfit_error.h"
|
||||
#include "gss.h"
|
||||
|
||||
/**
|
||||
* \def PHI
|
||||
*
|
||||
* The golden ratio, i.e. 1+sqrt(5)/2
|
||||
*/
|
||||
#define PHI 1.618033988749895
|
||||
|
||||
/**
|
||||
* \def RESPHI
|
||||
*
|
||||
* Constant defined as 2 - \c PHI
|
||||
*/
|
||||
#define RESPHI 0.3819660112501051
|
||||
|
||||
/**
|
||||
* \const _defparam
|
||||
*
|
||||
* Default parameters for the GSS algorithm.
|
||||
*/
|
||||
static const gss_parameter_t _defparam = {
|
||||
/* .epsilon = */ DBL_MIN,
|
||||
/* .on_error = */ GSS_ERROR_STOP
|
||||
};
|
||||
|
||||
/**
|
||||
* Stores whether the last optimization run triggered a warning or not.
|
||||
*/
|
||||
static unsigned short int gss_i_warning_flag = 0;
|
||||
|
||||
void gss_parameter_init(gss_parameter_t *param) {
|
||||
memcpy(param, &_defparam, sizeof(*param));
|
||||
}
|
||||
|
||||
unsigned short int gss_get_warning_flag(void) {
|
||||
return gss_i_warning_flag;
|
||||
}
|
||||
|
||||
#define TERMINATE { \
|
||||
if (_min) { \
|
||||
*(_min) = min; \
|
||||
} \
|
||||
if (_fmin) { \
|
||||
*(_fmin) = fmin; \
|
||||
} \
|
||||
}
|
||||
|
||||
#define EVALUATE(x, fx) { \
|
||||
fx = proc_evaluate(instance, x); \
|
||||
if (fmin > fx) { \
|
||||
min = x; \
|
||||
fmin = fx; \
|
||||
} \
|
||||
if (proc_progress) { \
|
||||
retval = proc_progress(instance, x, fx, min, fmin, \
|
||||
(a < b) ? a : b, (a < b) ? b : a, k); \
|
||||
if (retval) { \
|
||||
TERMINATE; \
|
||||
return PLFIT_SUCCESS; \
|
||||
} \
|
||||
} \
|
||||
}
|
||||
|
||||
int gss(double a, double b, double *_min, double *_fmin,
|
||||
gss_evaluate_t proc_evaluate, gss_progress_t proc_progress,
|
||||
void* instance, const gss_parameter_t *_param) {
|
||||
double c, d, min;
|
||||
double fa, fb, fc, fd, fmin;
|
||||
int k = 0;
|
||||
int retval;
|
||||
unsigned short int successful = 1;
|
||||
|
||||
gss_parameter_t param = _param ? (*_param) : _defparam;
|
||||
|
||||
gss_i_warning_flag = 0;
|
||||
|
||||
if (a > b) {
|
||||
c = a; a = b; b = c;
|
||||
}
|
||||
|
||||
min = a;
|
||||
fmin = proc_evaluate(instance, a);
|
||||
|
||||
c = a + RESPHI*(b-a);
|
||||
|
||||
EVALUATE(a, fa);
|
||||
EVALUATE(b, fb);
|
||||
EVALUATE(c, fc);
|
||||
|
||||
if (fc >= fa || fc >= fb) {
|
||||
if (param.on_error == GSS_ERROR_STOP) {
|
||||
return PLFIT_FAILURE;
|
||||
} else {
|
||||
gss_i_warning_flag = 1;
|
||||
}
|
||||
}
|
||||
|
||||
while (fabs(a-b) > param.epsilon) {
|
||||
k++;
|
||||
|
||||
d = c + RESPHI*(b-c);
|
||||
EVALUATE(d, fd);
|
||||
|
||||
if (fd >= fa || fd >= fb) {
|
||||
if (param.on_error == GSS_ERROR_STOP) {
|
||||
successful = 0;
|
||||
break;
|
||||
} else {
|
||||
gss_i_warning_flag = 1;
|
||||
}
|
||||
}
|
||||
|
||||
if (fc <= fd) {
|
||||
b = a; a = d;
|
||||
} else {
|
||||
a = c; c = d; fc = fd;
|
||||
}
|
||||
}
|
||||
|
||||
if (successful) {
|
||||
c = (a+b) / 2.0;
|
||||
k++;
|
||||
EVALUATE(c, fc);
|
||||
TERMINATE;
|
||||
}
|
||||
|
||||
return successful ? PLFIT_SUCCESS : PLFIT_FAILURE;
|
||||
}
|
||||
+138
@@ -0,0 +1,138 @@
|
||||
/* gss.h
|
||||
*
|
||||
* Copyright (C) 2012 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef __GSS_H__
|
||||
#define __GSS_H__
|
||||
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
/**
|
||||
* Enum specifying what the search should do when the function is not U-shaped.
|
||||
*/
|
||||
typedef enum {
|
||||
GSS_ERROR_STOP, /**< Stop and return an error code */
|
||||
GSS_ERROR_WARN /**< Continue and set the warning flag */
|
||||
} gss_error_handling_t;
|
||||
|
||||
/**
|
||||
* Parameter settings for a golden section search.
|
||||
*/
|
||||
typedef struct {
|
||||
double epsilon;
|
||||
gss_error_handling_t on_error;
|
||||
} gss_parameter_t;
|
||||
|
||||
/**
|
||||
* Callback interface to provide objective function evaluations for the golden
|
||||
* section search.
|
||||
*
|
||||
* The gss() function calls this function to obtain the values of the objective
|
||||
* function when needed. A client program must implement this function to evaluate
|
||||
* the value of the objective function, given the location.
|
||||
*
|
||||
* @param instance The user data sent for the gss() function by the client.
|
||||
* @param x The current value of the variable.
|
||||
* @retval double The value of the objective function for the current
|
||||
* variable.
|
||||
*/
|
||||
typedef double (*gss_evaluate_t)(void *instance, double x);
|
||||
|
||||
/**
|
||||
* Callback interface to receive the progress of the optimization process for
|
||||
* the golden section search.
|
||||
*
|
||||
* The gss() function calls this function for each iteration. Implementing
|
||||
* this function, a client program can store or display the current progress
|
||||
* of the optimization process.
|
||||
*
|
||||
* @param instance The user data sent for the gss() function by the client.
|
||||
* @param x The current value of the variable.
|
||||
* @param fx The value of the objective function at x.
|
||||
* @param min The location of the minimum value of the objective
|
||||
* function found so far.
|
||||
* @param fmin The minimum value of the objective function found so far.
|
||||
* @param left The left side of the current bracket.
|
||||
* @param right The right side of the current bracket.
|
||||
* @param k The index of the current iteration.
|
||||
* @retval int Zero to continue the optimization process. Returning a
|
||||
* non-zero value will cancel the optimization process.
|
||||
*/
|
||||
typedef int (*gss_progress_t)(void *instance, double x, double fx, double min,
|
||||
double fmin, double left, double right, int k);
|
||||
|
||||
/**
|
||||
* Start a golden section search optimization.
|
||||
*
|
||||
* @param a The left side of the bracket to start from
|
||||
* @param b The right side of the bracket to start from
|
||||
* @param min The pointer to the variable that receives the location of the
|
||||
* final value of the objective function. This argument can be set to
|
||||
* \c NULL if the location of the final value of the objective
|
||||
* function is unnecessary.
|
||||
* @param fmin The pointer to the variable that receives the final value of
|
||||
* the objective function. This argument can be st to \c NULL if the
|
||||
* final value of the objective function is unnecessary.
|
||||
* @param proc_evaluate The callback function to evaluate the objective
|
||||
* function at a given location.
|
||||
* @param proc_progress The callback function to receive the progress (the
|
||||
* last evaluated location, the value of the objective
|
||||
* function at that location, the width of the current
|
||||
* bracket, the minimum found so far and the step
|
||||
* count). This argument can be set to \c NULL if
|
||||
* a progress report is unnecessary.
|
||||
* @param instance A user data for the client program. The callback
|
||||
* functions will receive the value of this argument.
|
||||
* @param param The pointer to a structure representing parameters for
|
||||
* GSS algorithm. A client program can set this parameter
|
||||
* to \c NULL to use the default parameters.
|
||||
* Call the \ref gss_parameter_init() function to fill a
|
||||
* structure with the default values.
|
||||
* @retval int The status code. This function returns zero if the
|
||||
* minimization process terminates without an error. A
|
||||
* non-zero value indicates an error; in particular,
|
||||
* \c PLFIT_FAILURE means that the function is not
|
||||
* U-shaped.
|
||||
*/
|
||||
int gss(double a, double b, double *min, double *fmin,
|
||||
gss_evaluate_t proc_evaluate, gss_progress_t proc_progress,
|
||||
void* instance, const gss_parameter_t *_param);
|
||||
|
||||
/**
|
||||
* Return the state of the warning flag.
|
||||
*
|
||||
* The warning flag is 1 if the last optimization was run on a function that
|
||||
* was not U-shaped.
|
||||
*/
|
||||
unsigned short int gss_get_warning_flag(void);
|
||||
|
||||
/**
|
||||
* Initialize GSS parameters to the default values.
|
||||
*
|
||||
* Call this function to fill a parameter structure with the default values
|
||||
* and overwrite parameter values if necessary.
|
||||
*
|
||||
* @param param The pointer to the parameter structure.
|
||||
*/
|
||||
void gss_parameter_init(gss_parameter_t *param);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif /* __GSS_H__ */
|
||||
+672
@@ -0,0 +1,672 @@
|
||||
/* vim:set ts=4 sw=2 sts=2 et: */
|
||||
|
||||
/* This file was imported from a private scientific library
|
||||
* based on GSL coined Home Scientific Libray (HSL) by its author
|
||||
* Jerome Benoit; this very material is itself inspired from the
|
||||
* material written by G. Jungan and distributed by GSL.
|
||||
* Ultimately, some modifications were done in order to render the
|
||||
* imported material independent from the rest of GSL.
|
||||
*/
|
||||
|
||||
/* `hsl/specfunc/hzeta.c' C source file
|
||||
// HSL - Home Scientific Library
|
||||
// Copyright (C) 2017-2022 Jerome Benoit
|
||||
//
|
||||
// HSL 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.
|
||||
*/
|
||||
|
||||
/*
|
||||
// The material in this file is mainly inspired by the material written by
|
||||
// G. Jungan and distributed under GPLv2 by the GNU Scientific Library (GSL)
|
||||
// ( https://www.gnu.org/software/gsl/ [specfunc/zeta.c]), itself inspired by
|
||||
// the material written by Moshier and distributed in the Cephes Mathematical
|
||||
// Library ( http://www.moshier.net/ [zeta.c]).
|
||||
//
|
||||
// More specifically, hsl_sf_hzeta_e is a slightly modifed clone of
|
||||
// gsl_sf_hzeta_e as found in GSL 2.4; the remaining is `inspired by'.
|
||||
// [Sooner or later a _Working_Note_ may be deposited at ResearchGate
|
||||
// ( https://www.researchgate.net/profile/Jerome_Benoit )]
|
||||
*/
|
||||
|
||||
/* Author: Jerome G. Benoit < jgmbenoit _at_ rezozer _dot_ net > */
|
||||
|
||||
#ifdef _MSC_VER
|
||||
#define _USE_MATH_DEFINES
|
||||
#endif
|
||||
|
||||
/* Work around bug in some Windows SDK / MSVC versions where NAN is not a
|
||||
* constant expression, triggering an error in the definition of
|
||||
* hsl_sf_hzeta_eulermaclaurin_series_coeffs[] and
|
||||
* hsl_sf_hzeta_eulermaclaurin_series_majorantratios[] below.
|
||||
* We re-define NAN to the value it had in earlier MSVC versions.
|
||||
* See https://github.com/igraph/igraph/issues/2701
|
||||
* and https://developercommunity.visualstudio.com/t/NAN-is-no-longer-compile-time-constant-i/10688907
|
||||
*/
|
||||
#ifdef _MSC_VER
|
||||
#define _UCRT_NOISY_NAN
|
||||
#endif
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include "hzeta.h"
|
||||
#include "plfit_error.h"
|
||||
|
||||
/* imported from gsl_machine.h */
|
||||
|
||||
#define GSL_LOG_DBL_MIN (-7.0839641853226408e+02)
|
||||
#define GSL_LOG_DBL_MAX 7.0978271289338397e+02
|
||||
#define GSL_DBL_EPSILON 2.2204460492503131e-16
|
||||
|
||||
/* Math constants are not part of standard C.
|
||||
* The following are borrowed from igraph/src/core/math.h */
|
||||
|
||||
#ifndef M_LOG2E
|
||||
#define M_LOG2E 1.44269504088896340735992468100189214
|
||||
#endif
|
||||
|
||||
#ifndef M_LN2
|
||||
#define M_LN2 0.693147180559945309417232121458176568
|
||||
#endif
|
||||
|
||||
/* imported from gsl_sf_result.h */
|
||||
|
||||
struct gsl_sf_result_struct {
|
||||
double val;
|
||||
double err;
|
||||
};
|
||||
typedef struct gsl_sf_result_struct gsl_sf_result;
|
||||
|
||||
/* imported and adapted from hsl/specfunc/specfunc_def.h */
|
||||
|
||||
#define HSL_SF_EVAL_RESULT(FnE) \
|
||||
gsl_sf_result result; \
|
||||
FnE ; \
|
||||
return (result.val);
|
||||
|
||||
#define HSL_SF_EVAL_TUPLE_RESULT(FnET) \
|
||||
gsl_sf_result result0; \
|
||||
gsl_sf_result result1; \
|
||||
FnET ; \
|
||||
*tuple1=result1.val; \
|
||||
*tuple0=result0.val; \
|
||||
return (result0.val);
|
||||
|
||||
/* */
|
||||
|
||||
|
||||
#define HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT 10
|
||||
#define HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER 32
|
||||
|
||||
#define HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX 256
|
||||
|
||||
// B_{2j}/(2j)
|
||||
static
|
||||
double hsl_sf_hzeta_eulermaclaurin_series_coeffs[HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+2]={
|
||||
+1.0,
|
||||
+1.0/12.0,
|
||||
-1.0/720.0,
|
||||
+1.0/30240.0,
|
||||
-1.0/1209600.0,
|
||||
+1.0/47900160.0,
|
||||
-691.0/1307674368000.0,
|
||||
+1.0/74724249600.0,
|
||||
-3.38968029632258286683019539125e-13,
|
||||
+8.58606205627784456413590545043e-15,
|
||||
-2.17486869855806187304151642387e-16,
|
||||
+5.50900282836022951520265260890e-18,
|
||||
-1.39544646858125233407076862641e-19,
|
||||
+3.53470703962946747169322997780e-21,
|
||||
-8.95351742703754685040261131811e-23,
|
||||
+2.26795245233768306031095073887e-24,
|
||||
-5.74479066887220244526388198761e-26,
|
||||
+1.45517247561486490186626486727e-27,
|
||||
-3.68599494066531017818178247991e-29,
|
||||
+9.33673425709504467203255515279e-31,
|
||||
-2.36502241570062993455963519637e-32,
|
||||
+5.99067176248213430465991239682e-34,
|
||||
-1.51745488446829026171081313586e-35,
|
||||
+3.84375812545418823222944529099e-37,
|
||||
-9.73635307264669103526762127925e-39,
|
||||
+2.46624704420068095710640028029e-40,
|
||||
-6.24707674182074369314875679472e-42,
|
||||
+1.58240302446449142975108170683e-43,
|
||||
-4.00827368594893596853001219052e-45,
|
||||
+1.01530758555695563116307139454e-46,
|
||||
-2.57180415824187174992481940976e-48,
|
||||
+6.51445603523381493155843485864e-50,
|
||||
-1.65013099068965245550609878048e-51,
|
||||
NAN}; // hsl_sf_hzeta_eulermaclaurin_series_coeffs
|
||||
|
||||
// 4\zeta(2j)/(2\pi)^(2j)
|
||||
static
|
||||
double hsl_sf_hzeta_eulermaclaurin_series_majorantratios[HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+2]={
|
||||
-2.0,
|
||||
+1.0/6.0,
|
||||
+1.0/360.0,
|
||||
+1.0/15120.0,
|
||||
+1.0/604800.0,
|
||||
+1.0/23950080.0,
|
||||
+691.0/653837184000.0,
|
||||
+1.0/37362124800.0,
|
||||
+3617.0/5335311421440000.0,
|
||||
+1.71721241125556891282718109009e-14,
|
||||
+4.34973739711612374608303284773e-16,
|
||||
+1.10180056567204590304053052178e-17,
|
||||
+2.79089293716250466814153725281e-19,
|
||||
+7.06941407925893494338645995561e-21,
|
||||
+1.79070348540750937008052226362e-22,
|
||||
+4.53590490467536612062190147774e-24,
|
||||
+1.14895813377444048905277639752e-25,
|
||||
+2.91034495122972980373252973454e-27,
|
||||
+7.37198988133062035636356495982e-29,
|
||||
+1.86734685141900893440651103056e-30,
|
||||
+4.73004483140125986911927039274e-32,
|
||||
+1.19813435249642686093198247936e-33,
|
||||
+3.03490976893658052342162627173e-35,
|
||||
+7.68751625090837646445889058198e-37,
|
||||
+1.94727061452933820705352425585e-38,
|
||||
+4.93249408840136191421280056051e-40,
|
||||
+1.24941534836414873862975135893e-41,
|
||||
+3.16480604892898285950216341362e-43,
|
||||
+8.01654737189787193706002438098e-45,
|
||||
+2.03061517111391126232614278906e-46,
|
||||
+5.14360831648374349984963881946e-48,
|
||||
+1.30289120704676298631168697172e-49,
|
||||
+3.30026198137930491101219756091e-51,
|
||||
NAN}; // hsl_sf_hzeta_eulermaclaurin_series_majorantratios
|
||||
|
||||
|
||||
extern
|
||||
int hsl_sf_hzeta_e(const double s, const double q, gsl_sf_result * result) {
|
||||
|
||||
/* CHECK_POINTER(result) */
|
||||
|
||||
if ((s <= 1.0) || (q <= 0.0)) {
|
||||
PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
|
||||
}
|
||||
else {
|
||||
const double max_bits=54.0; // max_bits=\lceil{s}\rceil with \zeta(s,2)=\zeta(s)-1=GSL_DBL_EPSILON
|
||||
const double ln_term0=-s*log(q);
|
||||
if (ln_term0 < GSL_LOG_DBL_MIN+1.0) {
|
||||
PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
|
||||
}
|
||||
else if (GSL_LOG_DBL_MAX-1.0 < ln_term0) {
|
||||
PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
|
||||
}
|
||||
#if 1
|
||||
else if (((max_bits < s) && (q < 1.0)) || ((0.5*max_bits < s) && (q < 0.25))) {
|
||||
result->val=pow(q,-s);
|
||||
result->err=2.0*GSL_DBL_EPSILON*fabs(result->val);
|
||||
return (PLFIT_SUCCESS);
|
||||
}
|
||||
else if ((0.5*max_bits < s) && (q < 1.0)) {
|
||||
const double a0=pow(q,-s);
|
||||
const double p1=pow(q/(1.0+q),s);
|
||||
const double p2=pow(q/(2.0+q),s);
|
||||
const double ans=a0*(1.0+p1+p2);
|
||||
result->val=ans;
|
||||
result->err=GSL_DBL_EPSILON*(2.0+0.5*s)*fabs(result->val);
|
||||
return (PLFIT_SUCCESS);
|
||||
}
|
||||
#endif
|
||||
else { // Euler-Maclaurin summation formula
|
||||
const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
|
||||
const double inv_qshift=1.0/qshift;
|
||||
const double sqr_inv_qshift=inv_qshift*inv_qshift;
|
||||
const double inv_sm1=1.0/(s-1.0);
|
||||
const double pmax=pow(qshift,-s);
|
||||
double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
|
||||
double delta=NAN;
|
||||
double tscp=s;
|
||||
double scp=tscp;
|
||||
double pcp=pmax*inv_qshift;
|
||||
double ratio=scp*pcp;
|
||||
size_t n=0;
|
||||
size_t j=0;
|
||||
double ans=0.0;
|
||||
double mjr=NAN;
|
||||
|
||||
for(j=0;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++j) ans+=(terms[n++]=pow(j+q,-s));
|
||||
ans+=(terms[n++]=0.5*pmax);
|
||||
ans+=(terms[n++]=pmax*qshift*inv_sm1);
|
||||
for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
|
||||
delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
|
||||
ans+=(terms[n++]=delta);
|
||||
scp*=++tscp;
|
||||
scp*=++tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
ratio=scp*pcp;
|
||||
if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
|
||||
}
|
||||
if (HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER<j) PLFIT_ERROR("maximum iterations exceeded",PLFIT_EMAXITER);
|
||||
|
||||
ans=0.0; while (n) ans+=terms[--n];
|
||||
mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
|
||||
|
||||
result->val=+ans;
|
||||
result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
|
||||
return (PLFIT_SUCCESS);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
extern
|
||||
double hsl_sf_hzeta(const double s, const double q) {
|
||||
HSL_SF_EVAL_RESULT(hsl_sf_hzeta_e(s,q,&result)); }
|
||||
|
||||
extern
|
||||
int hsl_sf_hzeta_deriv_e(const double s, const double q, gsl_sf_result * result) {
|
||||
|
||||
/* CHECK_POINTER(result) */
|
||||
|
||||
if ((s <= 1.0) || (q <= 0.0)) {
|
||||
PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
|
||||
}
|
||||
else {
|
||||
const double ln_hz_term0=-s*log(q);
|
||||
if (ln_hz_term0 < GSL_LOG_DBL_MIN+1.0) {
|
||||
PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
|
||||
}
|
||||
else if (GSL_LOG_DBL_MAX-1.0 < ln_hz_term0) {
|
||||
PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
|
||||
}
|
||||
else { // Euler-Maclaurin summation formula
|
||||
const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
|
||||
const double inv_qshift=1.0/qshift;
|
||||
const double sqr_inv_qshift=inv_qshift*inv_qshift;
|
||||
const double inv_sm1=1.0/(s-1.0);
|
||||
const double pmax=pow(qshift,-s);
|
||||
const double lmax=log(qshift);
|
||||
double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
|
||||
double delta=NAN;
|
||||
double tscp=s;
|
||||
double scp=tscp;
|
||||
double pcp=pmax*inv_qshift;
|
||||
double lcp=lmax-1.0/s;
|
||||
double ratio=scp*pcp*lcp;
|
||||
double qs=NAN;
|
||||
size_t n=0;
|
||||
size_t j=0;
|
||||
double ans=0.0;
|
||||
double mjr=NAN;
|
||||
|
||||
for(j=0,qs=q;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++qs,++j) ans+=(terms[n++]=log(qs)*pow(qs,-s));
|
||||
ans+=(terms[n++]=0.5*lmax*pmax);
|
||||
ans+=(terms[n++]=pmax*qshift*inv_sm1*(lmax+inv_sm1));
|
||||
for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
|
||||
delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
|
||||
ans+=(terms[n++]=delta);
|
||||
scp*=++tscp; lcp-=1.0/tscp;
|
||||
scp*=++tscp; lcp-=1.0/tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
ratio=scp*pcp*lcp;
|
||||
if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
|
||||
}
|
||||
if (HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER<j) PLFIT_ERROR("maximum iterations exceeded",PLFIT_EMAXITER);
|
||||
|
||||
ans=0.0; while (n) ans+=terms[--n];
|
||||
mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
|
||||
|
||||
result->val=-ans;
|
||||
result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
|
||||
return (PLFIT_SUCCESS);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
extern
|
||||
double hsl_sf_hzeta_deriv(const double s, const double q) {
|
||||
HSL_SF_EVAL_RESULT(hsl_sf_hzeta_deriv_e(s,q,&result)); }
|
||||
|
||||
extern
|
||||
int hsl_sf_hzeta_deriv2_e(const double s, const double q, gsl_sf_result * result) {
|
||||
|
||||
/* CHECK_POINTER(result) */
|
||||
|
||||
if ((s <= 1.0) || (q <= 0.0)) {
|
||||
PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
|
||||
}
|
||||
else {
|
||||
const double ln_hz_term0=-s*log(q);
|
||||
if (ln_hz_term0 < GSL_LOG_DBL_MIN+1.0) {
|
||||
PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
|
||||
}
|
||||
else if (GSL_LOG_DBL_MAX-1.0 < ln_hz_term0) {
|
||||
PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
|
||||
}
|
||||
else { // Euler-Maclaurin summation formula
|
||||
const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
|
||||
const double inv_qshift=1.0/qshift;
|
||||
const double sqr_inv_qshift=inv_qshift*inv_qshift;
|
||||
const double inv_sm1=1.0/(s-1.0);
|
||||
const double pmax=pow(qshift,-s);
|
||||
const double lmax=log(qshift);
|
||||
const double lmax_p_inv_sm1=lmax+inv_sm1;
|
||||
const double sqr_inv_sm1=inv_sm1*inv_sm1;
|
||||
const double sqr_lmax=lmax*lmax;
|
||||
const double sqr_lmax_p_inv_sm1=lmax_p_inv_sm1*lmax_p_inv_sm1;
|
||||
double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
|
||||
double delta=NAN;
|
||||
double tscp=s;
|
||||
double slcp=NAN;
|
||||
double plcp=NAN;
|
||||
double scp=tscp;
|
||||
double pcp=pmax*inv_qshift;
|
||||
double lcp=1.0/s-lmax;
|
||||
double sqr_lcp=lmax*(lmax-2.0/s);
|
||||
double ratio=scp*pcp*sqr_lcp;
|
||||
double qs=NAN;
|
||||
double lqs=NAN;
|
||||
size_t n=0;
|
||||
size_t j=0;
|
||||
double ans=0.0;
|
||||
double mjr=NAN;
|
||||
|
||||
for(j=0,qs=q;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++qs,++j) {
|
||||
lqs=log(qs);
|
||||
ans+=(terms[n++]=lqs*lqs*pow(qs,-s));
|
||||
}
|
||||
ans+=(terms[n++]=0.5*sqr_lmax*pmax);
|
||||
ans+=(terms[n++]=pmax*qshift*inv_sm1*(sqr_lmax_p_inv_sm1+sqr_inv_sm1));
|
||||
for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
|
||||
delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
|
||||
ans+=(terms[n++]=delta);
|
||||
scp*=++tscp; slcp=plcp=1.0/tscp;
|
||||
scp*=++tscp; slcp+=1.0/tscp; plcp/=tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
sqr_lcp+=2.0*(plcp+slcp*lcp);
|
||||
ratio=scp*pcp*sqr_lcp;
|
||||
if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
|
||||
lcp+=slcp;
|
||||
}
|
||||
if (HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER<j) PLFIT_ERROR("maximum iterations exceeded",PLFIT_EMAXITER);
|
||||
|
||||
ans=0.0; while (n) ans+=terms[--n];
|
||||
mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
|
||||
|
||||
result->val=+ans;
|
||||
result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
|
||||
return (PLFIT_SUCCESS);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
extern
|
||||
double hsl_sf_hzeta_deriv2(const double s, const double q) {
|
||||
HSL_SF_EVAL_RESULT(hsl_sf_hzeta_deriv2_e(s,q,&result)); }
|
||||
|
||||
static inline
|
||||
double hsl_sf_hZeta0_zed(const double s, const double q) {
|
||||
#if 1
|
||||
const long double ld_q=(long double)(q);
|
||||
const long double ld_s=(long double)(s);
|
||||
const long double ld_log1prq=log1pl(1.0L/ld_q);
|
||||
const long double ld_epsilon=expm1l(-ld_s*ld_log1prq);
|
||||
const long double ld_z=ld_s+(ld_q+0.5L*ld_s+0.5L)*ld_epsilon;
|
||||
const double z=(double)(ld_z);
|
||||
#else
|
||||
double z=s+(q+0.5*s+0.5)*expm1(-s*log1p(1.0/q));
|
||||
#endif
|
||||
return (z); }
|
||||
|
||||
// Z_{0}(s,a) = a^s \left(\frac{1}{2}+\frac{a}{s-1}\right)^{-1} \zeta(s,a) - 1
|
||||
// Z_{0}(s,a) = O\left(\frac{(s-1)s}{6a^{2}}\right)
|
||||
static
|
||||
int hsl_sf_hZeta0(const double s, const double q, double * value, double * abserror) {
|
||||
const double criterion=ceil(10.0*s-q);
|
||||
const size_t shift=(criterion<0.0)?0:
|
||||
(criterion<HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX)?(size_t)(llrint(criterion)):
|
||||
HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX;
|
||||
const double qshift=(double)(shift)+q;
|
||||
const double inv_qshift=1.0/qshift;
|
||||
const double sqr_inv_qshift=inv_qshift*inv_qshift;
|
||||
const double sm1=s-1.0;
|
||||
double terms[HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
|
||||
double delta=NAN;
|
||||
double tscp=s;
|
||||
double scp=s*sm1;
|
||||
double pcp=inv_qshift/(2.0*qshift+sm1);
|
||||
double ratio=NAN;
|
||||
size_t n=0;
|
||||
size_t j=0;
|
||||
double ans=0.0;
|
||||
double mjr=NAN;
|
||||
|
||||
if (shift) {
|
||||
const double hsm1=0.5*sm1;
|
||||
const double inv_q=1.0/q;
|
||||
const double qphsm1=q+hsm1;
|
||||
const double inv_qphsm1=1.0/qphsm1;
|
||||
const double qshiftphsm1=qshift+hsm1;
|
||||
double qs=q;
|
||||
double a=1.0;
|
||||
for(j=0;j<shift;) {
|
||||
ans+=(terms[n++]=a*hsl_sf_hZeta0_zed(s,qs++)*inv_qphsm1);
|
||||
a=exp(-s*log1p((++j)*inv_q));
|
||||
}
|
||||
pcp*=a*qshiftphsm1*inv_qphsm1;
|
||||
}
|
||||
ratio=scp*pcp;
|
||||
ans+=terms[n++]=ratio/6.0;
|
||||
scp*=++tscp;
|
||||
scp*=++tscp;
|
||||
pcp*=2.0*sqr_inv_qshift;
|
||||
ratio=scp*pcp;
|
||||
for(j=2;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
|
||||
delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
|
||||
ans+=(terms[n++]=delta);
|
||||
scp*=++tscp;
|
||||
scp*=++tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
ratio=scp*pcp;
|
||||
if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
|
||||
}
|
||||
if (HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER<j) PLFIT_ERROR("maximum iterations exceeded",PLFIT_EMAXITER);
|
||||
|
||||
ans=0.0; while (n) ans+=terms[--n];
|
||||
mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
|
||||
|
||||
*value=ans;
|
||||
*abserror=2.0*((shift+1)*GSL_DBL_EPSILON*fabs(ans)+mjr);
|
||||
|
||||
return (PLFIT_SUCCESS); }
|
||||
|
||||
static inline
|
||||
double hsl_sf_hZeta1_zed(const double s, const double q) {
|
||||
#if 1
|
||||
const long double ld_q=(long double)(q);
|
||||
const long double ld_s=(long double)(s);
|
||||
const long double ld_sm1=ld_s-1.0L;
|
||||
const long double ld_logq=logl(ld_q);
|
||||
const long double ld_log1prq=log1pl(1.0L/ld_q);
|
||||
const long double ld_inv_logq=1.0L/ld_logq;
|
||||
const long double ld_logratiom1=ld_log1prq*ld_inv_logq;
|
||||
const long double ld_powratiom1=expm1l(-ld_s*ld_log1prq);
|
||||
const long double ld_varepsilon=expm1l(-ld_sm1*ld_log1prq);
|
||||
const long double ld_epsilon=ld_logratiom1+ld_powratiom1+ld_logratiom1*ld_powratiom1;
|
||||
const long double ld_z=ld_s+(ld_q+0.5L*ld_s+0.5L)*ld_epsilon+ld_q/ld_sm1*ld_inv_logq*ld_varepsilon;
|
||||
const double z=(double)(ld_z);
|
||||
#else
|
||||
const double sm1=s-1.0;
|
||||
const double inv_ln_q=1.0/log(q);
|
||||
const double log1prq=log1p(1.0/q);
|
||||
const double logratiom1=log1prq*inv_ln_q;
|
||||
const double powratiom1=expm1(-s*log1prq);
|
||||
const double epsilon=logratiom1+powratiom1+logratiom1*powratiom1;
|
||||
const double z=s+(q+0.5*s+0.5)*epsilon+q/sm1*inv_ln_q*expm1(-sm1*log1prq);
|
||||
#endif
|
||||
return (z); }
|
||||
|
||||
// Z_{1}(s,a) = -\frac{a^s}{\ln(a)} \left(\frac{1}{2}+\frac{a}{s-1}\,\left[1+\frac{1}{(s-1)\,\ln(a)}\right]\right)^{-1} \zeta^{\prime}(s,a) - 1
|
||||
// Z_{1}(s,a) = O\left(\frac{(s-1)s}{6a^{2}}\right)
|
||||
static
|
||||
int hsl_sf_hZeta1(const double s, const double q, const double ln_q, double * value, double * abserror, double * coeff1) {
|
||||
const double criterion=ceil(10.0*s-q);
|
||||
const size_t shift=(criterion<0.0)?0:
|
||||
(criterion<HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX)?(size_t)(llrint(criterion)):
|
||||
HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX;
|
||||
const double qshift=(double)(shift)+q;
|
||||
const double ln_qshift=log(qshift);
|
||||
const double inv_qshift=1.0/qshift;
|
||||
const double inv_ln_q=1.0/ln_q;
|
||||
const double inv_ln_qshift=1.0/ln_qshift;
|
||||
const double sqr_inv_qshift=inv_qshift*inv_qshift;
|
||||
const double sm1=s-1.0;
|
||||
const double hsm1=0.5*sm1;
|
||||
const double q_over_ln_q=q*inv_ln_q;
|
||||
const double qshift_over_ln_qshift=qshift*inv_ln_qshift;
|
||||
const double qphsm1=q+hsm1;
|
||||
const double qshiftphsm1=qshift+hsm1;
|
||||
double terms[HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
|
||||
double delta=NAN;
|
||||
double tscp=s;
|
||||
double scp=s*sm1;
|
||||
double pcp=inv_qshift*sm1/(qshift_over_ln_qshift+sm1*qshiftphsm1);
|
||||
double lcp=1.0-inv_ln_qshift/s;
|
||||
double ratio=NAN;
|
||||
size_t n=0;
|
||||
size_t j=0;
|
||||
double ans=0.0;
|
||||
double mjr=NAN;
|
||||
|
||||
if (shift) {
|
||||
const double inv_q=1.0/q;
|
||||
const double inv_sm1=1.0/sm1;
|
||||
const double w=1.0+inv_sm1*inv_ln_q;
|
||||
const double wshift=1.0+inv_sm1*inv_ln_qshift;
|
||||
const double qwphsm1=q*w+hsm1;
|
||||
const double inv_qwphsm1=1.0/qwphsm1;
|
||||
const double qshiftwshiftphsm1=qshift*wshift+hsm1;
|
||||
double qs=q;
|
||||
double a=1.0;
|
||||
for(j=0;j<shift;) {
|
||||
ans+=(terms[n++]=a*hsl_sf_hZeta1_zed(s,qs++)*inv_qwphsm1);
|
||||
a=log(qs)*inv_ln_q*exp(-s*log1p((++j)*inv_q));
|
||||
}
|
||||
pcp*=a*qshiftwshiftphsm1*inv_qwphsm1;
|
||||
}
|
||||
ratio=scp*pcp*lcp;
|
||||
ans+=terms[n++]=ratio/12.0;
|
||||
scp*=++tscp; lcp-=inv_ln_qshift/tscp;
|
||||
scp*=++tscp; lcp-=inv_ln_qshift/tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
ratio=scp*pcp*lcp;
|
||||
for(j=2;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
|
||||
delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
|
||||
ans+=(terms[n++]=delta);
|
||||
scp*=++tscp; lcp-=inv_ln_qshift/tscp;
|
||||
scp*=++tscp; lcp-=inv_ln_qshift/tscp;
|
||||
pcp*=sqr_inv_qshift;
|
||||
ratio=scp*pcp*lcp;
|
||||
if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
|
||||
}
|
||||
if (HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER<j) PLFIT_ERROR("maximum iterations exceeded",PLFIT_EMAXITER);
|
||||
|
||||
ans=0.0; while (n) ans+=terms[--n];
|
||||
mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
|
||||
|
||||
*value=ans;
|
||||
*abserror=2.0*((shift+1)*GSL_DBL_EPSILON*fabs(ans)+mjr);
|
||||
|
||||
if (coeff1) *coeff1=1.0+q_over_ln_q/qphsm1/sm1;
|
||||
|
||||
return (PLFIT_SUCCESS); }
|
||||
|
||||
extern
|
||||
int hsl_sf_lnhzeta_deriv_tuple_e(const double s, const double q, gsl_sf_result * result, gsl_sf_result * result_deriv) {
|
||||
|
||||
/* CHECK_POINTER(result) */
|
||||
|
||||
if ((s <= 1.0) || (q <= 0.0)) {
|
||||
PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
|
||||
}
|
||||
else if (q == 1.0) {
|
||||
const double inv_sm1=1.0/(s-1.0);
|
||||
const double inv_qsm1=4.0*inv_sm1;
|
||||
const double hz_coeff0=exp2(s+1.0);
|
||||
const double hz_coeff1=1.0+inv_qsm1;
|
||||
double hZeta0_value=NAN;
|
||||
double hZeta0_abserror=NAN;
|
||||
hsl_sf_hZeta0(s,2.0,&hZeta0_value,&hZeta0_abserror);
|
||||
hZeta0_value+=1.0;
|
||||
if (result) {
|
||||
const double ln_hz_coeff=hz_coeff1/hz_coeff0;
|
||||
const double ln_hZeta0_value=ln_hz_coeff*hZeta0_value;
|
||||
result->val=log1p(ln_hZeta0_value);
|
||||
result->err=(2.0*GSL_DBL_EPSILON*ln_hz_coeff+hZeta0_abserror)/(1.0+ln_hZeta0_value);
|
||||
}
|
||||
|
||||
if (result_deriv) {
|
||||
const double ld_hz_coeff2=1.0+inv_sm1*M_LOG2E;
|
||||
const double ld_hz_coeff1=1.0+inv_qsm1*ld_hz_coeff2;
|
||||
double hZeta1_value=NAN;
|
||||
double hZeta1_abserror=NAN;
|
||||
hsl_sf_hZeta1(s,2.0,M_LN2,&hZeta1_value,&hZeta1_abserror,NULL);
|
||||
hZeta0_value*=hz_coeff1;
|
||||
hZeta0_value+=hz_coeff0;
|
||||
hZeta1_value+=1.0;
|
||||
hZeta1_value*=-M_LN2*ld_hz_coeff1;
|
||||
result_deriv->val=hZeta1_value/hZeta0_value;
|
||||
result_deriv->err=2.0*GSL_DBL_EPSILON*fabs(result_deriv->val)+(hZeta0_abserror+hZeta1_abserror);
|
||||
}
|
||||
}
|
||||
else {
|
||||
const double ln_q=log(q);
|
||||
double hZeta0_value=NAN;
|
||||
double hZeta0_abserror=NAN;
|
||||
hsl_sf_hZeta0(s,q,&hZeta0_value,&hZeta0_abserror);
|
||||
if (result) {
|
||||
const double ln_hz_term0=-s*ln_q;
|
||||
const double ln_hz_term1=log(0.5+q/(s-1.0));
|
||||
result->val=ln_hz_term0+ln_hz_term1+log1p(hZeta0_value);
|
||||
result->err=2.0*GSL_DBL_EPSILON*(fabs(ln_hz_term0)+fabs(ln_hz_term1))+hZeta0_abserror/(1.0+hZeta0_value);
|
||||
}
|
||||
if (result_deriv) {
|
||||
double hZeta1_value=NAN;
|
||||
double hZeta1_abserror=NAN;
|
||||
double ld_hz_coeff1=NAN;
|
||||
hsl_sf_hZeta1(s,q,ln_q,&hZeta1_value,&hZeta1_abserror,&ld_hz_coeff1);
|
||||
result_deriv->val=-ln_q*ld_hz_coeff1*(1.0+hZeta1_value)/(1.0+hZeta0_value);
|
||||
result_deriv->err=2.0*GSL_DBL_EPSILON*fabs(result_deriv->val)+(hZeta0_abserror+hZeta1_abserror);
|
||||
}
|
||||
}
|
||||
|
||||
return (PLFIT_SUCCESS); }
|
||||
|
||||
extern
|
||||
double hsl_sf_lnhzeta_deriv_tuple(const double s, const double q, double * tuple0, double * tuple1) {
|
||||
HSL_SF_EVAL_TUPLE_RESULT(hsl_sf_lnhzeta_deriv_tuple_e(s,q,&result0,&result1)); }
|
||||
|
||||
extern
|
||||
int hsl_sf_lnhzeta_e(const double s, const double q, gsl_sf_result * result) {
|
||||
return (hsl_sf_lnhzeta_deriv_tuple_e(s,q,result,NULL)); }
|
||||
|
||||
extern
|
||||
double hsl_sf_lnhzeta(const double s, const double q) {
|
||||
HSL_SF_EVAL_RESULT(hsl_sf_lnhzeta_e(s,q,&result)); }
|
||||
|
||||
extern
|
||||
int hsl_sf_lnhzeta_deriv_e(const double s, const double q, gsl_sf_result * result) {
|
||||
return (hsl_sf_lnhzeta_deriv_tuple_e(s,q,NULL,result)); }
|
||||
|
||||
extern
|
||||
double hsl_sf_lnhzeta_deriv(const double s, const double q) {
|
||||
HSL_SF_EVAL_RESULT(hsl_sf_lnhzeta_deriv_e(s,q,&result)); }
|
||||
|
||||
//
|
||||
// End of file `hsl/specfunc/hzeta.c'.
|
||||
+88
@@ -0,0 +1,88 @@
|
||||
/* This file was imported from a private scientific library
|
||||
* based on GSL coined Home Scientific Libray (HSL) by its author
|
||||
* Jerome Benoit; this very material is itself inspired from the
|
||||
* material written by G. Jungan and distributed by GSL.
|
||||
* Ultimately, some modifications were done in order to render the
|
||||
* imported material independent from the rest of GSL.
|
||||
*/
|
||||
|
||||
/* `hsl/hsl_sf_zeta.h' C header file
|
||||
// HSL - Home Scientific Library
|
||||
// Copyright (C) 2005-2018 Jerome Benoit
|
||||
//
|
||||
// HSL 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.
|
||||
*/
|
||||
|
||||
/* For futher details, see its source conterpart src/hzeta.c */
|
||||
|
||||
/* Author: Jerome G. Benoit < jgmbenoit _at_ rezozer _dot_ net > */
|
||||
|
||||
#ifndef __HZETA_H__
|
||||
#define __HZETA_H__
|
||||
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
|
||||
/* Hurwitz Zeta Function
|
||||
* zeta(s,q) = Sum[ (k+q)^(-s), {k,0,Infinity} ]
|
||||
*
|
||||
* s > 1.0, q > 0.0
|
||||
*/
|
||||
double hsl_sf_hzeta(const double s, const double q);
|
||||
|
||||
/* First Derivative of Hurwitz Zeta Function
|
||||
* zeta'(s,q) = - Sum[ Ln(k+q)/(k+q)^(s), {k,0,Infinity} ]
|
||||
*
|
||||
* s > 1.0, q > 0.0
|
||||
*/
|
||||
double hsl_sf_hzeta_deriv(const double s, const double q);
|
||||
|
||||
/* Second Derivative of Hurwitz Zeta Function
|
||||
* zeta''(s,q) = + Sum[ Ln(k+q)^2/(k+q)^(s), {k,0,Infinity} ]
|
||||
*
|
||||
* s > 1.0, q > 0.0
|
||||
*/
|
||||
double hsl_sf_hzeta_deriv2(const double s, const double q);
|
||||
|
||||
/* Logarithm of Hurwitz Zeta Function
|
||||
* lnzeta(s,q) = ln(zeta(s,q))
|
||||
*
|
||||
* s > 1.0, q > 0.0 (and q >> 1)
|
||||
*/
|
||||
double hsl_sf_lnhzeta(const double s, const double q);
|
||||
|
||||
/* Logarithmic Derivative of Hurwitz Zeta Function
|
||||
* lnzeta'(s,q) = zeta'(s,q)/zeta(s,q)
|
||||
*
|
||||
* s > 1.0, q > 0.0 (and q >> 1)
|
||||
*/
|
||||
double hsl_sf_lnhzeta_deriv(const double s, const double q);
|
||||
|
||||
/* Logarithm and Logarithmic Derivative of Hurwitz Zeta Function:
|
||||
* nonredundant computation version:
|
||||
* - lnzeta(s,q) and lnzeta'(s,q) are stored in *deriv0 and *deriv1, respectively;
|
||||
* - the return value and the value stored in *deriv0 are the same;
|
||||
* - deriv0 and deriv1 must be effective pointers, that is, not the NULL pointer.
|
||||
*
|
||||
* s > 1.0, q > 0.0 (and q >> 1)
|
||||
*/
|
||||
double hsl_sf_lnhzeta_deriv_tuple(const double s, const double q, double * deriv0, double * deriv1);
|
||||
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif // __HZETA_H__
|
||||
@@ -0,0 +1,66 @@
|
||||
/* kolmogorov.c
|
||||
*
|
||||
* Copyright (C) 2010-2011 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include "kolmogorov.h"
|
||||
|
||||
double plfit_kolmogorov(double z) {
|
||||
const double fj[4] = { -2, -8, -18, -32 };
|
||||
const double w = 2.50662827;
|
||||
const double c1 = -1.2337005501361697; /* -pi^2 / 8 */
|
||||
const double c2 = -11.103304951225528; /* 9*c1 */
|
||||
const double c3 = -30.842513753404244; /* 25*c1 */
|
||||
|
||||
double u = fabs(z);
|
||||
double v;
|
||||
|
||||
if (u < 0.2)
|
||||
return 1;
|
||||
|
||||
if (u < 0.755) {
|
||||
v = 1.0 / (u*u);
|
||||
return 1 - w * (exp(c1*v) + exp(c2*v) + exp(c3*v)) / u;
|
||||
}
|
||||
|
||||
if (u < 6.8116) {
|
||||
double r[4] = { 0, 0, 0, 0 };
|
||||
long int maxj = (long int)(3.0 / u + 0.5);
|
||||
long int j;
|
||||
|
||||
if (maxj < 1)
|
||||
maxj = 1;
|
||||
|
||||
v = u*u;
|
||||
for (j = 0; j < maxj; j++) {
|
||||
r[j] = exp(fj[j] * v);
|
||||
}
|
||||
|
||||
return 2*(r[0] - r[1] + r[2] - r[3]);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
double plfit_ks_test_one_sample_p(double d, size_t n) {
|
||||
return plfit_kolmogorov(d * sqrt((double) n));
|
||||
}
|
||||
|
||||
double plfit_ks_test_two_sample_p(double d, size_t n1, size_t n2) {
|
||||
return plfit_kolmogorov(d * sqrt(n1*n2 / ((double)(n1+n2))));
|
||||
}
|
||||
@@ -0,0 +1,34 @@
|
||||
/* kolmogorov.h
|
||||
*
|
||||
* Copyright (C) 2010-2011 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef __KOLMOGOROV_H__
|
||||
#define __KOLMOGOROV_H__
|
||||
|
||||
#include <stdlib.h>
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
double plfit_kolmogorov(double z);
|
||||
double plfit_ks_test_one_sample_p(double d, size_t n);
|
||||
double plfit_ks_test_two_sample_p(double d, size_t n1, size_t n2);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif
|
||||
+1492
File diff suppressed because it is too large
Load Diff
+753
@@ -0,0 +1,753 @@
|
||||
/*
|
||||
* C library of Limited memory BFGS (L-BFGS).
|
||||
*
|
||||
* Copyright (c) 1990, Jorge Nocedal
|
||||
* Copyright (c) 2007-2010 Naoaki Okazaki
|
||||
* All rights reserved.
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
* of this software and associated documentation files (the "Software"), to deal
|
||||
* in the Software without restriction, including without limitation the rights
|
||||
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
* copies of the Software, and to permit persons to whom the Software is
|
||||
* furnished to do so, subject to the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be included in
|
||||
* all copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
* THE SOFTWARE.
|
||||
*/
|
||||
|
||||
/* $Id$ */
|
||||
|
||||
#ifndef __LBFGS_H__
|
||||
#define __LBFGS_H__
|
||||
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
/*
|
||||
* The default precision of floating point values is 64bit (double).
|
||||
*/
|
||||
#ifndef LBFGS_FLOAT
|
||||
#define LBFGS_FLOAT 64
|
||||
#endif/*LBFGS_FLOAT*/
|
||||
|
||||
/*
|
||||
* Activate optimization routines for IEEE754 floating point values.
|
||||
*/
|
||||
#ifndef LBFGS_IEEE_FLOAT
|
||||
#define LBFGS_IEEE_FLOAT 1
|
||||
#endif/*LBFGS_IEEE_FLOAT*/
|
||||
|
||||
#if LBFGS_FLOAT == 32
|
||||
typedef float lbfgsfloatval_t;
|
||||
|
||||
#elif LBFGS_FLOAT == 64
|
||||
typedef double lbfgsfloatval_t;
|
||||
|
||||
#else
|
||||
#error "libLBFGS supports single (float; LBFGS_FLOAT = 32) or double (double; LBFGS_FLOAT=64) precision only."
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
/**
|
||||
* \addtogroup liblbfgs_api libLBFGS API
|
||||
* @{
|
||||
*
|
||||
* The libLBFGS API.
|
||||
*/
|
||||
|
||||
/**
|
||||
* Return values of lbfgs().
|
||||
*
|
||||
* Roughly speaking, a negative value indicates an error.
|
||||
*/
|
||||
enum {
|
||||
/** L-BFGS reaches convergence. */
|
||||
LBFGS_SUCCESS = 0,
|
||||
LBFGS_CONVERGENCE = 0,
|
||||
LBFGS_STOP,
|
||||
/** The initial variables already minimize the objective function. */
|
||||
LBFGS_ALREADY_MINIMIZED,
|
||||
|
||||
/** Unknown error. */
|
||||
LBFGSERR_UNKNOWNERROR = -1024,
|
||||
/** Logic error. */
|
||||
LBFGSERR_LOGICERROR,
|
||||
/** Insufficient memory. */
|
||||
LBFGSERR_OUTOFMEMORY,
|
||||
/** The minimization process has been canceled. */
|
||||
LBFGSERR_CANCELED,
|
||||
/** Invalid number of variables specified. */
|
||||
LBFGSERR_INVALID_N,
|
||||
/** Invalid number of variables (for SSE) specified. */
|
||||
LBFGSERR_INVALID_N_SSE,
|
||||
/** The array x must be aligned to 16 (for SSE). */
|
||||
LBFGSERR_INVALID_X_SSE,
|
||||
/** Invalid parameter lbfgs_parameter_t::epsilon specified. */
|
||||
LBFGSERR_INVALID_EPSILON,
|
||||
/** Invalid parameter lbfgs_parameter_t::past specified. */
|
||||
LBFGSERR_INVALID_TESTPERIOD,
|
||||
/** Invalid parameter lbfgs_parameter_t::delta specified. */
|
||||
LBFGSERR_INVALID_DELTA,
|
||||
/** Invalid parameter lbfgs_parameter_t::linesearch specified. */
|
||||
LBFGSERR_INVALID_LINESEARCH,
|
||||
/** Invalid parameter lbfgs_parameter_t::max_step specified. */
|
||||
LBFGSERR_INVALID_MINSTEP,
|
||||
/** Invalid parameter lbfgs_parameter_t::max_step specified. */
|
||||
LBFGSERR_INVALID_MAXSTEP,
|
||||
/** Invalid parameter lbfgs_parameter_t::ftol specified. */
|
||||
LBFGSERR_INVALID_FTOL,
|
||||
/** Invalid parameter lbfgs_parameter_t::wolfe specified. */
|
||||
LBFGSERR_INVALID_WOLFE,
|
||||
/** Invalid parameter lbfgs_parameter_t::gtol specified. */
|
||||
LBFGSERR_INVALID_GTOL,
|
||||
/** Invalid parameter lbfgs_parameter_t::xtol specified. */
|
||||
LBFGSERR_INVALID_XTOL,
|
||||
/** Invalid parameter lbfgs_parameter_t::max_linesearch specified. */
|
||||
LBFGSERR_INVALID_MAXLINESEARCH,
|
||||
/** Invalid parameter lbfgs_parameter_t::orthantwise_c specified. */
|
||||
LBFGSERR_INVALID_ORTHANTWISE,
|
||||
/** Invalid parameter lbfgs_parameter_t::orthantwise_start specified. */
|
||||
LBFGSERR_INVALID_ORTHANTWISE_START,
|
||||
/** Invalid parameter lbfgs_parameter_t::orthantwise_end specified. */
|
||||
LBFGSERR_INVALID_ORTHANTWISE_END,
|
||||
/** The line-search step went out of the interval of uncertainty. */
|
||||
LBFGSERR_OUTOFINTERVAL,
|
||||
/** A logic error occurred; alternatively, the interval of uncertainty
|
||||
became too small. */
|
||||
LBFGSERR_INCORRECT_TMINMAX,
|
||||
/** A rounding error occurred; alternatively, no line-search step
|
||||
satisfies the sufficient decrease and curvature conditions. */
|
||||
LBFGSERR_ROUNDING_ERROR,
|
||||
/** The line-search step became smaller than lbfgs_parameter_t::min_step. */
|
||||
LBFGSERR_MINIMUMSTEP,
|
||||
/** The line-search step became larger than lbfgs_parameter_t::max_step. */
|
||||
LBFGSERR_MAXIMUMSTEP,
|
||||
/** The line-search routine reaches the maximum number of evaluations. */
|
||||
LBFGSERR_MAXIMUMLINESEARCH,
|
||||
/** The algorithm routine reaches the maximum number of iterations. */
|
||||
LBFGSERR_MAXIMUMITERATION,
|
||||
/** Relative width of the interval of uncertainty is at most
|
||||
lbfgs_parameter_t::xtol. */
|
||||
LBFGSERR_WIDTHTOOSMALL,
|
||||
/** A logic error (negative line-search step) occurred. */
|
||||
LBFGSERR_INVALIDPARAMETERS,
|
||||
/** The current search direction increases the objective function value. */
|
||||
LBFGSERR_INCREASEGRADIENT,
|
||||
};
|
||||
|
||||
/**
|
||||
* Line search algorithms.
|
||||
*/
|
||||
enum {
|
||||
/** The default algorithm (MoreThuente method). */
|
||||
LBFGS_LINESEARCH_DEFAULT = 0,
|
||||
/** MoreThuente method proposd by More and Thuente. */
|
||||
LBFGS_LINESEARCH_MORETHUENTE = 0,
|
||||
/**
|
||||
* Backtracking method with the Armijo condition.
|
||||
* The backtracking method finds the step length such that it satisfies
|
||||
* the sufficient decrease (Armijo) condition,
|
||||
* - f(x + a * d) <= f(x) + lbfgs_parameter_t::ftol * a * g(x)^T d,
|
||||
*
|
||||
* where x is the current point, d is the current search direction, and
|
||||
* a is the step length.
|
||||
*/
|
||||
LBFGS_LINESEARCH_BACKTRACKING_ARMIJO = 1,
|
||||
/** The backtracking method with the defualt (regular Wolfe) condition. */
|
||||
LBFGS_LINESEARCH_BACKTRACKING = 2,
|
||||
/**
|
||||
* Backtracking method with regular Wolfe condition.
|
||||
* The backtracking method finds the step length such that it satisfies
|
||||
* both the Armijo condition (LBFGS_LINESEARCH_BACKTRACKING_ARMIJO)
|
||||
* and the curvature condition,
|
||||
* - g(x + a * d)^T d >= lbfgs_parameter_t::wolfe * g(x)^T d,
|
||||
*
|
||||
* where x is the current point, d is the current search direction, and
|
||||
* a is the step length.
|
||||
*/
|
||||
LBFGS_LINESEARCH_BACKTRACKING_WOLFE = 2,
|
||||
/**
|
||||
* Backtracking method with strong Wolfe condition.
|
||||
* The backtracking method finds the step length such that it satisfies
|
||||
* both the Armijo condition (LBFGS_LINESEARCH_BACKTRACKING_ARMIJO)
|
||||
* and the following condition,
|
||||
* - |g(x + a * d)^T d| <= lbfgs_parameter_t::wolfe * |g(x)^T d|,
|
||||
*
|
||||
* where x is the current point, d is the current search direction, and
|
||||
* a is the step length.
|
||||
*/
|
||||
LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE = 3,
|
||||
};
|
||||
|
||||
/**
|
||||
* L-BFGS optimization parameters.
|
||||
* Call lbfgs_parameter_init() function to initialize parameters to the
|
||||
* default values.
|
||||
*/
|
||||
typedef struct {
|
||||
/**
|
||||
* The number of corrections to approximate the inverse hessian matrix.
|
||||
* The L-BFGS routine stores the computation results of previous \ref m
|
||||
* iterations to approximate the inverse hessian matrix of the current
|
||||
* iteration. This parameter controls the size of the limited memories
|
||||
* (corrections). The default value is \c 6. Values less than \c 3 are
|
||||
* not recommended. Large values will result in excessive computing time.
|
||||
*/
|
||||
int m;
|
||||
|
||||
/**
|
||||
* Epsilon for convergence test.
|
||||
* This parameter determines the accuracy with which the solution is to
|
||||
* be found. A minimization terminates when
|
||||
* ||g|| < \ref epsilon * max(1, ||x||),
|
||||
* where ||.|| denotes the Euclidean (L2) norm. The default value is
|
||||
* \c 1e-5.
|
||||
*/
|
||||
lbfgsfloatval_t epsilon;
|
||||
|
||||
/**
|
||||
* Distance for delta-based convergence test.
|
||||
* This parameter determines the distance, in iterations, to compute
|
||||
* the rate of decrease of the objective function. If the value of this
|
||||
* parameter is zero, the library does not perform the delta-based
|
||||
* convergence test. The default value is \c 0.
|
||||
*/
|
||||
int past;
|
||||
|
||||
/**
|
||||
* Delta for convergence test.
|
||||
* This parameter determines the minimum rate of decrease of the
|
||||
* objective function. The library stops iterations when the
|
||||
* following condition is met:
|
||||
* (f' - f) / f < \ref delta,
|
||||
* where f' is the objective value of \ref past iterations ago, and f is
|
||||
* the objective value of the current iteration.
|
||||
* The default value is \c 1e-5.
|
||||
*/
|
||||
lbfgsfloatval_t delta;
|
||||
|
||||
/**
|
||||
* The maximum number of iterations.
|
||||
* The lbfgs() function terminates an optimization process with
|
||||
* ::LBFGSERR_MAXIMUMITERATION status code when the iteration count
|
||||
* exceedes this parameter. Setting this parameter to zero continues an
|
||||
* optimization process until a convergence or error. The default value
|
||||
* is \c 0.
|
||||
*/
|
||||
int max_iterations;
|
||||
|
||||
/**
|
||||
* The line search algorithm.
|
||||
* This parameter specifies a line search algorithm to be used by the
|
||||
* L-BFGS routine.
|
||||
*/
|
||||
int linesearch;
|
||||
|
||||
/**
|
||||
* The maximum number of trials for the line search.
|
||||
* This parameter controls the number of function and gradients evaluations
|
||||
* per iteration for the line search routine. The default value is \c 40.
|
||||
*/
|
||||
int max_linesearch;
|
||||
|
||||
/**
|
||||
* The minimum step of the line search routine.
|
||||
* The default value is \c 1e-20. This value need not be modified unless
|
||||
* the exponents are too large for the machine being used, or unless the
|
||||
* problem is extremely badly scaled (in which case the exponents should
|
||||
* be increased).
|
||||
*/
|
||||
lbfgsfloatval_t min_step;
|
||||
|
||||
/**
|
||||
* The maximum step of the line search.
|
||||
* The default value is \c 1e+20. This value need not be modified unless
|
||||
* the exponents are too large for the machine being used, or unless the
|
||||
* problem is extremely badly scaled (in which case the exponents should
|
||||
* be increased).
|
||||
*/
|
||||
lbfgsfloatval_t max_step;
|
||||
|
||||
/**
|
||||
* A parameter to control the accuracy of the line search routine.
|
||||
* The default value is \c 1e-4. This parameter should be greater
|
||||
* than zero and smaller than \c 0.5.
|
||||
*/
|
||||
lbfgsfloatval_t ftol;
|
||||
|
||||
/**
|
||||
* A coefficient for the Wolfe condition.
|
||||
* This parameter is valid only when the backtracking line-search
|
||||
* algorithm is used with the Wolfe condition,
|
||||
* ::LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE or
|
||||
* ::LBFGS_LINESEARCH_BACKTRACKING_WOLFE .
|
||||
* The default value is \c 0.9. This parameter should be greater
|
||||
* the \ref ftol parameter and smaller than \c 1.0.
|
||||
*/
|
||||
lbfgsfloatval_t wolfe;
|
||||
|
||||
/**
|
||||
* A parameter to control the accuracy of the line search routine.
|
||||
* The default value is \c 0.9. If the function and gradient
|
||||
* evaluations are inexpensive with respect to the cost of the
|
||||
* iteration (which is sometimes the case when solving very large
|
||||
* problems) it may be advantageous to set this parameter to a small
|
||||
* value. A typical small value is \c 0.1. This parameter should be
|
||||
* greater than the \ref ftol parameter (\c 1e-4) and smaller than
|
||||
* \c 1.0.
|
||||
*/
|
||||
lbfgsfloatval_t gtol;
|
||||
|
||||
/**
|
||||
* The machine precision for floating-point values.
|
||||
* This parameter must be a positive value set by a client program to
|
||||
* estimate the machine precision. The line search routine will terminate
|
||||
* with the status code (::LBFGSERR_ROUNDING_ERROR) if the relative width
|
||||
* of the interval of uncertainty is less than this parameter.
|
||||
*/
|
||||
lbfgsfloatval_t xtol;
|
||||
|
||||
/**
|
||||
* Coeefficient for the L1 norm of variables.
|
||||
* This parameter should be set to zero for standard minimization
|
||||
* problems. Setting this parameter to a positive value activates
|
||||
* Orthant-Wise Limited-memory Quasi-Newton (OWL-QN) method, which
|
||||
* minimizes the objective function F(x) combined with the L1 norm |x|
|
||||
* of the variables, {F(x) + C |x|}. This parameter is the coeefficient
|
||||
* for the |x|, i.e., C. As the L1 norm |x| is not differentiable at
|
||||
* zero, the library modifies function and gradient evaluations from
|
||||
* a client program suitably; a client program thus have only to return
|
||||
* the function value F(x) and gradients G(x) as usual. The default value
|
||||
* is zero.
|
||||
*/
|
||||
lbfgsfloatval_t orthantwise_c;
|
||||
|
||||
/**
|
||||
* Start index for computing L1 norm of the variables.
|
||||
* This parameter is valid only for OWL-QN method
|
||||
* (i.e., \ref orthantwise_c != 0). This parameter b (0 <= b < N)
|
||||
* specifies the index number from which the library computes the
|
||||
* L1 norm of the variables x,
|
||||
* |x| := |x_{b}| + |x_{b+1}| + ... + |x_{N}| .
|
||||
* In other words, variables x_1, ..., x_{b-1} are not used for
|
||||
* computing the L1 norm. Setting b (0 < b < N), one can protect
|
||||
* variables, x_1, ..., x_{b-1} (e.g., a bias term of logistic
|
||||
* regression) from being regularized. The default value is zero.
|
||||
*/
|
||||
int orthantwise_start;
|
||||
|
||||
/**
|
||||
* End index for computing L1 norm of the variables.
|
||||
* This parameter is valid only for OWL-QN method
|
||||
* (i.e., \ref orthantwise_c != 0). This parameter e (0 < e <= N)
|
||||
* specifies the index number at which the library stops computing the
|
||||
* L1 norm of the variables x,
|
||||
*/
|
||||
int orthantwise_end;
|
||||
} lbfgs_parameter_t;
|
||||
|
||||
|
||||
/**
|
||||
* Callback interface to provide objective function and gradient evaluations.
|
||||
*
|
||||
* The lbfgs() function call this function to obtain the values of objective
|
||||
* function and its gradients when needed. A client program must implement
|
||||
* this function to evaluate the values of the objective function and its
|
||||
* gradients, given current values of variables.
|
||||
*
|
||||
* @param instance The user data sent for lbfgs() function by the client.
|
||||
* @param x The current values of variables.
|
||||
* @param g The gradient vector. The callback function must compute
|
||||
* the gradient values for the current variables.
|
||||
* @param n The number of variables.
|
||||
* @param step The current step of the line search routine.
|
||||
* @retval lbfgsfloatval_t The value of the objective function for the current
|
||||
* variables.
|
||||
*/
|
||||
typedef lbfgsfloatval_t (*lbfgs_evaluate_t)(
|
||||
void *instance,
|
||||
const lbfgsfloatval_t *x,
|
||||
lbfgsfloatval_t *g,
|
||||
const int n,
|
||||
const lbfgsfloatval_t step
|
||||
);
|
||||
|
||||
/**
|
||||
* Callback interface to receive the progress of the optimization process.
|
||||
*
|
||||
* The lbfgs() function call this function for each iteration. Implementing
|
||||
* this function, a client program can store or display the current progress
|
||||
* of the optimization process.
|
||||
*
|
||||
* @param instance The user data sent for lbfgs() function by the client.
|
||||
* @param x The current values of variables.
|
||||
* @param g The current gradient values of variables.
|
||||
* @param fx The current value of the objective function.
|
||||
* @param xnorm The Euclidean norm of the variables.
|
||||
* @param gnorm The Euclidean norm of the gradients.
|
||||
* @param step The line-search step used for this iteration.
|
||||
* @param n The number of variables.
|
||||
* @param k The iteration count.
|
||||
* @param ls The number of evaluations called for this iteration.
|
||||
* @retval int Zero to continue the optimization process. Returning a
|
||||
* non-zero value will cancel the optimization process.
|
||||
*/
|
||||
typedef int (*lbfgs_progress_t)(
|
||||
void *instance,
|
||||
const lbfgsfloatval_t *x,
|
||||
const lbfgsfloatval_t *g,
|
||||
const lbfgsfloatval_t fx,
|
||||
const lbfgsfloatval_t xnorm,
|
||||
const lbfgsfloatval_t gnorm,
|
||||
const lbfgsfloatval_t step,
|
||||
int n,
|
||||
int k,
|
||||
int ls
|
||||
);
|
||||
|
||||
/*
|
||||
A user must implement a function compatible with ::lbfgs_evaluate_t (evaluation
|
||||
callback) and pass the pointer to the callback function to lbfgs() arguments.
|
||||
Similarly, a user can implement a function compatible with ::lbfgs_progress_t
|
||||
(progress callback) to obtain the current progress (e.g., variables, function
|
||||
value, ||G||, etc) and to cancel the iteration process if necessary.
|
||||
Implementation of a progress callback is optional: a user can pass \c NULL if
|
||||
progress notification is not necessary.
|
||||
|
||||
In addition, a user must preserve two requirements:
|
||||
- The number of variables must be multiples of 16 (this is not 4).
|
||||
- The memory block of variable array ::x must be aligned to 16.
|
||||
|
||||
This algorithm terminates an optimization
|
||||
when:
|
||||
|
||||
||G|| < \epsilon \cdot \max(1, ||x||) .
|
||||
|
||||
In this formula, ||.|| denotes the Euclidean norm.
|
||||
*/
|
||||
|
||||
/**
|
||||
* Start a L-BFGS optimization.
|
||||
*
|
||||
* @param n The number of variables.
|
||||
* @param x The array of variables. A client program can set
|
||||
* default values for the optimization and receive the
|
||||
* optimization result through this array. This array
|
||||
* must be allocated by ::lbfgs_malloc function
|
||||
* for libLBFGS built with SSE/SSE2 optimization routine
|
||||
* enabled. The library built without SSE/SSE2
|
||||
* optimization does not have such a requirement.
|
||||
* @param ptr_fx The pointer to the variable that receives the final
|
||||
* value of the objective function for the variables.
|
||||
* This argument can be set to \c NULL if the final
|
||||
* value of the objective function is unnecessary.
|
||||
* @param proc_evaluate The callback function to provide function and
|
||||
* gradient evaluations given a current values of
|
||||
* variables. A client program must implement a
|
||||
* callback function compatible with \ref
|
||||
* lbfgs_evaluate_t and pass the pointer to the
|
||||
* callback function.
|
||||
* @param proc_progress The callback function to receive the progress
|
||||
* (the number of iterations, the current value of
|
||||
* the objective function) of the minimization
|
||||
* process. This argument can be set to \c NULL if
|
||||
* a progress report is unnecessary.
|
||||
* @param instance A user data for the client program. The callback
|
||||
* functions will receive the value of this argument.
|
||||
* @param param The pointer to a structure representing parameters for
|
||||
* L-BFGS optimization. A client program can set this
|
||||
* parameter to \c NULL to use the default parameters.
|
||||
* Call lbfgs_parameter_init() function to fill a
|
||||
* structure with the default values.
|
||||
* @retval int The status code. This function returns zero if the
|
||||
* minimization process terminates without an error. A
|
||||
* non-zero value indicates an error.
|
||||
*/
|
||||
int lbfgs(
|
||||
int n,
|
||||
lbfgsfloatval_t *x,
|
||||
lbfgsfloatval_t *ptr_fx,
|
||||
lbfgs_evaluate_t proc_evaluate,
|
||||
lbfgs_progress_t proc_progress,
|
||||
void *instance,
|
||||
lbfgs_parameter_t *param
|
||||
);
|
||||
|
||||
/**
|
||||
* Initialize L-BFGS parameters to the default values.
|
||||
*
|
||||
* Call this function to fill a parameter structure with the default values
|
||||
* and overwrite parameter values if necessary.
|
||||
*
|
||||
* @param param The pointer to the parameter structure.
|
||||
*/
|
||||
void lbfgs_parameter_init(lbfgs_parameter_t *param);
|
||||
|
||||
/**
|
||||
* Allocate an array for variables.
|
||||
*
|
||||
* This function allocates an array of variables for the convenience of
|
||||
* ::lbfgs function; the function has a requreiemt for a variable array
|
||||
* when libLBFGS is built with SSE/SSE2 optimization routines. A user does
|
||||
* not have to use this function for libLBFGS built without SSE/SSE2
|
||||
* optimization.
|
||||
*
|
||||
* @param n The number of variables.
|
||||
*/
|
||||
lbfgsfloatval_t* lbfgs_malloc(int n);
|
||||
|
||||
/**
|
||||
* Free an array of variables.
|
||||
*
|
||||
* @param x The array of variables allocated by ::lbfgs_malloc
|
||||
* function.
|
||||
*/
|
||||
void lbfgs_free(lbfgsfloatval_t *x);
|
||||
|
||||
/**
|
||||
* Get string description of an lbfgs() return code.
|
||||
*
|
||||
* @param err A value returned by lbfgs().
|
||||
*/
|
||||
const char* lbfgs_strerror(int err);
|
||||
|
||||
/** @} */
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
|
||||
|
||||
/**
|
||||
@mainpage libLBFGS: a library of Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS)
|
||||
|
||||
@section intro Introduction
|
||||
|
||||
This library is a C port of the implementation of Limited-memory
|
||||
Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method written by Jorge Nocedal.
|
||||
The original FORTRAN source code is available at:
|
||||
http://www.ece.northwestern.edu/~nocedal/lbfgs.html
|
||||
|
||||
The L-BFGS method solves the unconstrainted minimization problem,
|
||||
|
||||
<pre>
|
||||
minimize F(x), x = (x1, x2, ..., xN),
|
||||
</pre>
|
||||
|
||||
only if the objective function F(x) and its gradient G(x) are computable. The
|
||||
well-known Newton's method requires computation of the inverse of the hessian
|
||||
matrix of the objective function. However, the computational cost for the
|
||||
inverse hessian matrix is expensive especially when the objective function
|
||||
takes a large number of variables. The L-BFGS method iteratively finds a
|
||||
minimizer by approximating the inverse hessian matrix by information from last
|
||||
m iterations. This innovation saves the memory storage and computational time
|
||||
drastically for large-scaled problems.
|
||||
|
||||
Among the various ports of L-BFGS, this library provides several features:
|
||||
- <b>Optimization with L1-norm (Orthant-Wise Limited-memory Quasi-Newton
|
||||
(OWL-QN) method)</b>:
|
||||
In addition to standard minimization problems, the library can minimize
|
||||
a function F(x) combined with L1-norm |x| of the variables,
|
||||
{F(x) + C |x|}, where C is a constant scalar parameter. This feature is
|
||||
useful for estimating parameters of sparse log-linear models (e.g.,
|
||||
logistic regression and maximum entropy) with L1-regularization (or
|
||||
Laplacian prior).
|
||||
- <b>Clean C code</b>:
|
||||
Unlike C codes generated automatically by f2c (Fortran 77 into C converter),
|
||||
this port includes changes based on my interpretations, improvements,
|
||||
optimizations, and clean-ups so that the ported code would be well-suited
|
||||
for a C code. In addition to comments inherited from the original code,
|
||||
a number of comments were added through my interpretations.
|
||||
- <b>Callback interface</b>:
|
||||
The library receives function and gradient values via a callback interface.
|
||||
The library also notifies the progress of the optimization by invoking a
|
||||
callback function. In the original implementation, a user had to set
|
||||
function and gradient values every time the function returns for obtaining
|
||||
updated values.
|
||||
- <b>Thread safe</b>:
|
||||
The library is thread-safe, which is the secondary gain from the callback
|
||||
interface.
|
||||
- <b>Cross platform.</b> The source code can be compiled on Microsoft Visual
|
||||
Studio 2010, GNU C Compiler (gcc), etc.
|
||||
- <b>Configurable precision</b>: A user can choose single-precision (float)
|
||||
or double-precision (double) accuracy by changing ::LBFGS_FLOAT macro.
|
||||
- <b>SSE/SSE2 optimization</b>:
|
||||
This library includes SSE/SSE2 optimization (written in compiler intrinsics)
|
||||
for vector arithmetic operations on Intel/AMD processors. The library uses
|
||||
SSE for float values and SSE2 for double values. The SSE/SSE2 optimization
|
||||
routine is disabled by default.
|
||||
|
||||
This library is used by:
|
||||
- <a href="http://www.chokkan.org/software/crfsuite/">CRFsuite: A fast implementation of Conditional Random Fields (CRFs)</a>
|
||||
- <a href="http://www.chokkan.org/software/classias/">Classias: A collection of machine-learning algorithms for classification</a>
|
||||
- <a href="http://www.public.iastate.edu/~gdancik/mlegp/">mlegp: an R package for maximum likelihood estimates for Gaussian processes</a>
|
||||
- <a href="http://infmath.uibk.ac.at/~matthiasf/imaging2/">imaging2: the imaging2 class library</a>
|
||||
|
||||
@section download Download
|
||||
|
||||
- <a href="https://github.com/downloads/chokkan/liblbfgs/liblbfgs-1.10.tar.gz">Source code</a>
|
||||
- <a href="https://github.com/chokkan/liblbfgs">GitHub repository</a>
|
||||
|
||||
libLBFGS is distributed under the term of the
|
||||
<a href="http://opensource.org/licenses/mit-license.php">MIT license</a>.
|
||||
|
||||
@section modules Third-party modules
|
||||
- <a href="http://cran.r-project.org/web/packages/lbfgs/index.html">lbfgs: Limited-memory BFGS Optimization (a wrapper for R)</a> maintained by Antonio Coppola.
|
||||
- <a href="http://search.cpan.org/~laye/Algorithm-LBFGS-0.16/">Algorithm::LBFGS - Perl extension for L-BFGS</a> maintained by Lei Sun.
|
||||
- <a href="http://www.cs.kuleuven.be/~bernd/yap-lbfgs/">YAP-LBFGS (an interface to call libLBFGS from YAP Prolog)</a> maintained by Bernd Gutmann.
|
||||
|
||||
@section changelog History
|
||||
- Version 1.10 (2010-12-22):
|
||||
- Fixed compiling errors on Mac OS X; this patch was kindly submitted by
|
||||
Nic Schraudolph.
|
||||
- Reduced compiling warnings on Mac OS X; this patch was kindly submitted
|
||||
by Tamas Nepusz.
|
||||
- Replaced memalign() with posix_memalign().
|
||||
- Updated solution and project files for Microsoft Visual Studio 2010.
|
||||
- Version 1.9 (2010-01-29):
|
||||
- Fixed a mistake in checking the validity of the parameters "ftol" and
|
||||
"wolfe"; this was discovered by Kevin S. Van Horn.
|
||||
- Version 1.8 (2009-07-13):
|
||||
- Accepted the patch submitted by Takashi Imamichi;
|
||||
the backtracking method now has three criteria for choosing the step
|
||||
length:
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING_ARMIJO: sufficient decrease (Armijo)
|
||||
condition only
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING_WOLFE: regular Wolfe condition
|
||||
(sufficient decrease condition + curvature condition)
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE: strong Wolfe condition
|
||||
- Updated the documentation to explain the above three criteria.
|
||||
- Version 1.7 (2009-02-28):
|
||||
- Improved OWL-QN routines for stability.
|
||||
- Removed the support of OWL-QN method in MoreThuente algorithm because
|
||||
it accidentally fails in early stages of iterations for some objectives.
|
||||
Because of this change, <b>the OW-LQN method must be used with the
|
||||
backtracking algorithm (::LBFGS_LINESEARCH_BACKTRACKING)</b>, or the
|
||||
library returns ::LBFGSERR_INVALID_LINESEARCH.
|
||||
- Renamed line search algorithms as follows:
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING: regular Wolfe condition.
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING_LOOSE: regular Wolfe condition.
|
||||
- ::LBFGS_LINESEARCH_BACKTRACKING_STRONG: strong Wolfe condition.
|
||||
- Source code clean-up.
|
||||
- Version 1.6 (2008-11-02):
|
||||
- Improved line-search algorithm with strong Wolfe condition, which was
|
||||
contributed by Takashi Imamichi. This routine is now default for
|
||||
::LBFGS_LINESEARCH_BACKTRACKING. The previous line search algorithm
|
||||
with regular Wolfe condition is still available as
|
||||
::LBFGS_LINESEARCH_BACKTRACKING_LOOSE.
|
||||
- Configurable stop index for L1-norm computation. A member variable
|
||||
::lbfgs_parameter_t::orthantwise_end was added to specify the index
|
||||
number at which the library stops computing the L1 norm of the
|
||||
variables. This is useful to prevent some variables from being
|
||||
regularized by the OW-LQN method.
|
||||
- A sample program written in C++ (sample/sample.cpp).
|
||||
- Version 1.5 (2008-07-10):
|
||||
- Configurable starting index for L1-norm computation. A member variable
|
||||
::lbfgs_parameter_t::orthantwise_start was added to specify the index
|
||||
number from which the library computes the L1 norm of the variables.
|
||||
This is useful to prevent some variables from being regularized by the
|
||||
OWL-QN method.
|
||||
- Fixed a zero-division error when the initial variables have already
|
||||
been a minimizer (reported by Takashi Imamichi). In this case, the
|
||||
library returns ::LBFGS_ALREADY_MINIMIZED status code.
|
||||
- Defined ::LBFGS_SUCCESS status code as zero; removed unused constants,
|
||||
LBFGSFALSE and LBFGSTRUE.
|
||||
- Fixed a compile error in an implicit down-cast.
|
||||
- Version 1.4 (2008-04-25):
|
||||
- Configurable line search algorithms. A member variable
|
||||
::lbfgs_parameter_t::linesearch was added to choose either MoreThuente
|
||||
method (::LBFGS_LINESEARCH_MORETHUENTE) or backtracking algorithm
|
||||
(::LBFGS_LINESEARCH_BACKTRACKING).
|
||||
- Fixed a bug: the previous version did not compute psuedo-gradients
|
||||
properly in the line search routines for OWL-QN. This bug might quit
|
||||
an iteration process too early when the OWL-QN routine was activated
|
||||
(0 < ::lbfgs_parameter_t::orthantwise_c).
|
||||
- Configure script for POSIX environments.
|
||||
- SSE/SSE2 optimizations with GCC.
|
||||
- New functions ::lbfgs_malloc and ::lbfgs_free to use SSE/SSE2 routines
|
||||
transparently. It is uncessary to use these functions for libLBFGS built
|
||||
without SSE/SSE2 routines; you can still use any memory allocators if
|
||||
SSE/SSE2 routines are disabled in libLBFGS.
|
||||
- Version 1.3 (2007-12-16):
|
||||
- An API change. An argument was added to lbfgs() function to receive the
|
||||
final value of the objective function. This argument can be set to
|
||||
\c NULL if the final value is unnecessary.
|
||||
- Fixed a null-pointer bug in the sample code (reported by Takashi Imamichi).
|
||||
- Added build scripts for Microsoft Visual Studio 2005 and GCC.
|
||||
- Added README file.
|
||||
- Version 1.2 (2007-12-13):
|
||||
- Fixed a serious bug in orthant-wise L-BFGS.
|
||||
An important variable was used without initialization.
|
||||
- Version 1.1 (2007-12-01):
|
||||
- Implemented orthant-wise L-BFGS.
|
||||
- Implemented lbfgs_parameter_init() function.
|
||||
- Fixed several bugs.
|
||||
- API documentation.
|
||||
- Version 1.0 (2007-09-20):
|
||||
- Initial release.
|
||||
|
||||
@section api Documentation
|
||||
|
||||
- @ref liblbfgs_api "libLBFGS API"
|
||||
|
||||
@section sample Sample code
|
||||
|
||||
@include sample.c
|
||||
|
||||
@section ack Acknowledgements
|
||||
|
||||
The L-BFGS algorithm is described in:
|
||||
- Jorge Nocedal.
|
||||
Updating Quasi-Newton Matrices with Limited Storage.
|
||||
<i>Mathematics of Computation</i>, Vol. 35, No. 151, pp. 773--782, 1980.
|
||||
- Dong C. Liu and Jorge Nocedal.
|
||||
On the limited memory BFGS method for large scale optimization.
|
||||
<i>Mathematical Programming</i> B, Vol. 45, No. 3, pp. 503-528, 1989.
|
||||
|
||||
The line search algorithms used in this implementation are described in:
|
||||
- John E. Dennis and Robert B. Schnabel.
|
||||
<i>Numerical Methods for Unconstrained Optimization and Nonlinear
|
||||
Equations</i>, Englewood Cliffs, 1983.
|
||||
- Jorge J. More and David J. Thuente.
|
||||
Line search algorithm with guaranteed sufficient decrease.
|
||||
<i>ACM Transactions on Mathematical Software (TOMS)</i>, Vol. 20, No. 3,
|
||||
pp. 286-307, 1994.
|
||||
|
||||
This library also implements Orthant-Wise Limited-memory Quasi-Newton (OWL-QN)
|
||||
method presented in:
|
||||
- Galen Andrew and Jianfeng Gao.
|
||||
Scalable training of L1-regularized log-linear models.
|
||||
In <i>Proceedings of the 24th International Conference on Machine
|
||||
Learning (ICML 2007)</i>, pp. 33-40, 2007.
|
||||
|
||||
Special thanks go to:
|
||||
- Yoshimasa Tsuruoka and Daisuke Okanohara for technical information about
|
||||
OWL-QN
|
||||
- Takashi Imamichi for the useful enhancements of the backtracking method
|
||||
- Kevin S. Van Horn, Nic Schraudolph, and Tamas Nepusz for bug fixes
|
||||
|
||||
Finally I would like to thank the original author, Jorge Nocedal, who has been
|
||||
distributing the effieicnt and explanatory implementation in an open source
|
||||
licence.
|
||||
|
||||
@section reference Reference
|
||||
|
||||
- <a href="http://www.ece.northwestern.edu/~nocedal/lbfgs.html">L-BFGS</a> by Jorge Nocedal.
|
||||
- <a href="http://research.microsoft.com/en-us/downloads/b1eb1016-1738-4bd5-83a9-370c9d498a03/default.aspx">Orthant-Wise Limited-memory Quasi-Newton Optimizer for L1-regularized Objectives</a> by Galen Andrew.
|
||||
- <a href="http://chasen.org/~taku/software/misc/lbfgs/">C port (via f2c)</a> by Taku Kudo.
|
||||
- <a href="http://www.alglib.net/optimization/lbfgs.php">C#/C++/Delphi/VisualBasic6 port</a> in ALGLIB.
|
||||
- <a href="http://cctbx.sourceforge.net/">Computational Crystallography Toolbox</a> includes
|
||||
<a href="http://cctbx.sourceforge.net/current_cvs/c_plus_plus/namespacescitbx_1_1lbfgs.html">scitbx::lbfgs</a>.
|
||||
*/
|
||||
|
||||
#endif/*__LBFGS_H__*/
|
||||
+91
@@ -0,0 +1,91 @@
|
||||
/* mt.c
|
||||
*
|
||||
* Mersenne Twister random number generator, based on the implementation of
|
||||
* Michael Brundage (which has been placed in the public domain).
|
||||
*
|
||||
* Author: Tamas Nepusz (original by Michael Brundage)
|
||||
*
|
||||
* See the following URL for the original implementation:
|
||||
* http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html
|
||||
*
|
||||
* This file has been placed in the public domain.
|
||||
*/
|
||||
|
||||
#include "igraph_random.h"
|
||||
|
||||
#include "plfit_mt.h"
|
||||
|
||||
static uint16_t get_random_uint16(void) {
|
||||
return RNG_INTEGER(0, 0xffff);
|
||||
}
|
||||
|
||||
void plfit_mt_init(plfit_mt_rng_t* rng) {
|
||||
plfit_mt_init_from_rng(rng, 0);
|
||||
}
|
||||
|
||||
void plfit_mt_init_from_rng(plfit_mt_rng_t* rng, plfit_mt_rng_t* seeder) {
|
||||
int i;
|
||||
|
||||
if (seeder == 0) {
|
||||
for (i = 0; i < PLFIT_MT_LEN; i++) {
|
||||
/* RAND_MAX is guaranteed to be at least 32767, so we can use two
|
||||
* calls to rand() to produce a random 32-bit number */
|
||||
rng->mt_buffer[i] = (((uint32_t) get_random_uint16()) << 16) + get_random_uint16();
|
||||
}
|
||||
} else {
|
||||
for (i = 0; i < PLFIT_MT_LEN; i++) {
|
||||
rng->mt_buffer[i] = plfit_mt_random(seeder);
|
||||
}
|
||||
}
|
||||
|
||||
rng->mt_index = 0;
|
||||
}
|
||||
|
||||
#define MT_IA 397
|
||||
#define MT_IB (PLFIT_MT_LEN - MT_IA)
|
||||
#define UPPER_MASK 0x80000000
|
||||
#define LOWER_MASK 0x7FFFFFFF
|
||||
#define MATRIX_A 0x9908B0DF
|
||||
#define TWIST(b,i,j) ((b)[i] & UPPER_MASK) | ((b)[j] & LOWER_MASK)
|
||||
#define MAGIC(s) (((s)&1)*MATRIX_A)
|
||||
|
||||
uint32_t plfit_mt_random(plfit_mt_rng_t* rng) {
|
||||
uint32_t * b = rng->mt_buffer;
|
||||
int idx = rng->mt_index;
|
||||
uint32_t s;
|
||||
int i;
|
||||
|
||||
if (idx == PLFIT_MT_LEN * sizeof(uint32_t)) {
|
||||
idx = 0;
|
||||
i = 0;
|
||||
for (; i < MT_IB; i++) {
|
||||
s = TWIST(b, i, i+1);
|
||||
b[i] = b[i + MT_IA] ^ (s >> 1) ^ MAGIC(s);
|
||||
}
|
||||
for (; i < PLFIT_MT_LEN-1; i++) {
|
||||
s = TWIST(b, i, i+1);
|
||||
b[i] = b[i - MT_IB] ^ (s >> 1) ^ MAGIC(s);
|
||||
}
|
||||
|
||||
s = TWIST(b, PLFIT_MT_LEN-1, 0);
|
||||
b[PLFIT_MT_LEN-1] = b[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
|
||||
}
|
||||
|
||||
rng->mt_index = idx + sizeof(uint32_t);
|
||||
return *(uint32_t *)((unsigned char *)b + idx);
|
||||
/*
|
||||
Matsumoto and Nishimura additionally confound the bits returned to the caller
|
||||
but this doesn't increase the randomness, and slows down the generator by
|
||||
as much as 25%. So I omit these operations here.
|
||||
|
||||
r ^= (r >> 11);
|
||||
r ^= (r << 7) & 0x9D2C5680;
|
||||
r ^= (r << 15) & 0xEFC60000;
|
||||
r ^= (r >> 18);
|
||||
*/
|
||||
}
|
||||
|
||||
|
||||
double plfit_mt_uniform_01(plfit_mt_rng_t* rng) {
|
||||
return ((double)plfit_mt_random(rng)) / PLFIT_MT_RAND_MAX;
|
||||
}
|
||||
+52
@@ -0,0 +1,52 @@
|
||||
/* options.c
|
||||
*
|
||||
* Copyright (C) 2012 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#include "plfit_error.h"
|
||||
#include "plfit.h"
|
||||
|
||||
const plfit_continuous_options_t plfit_continuous_default_options = {
|
||||
/* .finite_size_correction = */ 0,
|
||||
/* .xmin_method = */ PLFIT_DEFAULT_CONTINUOUS_METHOD,
|
||||
/* .p_value_method = */ PLFIT_DEFAULT_P_VALUE_METHOD,
|
||||
/* .p_value_precision = */ 0.01,
|
||||
/* .rng = */ 0
|
||||
};
|
||||
|
||||
const plfit_discrete_options_t plfit_discrete_default_options = {
|
||||
/* .finite_size_correction = */ 0,
|
||||
/* .alpha_method = */ PLFIT_DEFAULT_DISCRETE_METHOD,
|
||||
/* .alpha = */ {
|
||||
/* .min = */ 1.01,
|
||||
/* .max = */ 5,
|
||||
/* .step = */ 0.01
|
||||
},
|
||||
/* .p_value_method = */ PLFIT_DEFAULT_P_VALUE_METHOD,
|
||||
/* .p_value_precision = */ 0.01,
|
||||
/* .rng = */ 0
|
||||
};
|
||||
|
||||
int plfit_continuous_options_init(plfit_continuous_options_t* options) {
|
||||
*options = plfit_continuous_default_options;
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
|
||||
int plfit_discrete_options_init(plfit_discrete_options_t* options) {
|
||||
*options = plfit_discrete_default_options;
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
+1376
File diff suppressed because it is too large
Load Diff
+132
@@ -0,0 +1,132 @@
|
||||
/* plfit.h
|
||||
*
|
||||
* Copyright (C) 2010-2011 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_H
|
||||
#define PLFIT_H
|
||||
|
||||
#include <stdlib.h>
|
||||
#include "plfit_decls.h"
|
||||
#include "plfit_error.h"
|
||||
#include "plfit_mt.h"
|
||||
#include "plfit_sampling.h"
|
||||
#include "plfit_version.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
typedef unsigned short int plfit_bool_t;
|
||||
|
||||
typedef enum {
|
||||
PLFIT_LINEAR_ONLY,
|
||||
PLFIT_STRATIFIED_SAMPLING,
|
||||
PLFIT_GSS_OR_LINEAR,
|
||||
PLFIT_DEFAULT_CONTINUOUS_METHOD = PLFIT_STRATIFIED_SAMPLING
|
||||
} plfit_continuous_method_t;
|
||||
|
||||
typedef enum {
|
||||
PLFIT_LBFGS,
|
||||
PLFIT_LINEAR_SCAN,
|
||||
PLFIT_PRETEND_CONTINUOUS,
|
||||
PLFIT_DEFAULT_DISCRETE_METHOD = PLFIT_LBFGS
|
||||
} plfit_discrete_method_t;
|
||||
|
||||
typedef enum {
|
||||
PLFIT_P_VALUE_SKIP,
|
||||
PLFIT_P_VALUE_APPROXIMATE,
|
||||
PLFIT_P_VALUE_EXACT,
|
||||
PLFIT_DEFAULT_P_VALUE_METHOD = PLFIT_P_VALUE_EXACT
|
||||
} plfit_p_value_method_t;
|
||||
|
||||
typedef struct _plfit_result_t {
|
||||
double alpha; /* fitted power-law exponent */
|
||||
double xmin; /* cutoff where the power-law behaviour kicks in */
|
||||
double L; /* log-likelihood of the sample */
|
||||
double D; /* test statistic for the KS test */
|
||||
double p; /* p-value of the KS test */
|
||||
} plfit_result_t;
|
||||
|
||||
/********** structure that holds the options of plfit **********/
|
||||
|
||||
typedef struct _plfit_continuous_options_t {
|
||||
plfit_bool_t finite_size_correction;
|
||||
plfit_continuous_method_t xmin_method;
|
||||
plfit_p_value_method_t p_value_method;
|
||||
double p_value_precision;
|
||||
plfit_mt_rng_t* rng;
|
||||
} plfit_continuous_options_t;
|
||||
|
||||
typedef struct _plfit_discrete_options_t {
|
||||
plfit_bool_t finite_size_correction;
|
||||
plfit_discrete_method_t alpha_method;
|
||||
struct {
|
||||
double min;
|
||||
double max;
|
||||
double step;
|
||||
} alpha;
|
||||
plfit_p_value_method_t p_value_method;
|
||||
double p_value_precision;
|
||||
plfit_mt_rng_t* rng;
|
||||
} plfit_discrete_options_t;
|
||||
|
||||
PLFIT_EXPORT int plfit_continuous_options_init(plfit_continuous_options_t* options);
|
||||
PLFIT_EXPORT int plfit_discrete_options_init(plfit_discrete_options_t* options);
|
||||
|
||||
PLFIT_EXPORT extern const plfit_continuous_options_t plfit_continuous_default_options;
|
||||
PLFIT_EXPORT extern const plfit_discrete_options_t plfit_discrete_default_options;
|
||||
|
||||
/********** continuous power law distribution fitting **********/
|
||||
|
||||
PLFIT_EXPORT int plfit_log_likelihood_continuous(const double* xs, size_t n, double alpha,
|
||||
double xmin, double* l);
|
||||
PLFIT_EXPORT int plfit_estimate_alpha_continuous(const double* xs, size_t n, double xmin,
|
||||
const plfit_continuous_options_t* options, plfit_result_t* result);
|
||||
PLFIT_EXPORT int plfit_continuous(const double* xs, size_t n,
|
||||
const plfit_continuous_options_t* options, plfit_result_t* result);
|
||||
|
||||
/*********** discrete power law distribution fitting ***********/
|
||||
|
||||
PLFIT_EXPORT int plfit_estimate_alpha_discrete(const double* xs, size_t n, double xmin,
|
||||
const plfit_discrete_options_t* options, plfit_result_t *result);
|
||||
PLFIT_EXPORT int plfit_log_likelihood_discrete(const double* xs, size_t n, double alpha, double xmin, double* l);
|
||||
PLFIT_EXPORT int plfit_discrete(const double* xs, size_t n, const plfit_discrete_options_t* options,
|
||||
plfit_result_t* result);
|
||||
|
||||
/***** resampling routines to generate synthetic replicates ****/
|
||||
|
||||
PLFIT_EXPORT int plfit_resample_continuous(const double* xs, size_t n, double alpha, double xmin,
|
||||
size_t num_samples, plfit_mt_rng_t* rng, double* result);
|
||||
PLFIT_EXPORT int plfit_resample_discrete(const double* xs, size_t n, double alpha, double xmin,
|
||||
size_t num_samples, plfit_mt_rng_t* rng, double* result);
|
||||
|
||||
/******** calculating the p-value of a fitted model only *******/
|
||||
|
||||
PLFIT_EXPORT int plfit_calculate_p_value_continuous(const double* xs, size_t n,
|
||||
const plfit_continuous_options_t* options, plfit_bool_t xmin_fixed,
|
||||
plfit_result_t *result);
|
||||
PLFIT_EXPORT int plfit_calculate_p_value_discrete(const double* xs, size_t n,
|
||||
const plfit_discrete_options_t* options, plfit_bool_t xmin_fixed,
|
||||
plfit_result_t *result);
|
||||
|
||||
/************* calculating descriptive statistics **************/
|
||||
|
||||
PLFIT_EXPORT int plfit_moments(const double* data, size_t n, double* mean, double* variance,
|
||||
double* skewness, double* kurtosis);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif /* PLFIT_H */
|
||||
@@ -0,0 +1,35 @@
|
||||
/* plfit_decls.h
|
||||
*
|
||||
* Copyright (C) 2024 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_DECLS_H
|
||||
#define PLFIT_DECLS_H
|
||||
|
||||
#undef PLFIT_BEGIN_C_DECLS
|
||||
#undef PLFIT_END_C_DECLS
|
||||
#ifdef __cplusplus
|
||||
#define PLFIT_BEGIN_C_DECLS extern "C" {
|
||||
#define PLFIT_END_C_DECLS }
|
||||
#else
|
||||
#define PLFIT_BEGIN_C_DECLS /* empty */
|
||||
#define PLFIT_END_C_DECLS /* empty */
|
||||
#endif
|
||||
|
||||
#define PLFIT_EXPORT /* empty */
|
||||
|
||||
#endif /* PLFIT_DECLS_H */
|
||||
@@ -0,0 +1,66 @@
|
||||
/* error.c
|
||||
*
|
||||
* Copyright (C) 2010-2011 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include "plfit_error.h"
|
||||
|
||||
static const char *plfit_i_error_strings[] = {
|
||||
"No error",
|
||||
"Failed",
|
||||
"Invalid value",
|
||||
"Underflow",
|
||||
"Overflow",
|
||||
"Not enough memory",
|
||||
"Maximum number of iterations exceeded"
|
||||
};
|
||||
|
||||
#ifndef USING_R
|
||||
static plfit_error_handler_t* plfit_error_handler = plfit_error_handler_printignore;
|
||||
#else
|
||||
/* This is overwritten, anyway */
|
||||
static plfit_error_handler_t* plfit_error_handler = plfit_error_handler_ignore;
|
||||
#endif
|
||||
|
||||
const char* plfit_strerror(const int plfit_errno) {
|
||||
return plfit_i_error_strings[plfit_errno];
|
||||
}
|
||||
|
||||
plfit_error_handler_t* plfit_set_error_handler(plfit_error_handler_t* new_handler) {
|
||||
plfit_error_handler_t* old_handler = plfit_error_handler;
|
||||
plfit_error_handler = new_handler;
|
||||
return old_handler;
|
||||
}
|
||||
|
||||
void plfit_error(const char *reason, const char *file, int line,
|
||||
int plfit_errno) {
|
||||
plfit_error_handler(reason, file, line, plfit_errno);
|
||||
}
|
||||
|
||||
#ifndef USING_R
|
||||
void plfit_error_handler_printignore(const char *reason, const char *file, int line,
|
||||
int plfit_errno) {
|
||||
fprintf(stderr, "Error at %s:%i : %s, %s\n", file, line, reason,
|
||||
plfit_strerror(plfit_errno));
|
||||
}
|
||||
#endif
|
||||
|
||||
void plfit_error_handler_ignore(const char* reason, const char* file, int line,
|
||||
int plfit_errno) {
|
||||
}
|
||||
@@ -0,0 +1,72 @@
|
||||
/* plfit_error.h
|
||||
*
|
||||
* Copyright (C) 2010-2011 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_ERROR_H
|
||||
#define PLFIT_ERROR_H
|
||||
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
enum {
|
||||
PLFIT_SUCCESS = 0,
|
||||
PLFIT_FAILURE = 1,
|
||||
PLFIT_EINVAL = 2,
|
||||
PLFIT_UNDRFLOW = 3,
|
||||
PLFIT_OVERFLOW = 4,
|
||||
PLFIT_ENOMEM = 5,
|
||||
PLFIT_EMAXITER = 6
|
||||
};
|
||||
|
||||
#if (defined(__GNUC__) && GCC_VERSION_MAJOR >= 3)
|
||||
# define PLFIT_UNLIKELY(a) __builtin_expect((a), 0)
|
||||
# define PLFIT_LIKELY(a) __builtin_expect((a), 1)
|
||||
#else
|
||||
# define PLFIT_UNLIKELY(a) a
|
||||
# define PLFIT_LIKELY(a) a
|
||||
#endif
|
||||
|
||||
#define PLFIT_CHECK(a) \
|
||||
do {\
|
||||
int plfit_i_ret=(a); \
|
||||
if (PLFIT_UNLIKELY(plfit_i_ret != PLFIT_SUCCESS)) {\
|
||||
return plfit_i_ret; \
|
||||
} \
|
||||
} while (0)
|
||||
|
||||
#define PLFIT_ERROR(reason,plfit_errno) \
|
||||
do {\
|
||||
plfit_error (reason, __FILE__, __LINE__, plfit_errno) ; \
|
||||
return plfit_errno ; \
|
||||
} while (0)
|
||||
|
||||
typedef void plfit_error_handler_t(const char*, const char*, int, int);
|
||||
|
||||
PLFIT_EXPORT extern plfit_error_handler_t plfit_error_handler_abort;
|
||||
PLFIT_EXPORT extern plfit_error_handler_t plfit_error_handler_ignore;
|
||||
PLFIT_EXPORT extern plfit_error_handler_t plfit_error_handler_printignore;
|
||||
|
||||
PLFIT_EXPORT plfit_error_handler_t* plfit_set_error_handler(plfit_error_handler_t* new_handler);
|
||||
|
||||
PLFIT_EXPORT void plfit_error(const char *reason, const char *file, int line, int plfit_errno);
|
||||
PLFIT_EXPORT const char* plfit_strerror(const int plfit_errno);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif /* PLFIT_ERROR_H */
|
||||
+87
@@ -0,0 +1,87 @@
|
||||
/* plfit_mt.h
|
||||
*
|
||||
* Mersenne Twister random number generator, based on the implementation of
|
||||
* Michael Brundage (which has been placed in the public domain).
|
||||
*
|
||||
* Author: Tamas Nepusz (original by Michael Brundage)
|
||||
*
|
||||
* See the following URL for the original implementation:
|
||||
* http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html
|
||||
*
|
||||
* This file has been placed in the public domain.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_MT_H
|
||||
#define PLFIT_MT_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include "plfit_decls.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
#define PLFIT_MT_LEN 624
|
||||
|
||||
/**
|
||||
* \def PLFIT_MT_RAND_MAX
|
||||
*
|
||||
* The maximum random number that \c plfit_mt_random() can generate.
|
||||
*/
|
||||
#define PLFIT_MT_RAND_MAX 0xFFFFFFFF
|
||||
|
||||
/**
|
||||
* Struct that stores the internal state of a Mersenne Twister random number
|
||||
* generator.
|
||||
*/
|
||||
typedef struct {
|
||||
int mt_index;
|
||||
uint32_t mt_buffer[PLFIT_MT_LEN];
|
||||
} plfit_mt_rng_t;
|
||||
|
||||
/**
|
||||
* \brief Initializes a Mersenne Twister random number generator.
|
||||
*
|
||||
* The random number generator is seeded with random 32-bit numbers obtained
|
||||
* from the \em built-in random number generator using consecutive calls to
|
||||
* \c rand().
|
||||
*
|
||||
* \param rng the random number generator to initialize
|
||||
*/
|
||||
PLFIT_EXPORT void plfit_mt_init(plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* \brief Initializes a Mersenne Twister random number generator, seeding it
|
||||
* from another one.
|
||||
*
|
||||
* The random number generator is seeded with random 32-bit numbers obtained
|
||||
* from another, initialized Mersenne Twister random number generator.
|
||||
*
|
||||
* \param rng the random number generator to initialize
|
||||
* \param seeder the random number generator that will seed the one being
|
||||
* initialized. When null, the random number generator will
|
||||
* be initialized from the built-in RNG as if \ref plfit_mt_init()
|
||||
* was called.
|
||||
*/
|
||||
PLFIT_EXPORT void plfit_mt_init_from_rng(plfit_mt_rng_t* rng, plfit_mt_rng_t* seeder);
|
||||
|
||||
/**
|
||||
* \brief Returns the next 32-bit random number from the given Mersenne Twister
|
||||
* random number generator.
|
||||
*
|
||||
* \param rng the random number generator to use
|
||||
* \return the next 32-bit random number from the generator
|
||||
*/
|
||||
PLFIT_EXPORT uint32_t plfit_mt_random(plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* \brief Returns a uniformly distributed double from the interval [0;1)
|
||||
* based on the next value of the given Mersenne Twister random number
|
||||
* generator.
|
||||
*
|
||||
* \param rng the random number generator to use
|
||||
* \return a uniformly distributed random number from the interval [0;1)
|
||||
*/
|
||||
PLFIT_EXPORT double plfit_mt_uniform_01(plfit_mt_rng_t* rng);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif /* PLFIT_MT_H */
|
||||
@@ -0,0 +1,168 @@
|
||||
/* plfit_sampling.h
|
||||
*
|
||||
* Copyright (C) 2012 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_SAMPLING_H
|
||||
#define PLFIT_SAMPLING_H
|
||||
|
||||
#include <stdlib.h>
|
||||
#include "plfit_decls.h"
|
||||
#include "plfit_mt.h"
|
||||
|
||||
PLFIT_BEGIN_C_DECLS
|
||||
|
||||
/**
|
||||
* Draws a sample from a binomial distribution with the given count and
|
||||
* probability values.
|
||||
*
|
||||
* This function is borrowed from R; see the corresponding license in
|
||||
* \c rbinom.c. The return value is always an integer.
|
||||
*
|
||||
* The function is \em not thread-safe.
|
||||
*
|
||||
* \param n the number of trials
|
||||
* \param p the success probability of each trial
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \return the value drawn from the given binomial distribution.
|
||||
*/
|
||||
PLFIT_EXPORT double plfit_rbinom(double n, double p, plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* Draws a sample from a Pareto distribution with the given minimum value and
|
||||
* power-law exponent.
|
||||
*
|
||||
* \param xmin the minimum value of the distribution. Must be positive.
|
||||
* \param alpha the exponent. Must be positive
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
*
|
||||
* \return the sample or NaN if one of the parameters is invalid
|
||||
*/
|
||||
PLFIT_EXPORT extern double plfit_rpareto(double xmin, double alpha, plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* Draws a given number of samples from a Pareto distribution with the given
|
||||
* minimum value and power-law exponent.
|
||||
*
|
||||
* \param xmin the minimum value of the distribution. Must be positive.
|
||||
* \param alpha the exponent. Must be positive
|
||||
* \param n the number of samples to draw
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \param result the array where the result should be written. It must
|
||||
* have enough space to store n items
|
||||
*
|
||||
* \return \c PLFIT_EINVAL if one of the parameters is invalid, zero otherwise
|
||||
*/
|
||||
PLFIT_EXPORT int plfit_rpareto_array(double xmin, double alpha, size_t n, plfit_mt_rng_t* rng,
|
||||
double* result);
|
||||
|
||||
/**
|
||||
* Draws a sample from a zeta distribution with the given minimum value and
|
||||
* power-law exponent.
|
||||
*
|
||||
* \param xmin the minimum value of the distribution. Must be positive.
|
||||
* \param alpha the exponent. Must be positive
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
*
|
||||
* \return the sample or NaN if one of the parameters is invalid
|
||||
*/
|
||||
PLFIT_EXPORT extern double plfit_rzeta(long int xmin, double alpha, plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* Draws a given number of samples from a zeta distribution with the given
|
||||
* minimum value and power-law exponent.
|
||||
*
|
||||
* \param xmin the minimum value of the distribution. Must be positive.
|
||||
* \param alpha the exponent. Must be positive
|
||||
* \param n the number of samples to draw
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \param result the array where the result should be written. It must
|
||||
* have enough space to store n items
|
||||
*
|
||||
* \return \c PLFIT_EINVAL if one of the parameters is invalid, zero otherwise
|
||||
*/
|
||||
PLFIT_EXPORT int plfit_rzeta_array(long int xmin, double alpha, size_t n, plfit_mt_rng_t* rng,
|
||||
double* result);
|
||||
|
||||
/**
|
||||
* Draws a sample from a uniform distribution with the given lower and
|
||||
* upper bounds.
|
||||
*
|
||||
* The lower bound is inclusive, the uppoer bound is not.
|
||||
*
|
||||
* \param lo the lower bound
|
||||
* \param hi the upper bound
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \return the value drawn from the given uniform distribution.
|
||||
*/
|
||||
PLFIT_EXPORT extern double plfit_runif(double lo, double hi, plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* Draws a sample from a uniform distribution over the [0; 1) interval.
|
||||
*
|
||||
* The interval is closed from the left and open from the right.
|
||||
*
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \return the value drawn from the given uniform distribution.
|
||||
*/
|
||||
PLFIT_EXPORT extern double plfit_runif_01(plfit_mt_rng_t* rng);
|
||||
|
||||
/**
|
||||
* Random sampler using Walker's alias method.
|
||||
*/
|
||||
typedef struct {
|
||||
long int num_bins; /**< Number of bins */
|
||||
long int* indexes; /**< Index of the "other" element in each bin */
|
||||
double* probs; /**< Probability of drawing the "own" element from a bin */
|
||||
} plfit_walker_alias_sampler_t;
|
||||
|
||||
/**
|
||||
* \brief Initializes the sampler with item probabilities.
|
||||
*
|
||||
* \param sampler the sampler to initialize
|
||||
* \param ps pointer to an array containing a value proportional to the
|
||||
* sampling probability of each item in the set being sampled.
|
||||
* \param n the number of items in the array
|
||||
* \return error code
|
||||
*/
|
||||
PLFIT_EXPORT int plfit_walker_alias_sampler_init(plfit_walker_alias_sampler_t* sampler,
|
||||
double* ps, size_t n);
|
||||
|
||||
/**
|
||||
* \brief Destroys an initialized sampler and frees the allocated memory.
|
||||
*
|
||||
* \param sampler the sampler to destroy
|
||||
*/
|
||||
PLFIT_EXPORT void plfit_walker_alias_sampler_destroy(plfit_walker_alias_sampler_t* sampler);
|
||||
|
||||
/**
|
||||
* \brief Draws a given number of samples from the sampler and writes them
|
||||
* to a given array.
|
||||
*
|
||||
* \param sampler the sampler to use
|
||||
* \param xs pointer to an array where the sampled items should be
|
||||
* written
|
||||
* \param n the number of samples to draw
|
||||
* \param rng the Mersenne Twister random number generator to use
|
||||
* \return error code
|
||||
*/
|
||||
PLFIT_EXPORT int plfit_walker_alias_sampler_sample(const plfit_walker_alias_sampler_t* sampler,
|
||||
long int* xs, size_t n, plfit_mt_rng_t* rng);
|
||||
|
||||
PLFIT_END_C_DECLS
|
||||
|
||||
#endif /* PLFIT_SAMPLING_H */
|
||||
@@ -0,0 +1,28 @@
|
||||
/* plfit_version.h
|
||||
*
|
||||
* Copyright (C) 2021 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef PLFIT_VERSION_H
|
||||
#define PLFIT_VERSION_H
|
||||
|
||||
#define PLFIT_VERSION_MAJOR 1
|
||||
#define PLFIT_VERSION_MINOR 0
|
||||
#define PLFIT_VERSION_PATCH 0
|
||||
#define PLFIT_VERSION_STRING "1.0.0"
|
||||
|
||||
#endif
|
||||
+206
@@ -0,0 +1,206 @@
|
||||
/*
|
||||
* Mathlib : A C Library of Special Functions
|
||||
* Copyright (C) 1998 Ross Ihaka
|
||||
* Copyright (C) 2000-2002 The R Core Team
|
||||
* Copyright (C) 2007 The R Foundation
|
||||
*
|
||||
* 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, a copy is available at
|
||||
* http://www.r-project.org/Licenses/
|
||||
*
|
||||
* SYNOPSIS
|
||||
*
|
||||
* #include <Rmath.h>
|
||||
* double rbinom(double nin, double pp)
|
||||
*
|
||||
* DESCRIPTION
|
||||
*
|
||||
* Random variates from the binomial distribution.
|
||||
*
|
||||
* REFERENCE
|
||||
*
|
||||
* Kachitvichyanukul, V. and Schmeiser, B. W. (1988).
|
||||
* Binomial random variate generation.
|
||||
* Communications of the ACM 31, 216-222.
|
||||
* (Algorithm BTPEC).
|
||||
*/
|
||||
|
||||
/*
|
||||
* Modifications for this file were performed by Tamas Nepusz to make it fit
|
||||
* better with plfit. The license of the original file applies to the
|
||||
* modifications as well.
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include "plfit_sampling.h"
|
||||
|
||||
#define repeat for(;;)
|
||||
|
||||
double plfit_rbinom(double nin, double pp, plfit_mt_rng_t* rng)
|
||||
{
|
||||
/* FIXME: These should become THREAD_specific globals : */
|
||||
|
||||
static double c, fm, npq, p1, p2, p3, p4, qn;
|
||||
static double xl, xll, xlr, xm, xr;
|
||||
|
||||
static double psave = -1.0;
|
||||
static int nsave = -1;
|
||||
static int m;
|
||||
|
||||
double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
|
||||
double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
|
||||
int i, ix, k, n;
|
||||
|
||||
if (!isfinite(nin)) return NAN;
|
||||
r = floor(nin + 0.5);
|
||||
if (r != nin) return NAN;
|
||||
if (!isfinite(pp) ||
|
||||
/* n=0, p=0, p=1 are not errors <TSL>*/
|
||||
r < 0 || pp < 0. || pp > 1.) return NAN;
|
||||
|
||||
if (r == 0 || pp == 0.) return 0;
|
||||
if (pp == 1.) return r;
|
||||
|
||||
n = (int) r;
|
||||
|
||||
p = fmin(pp, 1. - pp);
|
||||
q = 1. - p;
|
||||
np = n * p;
|
||||
r = p / q;
|
||||
g = r * (n + 1);
|
||||
|
||||
/* Setup, perform only when parameters change [using static (globals): */
|
||||
|
||||
/* FIXING: Want this thread safe
|
||||
-- use as little (thread globals) as possible
|
||||
*/
|
||||
if (pp != psave || n != nsave) {
|
||||
psave = pp;
|
||||
nsave = n;
|
||||
if (np < 30.0) {
|
||||
/* inverse cdf logic for mean less than 30 */
|
||||
qn = pow(q, (double) n);
|
||||
goto L_np_small;
|
||||
} else {
|
||||
ffm = np + p;
|
||||
m = (int) ffm;
|
||||
fm = m;
|
||||
npq = np * q;
|
||||
p1 = (int)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
|
||||
xm = fm + 0.5;
|
||||
xl = xm - p1;
|
||||
xr = xm + p1;
|
||||
c = 0.134 + 20.5 / (15.3 + fm);
|
||||
al = (ffm - xl) / (ffm - xl * p);
|
||||
xll = al * (1.0 + 0.5 * al);
|
||||
al = (xr - ffm) / (xr * q);
|
||||
xlr = al * (1.0 + 0.5 * al);
|
||||
p2 = p1 * (1.0 + c + c);
|
||||
p3 = p2 + c / xll;
|
||||
p4 = p3 + c / xlr;
|
||||
}
|
||||
} else if (n == nsave) {
|
||||
if (np < 30.0)
|
||||
goto L_np_small;
|
||||
}
|
||||
|
||||
/*-------------------------- np = n*p >= 30 : ------------------- */
|
||||
repeat {
|
||||
u = plfit_runif_01(rng) * p4;
|
||||
v = plfit_runif_01(rng);
|
||||
/* triangular region */
|
||||
if (u <= p1) {
|
||||
ix = (int)(xm - p1 * v + u);
|
||||
goto finis;
|
||||
}
|
||||
/* parallelogram region */
|
||||
if (u <= p2) {
|
||||
x = xl + (u - p1) / c;
|
||||
v = v * c + 1.0 - fabs(xm - x) / p1;
|
||||
if (v > 1.0 || v <= 0.)
|
||||
continue;
|
||||
ix = (int) x;
|
||||
} else {
|
||||
if (u > p3) { /* right tail */
|
||||
ix = (int)(xr - log(v) / xlr);
|
||||
if (ix > n)
|
||||
continue;
|
||||
v = v * (u - p3) * xlr;
|
||||
} else {/* left tail */
|
||||
ix = (int)(xl + log(v) / xll);
|
||||
if (ix < 0)
|
||||
continue;
|
||||
v = v * (u - p2) * xll;
|
||||
}
|
||||
}
|
||||
/* determine appropriate way to perform accept/reject test */
|
||||
k = abs(ix - m);
|
||||
if (k <= 20 || k >= npq / 2 - 1) {
|
||||
/* explicit evaluation */
|
||||
f = 1.0;
|
||||
if (m < ix) {
|
||||
for (i = m + 1; i <= ix; i++)
|
||||
f *= (g / i - r);
|
||||
} else if (m != ix) {
|
||||
for (i = ix + 1; i <= m; i++)
|
||||
f /= (g / i - r);
|
||||
}
|
||||
if (v <= f)
|
||||
goto finis;
|
||||
} else {
|
||||
/* squeezing using upper and lower bounds on log(f(x)) */
|
||||
amaxp = (k / npq) * ((k * (k / 3. + 0.625) + (1.0 / 6.0)) / npq + 0.5);
|
||||
ynorm = -k * k / (2.0 * npq);
|
||||
alv = log(v);
|
||||
if (alv < ynorm - amaxp)
|
||||
goto finis;
|
||||
if (alv <= ynorm + amaxp) {
|
||||
/* stirling's formula to machine accuracy */
|
||||
/* for the final acceptance/rejection test */
|
||||
x1 = ix + 1;
|
||||
f1 = fm + 1.0;
|
||||
z = n + 1 - fm;
|
||||
w = n - ix + 1.0;
|
||||
z2 = z * z;
|
||||
x2 = x1 * x1;
|
||||
f2 = f1 * f1;
|
||||
w2 = w * w;
|
||||
if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.)
|
||||
goto finis;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
L_np_small:
|
||||
/*---------------------- np = n*p < 30 : ------------------------- */
|
||||
|
||||
repeat {
|
||||
ix = 0;
|
||||
f = qn;
|
||||
u = plfit_runif_01(rng);
|
||||
repeat {
|
||||
if (u < f)
|
||||
goto finis;
|
||||
if (ix > 110)
|
||||
break;
|
||||
u -= f;
|
||||
ix++;
|
||||
f *= (g / ix - r);
|
||||
}
|
||||
}
|
||||
finis:
|
||||
if (psave > 0.5)
|
||||
ix = n - ix;
|
||||
return (double)ix;
|
||||
}
|
||||
+311
@@ -0,0 +1,311 @@
|
||||
/* sampling.c
|
||||
*
|
||||
* Copyright (C) 2012 Tamas Nepusz
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
|
||||
#include "igraph_random.h"
|
||||
|
||||
#include "plfit_error.h"
|
||||
#include "plfit_sampling.h"
|
||||
|
||||
inline double plfit_runif(double lo, double hi, plfit_mt_rng_t* rng) {
|
||||
if (rng == 0) {
|
||||
return RNG_UNIF(lo, hi);
|
||||
}
|
||||
return lo + plfit_mt_uniform_01(rng) * (hi-lo);
|
||||
}
|
||||
|
||||
inline double plfit_runif_01(plfit_mt_rng_t* rng) {
|
||||
if (rng == 0) {
|
||||
return RNG_UNIF01();
|
||||
}
|
||||
return plfit_mt_uniform_01(rng);
|
||||
}
|
||||
|
||||
inline double plfit_rpareto(double xmin, double alpha, plfit_mt_rng_t* rng) {
|
||||
if (alpha <= 0 || xmin <= 0)
|
||||
return NAN;
|
||||
|
||||
/* 1-u is used in the base here because we want to avoid the case of
|
||||
* sampling zero */
|
||||
return pow(1-plfit_runif_01(rng), -1.0 / alpha) * xmin;
|
||||
}
|
||||
|
||||
int plfit_rpareto_array(double xmin, double alpha, size_t n, plfit_mt_rng_t* rng,
|
||||
double* result) {
|
||||
double gamma;
|
||||
|
||||
if (alpha <= 0 || xmin <= 0)
|
||||
return PLFIT_EINVAL;
|
||||
|
||||
if (result == 0 || n == 0)
|
||||
return PLFIT_SUCCESS;
|
||||
|
||||
gamma = -1.0 / alpha;
|
||||
while (n > 0) {
|
||||
/* 1-u is used in the base here because we want to avoid the case of
|
||||
* sampling zero */
|
||||
*result = pow(1-plfit_runif_01(rng), gamma) * xmin;
|
||||
result++; n--;
|
||||
}
|
||||
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
|
||||
inline double plfit_rzeta(long int xmin, double alpha, plfit_mt_rng_t* rng) {
|
||||
double u, v, t;
|
||||
long int x;
|
||||
double alpha_minus_1 = alpha-1;
|
||||
double minus_1_over_alpha_minus_1 = -1.0 / (alpha-1);
|
||||
double b;
|
||||
double one_over_b_minus_1;
|
||||
|
||||
if (alpha <= 0 || xmin < 1)
|
||||
return NAN;
|
||||
|
||||
xmin = (long int) round(xmin);
|
||||
|
||||
/* Rejection sampling for the win. We use Y=floor(U^{-1/alpha} * xmin) as the
|
||||
* envelope distribution, similarly to Chapter X.6 of Luc Devroye's book
|
||||
* (where xmin is assumed to be 1): http://luc.devroye.org/chapter_ten.pdf
|
||||
*
|
||||
* Some notes that should help me recover what I was doing:
|
||||
*
|
||||
* p_i = 1/zeta(alpha, xmin) * i^-alpha
|
||||
* q_i = (xmin/i)^{alpha-1} - (xmin/(i+1))^{alpha-1}
|
||||
* = (i/xmin)^{1-alpha} - ((i+1)/xmin)^{1-alpha}
|
||||
* = [i^{1-alpha} - (i+1)^{1-alpha}] / xmin^{1-alpha}
|
||||
*
|
||||
* p_i / q_i attains its maximum at xmin=i, so the rejection constant is:
|
||||
*
|
||||
* c = p_xmin / q_xmin
|
||||
*
|
||||
* We have to accept the sample if V <= (p_i / q_i) * (q_xmin / p_xmin) =
|
||||
* (i/xmin)^-alpha * [xmin^{1-alpha} - (xmin+1)^{1-alpha}] / [i^{1-alpha} - (i+1)^{1-alpha}] =
|
||||
* [xmin - xmin^alpha / (xmin+1)^{alpha-1}] / [i - i^alpha / (i+1)^{alpha-1}] =
|
||||
* xmin/i * [1-(xmin/(xmin+1))^{alpha-1}]/[1-(i/(i+1))^{alpha-1}]
|
||||
*
|
||||
* In other words (and substituting i with X, which is the same),
|
||||
*
|
||||
* V * (X/xmin) <= [1 - (1+1/xmin)^{1-alpha}] / [1 - (1+1/i)^{1-alpha}]
|
||||
*
|
||||
* Let b := (1+1/xmin)^{alpha-1} and let T := (1+1/i)^{alpha-1}. Then:
|
||||
*
|
||||
* V * (X/xmin) <= [(b-1)/b] / [(T-1)/T]
|
||||
* V * (X/xmin) * (T-1) / (b-1) <= T / b
|
||||
*
|
||||
* which is the same as in Devroye's book, except for the X/xmin term, and
|
||||
* the definition of b.
|
||||
*/
|
||||
b = pow(1 + 1.0/xmin, alpha_minus_1);
|
||||
one_over_b_minus_1 = 1.0/(b-1);
|
||||
do {
|
||||
do {
|
||||
u = plfit_runif_01(rng);
|
||||
v = plfit_runif_01(rng);
|
||||
/* 1-u is used in the base here because we want to avoid the case of
|
||||
* having zero in x */
|
||||
x = (long int) floor(pow(1-u, minus_1_over_alpha_minus_1) * xmin);
|
||||
} while (x < xmin);
|
||||
t = pow((x+1.0)/x, alpha_minus_1);
|
||||
} while (v*x*(t-1)*one_over_b_minus_1*b > t*xmin);
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
int plfit_rzeta_array(long int xmin, double alpha, size_t n, plfit_mt_rng_t* rng,
|
||||
double* result) {
|
||||
double u, v, t;
|
||||
long int x;
|
||||
double alpha_minus_1 = alpha-1;
|
||||
double minus_1_over_alpha_minus_1 = -1.0 / (alpha-1);
|
||||
double b, one_over_b_minus_1;
|
||||
|
||||
if (alpha <= 0 || xmin < 1)
|
||||
return PLFIT_EINVAL;
|
||||
|
||||
if (result == 0 || n == 0)
|
||||
return PLFIT_SUCCESS;
|
||||
|
||||
/* See the comments in plfit_rzeta for an explanation of the algorithm
|
||||
* below. */
|
||||
xmin = (long int) round(xmin);
|
||||
b = pow(1 + 1.0/xmin, alpha_minus_1);
|
||||
one_over_b_minus_1 = 1.0/(b-1);
|
||||
|
||||
while (n > 0) {
|
||||
do {
|
||||
do {
|
||||
u = plfit_runif_01(rng);
|
||||
v = plfit_runif_01(rng);
|
||||
/* 1-u is used in the base here because we want to avoid the case of
|
||||
* having zero in x */
|
||||
x = (long int) floor(pow(1-u, minus_1_over_alpha_minus_1) * xmin);
|
||||
} while (x < xmin); /* handles overflow as well */
|
||||
t = pow((x+1.0)/x, alpha_minus_1);
|
||||
} while (v*x*(t-1)*one_over_b_minus_1*b > t*xmin);
|
||||
*result = x;
|
||||
if (x < 0) return PLFIT_EINVAL;
|
||||
result++; n--;
|
||||
}
|
||||
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
|
||||
int plfit_walker_alias_sampler_init(plfit_walker_alias_sampler_t* sampler,
|
||||
double* ps, size_t n) {
|
||||
double *p, *p2, *ps_end;
|
||||
double sum;
|
||||
long int *short_sticks, *long_sticks;
|
||||
long int num_short_sticks, num_long_sticks;
|
||||
long int i;
|
||||
|
||||
if (n > LONG_MAX) {
|
||||
return PLFIT_EINVAL;
|
||||
}
|
||||
|
||||
sampler->num_bins = (long int) n;
|
||||
|
||||
ps_end = ps + n;
|
||||
|
||||
/* Initialize indexes and probs */
|
||||
sampler->indexes = (long int*)calloc(n > 0 ? n : 1, sizeof(long int));
|
||||
if (sampler->indexes == NULL) {
|
||||
return PLFIT_ENOMEM;
|
||||
}
|
||||
sampler->probs = (double*)calloc(n > 0 ? n : 1, sizeof(double));
|
||||
if (sampler->probs == NULL) {
|
||||
free(sampler->indexes);
|
||||
return PLFIT_ENOMEM;
|
||||
}
|
||||
|
||||
/* Normalize the probability vector; count how many short and long sticks
|
||||
* are there initially */
|
||||
for (sum = 0.0, p = ps; p != ps_end; p++) {
|
||||
sum += *p;
|
||||
}
|
||||
sum = n / sum;
|
||||
|
||||
num_short_sticks = num_long_sticks = 0;
|
||||
for (p = ps, p2 = sampler->probs; p != ps_end; p++, p2++) {
|
||||
*p2 = *p * sum;
|
||||
if (*p2 < 1) {
|
||||
num_short_sticks++;
|
||||
} else if (*p2 > 1) {
|
||||
num_long_sticks++;
|
||||
}
|
||||
}
|
||||
|
||||
/* Allocate space for short & long stick indexes */
|
||||
long_sticks = (long int*)calloc(num_long_sticks > 0 ? num_long_sticks : 1, sizeof(long int));
|
||||
if (long_sticks == NULL) {
|
||||
free(sampler->probs);
|
||||
free(sampler->indexes);
|
||||
return PLFIT_ENOMEM;
|
||||
}
|
||||
short_sticks = (long int*)calloc(num_short_sticks > 0 ? num_short_sticks : 1, sizeof(long int));
|
||||
if (short_sticks == NULL) {
|
||||
free(sampler->probs);
|
||||
free(sampler->indexes);
|
||||
free(long_sticks);
|
||||
return PLFIT_ENOMEM;
|
||||
}
|
||||
|
||||
/* Initialize short_sticks and long_sticks */
|
||||
num_short_sticks = num_long_sticks = 0;
|
||||
for (i = 0, p = sampler->probs; i < n; i++, p++) {
|
||||
if (*p < 1) {
|
||||
short_sticks[num_short_sticks++] = i;
|
||||
} else if (*p > 1) {
|
||||
long_sticks[num_long_sticks++] = i;
|
||||
}
|
||||
}
|
||||
|
||||
/* Prepare the index table */
|
||||
while (num_short_sticks && num_long_sticks) {
|
||||
long int short_index, long_index;
|
||||
short_index = short_sticks[--num_short_sticks];
|
||||
long_index = long_sticks[num_long_sticks-1];
|
||||
sampler->indexes[short_index] = long_index;
|
||||
sampler->probs[long_index] = /* numerical stability */
|
||||
(sampler->probs[long_index] + sampler->probs[short_index]) - 1;
|
||||
if (sampler->probs[long_index] < 1) {
|
||||
short_sticks[num_short_sticks++] = long_index;
|
||||
num_long_sticks--;
|
||||
}
|
||||
}
|
||||
|
||||
/* Fix numerical stability issues */
|
||||
while (num_long_sticks) {
|
||||
i = long_sticks[--num_long_sticks];
|
||||
sampler->probs[i] = 1;
|
||||
}
|
||||
while (num_short_sticks) {
|
||||
i = short_sticks[--num_short_sticks];
|
||||
sampler->probs[i] = 1;
|
||||
}
|
||||
|
||||
free(short_sticks);
|
||||
free(long_sticks);
|
||||
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
||||
void plfit_walker_alias_sampler_destroy(plfit_walker_alias_sampler_t* sampler) {
|
||||
if (sampler->indexes) {
|
||||
free(sampler->indexes);
|
||||
sampler->indexes = 0;
|
||||
}
|
||||
if (sampler->probs) {
|
||||
free(sampler->probs);
|
||||
sampler->probs = 0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int plfit_walker_alias_sampler_sample(const plfit_walker_alias_sampler_t* sampler,
|
||||
long int *xs, size_t n, plfit_mt_rng_t* rng) {
|
||||
double u;
|
||||
long int j;
|
||||
long int *x;
|
||||
|
||||
x = xs;
|
||||
|
||||
if (rng == 0) {
|
||||
/* Using built-in RNG */
|
||||
while (n > 0) {
|
||||
u = RNG_UNIF01();
|
||||
j = RNG_INTEGER(0, sampler->num_bins - 1);
|
||||
*x = (u < sampler->probs[j]) ? j : sampler->indexes[j];
|
||||
n--; x++;
|
||||
}
|
||||
} else {
|
||||
/* Using Mersenne Twister */
|
||||
while (n > 0) {
|
||||
u = plfit_mt_uniform_01(rng);
|
||||
j = plfit_mt_random(rng) % sampler->num_bins;
|
||||
*x = (u < sampler->probs[j]) ? j : sampler->indexes[j];
|
||||
n--; x++;
|
||||
}
|
||||
}
|
||||
|
||||
return PLFIT_SUCCESS;
|
||||
}
|
||||
Reference in New Issue
Block a user