N_arrays_calc.c

Go to the documentation of this file.

/*****************************************************************************
*
* MODULE:       Grass PDE Numerical Library
* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
*               soerengebbert <at> gmx <dot> de
*               
* PURPOSE:      Higher level array managment functions 
*               part of the gpde library
*
* COPYRIGHT:    (C) 2000 by the GRASS Development Team
*
*               This program is free software under the GNU General Public
*               License (>=v2). Read the file COPYING that comes with GRASS
*               for details.
*
*****************************************************************************/

#include "grass/N_pde.h"
#include "grass/glocale.h"
#include <math.h>

/* ******************** 2D ARRAY FUNCTIONS *********************** */

void N_copy_array_2d(N_array_2d * source, N_array_2d * target)
{
    int i;
    int null = 0;

#pragma omp single
    {
        if (source->cols_intern != target->cols_intern)
            G_fatal_error
                ("N_copy_array_2d: the arrays are not of equal size");

        if (source->rows_intern != target->rows_intern)
            G_fatal_error
                ("N_copy_array_2d: the arrays are not of equal size");

        G_debug(3,
                "N_copy_array_2d: copy source array to target array size %i",
                source->cols_intern * source->rows_intern);
    }

#pragma omp for
    for (i = 0; i < source->cols_intern * source->rows_intern; i++) {
        null = 0;
        if (source->type == CELL_TYPE) {
            if (G_is_c_null_value((void *)&source->cell_array[i]))
                null = 1;

            if (target->type == CELL_TYPE) {
                target->cell_array[i] = source->cell_array[i];
            }
            if (target->type == FCELL_TYPE) {
                if (null)
                    G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
                else
                    target->fcell_array[i] = (FCELL) source->cell_array[i];
            }
            if (target->type == DCELL_TYPE) {
                if (null)
                    G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
                else
                    target->dcell_array[i] = (DCELL) source->cell_array[i];
            }

        }
        if (source->type == FCELL_TYPE) {
            if (G_is_f_null_value((void *)&source->fcell_array[i]))
                null = 1;

            if (target->type == CELL_TYPE) {
                if (null)
                    G_set_c_null_value((void *)&(target->cell_array[i]), 1);
                else
                    target->cell_array[i] = (CELL) source->fcell_array[i];
            }
            if (target->type == FCELL_TYPE) {
                target->fcell_array[i] = source->fcell_array[i];
            }
            if (target->type == DCELL_TYPE) {
                if (null)
                    G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
                else
                    target->dcell_array[i] = (DCELL) source->fcell_array[i];
            }
        }
        if (source->type == DCELL_TYPE) {
            if (G_is_d_null_value((void *)&source->dcell_array[i]))
                null = 1;

            if (target->type == CELL_TYPE) {
                if (null)
                    G_set_c_null_value((void *)&(target->cell_array[i]), 1);
                else
                    target->cell_array[i] = (CELL) source->dcell_array[i];
            }
            if (target->type == FCELL_TYPE) {
                if (null)
                    G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
                else
                    target->fcell_array[i] = (FCELL) source->dcell_array[i];
            }
            if (target->type == DCELL_TYPE) {
                target->dcell_array[i] = source->dcell_array[i];
            }
        }
    }

    return;
}

