The R Project SVN R

/* Algorithm AS 51 Appl. Statist. (1972), vol. 21, p. 218
   original (C) Royal Statistical Society 1972

   Performs an iterative proportional fit of the marginal totals of a
   contingency table.
*/

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include <math.h>

#include <stdio.h>
#include <R_ext/Memory.h>
#include <R_ext/Applic.h>

#undef max
#undef min
#undef abs
#define max(a, b)       ((a) < (b) ? (b) : (a))
#define min(a, b)       ((a) > (b) ? (b) : (a))
#define abs(x)          ((x) >= 0 ? (x) : -(x))

static void collap(int *nvar, double *x, double *y, int *locy,
           int *nx, int *ny, int *dim, int *config);
static void adjust(int *nvar, double *x, double *y, double *z,
           int *locz, int *nx, int *ny, int *nz, int *dim,
           int *config, double *d);
static int *lvector(int n);

/* Table of constant values */

static int c__1 = 1;

void loglin(int *nvar, int *dim, int *ncon, int *config, int *ntab,
        double *table, double *fit, int *locmar, int *nmar, double *marg,
        int *nu, double *u, double *maxdev, int *maxit,
        double *dev, int *nlast, int *ifault)
{
    int i, j, k, n, point, size, *icon, *check;
    double x, y, xmax;

    check = lvector(*nvar);
    icon = lvector(*nvar);

    /* Parameter adjustments */
    --dim;
    --locmar;
    config -= *nvar + 1;
    --fit;
    --table;
    --marg;
    --u;
    --dev;

    /* Function body */

    *ifault = 0;
    *nlast = 0;

    /* Check validity of NVAR, the number of variables, and of maxit,
       the maximum number of iterations */

    if (*nvar > 0 && *maxit > 0) {
    goto L10;
    }
L5:
    *ifault = 4;
    return;

    /* Look at table and fit constants */

L10:
    size = 1;
    for (j = 1; j <= *nvar; ++j) {
    if (dim[j] <= 0) {
        goto L5;
    }
    size *= dim[j];
    }
    if (size <= *ntab) {
    goto L40;
    }
L35:
    *ifault = 2;
    return;
L40:
    x = 0.;
    y = 0.;
    for (i = 1; i <= size; ++i) {
    if (table[i] < 0. || fit[i] < 0.) {
        goto L5;
    }
    x += table[i];
    y += fit[i];
    }

    /* Make a preliminary adjustment to obtain the fit to an empty
       configuration list */

    if (y == 0.) {
    goto L5;
    }
    x /= y;
    for (i = 1; i <= size; ++i) {
    fit[i] = x * fit[i];
    }
    if (*ncon <= 0 || config[*nvar + 1] == 0) {
    return;
    }

    /* Allocate marginal tables */

    point = 1;
    for (i = 1; i <= *ncon; ++i) {
    /* A zero beginning a configuration indicates that the list is
       completed */
    if (config[i * *nvar + 1] == 0) {
        goto L160;
    }
    /* Get marginal table size.  While doing this task, see if the
       configuration list contains duplications or elements out of
       range. */
    size = 1;
    for (j = 0; j < *nvar; ++j) {
        check[j] = 0;
    }
    for (j = 1; j <= *nvar; ++j) {
        k = config[j + i * *nvar];
        /* A zero indicates the end of the string. */
        if (k == 0) {
        goto L130;
        }
        /* See if element is valid. */
        if (k >= 0 && k <= *nvar) {
        goto L100;
        }
L95:
        *ifault = 1;
        return;
        /* Check for duplication */
L100:
        if (check[k - 1]) {
        goto L95;
        }
        check[k - 1] = (1);
        /* Get size */
        size *= dim[k];
    }

    /* Since U is used to store fitted marginals, size must not
       exceed NU */
L130:
    if (size > *nu) {
        goto L35;
    }

    /* LOCMAR points to marginal tables to be placed in MARG */
    locmar[i] = point;
    point += size;
    }

    /* Get N, number of valid configurations */

    i = *ncon + 1;
L160:
    n = i - 1;

    /* See if MARG can hold all marginal tables */

    if (point > *nmar + 1) {
    goto L35;
    }

    /* Obtain marginal tables */

    for (i = 1; i <= n; ++i) {
    for (j = 1; j <= *nvar; ++j) {
        icon[j - 1] = config[j + i * *nvar];
    }
    collap(nvar, &table[1], &marg[1], &locmar[i], ntab, nmar,
           &dim[1], icon);
    }

    /* Perform iterations */

    for (k = 1; k <= *maxit; ++k) {
    /* XMAX is maximum deviation observed between fitted and true
       marginal during a cycle */
    xmax = 0.;
    for (i = 1; i <= n; ++i) {
        for (j = 1; j <= *nvar; ++j) {
        icon[j - 1] = config[j + i * *nvar];
        }
        collap(nvar, &fit[1], &u[1], &c__1, ntab, nu, &dim[1],
           icon);
        adjust(nvar, &fit[1], &u[1], &marg[1], &locmar[i], ntab, nu,
           nmar, &dim[1], icon, &xmax);
    }
    /* Test convergence */
    dev[k] = xmax;
    if (xmax < *maxdev) {
        goto L240;
    }
    }
    if (*maxit > 1) {
    goto L230;
    }
    *nlast = 1;
    return;

    /* No convergence */
L230:
    *ifault = 3;
    *nlast = *maxit;
    return;

    /* Normal termination */
L240:
    *nlast = k;

    return;
}

