logmath.c
/* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
/* ====================================================================
* Copyright (c) 1999-2007 Carnegie Mellon University. All rights
* reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in
* the documentation and/or other materials provided with the
* distribution.
*
* This work was supported in part by funding from the Defense Advanced
* Research Projects Agency and the National Science Foundation of the
* United States of America, and the CMU Sphinx Speech Consortium.
*
* THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
* ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
* THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
* NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* ====================================================================
*
*/
#include <math.h>
#include <string.h>
#include <assert.h>
#include "sphinxbase/logmath.h"
#include "sphinxbase/err.h"
#include "sphinxbase/ckd_alloc.h"
#include "sphinxbase/mmio.h"
#include "sphinxbase/bio.h"
#include "sphinxbase/strfuncs.h"
struct logmath_s {
logadd_t t;
int refcount;
mmio_file_t *filemap;
float64 base;
float64 log_of_base;
float64 log10_of_base;
float64 inv_log_of_base;
float64 inv_log10_of_base;
int32 zero;
};
logmath_t *
logmath_init(float64 base, int shift, int use_table)
{
logmath_t *lmath;
uint32 maxyx, i;
float64 byx;
int width;
/* Check that the base is correct. */
if (base <= 1.0) {
E_ERROR("Base must be greater than 1.0\n");
return NULL;
}
/* Set up various necessary constants. */
lmath = ckd_calloc(1, sizeof(*lmath));
lmath->refcount = 1;
lmath->base = base;
lmath->log_of_base = log(base);
lmath->log10_of_base = log10(base);
lmath->inv_log_of_base = 1.0/lmath->log_of_base;
lmath->inv_log10_of_base = 1.0/lmath->log10_of_base;
lmath->t.shift = shift;
/* Shift this sufficiently that overflows can be avoided. */
lmath->zero = MAX_NEG_INT32 >> (shift + 2);
if (!use_table)
return lmath;
/* Create a logadd table with the appropriate width */
maxyx = (uint32) (log(2.0) / log(base) + 0.5) >> shift;
/* Poor man's log2 */
if (maxyx < 256) width = 1;
else if (maxyx < 65536) width = 2;
else width = 4;
lmath->t.width = width;
/* Figure out size of add table required. */
byx = 1.0; /* Maximum possible base^{y-x} value - note that this implies that y-x == 0 */
for (i = 0;; ++i) {
float64 lobyx = log(1.0 + byx) * lmath->inv_log_of_base; /* log_{base}(1 + base^{y-x}); */
int32 k = (int32) (lobyx + 0.5 * (1<<shift)) >> shift; /* Round to shift */
/* base^{y-x} has reached the smallest representable value. */
if (k <= 0)
break;
/* This table is indexed by -(y-x), so we multiply byx by
* base^{-1} here which is equivalent to subtracting one from
* (y-x). */
byx /= base;
}
i >>= shift;
/* Never produce a table smaller than 256 entries. */
if (i < 255) i = 255;
lmath->t.table = ckd_calloc(i+1, width);
lmath->t.table_size = i + 1;
/* Create the add table (see above). */
byx = 1.0;
for (i = 0;; ++i) {
float64 lobyx = log(1.0 + byx) * lmath->inv_log_of_base;
int32 k = (int32) (lobyx + 0.5 * (1<<shift)) >> shift; /* Round to shift */
uint32 prev = 0;
/* Check any previous value - if there is a shift, we want to
* only store the highest one. */
switch (width) {
case 1:
prev = ((uint8 *)lmath->t.table)[i >> shift];
break;
case 2:
prev = ((uint16 *)lmath->t.table)[i >> shift];
break;
case 4:
prev = ((uint32 *)lmath->t.table)[i >> shift];
break;
}
if (prev == 0) {
switch (width) {
case 1:
((uint8 *)lmath->t.table)[i >> shift] = (uint8) k;
break;
case 2:
((uint16 *)lmath->t.table)[i >> shift] = (uint16) k;
break;
case 4:
((uint32 *)lmath->t.table)[i >> shift] = (uint32) k;
break;
}
}
if (k <= 0)
break;
/* Decay base^{y-x} exponentially according to base. */
byx /= base;
}
return lmath;
}
logmath_t *
logmath_read(const char *file_name)
{
logmath_t *lmath;
char **argname, **argval;
int32 byteswap, i;
int chksum_present, do_mmap;
uint32 chksum;
long pos;
FILE *fp;
E_INFO("Reading log table file '%s'\n", file_name);
if ((fp = fopen(file_name, "rb")) == NULL) {
E_ERROR_SYSTEM("Failed to open log table file '%s' for reading", file_name);
return NULL;
}
/* Read header, including argument-value info and 32-bit byteorder magic */
if (bio_readhdr(fp, &argname, &argval, &byteswap) < 0) {
E_ERROR("Failed to read the header from the file '%s'\n", file_name);
fclose(fp);
return NULL;
}
lmath = ckd_calloc(1, sizeof(*lmath));
/* Default values. */
lmath->t.shift = 0;
lmath->t.width = 2;
lmath->base = 1.0001;
/* Parse argument-value list */
chksum_present = 0;
for (i = 0; argname[i]; i++) {
if (strcmp(argname[i], "version") == 0) {
}
else if (strcmp(argname[i], "chksum0") == 0) {
if (strcmp(argval[i], "yes") == 0)
chksum_present = 1;
}
else if (strcmp(argname[i], "width") == 0) {
lmath->t.width = atoi(argval[i]);
}
else if (strcmp(argname[i], "shift") == 0) {
lmath->t.shift = atoi(argval[i]);
}
else if (strcmp(argname[i], "logbase") == 0) {
lmath->base = atof_c(argval[i]);
}
}
bio_hdrarg_free(argname, argval);
chksum = 0;
/* Set up various necessary constants. */
lmath->log_of_base = log(lmath->base);
lmath->log10_of_base = log10(lmath->base);
lmath->inv_log_of_base = 1.0/lmath->log_of_base;
lmath->inv_log10_of_base = 1.0/lmath->log10_of_base;
/* Shift this sufficiently that overflows can be avoided. */
lmath->zero = MAX_NEG_INT32 >> (lmath->t.shift + 2);
/* #Values to follow */
if (bio_fread(&lmath->t.table_size, sizeof(int32), 1, fp, byteswap, &chksum) != 1) {
E_ERROR("Failed to read values from the file '%s'", file_name);
goto error_out;
}
/* Check alignment constraints for memory mapping */
do_mmap = 1;
pos = ftell(fp);
if (pos & ((long)lmath->t.width - 1)) {
E_WARN("%s: Data start %ld is not aligned on %d-byte boundary, will not memory map\n",
file_name, pos, lmath->t.width);
do_mmap = 0;
}
/* Check byte order for memory mapping */
if (byteswap) {
E_WARN("%s: Data is wrong-endian, will not memory map\n", file_name);
do_mmap = 0;
}
if (do_mmap) {
lmath->filemap = mmio_file_read(file_name);
lmath->t.table = (char *)mmio_file_ptr(lmath->filemap) + pos;
}
else {
lmath->t.table = ckd_calloc(lmath->t.table_size, lmath->t.width);
if (bio_fread(lmath->t.table, lmath->t.width, lmath->t.table_size,
fp, byteswap, &chksum) != lmath->t.table_size) {
E_ERROR("Failed to read data (%d x %d bytes) from the file '%s' failed",
lmath->t.table_size, lmath->t.width, file_name);
goto error_out;
}
if (chksum_present)
bio_verify_chksum(fp, byteswap, chksum);
if (fread(&i, 1, 1, fp) == 1) {
E_ERROR("%s: More data than expected\n", file_name);
goto error_out;
}
}
fclose(fp);
return lmath;
error_out:
logmath_free(lmath);
return NULL;
}
int32
logmath_write(logmath_t *lmath, const char *file_name)
{
FILE *fp;
long pos;
uint32 chksum;
if (lmath->t.table == NULL) {
E_ERROR("No log table to write!\n");
return -1;
}
E_INFO("Writing log table file '%s'\n", file_name);
if ((fp = fopen(file_name, "wb")) == NULL) {
E_ERROR_SYSTEM("Failed to open logtable file '%s' for writing", file_name);
return -1;
}
/* For whatever reason, we have to do this manually at the
* moment. */
fprintf(fp, "s3\nversion 1.0\nchksum0 yes\n");
fprintf(fp, "width %d\n", lmath->t.width);
fprintf(fp, "shift %d\n", lmath->t.shift);
fprintf(fp, "logbase %f\n", lmath->base);
/* Pad it out to ensure alignment. */
pos = ftell(fp) + strlen("endhdr\n");
if (pos & ((long)lmath->t.width - 1)) {
size_t align = lmath->t.width - (pos & ((long)lmath->t.width - 1));
assert(lmath->t.width <= 8);
fwrite(" " /* 8 spaces */, 1, align, fp);
}
fprintf(fp, "endhdr\n");
/* Now write the binary data. */
chksum = (uint32)BYTE_ORDER_MAGIC;
fwrite(&chksum, sizeof(uint32), 1, fp);
chksum = 0;
/* #Values to follow */
if (bio_fwrite(&lmath->t.table_size, sizeof(uint32),
1, fp, 0, &chksum) != 1) {
E_ERROR("Failed to write data to a file '%s'", file_name);
goto error_out;
}
if (bio_fwrite(lmath->t.table, lmath->t.width, lmath->t.table_size,
fp, 0, &chksum) != lmath->t.table_size) {
E_ERROR("Failed to write data (%d x %d bytes) to the file '%s'",
lmath->t.table_size, lmath->t.width, file_name);
goto error_out;
}
if (bio_fwrite(&chksum, sizeof(uint32), 1, fp, 0, NULL) != 1) {
E_ERROR("Failed to write checksum to the file '%s'", file_name);
goto error_out;
}
fclose(fp);
return 0;
error_out:
fclose(fp);
return -1;
}
logmath_t *
logmath_retain(logmath_t *lmath)
{
++lmath->refcount;
return lmath;
}
int
logmath_free(logmath_t *lmath)
{
if (lmath == NULL)
return 0;
if (--lmath->refcount > 0)
return lmath->refcount;
if (lmath->filemap)
mmio_file_unmap(lmath->filemap);
else
ckd_free(lmath->t.table);
ckd_free(lmath);
return 0;
}
int32
logmath_get_table_shape(logmath_t *lmath, uint32 *out_size,
uint32 *out_width, uint32 *out_shift)
{
if (out_size) *out_size = lmath->t.table_size;
if (out_width) *out_width = lmath->t.width;
if (out_shift) *out_shift = lmath->t.shift;
return lmath->t.table_size * lmath->t.width;
}
float64
logmath_get_base(logmath_t *lmath)
{
return lmath->base;
}
int
logmath_get_zero(logmath_t *lmath)
{
return lmath->zero;
}
int
logmath_get_width(logmath_t *lmath)
{
return lmath->t.width;
}
int
logmath_get_shift(logmath_t *lmath)
{
return lmath->t.shift;
}
int
logmath_add(logmath_t *lmath, int logb_x, int logb_y)
{
logadd_t *t = LOGMATH_TABLE(lmath);
int d, r;
/* handle 0 + x = x case. */
if (logb_x <= lmath->zero)
return logb_y;
if (logb_y <= lmath->zero)
return logb_x;
if (t->table == NULL)
return logmath_add_exact(lmath, logb_x, logb_y);
/* d must be positive, obviously. */
if (logb_x > logb_y) {
d = (logb_x - logb_y);
r = logb_x;
}
else {
d = (logb_y - logb_x);
r = logb_y;
}
if (d < 0) {
/* Some kind of overflow has occurred, fail gracefully. */
return r;
}
if ((size_t)d >= t->table_size) {
/* If this happens, it's not actually an error, because the
* last entry in the logadd table is guaranteed to be zero.
* Therefore we just return the larger of the two values. */
return r;
}
switch (t->width) {
case 1:
return r + (((uint8 *)t->table)[d]);
case 2:
return r + (((uint16 *)t->table)[d]);
case 4:
return r + (((uint32 *)t->table)[d]);
}
return r;
}
int
logmath_add_exact(logmath_t *lmath, int logb_p, int logb_q)
{
return logmath_log(lmath,
logmath_exp(lmath, logb_p)
+ logmath_exp(lmath, logb_q));
}
int
logmath_log(logmath_t *lmath, float64 p)
{
if (p <= 0) {
return lmath->zero;
}
return (int)(log(p) * lmath->inv_log_of_base) >> lmath->t.shift;
}
float64
logmath_exp(logmath_t *lmath, int logb_p)
{
return pow(lmath->base, (float64)(logb_p << lmath->t.shift));
}
int
logmath_ln_to_log(logmath_t *lmath, float64 log_p)
{
return (int)(log_p * lmath->inv_log_of_base) >> lmath->t.shift;
}
float64
logmath_log_to_ln(logmath_t *lmath, int logb_p)
{
return (float64)(logb_p << lmath->t.shift) * lmath->log_of_base;
}
int
logmath_log10_to_log(logmath_t *lmath, float64 log_p)
{
return (int)(log_p * lmath->inv_log10_of_base) >> lmath->t.shift;
}
float
logmath_log10_to_log_float(logmath_t *lmath, float64 log_p)
{
int i;
float res = (float)(log_p * lmath->inv_log10_of_base);
for (i = 0; i < lmath->t.shift; i++)
res /= 2.0f;
return res;
}
float64
logmath_log_to_log10(logmath_t *lmath, int logb_p)
{
return (float64)(logb_p << lmath->t.shift) * lmath->log10_of_base;
}
float64
logmath_log_float_to_log10(logmath_t *lmath, float log_p)
{
int i;
for (i = 0; i < lmath->t.shift; i++) {
log_p *= 2;
}
return log_p * lmath->log10_of_base;
}
bio.h
Cross platform binary IO to process files in sphinx3 format.
bio_fwrite
SPHINXBASE_EXPORT int32 bio_fwrite(const void *buf, int32 el_sz, int32 n_el, FILE *fp, int32 swap, uint32 *chksum)
Like fwrite but perform byteswapping and accumulate checksum (the 2 extra arguments).
Definition bio.c:342
bio_fread
SPHINXBASE_EXPORT int32 bio_fread(void *buf, int32 el_sz, int32 n_el, FILE *fp, int32 swap, uint32 *chksum)
Like fread but perform byteswapping and accumulate checksum (the 2 extra arguments).
Definition bio.c:326
bio_readhdr
SPHINXBASE_EXPORT int32 bio_readhdr(FILE *fp, char ***name, char ***val, int32 *swap)
Read binary file format header: has the following format.
Definition bio.c:187
bio_verify_chksum
SPHINXBASE_EXPORT void bio_verify_chksum(FILE *fp, int32 byteswap, uint32 chksum)
Read and verify checksum at the end of binary file.
Definition bio.c:492
bio_hdrarg_free
SPHINXBASE_EXPORT void bio_hdrarg_free(char **name, char **val)
Free name and value strings previously allocated and returned by bio_readhdr.
Definition bio.c:121
ckd_alloc.h
Sphinx's memory allocation/deallocation routines.
ckd_free
SPHINXBASE_EXPORT void ckd_free(void *ptr)
Test and free a 1-D array.
Definition ckd_alloc.c:244
ckd_calloc
#define ckd_calloc(n, sz)
Macros to simplify the use of above functions.
Definition ckd_alloc.h:248
err.h
Implementation of logging routines.
E_ERROR
#define E_ERROR(...)
Print error message to error log.
Definition err.h:104
E_INFO
#define E_INFO(...)
Print logging information to standard error stream.
Definition err.h:114
E_ERROR_SYSTEM
#define E_ERROR_SYSTEM(...)
Print error text; Call perror("");.
Definition err.h:99
E_WARN
#define E_WARN(...)
Print warning message to error log.
Definition err.h:109
logmath.h
Fast integer logarithmic addition operations.
logmath_get_width
SPHINXBASE_EXPORT int logmath_get_width(logmath_t *lmath)
Get the width of the values in a log table.
Definition logmath.c:380
logmath_get_zero
SPHINXBASE_EXPORT int logmath_get_zero(logmath_t *lmath)
Get the smallest possible value represented in this base.
Definition logmath.c:374
logmath_log10_to_log_float
SPHINXBASE_EXPORT float logmath_log10_to_log_float(logmath_t *lmath, float64 log_p)
Convert base 10 log (in floating point) to float log in base B.
Definition logmath.c:480
logmath_log_float_to_log10
SPHINXBASE_EXPORT float64 logmath_log_float_to_log10(logmath_t *lmath, float log_p)
Convert float log in base B to base 10 log.
Definition logmath.c:496
logmath_get_table_shape
SPHINXBASE_EXPORT int32 logmath_get_table_shape(logmath_t *lmath, uint32 *out_size, uint32 *out_width, uint32 *out_shift)
Get the log table size and dimensions.
Definition logmath.c:357
logmath_ln_to_log
SPHINXBASE_EXPORT int logmath_ln_to_log(logmath_t *lmath, float64 log_p)
Convert natural log (in floating point) to integer log in base B.
Definition logmath.c:462
logmath_read
SPHINXBASE_EXPORT logmath_t * logmath_read(const char *filename)
Memory-map (or read) a log table from a file.
Definition logmath.c:164
logmath_add
SPHINXBASE_EXPORT int logmath_add(logmath_t *lmath, int logb_p, int logb_q)
Add two values in log space (i.e.
Definition logmath.c:392
logmath_get_base
SPHINXBASE_EXPORT float64 logmath_get_base(logmath_t *lmath)
Get the log base.
Definition logmath.c:368
logmath_add_exact
SPHINXBASE_EXPORT int logmath_add_exact(logmath_t *lmath, int logb_p, int logb_q)
Add two values in log space exactly and slowly (without using add table).
Definition logmath.c:439
logmath_write
SPHINXBASE_EXPORT int32 logmath_write(logmath_t *lmath, const char *filename)
Write a log table to a file.
Definition logmath.c:272
logmath_log_to_log10
SPHINXBASE_EXPORT float64 logmath_log_to_log10(logmath_t *lmath, int logb_p)
Convert integer log in base B to base 10 log (in floating point).
Definition logmath.c:490
logmath_log_to_ln
SPHINXBASE_EXPORT float64 logmath_log_to_ln(logmath_t *lmath, int logb_p)
Convert integer log in base B to natural log (in floating point).
Definition logmath.c:468
logadd_t
struct logadd_s logadd_t
Integer log math computation table.
Definition logmath.h:93
logmath_free
SPHINXBASE_EXPORT int logmath_free(logmath_t *lmath)
Free a log table.
Definition logmath.c:342
logmath_log10_to_log
SPHINXBASE_EXPORT int logmath_log10_to_log(logmath_t *lmath, float64 log_p)
Convert base 10 log (in floating point) to integer log in base B.
Definition logmath.c:474
logmath_init
SPHINXBASE_EXPORT logmath_t * logmath_init(float64 base, int shift, int use_table)
Initialize a log math computation table.
Definition logmath.c:62
logmath_retain
SPHINXBASE_EXPORT logmath_t * logmath_retain(logmath_t *lmath)
Retain ownership of a log table.
Definition logmath.c:335
LOGMATH_TABLE
#define LOGMATH_TABLE(lm)
Obtain the log-add table from a logmath_t *.
Definition logmath.h:113
logmath_t
struct logmath_s logmath_t
Integer log math computation class.
Definition logmath.h:108
logmath_exp
SPHINXBASE_EXPORT float64 logmath_exp(logmath_t *lmath, int logb_p)
Convert integer log in base B to linear floating point.
Definition logmath.c:456
logmath_log
SPHINXBASE_EXPORT int logmath_log(logmath_t *lmath, float64 p)
Convert linear floating point number to integer log in base B.
Definition logmath.c:447
logmath_get_shift
SPHINXBASE_EXPORT int logmath_get_shift(logmath_t *lmath)
Get the shift of the values in a log table.
Definition logmath.c:386
mmio.h
Memory-mapped I/O wrappers for files.
mmio_file_unmap
SPHINXBASE_EXPORT void mmio_file_unmap(mmio_file_t *mf)
Unmap a file, releasing memory associated with it.
Definition mmio.c:241
mmio_file_read
SPHINXBASE_EXPORT mmio_file_t * mmio_file_read(const char *filename)
Memory-map a file for reading.
Definition mmio.c:207
mmio_file_ptr
SPHINXBASE_EXPORT void * mmio_file_ptr(mmio_file_t *mf)
Get a pointer to the memory mapped for a file.
Definition mmio.c:252
mmio_file_t
struct mmio_file_s mmio_file_t
Abstract structure representing a memory-mapped file.
Definition mmio.h:59
strfuncs.h
Miscellaneous useful string functions.
atof_c
SPHINXBASE_EXPORT double atof_c(char const *str)
Locale independent version of atof().
Definition strfuncs.c:55
logadd_s::table
void * table
Table, in unsigned integers of (width) bytes.
Definition logmath.h:96
logadd_s::shift
int8 shift
Right shift applied to elements in (table).
Definition logmath.h:102
logadd_s::table_size
uint32 table_size
Number of elements in (table).
Definition logmath.h:98
logadd_s::width
uint8 width
Width of elements of (table).
Definition logmath.h:100
logmath_s
Definition logmath.c:49