double N_norm_array_2d(N_array_2d * a, N_array_2d * b, int type)
{
    int i = 0;
    double norm = 0.0, tmp = 0.0;
    double v1 = 0.0, v2 = 0.0;

    if (a->cols_intern != b->cols_intern)
        G_fatal_error("N_norm_array_2d: the arrays are not of equal size");

    if (a->rows_intern != b->rows_intern)
        G_fatal_error("N_norm_array_2d: the arrays are not of equal size");

    G_debug(3, "N_norm_array_2d: norm of a and b size %i",
            a->cols_intern * a->rows_intern);

    for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
        v1 = 0.0;
        v2 = 0.0;

        if (a->type == CELL_TYPE) {
            if (!G_is_f_null_value((void *)&(a->cell_array[i])))
                v1 = (double)a->cell_array[i];
        }
        if (a->type == FCELL_TYPE) {
            if (!G_is_f_null_value((void *)&(a->fcell_array[i])))
                v1 = (double)a->fcell_array[i];
        }
        if (a->type == DCELL_TYPE) {
            if (!G_is_f_null_value((void *)&(a->dcell_array[i])))
                v1 = (double)a->dcell_array[i];
        }
        if (b->type == CELL_TYPE) {
            if (!G_is_f_null_value((void *)&(b->cell_array[i])))
                v2 = (double)b->cell_array[i];
        }
        if (b->type == FCELL_TYPE) {
            if (!G_is_f_null_value((void *)&(b->fcell_array[i])))
                v2 = (double)b->fcell_array[i];
        }
        if (b->type == DCELL_TYPE) {
            if (!G_is_f_null_value((void *)&(b->dcell_array[i])))
                v2 = (double)b->dcell_array[i];
        }

        if (type == N_MAXIMUM_NORM) {
            tmp = fabs(v2 - v1);
            if ((tmp > norm))
                norm = tmp;
        }
        if (type == N_EUKLID_NORM) {
            norm += fabs(v2 - v1);
        }
    }

    return norm;
}

void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
                           double *sum, int *nonull, int withoffset)
{
    int i, j;
    double val;

    *sum = 0.0;
    *nonull = 0;

    if (withoffset == 1) {

        *min =
            (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
        *max =
            (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);

        for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
            for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
                if (!N_is_array_2d_value_null(a, i, j)) {
                    val = (double)N_get_array_2d_d_value(a, i, j);
                    if (*min > val)
                        *min = val;
                    if (*max < val)
                        *max = val;
                    *sum += val;
                    (*nonull)++;
                }
            }
        }
    }
    else {

        *min = (double)N_get_array_2d_d_value(a, 0, 0);
        *max = (double)N_get_array_2d_d_value(a, 0, 0);


        for (j = 0; j < a->rows; j++) {
            for (i = 0; i < a->cols; i++) {
                if (!N_is_array_2d_value_null(a, i, j)) {
                    val = (double)N_get_array_2d_d_value(a, i, j);
                    if (*min > val)
                        *min = val;
                    if (*max < val)
                        *max = val;
                    *sum += val;
                    (*nonull)++;
                }
            }
        }
    }

    G_debug(3,
            "N_calc_array_2d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
            *min, *max, *sum, *nonull);
    return;
}


N_array_2d *N_math_array_2d(N_array_2d * a, N_array_2d * b,
                            N_array_2d * result, int type)
{
    N_array_2d *c;
    int i, j, setnull = 0;
    double va = 0.0, vb = 0.0, vc = 0.0;        /*variables used for calculation */

    /*Set the pointer */
    c = result;

#pragma omp single
    {
        /*Check the array sizes */
        if (a->cols_intern != b->cols_intern)
            G_fatal_error
                ("N_math_array_2d: the arrays are not of equal size");
        if (a->rows_intern != b->rows_intern)
            G_fatal_error
                ("N_math_array_2d: the arrays are not of equal size");
        if (a->offset != b->offset)
            G_fatal_error
                ("N_math_array_2d: the arrays have different offsets");

        G_debug(3, "N_math_array_2d: mathematical calculations, size: %i",
                a->cols_intern * a->rows_intern);

        /*if the result array is null, allocate a new one, use the 
         * largest data type of the input arrays*/
        if (c == NULL) {
            if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
                c = N_alloc_array_2d(a->cols, a->rows, a->offset, DCELL_TYPE);
                G_debug(3,
                        "N_math_array_2d: array of type DCELL_TYPE created");
            }
            else if (a->type == FCELL_TYPE || b->type == FCELL_TYPE) {
                c = N_alloc_array_2d(a->cols, a->rows, a->offset, FCELL_TYPE);
                G_debug(3,
                        "N_math_array_2d: array of type FCELL_TYPE created");
            }
            else {
                c = N_alloc_array_2d(a->cols, a->rows, a->offset, CELL_TYPE);
                G_debug(3,
                        "N_math_array_2d: array of type CELL_TYPE created");
            }
        }
        else {
            /*Check the array sizes */
            if (a->cols_intern != c->cols_intern)
                G_fatal_error
                    ("N_math_array_2d: the arrays are not of equal size");
            if (a->rows_intern != c->rows_intern)
                G_fatal_error
                    ("N_math_array_2d: the arrays are not of equal size");
            if (a->offset != c->offset)
                G_fatal_error
                    ("N_math_array_2d: the arrays have different offsets");
        }
    }