/* Algorithm AS 51.1 Appl. Statist. (1972), vol. 21, p. 218

   Computes a marginal table from a complete table.
   All parameters are assumed valid without test.

   The larger table is X and the smaller one is Y.
*/

void collap(int *nvar, double *x, double *y, int *locy, int *nx, int
       *ny, int *dim, int *config)
{
    int i, j, k, l, n, locu, *coord, *size;

    size = lvector(*nvar + 1);
    coord = lvector(*nvar);

    /* Parameter adjustments */
    --config;
    --dim;
    --x;
    --y;

    /* Initialize arrays */

    size[0] = 1;
    for (k = 1; k <= *nvar; ++k) {
    l = config[k];
    if (l == 0) {
        goto L20;
    }
    size[k] = size[k - 1] * dim[l];
    }

    /* Find number of variables in configuration */

    k = *nvar + 1;
L20:
    n = k - 1;

    /* Initialize Y.  First cell of marginal table is at Y(LOCY) and
       table has SIZE(K) elements */

    locu = *locy + size[k - 1] - 1;
    for (j = *locy; j <= locu; ++j) {
    y[j] = 0.;
    }

    /* Initialize coordinates */

    for (k = 0; k < *nvar; ++k) {
    coord[k] = 0;
    }

    /* Find locations in tables */
    i = 1;
L60:
    j = *locy;
    for (k = 1; k <= n; ++k) {
    l = config[k];
    j += coord[l - 1] * size[k - 1];
    }
    y[j] += x[i];

    /* Update coordinates */

    ++i;
    for (k = 1; k <= *nvar; ++k) {
    ++coord[k - 1];
    if (coord[k - 1] < dim[k]) {
        goto L60;
    }
    coord[k - 1] = 0;
    }

    return;
}


/* Algorithm AS 51.2 Appl. Statist. (1972), vol. 21, p. 218

   Makes proportional adjustment corresponding to CONFIG.
   All parameters are assumed valid without test.
   */

void adjust(int *nvar, double *x, double *y, double *z, int *locz, int
       *nx, int *ny, int *nz, int *dim, int *config, double *d)
{
    int i, j, k, l, n, *coord, *size;
    double e;

    size = lvector(*nvar + 1);
    coord = lvector(*nvar);

    /* Parameter adjustments */
    --config;
    --dim;
    --x;
    --y;
    --z;

    /* Set size array */

    size[0] = 1;
    for (k = 1; k <= *nvar; ++k) {
    l = config[k];
    if (l == 0) {
        goto L20;
    }
    size[k] = size[k - 1] * dim[l];
    }

    /* Find number of variables in configuration */

    k = *nvar + 1;
L20:
    n = k - 1;

    /* Test size of deviation */

    l = size[k - 1];
    j = 1;
    k = *locz;
    for (i = 1; i <= l; ++i) {
    e = abs(z[k] - y[j]);
    if (e > *d) {
        *d = e;
    }
    ++j;
    ++k;
    }

    /* Initialize coordinates */

    for (k = 0; k < *nvar; ++k) {
    coord[k] = 0;
    }
    i = 1;

    /* Perform adjustment */

L50:
    j = 0;
    for (k = 1; k <= n; ++k) {
    l = config[k];
    j += coord[l - 1] * size[k - 1];
    }
    k = j + *locz;
    ++j;

    /* Note that Y(J) should be non-negative */

    if (y[j] <= 0.) {
    x[i] = 0.;
    }
    if (y[j] > 0.) {
    x[i] = x[i] * z[k] / y[j];
    }

    /* Update coordinates */

    ++i;
    for (k = 1; k <= *nvar; ++k) {
    ++coord[k - 1];
    if (coord[k - 1] < dim[k]) {
        goto L50;
    }
    coord[k - 1] = 0;
    }

    return;
}

/* Auxiliary routine to get rid of limitations on the number of factors
   in the model. 

   Changed to use R_alloc to avoid memory leak if routine was
   interrupted.
*/

static int *lvector(int n) {
    return (int *) R_alloc(n, sizeof(int));
}