LRlspline.c

Go to the documentation of this file.

/*
 * Copyright    2019    R.K. Owen, Ph.D.
 * License      see lgpl.md (Gnu Lesser General Public License)
 */
#ifdef __cplusplus
extern "C" {
#endif

#include <stdlib.h>
#include <math.h>       /* sqrt, sqrtf, NAN */
#include "libran.h"
#include <stdio.h>

int LR_lspl_new(LR_obj *o, int n) {
        return LR_pcs_new(o,n);
}

int LR_lspl_rm(LR_obj *o) {
        return LR_pcs_rm(o);
}

int LR_lspl_set(LR_obj *o, double x, double p) {
        int ret = LR_pcs_set(o,x,p);
        /* must have at least one good value set */
        if ((!ret) && (p > 0))
                ((LR_pcs *) o->aux)->flags |= LR_AUX_SET;
        return ret;
}

int LR_lspl_norm(LR_obj *o) {
        double zero = 0.0, half = 0.5, one = 1.0, delta = .000001;
        double v = zero;
        LR_pcs *aux = (LR_pcs *) o->aux;

        /* must have at least one good value set else why bother? */
        if (!(aux->flags & LR_AUX_SET))
                return o->errno = LRerr_UnmetPreconditions;

        /* move up the set of values */
        for (int i = aux->nn - 1, i1 = aux->nn - 2; i >= 1; i--, i1--) {
                aux->bdrs[i] = aux->bdrs[i1];
                aux->c[i] = aux->c[i1];
        }
        if (o->d == LR_double) {
                aux->bdrs[0] = o->a.d;
                aux->bdrs[aux->nn] = o->b.d;
                aux->c[0] = o->x.d;
        } else if (o->d == LR_float) {
                aux->bdrs[0] = (double) o->a.f;
                aux->bdrs[aux->nn] = (double) o->b.f;
                aux->c[0] = (double) o->x.f;
        }

        /* integrate over each interval */
        aux->norm = zero;
        for (int i = 0, i1 = 1; i < aux->nn; i++, i1++) {
                aux->norm += 
                aux->sc[i] = half * (aux->c[i1] + aux->c[i])
                        * (aux->bdrs[i1] - aux->bdrs[i]);
        }

        /* collect and rescale the CDF */
        aux->norm = one/aux->norm;
        for (int i = 0; i < aux->nn; i++) {
                v += aux->sc[i];
                aux->sc[i] = v * aux->norm;
        }
        /* shift the CDF set up */
        if (aux->sc[aux->nn-1] < one - delta
        ||  one + delta < aux->sc[aux->nn-1]) {
                /* the last value should be 1.0 */
                return o->errno = LRerr_SuspiciousValues;
        }
        for (int i = aux->nn-1; i > 0; i--) {
                aux->sc[i] = aux->sc[i-1];
        }
        aux->sc[0] = zero;
        aux->sc[aux->nn] = one;
        /* norm success */
        ((LR_pcs *) o->aux)->flags |= LR_AUX_NORM;

        return LRerr_OK;
}

double LRd_lspline_RAN(LR_obj *o) {
        LR_pcs *aux = (LR_pcs *) o->aux;
        double x, y, dy, slope, zero = 0.0, two = 2.0;
        int i = 0;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        x = o->ud(o);
        /* find interval */
        while (x >= aux->sc[i]) i++;

        /* use inverse method to return variate */
        slope = ((aux->c[i] - aux->c[i-1]) * aux->norm) /
                 (aux->bdrs[i] - aux->bdrs[i-1]);

        y = aux->c[i-1] * aux->norm;
        x -= aux->sc[i-1];
        dy = (sqrt(y*y + two*slope*(x)) - y)/slope;

        y = aux->bdrs[i-1] + dy;
/*
        printf(">>>>>% 8.5f % 8.5f % 8.5f\n",
                y, aux->sc[i-1] + x, LRd_CDF(o, y));
*/
        return y;
}

double LRd_lspline_PDF(LR_obj *o, double x) {
        LR_pcs *aux = (LR_pcs *) o->aux;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        if (x <= o->a.d || x >= o->b.d) {
                return 0.0;
        } else {
                /* find interval */
                int i = 0;
                double slope, y;
                while (x > aux->bdrs[i])        i++;
                slope = (aux->c[i] - aux->c[i-1]) /
                        (aux->bdrs[i] - aux->bdrs[i-1]);
                y =     (aux->c[i-1] + (slope * (x - aux->bdrs[i-1])))
                        *aux->norm;
/*
        printf("lspline_PDF>>>>>% 8.5f % 8.5f\n", x, y);
*/
                return y;
        }
}

double LRd_lspline_CDF(LR_obj *o, double x) {
        LR_pcs *aux = (LR_pcs *) o->aux;
        double zero = 0.0, one = 1.0, half = 0.5, y;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        if (x <= o->a.d) {
                return zero;
        } else if (x >= o->b.d) {
                return one;
        } else {
                /* find interval */
                int i = 0;
                double diff, slope;
                while (x > aux->bdrs[i])        i++;
                diff =  (x - aux->bdrs[i-1]);
                slope = (aux->c[i] - aux->c[i-1]) /
                        (aux->bdrs[i] - aux->bdrs[i-1]);
                y =     (aux->c[i-1] + (slope * diff));
                y = aux->sc[i-1] + half * (y + aux->c[i-1])*diff*aux->norm;
/*
                printf("lspline_CDF>>>>>% 8.5f % 8.5f\n", x, y);
*/
                return y;
        }
}

float LRf_lspline_RAN(LR_obj *o) {
        LR_pcs *aux = (LR_pcs *) o->aux;
        float x, y, dy, slope, zero = 0.0, two = 2.0;
        int i = 0;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        x = o->uf(o);
        /* find interval */
        while (x >= aux->sc[i]) i++;

        /* use inverse method to return variate */
        slope = ((aux->c[i] - aux->c[i-1]) * aux->norm) /
                 (aux->bdrs[i] - aux->bdrs[i-1]);

        y = aux->c[i-1] * aux->norm;
        x -= aux->sc[i-1];
        dy = (sqrtf(y*y + two*slope*(x)) - y)/slope;

        y = aux->bdrs[i-1] + dy;
/*
        printf(">>>>>% 8.5f % 8.5f % 8.5f\n",
                y, aux->sc[i-1] + x, LRf_CDF(o, y));
*/
        return y;
}

float LRf_lspline_PDF(LR_obj *o, float x) {
        LR_pcs *aux = (LR_pcs *) o->aux;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        if (x <= o->a.f || x >= o->b.f) {
                return 0.0;
        } else {
                /* find interval */
                int i = 0;
                float slope, y;
                while (x > aux->bdrs[i])        i++;
                slope = (aux->c[i] - aux->c[i-1]) /
                        (aux->bdrs[i] - aux->bdrs[i-1]);
                y =     (aux->c[i-1] + (slope * (x - aux->bdrs[i-1])))
                        *aux->norm;
/*
        printf("lspline_PDF>>>>>% 8.5f % 8.5f\n", x, y);
*/
                return y;
        }
}

float LRf_lspline_CDF(LR_obj *o, float x) {
        LR_pcs *aux = (LR_pcs *) o->aux;
        float zero = 0.0, one = 1.0, half = 0.5, y;

        /* must have successfully normalized */
        if (!(aux->flags & LR_AUX_NORM)) {
                o->errno = LRerr_NoAuxNormalizeDone;
                return NAN;
        }

        if (x <= o->a.f) {
                return zero;
        } else if (x >= o->b.f) {
                return one;
        } else {
                /* find interval */
                int i = 0;
                float diff, slope;
                while (x > aux->bdrs[i])        i++;
                diff =  (x - aux->bdrs[i-1]);
                slope = (aux->c[i] - aux->c[i-1]) /
                        (aux->bdrs[i] - aux->bdrs[i-1]);
                y =     (aux->c[i-1] + (slope * diff));
                y = aux->sc[i-1] + half * (y + aux->c[i-1])*diff*aux->norm;
/*
                printf("lspline_CDF>>>>>% 8.5f % 8.5f\n", x, y);
*/
                return y;
        }
}

#ifdef __cplusplus
}
#endif