#pragma omp for private(va, vb, vc, setnull)
    for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
        for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
            if (!N_is_array_2d_value_null(a, i, j) &&
                !N_is_array_2d_value_null(b, i, j)) {
                /*we always calulate internally with double values */
                va = (double)N_get_array_2d_d_value(a, i, j);
                vb = (double)N_get_array_2d_d_value(b, i, j);
                vc = 0;
                setnull = 0;

                switch (type) {
                case N_ARRAY_SUM:
                    vc = va + vb;
                    break;
                case N_ARRAY_DIF:
                    vc = va - vb;
                    break;
                case N_ARRAY_MUL:
                    vc = va * vb;
                    break;
                case N_ARRAY_DIV:
                    if (vb != 0)
                        vc = va / vb;
                    else
                        setnull = 1;
                    break;
                }

                if (c->type == CELL_TYPE) {
                    if (setnull)
                        N_put_array_2d_value_null(c, i, j);
                    else
                        N_put_array_2d_c_value(c, i, j, (CELL) vc);
                }
                if (c->type == FCELL_TYPE) {
                    if (setnull)
                        N_put_array_2d_value_null(c, i, j);
                    else
                        N_put_array_2d_f_value(c, i, j, (FCELL) vc);
                }
                if (c->type == DCELL_TYPE) {
                    if (setnull)
                        N_put_array_2d_value_null(c, i, j);
                    else
                        N_put_array_2d_d_value(c, i, j, (DCELL) vc);
                }

            }
            else {
                N_put_array_2d_value_null(c, i, j);
            }
        }
    }

    return c;
}

int N_convert_array_2d_null_to_zero(N_array_2d * a)
{
    int i = 0, count = 0;

    G_debug(3, "N_convert_array_2d_null_to_zero: convert array of size %i",
            a->cols_intern * a->rows_intern);

    if (a->type == CELL_TYPE)
        for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
            if (G_is_c_null_value((void *)&(a->cell_array[i]))) {
                a->cell_array[i] = 0;
                count++;
            }
        }

    if (a->type == FCELL_TYPE)
        for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
            if (G_is_f_null_value((void *)&(a->fcell_array[i]))) {
                a->fcell_array[i] = 0.0;
                count++;
            }
        }


    if (a->type == DCELL_TYPE)
        for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
            if (G_is_d_null_value((void *)&(a->dcell_array[i]))) {
                a->dcell_array[i] = 0.0;
                count++;
            }
        }


    if (a->type == CELL_TYPE)
        G_debug(2,
                "N_convert_array_2d_null_to_zero: %i values of type CELL_TYPE are converted",
                count);
    if (a->type == FCELL_TYPE)
        G_debug(2,
                "N_convert_array_2d_null_to_zero: %i valuess of type FCELL_TYPE are converted",
                count);
    if (a->type == DCELL_TYPE)
        G_debug(2,
                "N_convert_array_2d_null_to_zero: %i valuess of type DCELL_TYPE are converted",
                count);

    return count;
}

/* ******************** 3D ARRAY FUNCTIONS *********************** */

