504 lines
14 KiB
C
504 lines
14 KiB
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;
|
|
}
|