src/libsphinxbase/util/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 "logmath.h"
#include "err.h"
#include "ckd_alloc.h"
#include "mmio.h"
#include "bio.h"
#include "strfuncs.h"

#include <math.h>
#include <string.h>
#include <assert.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;
    lmath->zero = logmath_log(lmath, 1e-300);

    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("fopen(%s,rb) failed\n", 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("bio_readhdr(%s) failed\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;


    /* #Values to follow */
    if (bio_fread(&lmath->t.table_size, sizeof(int32), 1, fp, byteswap, &chksum) != 1) {
        E_ERROR("fread(%s) (total #values) failed\n", 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("fread(%s) (%d x %d bytes) failed\n",
                    file_name, lmath->t.table_size, lmath->t.width);
            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("fopen(%s,wb) failed\n", 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("fwrite(%s) (total #values) failed\n", 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("fwrite(%s) (%d x %d bytes) failed\n",
                file_name, lmath->t.table_size, lmath->t.width);
        goto error_out;
    }
    if (bio_fwrite(&chksum, sizeof(uint32), 1, fp, 0, NULL) != 1) {
        E_ERROR("fwrite(%s) checksum failed\n", 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;

    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;
}

float64
logmath_log_to_log10(logmath_t *lmath, int logb_p)
{
    return (float64)(logb_p << lmath->t.shift) * lmath->log10_of_base;
}

---