void N_copy_array_3d(N_array_3d * source, N_array_3d * target)
{
    int i;
    int null;

    if (source->cols_intern != target->cols_intern)
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");

    if (source->rows_intern != target->rows_intern)
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");

    if (source->depths_intern != target->depths_intern)
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");


    G_debug(3, "N_copy_array_3d: copy source array to target array size %i",
            source->cols_intern * source->rows_intern *
            source->depths_intern);

    for (i = 0;
         i <
         source->cols_intern * source->rows_intern * source->depths_intern;
         i++) {
        null = 0;
        if (source->type == FCELL_TYPE) {
            if (G3d_isNullValueNum
                ((void *)&(source->fcell_array[i]), FCELL_TYPE))
                null = 1;

            if (target->type == FCELL_TYPE) {
                target->fcell_array[i] = source->fcell_array[i];
            }
            if (target->type == DCELL_TYPE) {
                if (null)
                    G3d_setNullValue((void *)&(target->dcell_array[i]), 1,
                                     DCELL_TYPE);
                else
                    target->dcell_array[i] = (double)source->fcell_array[i];
            }

        }
        if (source->type == DCELL_TYPE) {
            if (G3d_isNullValueNum
                ((void *)&(source->dcell_array[i]), DCELL_TYPE))
                null = 1;

            if (target->type == FCELL_TYPE) {
                if (null)
                    G3d_setNullValue((void *)&(target->fcell_array[i]), 1,
                                     FCELL_TYPE);
                else
                    target->fcell_array[i] = (float)source->dcell_array[i];
            }
            if (target->type == DCELL_TYPE) {
                target->dcell_array[i] = source->dcell_array[i];
            }
        }
    }

    return;
}


double N_norm_array_3d(N_array_3d * a, N_array_3d * b, int type)
{
    int i = 0;
    double norm = 0.0, tmp = 0.0;
    double v1 = 0.0, v2 = 0.0;

    if (a->cols_intern != b->cols_intern)
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");

    if (a->rows_intern != b->rows_intern)
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");

    if (a->depths_intern != b->depths_intern)
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");

    G_debug(3, "N_norm_array_3d: norm of a and b size %i",
            a->cols_intern * a->rows_intern * a->depths_intern);

    for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern; i++) {
        v1 = 0.0;
        v2 = 0.0;

        if (a->type == FCELL_TYPE) {
            if (!G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE))
                v1 = (double)a->fcell_array[i];
        }
        if (a->type == DCELL_TYPE) {
            if (!G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE))
                v1 = (double)a->dcell_array[i];
        }
        if (b->type == FCELL_TYPE) {
            if (!G3d_isNullValueNum((void *)&(b->fcell_array[i]), FCELL_TYPE))
                v2 = (double)b->fcell_array[i];
        }
        if (b->type == DCELL_TYPE) {
            if (!G3d_isNullValueNum((void *)&(b->dcell_array[i]), DCELL_TYPE))
                v2 = (double)b->dcell_array[i];
        }

        if (type == N_MAXIMUM_NORM) {
            tmp = fabs(v2 - v1);
            if ((tmp > norm))
                norm = tmp;
        }
        if (type == N_EUKLID_NORM) {
            norm += fabs(v2 - v1);
        }
    }

    return norm;
}

void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max,
                           double *sum, int *nonull, int withoffset)
{
    int i, j, k;
    double val;

    *sum = 0.0;
    *nonull = 0;

    if (withoffset == 1) {

        *min =
            (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
                                           0 - a->offset);
        *max =
            (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
                                           0 - a->offset);

        for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
            for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
                for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
                    if (!N_is_array_3d_value_null(a, i, j, k)) {
                        val = (double)N_get_array_3d_d_value(a, i, j, k);
                        if (*min > val)
                            *min = val;
                        if (*max < val)
                            *max = val;
                        *sum += val;
                        (*nonull)++;
                    }
                }
            }
        }
    }
    else {

        *min = (double)N_get_array_3d_d_value(a, 0, 0, 0);
        *max = (double)N_get_array_3d_d_value(a, 0, 0, 0);

        for (k = 0; k < a->depths; k++) {
            for (j = 0; j < a->rows; j++) {
                for (i = 0; i < a->cols; i++) {
                    if (!N_is_array_3d_value_null(a, i, j, k)) {
                        val = (double)N_get_array_3d_d_value(a, i, j, k);
                        if (*min > val)
                            *min = val;
                        if (*max < val)
                            *max = val;
                        *sum += val;
                        (*nonull)++;
                    }
                }
            }
        }
    }

    G_debug(3,
            "N_calc_array_3d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
            *min, *max, *sum, *nonull);

    return;
}

