pmsco-public/pmsco/loess/loess.i

%module loess
%include "typemaps.i"

%{
#define SWIG_FILE_WITH_INIT
#include <errno.h>
#define EARRLEN 1000
#include "loess.h"

extern void loess(struct loess_struct *lo);

extern void loess_summary(struct loess_struct *lo);

// not implemented
// extern void predict(double *eval, int m, struct loess_struct *lo, struct pred_struct *pre, int se);
%}

%include "numpy.i"

%init %{
import_array();
%}

%apply (double *IN_ARRAY1, int DIM1) {(double *v, int n)};
%apply (int *IN_ARRAY1, int DIM1) {(int *v, int n)};
%apply (double **ARGOUTVIEWM_ARRAY1, int *DIM1) {(double **w, int *n)};
%apply (int **ARGOUTVIEWM_ARRAY1, int *DIM1) {(int **w, int *n)};

%include "loess.h"

extern void loess(struct loess_struct *lo);

extern void loess_summary(struct loess_struct *lo);

// not implemented
// extern void predict(double *eval, int m, struct loess_struct *lo, struct pred_struct *pre, int se);

%exception {
    errno = 0;
    $action
    
    if (errno != 0) {
        switch(errno) {
            case ENOMEM:
                PyErr_Format(PyExc_MemoryError, "memory allocation failed.");
                break;
            case EARRLEN:
                PyErr_Format(PyExc_ValueError, "unexpected array length.");
                break;
            default:
                PyErr_Format(PyExc_Exception, "unknown exception.");
        }
        SWIG_fail;
    }
}

%extend loess_struct {
    //// constructor of a loess_struct
    //
    //   @param n: number of data points.
    //
    //   @param p: number of factors (independent variables). maximum 8.
    
	loess_struct(int n, int p) {
		struct loess_struct *lo;
		lo = (struct loess_struct *) malloc(sizeof(loess_struct));

        int	i, max_kd;
        max_kd = n > 200 ? n : 200;
    
        lo->in.y = (double *) malloc(n * sizeof(double));
        lo->in.x = (double *) malloc(n * p * sizeof(double));
        lo->in.weights = (double *) malloc(n * sizeof(double));
        for(i = 0; i < (n * p); i++)
            lo->in.x[i] = 0.0;
        for(i = 0; i < n; i++) {
            lo->in.y[i] = 0.0;
            lo->in.weights[i] = 1.0;
        }
        lo->in.n = n;
        lo->in.p = p;
        lo->model.span = 0.75;
        lo->model.degree = 2;
        lo->model.normalize = TRUE;
        for(i = 0; i < 8; i++)
            lo->model.parametric[i] = lo->model.drop_square[i] = FALSE;
        lo->model.family = "gaussian";
        lo->control.surface = "interpolate";
        lo->control.statistics = "approximate";
        lo->control.cell = 0.2;
        lo->control.trace_hat = "wait.to.decide";
        lo->control.iterations = 4;
    
        lo->out.fitted_values = (double *) malloc(n * sizeof(double));
        lo->out.fitted_residuals = (double *) malloc(n * sizeof(double));
        lo->out.pseudovalues = (double *) malloc(n * sizeof(double));
        lo->out.diagonal = (double *) malloc(n * sizeof(double));
        lo->out.robust = (double *) malloc(n * sizeof(double));
        lo->out.divisor = (double *) malloc(p * sizeof(double));
    
        lo->kd_tree.parameter = (int *) malloc(7 * sizeof(int));
        lo->kd_tree.a = (int *) malloc(max_kd * sizeof(int));
        lo->kd_tree.xi = (double *) malloc(max_kd * sizeof(double));
        lo->kd_tree.vert = (double *) malloc(p * 2 * sizeof(double));
        lo->kd_tree.vval = (double *) malloc((p + 1) * max_kd * sizeof(double));

		return lo;
	}
    
	~loess_struct() {
        free($self->in.x);
        free($self->in.y);
        free($self->in.weights);
        free($self->out.fitted_values);
        free($self->out.fitted_residuals);
        free($self->out.pseudovalues);
        free($self->out.diagonal);
        free($self->out.robust);
        free($self->out.divisor);
        free($self->kd_tree.parameter);
        free($self->kd_tree.a);
        free($self->kd_tree.xi);
        free($self->kd_tree.vert);
        free($self->kd_tree.vval);
		free($self);
	}
    
    void set_x(double *v, int n) {
        int n_exp = $self->in.n * $self->in.p;
        if (n == n_exp) {
            int i;
            for(i = 0; i < n; i++)
                $self->in.x[i] = v[i];
        } else {
            errno = EARRLEN;
        }
    }

    void set_y(double *v, int n) {
        int n_exp = $self->in.n;
        if (n == n_exp) {
            int i;
            for(i = 0; i < n; i++)
                $self->in.y[i] = v[i];
        } else {
            errno = EARRLEN;
        }
    }

    void set_parametric(int *v, int n) {
        int n_exp = $self->in.p;
        if (n == n_exp) {
            int i;
            for(i = 0; i < n; i++)
                $self->model.parametric[i] = v[i];
        } else {
            errno = EARRLEN;
        }
    }

    void set_drop_square(int *v, int n) {
        int n_exp = $self->in.p;
        if (n == n_exp) {
            int i;
            for(i = 0; i < n; i++)
                $self->model.drop_square[i] = v[i];
        } else {
            errno = EARRLEN;
        }
    }

    void get_x(double **w, int *n) {
        int ni = $self->in.n * $self->in.p;
        
        double *temp;
        temp = (double *)malloc(ni * sizeof(double));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->in.x[i];
        
        *w = temp;
        *n = ni;
    }

    void get_y(double **w, int *n) {
        int ni = $self->in.n;
        
        double *temp;
        temp = (double *)malloc(ni * sizeof(double));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->in.y[i];
        
        *w = temp;
        *n = ni;
    }

    void get_weights(double **w, int *n) {
        int ni = $self->in.n;
        
        double *temp;
        temp = (double *)malloc(ni * sizeof(double));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->in.weights[i];
        
        *w = temp;
        *n = ni;
    }

    void get_fitted_values(double **w, int *n) {
        int ni = $self->in.n;
        
        double *temp;
        temp = (double *)malloc(ni * sizeof(double));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->out.fitted_values[i];
        
        *w = temp;
        *n = ni;
    }

    void get_fitted_residuals(double **w, int *n) {
        int ni = $self->in.n;
        
        double *temp;
        temp = (double *)malloc(ni * sizeof(double));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->out.fitted_residuals[i];
        
        *w = temp;
        *n = ni;
    }

    void get_parametric(int **w, int *n) {
        int ni = $self->in.p;
        
        int *temp;
        temp = (int *)malloc(ni * sizeof(int));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->model.parametric[i];
        
        *w = temp;
        *n = ni;
    }

    void get_drop_square(int **w, int *n) {
        int ni = $self->in.p;
        
        int *temp;
        temp = (int *)malloc(ni * sizeof(int));
        if (temp == NULL)
            errno = ENOMEM;

        int i;
        for(i = 0; i < ni; i++)
            temp[i] = $self->model.drop_square[i];
        
        *w = temp;
        *n = ni;
    }

};