N_array_3d *N_math_array_3d(N_array_3d * a, N_array_3d * b,
                            N_array_3d * result, int type)
{
    N_array_3d *c;
    int i, j, k, setnull = 0;
    double va = 0.0, vb = 0.0, vc = 0.0;        /*variables used for calculation */

    /*Set the pointer */
    c = result;

    /*Check the array sizes */
    if (a->cols_intern != b->cols_intern)
        G_fatal_error("N_math_array_3d: the arrays are not of equal size");
    if (a->rows_intern != b->rows_intern)
        G_fatal_error("N_math_array_3d: the arrays are not of equal size");
    if (a->depths_intern != b->depths_intern)
        G_fatal_error("N_math_array_3d: the arrays are not of equal size");
    if (a->offset != b->offset)
        G_fatal_error("N_math_array_3d: the arrays have different offsets");

    G_debug(3, "N_math_array_3d: mathematical calculations, size: %i",
            a->cols_intern * a->rows_intern * a->depths_intern);

    /*if the result array is null, allocate a new one, use the 
     * largest data type of the input arrays*/
    if (c == NULL) {
        if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
            c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
                                 DCELL_TYPE);
            G_debug(3, "N_math_array_3d: array of type DCELL_TYPE created");
        }
        else {
            c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
                                 FCELL_TYPE);
            G_debug(3, "N_math_array_3d: array of type FCELL_TYPE created");
        }
    }
    else {
        /*Check the array sizes */
        if (a->cols_intern != c->cols_intern)
            G_fatal_error
                ("N_math_array_3d: the arrays are not of equal size");
        if (a->rows_intern != c->rows_intern)
            G_fatal_error
                ("N_math_array_3d: the arrays are not of equal size");
        if (a->depths_intern != c->depths_intern)
            G_fatal_error
                ("N_math_array_3d: the arrays are not of equal size");
        if (a->offset != c->offset)
            G_fatal_error
                ("N_math_array_3d: the arrays have different offsets");
    }

    for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
        for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
            for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
                if (!N_is_array_3d_value_null(a, i, j, k) &&
                    !N_is_array_3d_value_null(a, i, j, k)) {
                    /*we always calulate internally with double values */
                    va = (double)N_get_array_3d_d_value(a, i, j, k);
                    vb = (double)N_get_array_3d_d_value(b, i, j, k);
                    vc = 0;
                    setnull = 0;

                    switch (type) {
                    case N_ARRAY_SUM:
                        vc = va + vb;
                        break;
                    case N_ARRAY_DIF:
                        vc = va - vb;
                        break;
                    case N_ARRAY_MUL:
                        vc = va * vb;
                        break;
                    case N_ARRAY_DIV:
                        if (vb != 0)
                            vc = va / vb;
                        else
                            setnull = 1;
                        break;
                    }

                    if (c->type == FCELL_TYPE) {
                        if (setnull)
                            N_put_array_3d_value_null(c, i, j, k);
                        else
                            N_put_array_3d_f_value(c, i, j, k, (float)vc);
                    }
                    if (c->type == DCELL_TYPE) {
                        if (setnull)
                            N_put_array_3d_value_null(c, i, j, k);
                        else
                            N_put_array_3d_d_value(c, i, j, k, vc);
                    }
                }
                else {
                    N_put_array_3d_value_null(c, i, j, k);
                }
            }
        }
    }

    return c;
}

int N_convert_array_3d_null_to_zero(N_array_3d * a)
{
    int i = 0, count = 0;

    G_debug(3, "N_convert_array_3d_null_to_zero: convert array of size %i",
            a->cols_intern * a->rows_intern * a->depths_intern);

    if (a->type == FCELL_TYPE)
        for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
             i++) {
            if (G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE)) {
                a->fcell_array[i] = 0.0;
                count++;
            }
        }

    if (a->type == DCELL_TYPE)
        for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
             i++) {
            if (G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE)) {
                a->dcell_array[i] = 0.0;
                count++;
            }
        }


    if (a->type == FCELL_TYPE)
        G_debug(3,
                "N_convert_array_3d_null_to_zero: %i values of type FCELL_TYPE are converted",
                count);

    if (a->type == DCELL_TYPE)
        G_debug(3,
                "N_convert_array_3d_null_to_zero: %i values of type DCELL_TYPE are converted",
                count);

    return count;
}