johnpye/datareader/dr.c

/*  ASCEND modelling environment
    Copyright (C) 2006 Carnegie Mellon University

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*/

#include <errno.h>

#include "dr.h"
#include "tmy.h"
#include "acdb.h"

#include <utilities/config.h>
#include <general/ospath.h>
#include <utilities/ascMalloc.h>
#include <utilities/error.h>
#include <utilities/ascPanic.h>
#include <utilities/ascEnvVar.h>

/*------------------------------------------------------------------------------
  FORWARD DECLARATIONS
*/

/*
    declare the possible formats that we will accept ('format' child in the
    DATA instance of the external relation
*/

#define FMTS(D,X) D(TMY2) X D(ACDB) X D(CSV) X D(TDV)

#define ENUM(F_) DATAREADER_FORMAT_##F_
#define COMMA ,
typedef enum{
    FMTS(ENUM,COMMA),
    DATAREADER_INVALID_FORMAT
} DataReaderFormat;
#undef ENUM
#undef COMMA


/*------------------------------------------------------------------------------
  API FUNCTIONS
*/

/**
    Create a data reader object, with the filename specified. The filename
    will be searched for in a specified path, eg ASCENDLIBRARY.
*/
DataReader *datareader_new(const char *fn){
    DataReader *d;

    d = ASC_NEW(DataReader);
    d->fn = fn;
    d->fp = NULL;
    d->f = NULL;
    d->noutputs = 0;

    d->datafn=NULL;
    d->headerfn=NULL;
    d->eoffn=NULL;

    CONSOLE_DEBUG("Datareader created...");
    return d;
}

/**
    Set data file format
    @return 0 on success
*/
int datareader_set_format(DataReader *d, const char *format){

#define STR(N) #N
#define COMMA ,
    const char *fmts[]={ FMTS(STR,COMMA) };
#undef STR
#undef COMMA

    int i;

    CONSOLE_DEBUG("FORMAT '%s'",format);

    DataReaderFormat found = DATAREADER_INVALID_FORMAT;
    for(i=0; i < DATAREADER_INVALID_FORMAT; ++i){
        if(strcmp(format,fmts[i])==0){
            found = (DataReaderFormat)i;
            break;
        }
    }

    CONSOLE_DEBUG("FOUND DATA FORMAT %d",found);

    switch(found){
        case DATAREADER_FORMAT_TMY2:
            d->headerfn=&datareader_tmy2_header;
            d->datafn=&datareader_tmy2_data;
            d->eoffn=&datareader_tmy2_eof;
            d->indepfn=&datareader_tmy2_time;
            d->valfn=&datareader_tmy2_vals;
            break;
        case DATAREADER_FORMAT_ACDB:
            d->headerfn=&datareader_acdb_header;
            d->datafn=&datareader_acdb_data;
            d->eoffn=&datareader_acdb_eof;
            d->indepfn=&datareader_acdb_time;
            d->valfn=&datareader_acdb_vals;
            break;
        case DATAREADER_FORMAT_TDV:
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"Tab delimited values (TDV) format not yet implemenented.");
            return 1;
        case DATAREADER_FORMAT_CSV:
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"Comma-separated values (TDV) format not yet implemenented.");
            return 1;
        default:
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unknown file format '%s' specified",format);
            return 1;
    }

    return 0;
}

typedef struct DataFileSearchData_struct{
    struct FilePath *fp; /**< the relative path we're searching for */
    ospath_stat_t buf; /**< saves memory allocation in the 'test' fn */
    int error;           /**< error code (in case stat or fopen failed) */
    struct FilePath *fp_found; /**< the full path we found */
} DataFileSearchData;

FilePathTestFn datareader_searchpath_test;

/**
    @return 1 on success
*/
int datareader_searchpath_test(struct FilePath *path,void *searchdata){
    struct FilePath *fp1;
    DataFileSearchData *sd;

    sd = (DataFileSearchData *)searchdata;
    assert(sd!=NULL);
    assert(sd->fp!=NULL);

    fp1 = ospath_concat(path,sd->fp);
    if(fp1==NULL){
        CONSOLE_DEBUG("Couldn't concatenate path");
        return 0;
    }

    if(ospath_stat(fp1,&sd->buf)){
        sd->error = errno;
        ospath_free(fp1);
        return 0;
    }

    sd->fp_found = fp1;
    return 1;
};


/**
    Initialise the datareader: open the file, check the number of columns, etc.
    @return 0 on success

    @TODO search for the file in the ASCENDLIBRARY if not found immediately
*/
int datareader_init(DataReader *d){
    ospath_stat_t s;
    char *tmp;
    struct FilePath **sp1, *fp2;
    DataFileSearchData sd;

    d->fp = ospath_new(d->fn);
    if(d->fp==NULL){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid filepath");
        return 1;
    }

    if(ospath_stat(d->fp,&s)){
        if(errno==ENOENT){
            /* file doesn't exist: check the search path instead */
            tmp = Asc_GetEnv(ASC_ENV_LIBRARY);
            if(tmp==NULL){
                ERROR_REPORTER_HERE(ASC_PROG_ERROR,"No paths to search (is env var '%s' set?)",ASC_ENV_LIBRARY);
                return 1;
            }

            sp1 = ospath_searchpath_new(tmp);
            if(sp1==NULL){
                ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Unable to process %s value '%s'",ASC_ENV_LIBRARY,tmp);
                /* memory error */
                ascfree(tmp);
                return -3;
            }
            ascfree(tmp);

            sd.fp = d->fp;

            fp2 = ospath_searchpath_iterate(sp1, &datareader_searchpath_test, &sd);

            if(fp2==NULL){
                ERROR_REPORTER_HERE(ASC_USER_ERROR,"File '%s' not found in search path (error %d)",d->fn,sd.error);
                ospath_searchpath_free(sp1);
                return -1;
            }

            ospath_searchpath_free(sp1);

            /* replace our relative path with an absolute one */
            ospath_free(d->fp);
            d->fp = sd.fp_found;

        }else{
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"The file '%s' cannot be accessed.\n"
                "Check the file privileges, or try specifying an absolute path.",d->fn
            );
            return 1;
        }
    }

    CONSOLE_DEBUG("About to open the data file");

    d->f = ospath_fopen(d->fp,"r");
    if(d->f == NULL){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unable to open file '%s' for reading.",d->fn);
        return 1;
    }

    CONSOLE_DEBUG("Data file open ok");

    asc_assert(d->headerfn);
    asc_assert(d->eoffn);
    asc_assert(d->datafn);

    if((*d->headerfn)(d)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error processing file header in '%s'",d->fn);
        fclose(d->f);
        return 1;
    }

    while(! (*d->eoffn)(d)){
        if((*d->datafn)(d)){
            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error reading file data in '%s'",d->fn);
            fclose(d->f);
            return 1;
        }
    }
    fclose(d->f);

    d->i = 0; /* set current position to zero */

    return 0;
}

/**
    Shut down the data reader and deallocate any memory owned by it, then
    free the memory at d.
*/
int datareader_delete(DataReader *d){
    if(d->fp){
        ospath_free(d->fp);
        d->fp = NULL;
    }
    if(d->f){
        fclose(d->f);
        d->f = NULL;
    }
    ASC_FREE(d);
    return 0;
}

/**
    Return the number of inputs (independent variables) supplied in the
    DataReader's current file. Can only be 1 at this stage.
*/
int datareader_num_inputs(const DataReader *d){
    return d->ninputs;
}

/**
    Return the number of outputs (dependent variables) found in the DataReader's
    current file. Should be one or more.
*/
int datareader_num_outputs(const DataReader *d){
    return d->noutputs;
}


int datareader_locate(DataReader *d, double t, double *t1, double *t2){
    (*d->indepfn)(d,t1);
    if(*t1 > t && d->i > 0){
        /* start of current interval is too late */
        do{
            --(d->i);
            (*d->indepfn)(d,t1);
            /*CONSOLE_DEBUG("STEPPING BACK (d->i = %d currently), t1=%lf",d->i, *t1); */
        }while(*t1 > t && d->i > 0);
    }
    /* now either d->i==0 or t1 < t*/
    /* CONSOLE_DEBUG("d->i==%d, t1=%lf",d->i,*t1); */
    ++d->i;
    (*d->indepfn)(d,t2);
    if(*t2 <= t){
        /* end of current interface is too early */
        do{
            /* CONSOLE_DEBUG("STEPPING FORWARD (d->i = %d currently), t1=%lf, t2=%lf",d->i,*t1,*t2); */
            (*d->indepfn)(d,t1);
            ++(d->i);
            (*d->indepfn)(d,t2);
        }while(*t2 < t && d->i < d->ndata);
    }
    CONSOLE_DEBUG("d->i==%d, t1[0] = %lf, t2[0] = %lf",d->i,t1[0],t2[0]);

    if(d->i == d->ndata || d->i == 0){
        return 1;
    }

    /* CONSOLE_DEBUG("INTERVAL OK"); */
    return 0;
}

/**
    Return an interpolated set of output values for the given input values.
    This should be computed such that the output values are smooth in their
    first derivatives.

    The required memory for the inputs and outputs must be allocated by the
    caller, and indicated by the pointers 'inputs' and 'outputs'.

    @see datareader_deriv
    @TODO implement this
*/
int datareader_func(DataReader *d, double *inputs, double *outputs){
    int i;
    double t1[1], t2[1];
    double v1[5], v2[5]; /** @TODO this is a fixed size... not good */

    double t,g,dt;
    t = inputs[0];

    CONSOLE_DEBUG("EVALUATING AT t = %lf",inputs[0]);

    asc_assert(d->indepfn);

    if(datareader_locate(d,t,t1,t2)){
        CONSOLE_DEBUG("LOCATION ERROR");
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Time value t=%f is out of range",t);
        return 1;
    }
    CONSOLE_DEBUG("LOCATED AT t1 = %lf, t2 = %lf",t1[0], t2[0]);

    (*d->valfn)(d,v2);
    --d->i;
    (*d->valfn)(d,v1);

    CONSOLE_DEBUG("LOCATED OK, d->i = %d, t1 = %lf, t2 = %lf, v1=%lf, v2=%lf", d->i, t1[0], t2[0],v1[0],v2[0]);

    for(i=0;i<d->noutputs;++i){
        dt = t2[0] - t1[0];
        g = (v2[i]-v1[i])/dt;
        outputs[i]=v1[i]+g*(t-(*t1));
        CONSOLE_DEBUG("[%d]: DT = %lf, START = %lf, GRADIENT = %lf, VALUE=%lf",i,dt, v1[i],g,outputs[i]);
    }

    return 0;
}

/**
    Return an interpolated set of output derivatives for the given input
    values. These should be smooth.
    @see datareader_func
    @TODO implement this
*/
int datareader_deriv(DataReader *d, double *inputs, double *jacobian){
    int i;
    double t1[1], t2[1];
    double v1[2], v2[2];

    double t;
    t = inputs[0];

    /* CONSOLE_DEBUG("EVALUATING"); */

    asc_assert(d->indepfn);

    if(datareader_locate(d,t,t1,t2)){
        CONSOLE_DEBUG("LOCATION ERROR");
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Time value t=%f is out of range",t);
        return 1;
    }

    (*d->valfn)(d,v2);
    --d->i;
    (*d->valfn)(d,v1);

    for(i=0;i<d->noutputs;++i){
        jacobian[i]=(v2[i]-v1[i])/(*t2-*t1);
    }

    return 0;
}


1	johnpye	811	/* ASCEND modelling environment
2			Copyright (C) 2006 Carnegie Mellon University
3	johnpye	801
4	johnpye	811	This program is free software; you can redistribute it and/or modify
5			it under the terms of the GNU General Public License as published by
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8	johnpye	801
9	johnpye	811	This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13	johnpye	801
14	johnpye	811	You should have received a copy of the GNU General Public License
15			along with this program; if not, write to the Free Software
16			Foundation, Inc., 59 Temple Place - Suite 330,
17			Boston, MA 02111-1307, USA.
18	johnpye	801	*/
19
20	johnpye	812	#include <errno.h>
21
22	johnpye	811	#include "dr.h"
23			#include "tmy.h"
24	jpye	1773	#include "acdb.h"
25	johnpye	801
26	johnpye	811	#include <utilities/config.h>
27			#include <general/ospath.h>
28			#include <utilities/ascMalloc.h>
29			#include <utilities/error.h>
30			#include <utilities/ascPanic.h>
31	johnpye	819	#include <utilities/ascEnvVar.h>
32	johnpye	811
33	johnpye	809	/*------------------------------------------------------------------------------
34			FORWARD DECLARATIONS
35			*/
36
37	johnpye	812	/*
38			declare the possible formats that we will accept ('format' child in the
39			DATA instance of the external relation
40			*/
41	johnpye	809
42	jpye	1773	#define FMTS(D,X) D(TMY2) X D(ACDB) X D(CSV) X D(TDV)
43	johnpye	812
44			#define ENUM(F_) DATAREADER_FORMAT_##F_
45			#define COMMA ,
46			typedef enum{
47			FMTS(ENUM,COMMA),
48			DATAREADER_INVALID_FORMAT
49			} DataReaderFormat;
50			#undef ENUM
51			#undef COMMA
52
53
54	johnpye	809	/*------------------------------------------------------------------------------
55			API FUNCTIONS
56			*/
57
58	johnpye	801	/**
59			Create a data reader object, with the filename specified. The filename
60			will be searched for in a specified path, eg ASCENDLIBRARY.
61	johnpye	824	*/
62	johnpye	811	DataReader datareader_new(const char fn){
63	johnpye	801	DataReader *d;
64	johnpye	824
65	johnpye	801	d = ASC_NEW(DataReader);
66			d->fn = fn;
67			d->fp = NULL;
68	johnpye	809	d->f = NULL;
69	johnpye	801	d->noutputs = 0;
70
71	johnpye	812	d->datafn=NULL;
72			d->headerfn=NULL;
73			d->eoffn=NULL;
74
75			CONSOLE_DEBUG("Datareader created...");
76	johnpye	801	return d;
77			}
78
79			/**
80	johnpye	809	Set data file format
81			@return 0 on success
82			*/
83	johnpye	812	int datareader_set_format(DataReader d, const char format){
84
85			#define STR(N) #N
86			#define COMMA ,
87			const char *fmts[]={ FMTS(STR,COMMA) };
88			#undef STR
89			#undef COMMA
90
91			int i;
92
93			CONSOLE_DEBUG("FORMAT '%s'",format);
94
95			DataReaderFormat found = DATAREADER_INVALID_FORMAT;
96			for(i=0; i < DATAREADER_INVALID_FORMAT; ++i){
97			if(strcmp(format,fmts[i])==0){
98			found = (DataReaderFormat)i;
99			break;
100			}
101			}
102
103			CONSOLE_DEBUG("FOUND DATA FORMAT %d",found);
104
105			switch(found){
106	johnpye	809	case DATAREADER_FORMAT_TMY2:
107			d->headerfn=&datareader_tmy2_header;
108	johnpye	811	d->datafn=&datareader_tmy2_data;
109	johnpye	812	d->eoffn=&datareader_tmy2_eof;
110			d->indepfn=&datareader_tmy2_time;
111			d->valfn=&datareader_tmy2_vals;
112	johnpye	809	break;
113	jpye	1773	case DATAREADER_FORMAT_ACDB:
114			d->headerfn=&datareader_acdb_header;
115			d->datafn=&datareader_acdb_data;
116			d->eoffn=&datareader_acdb_eof;
117			d->indepfn=&datareader_acdb_time;
118			d->valfn=&datareader_acdb_vals;
119			break;
120	johnpye	812	case DATAREADER_FORMAT_TDV:
121			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Tab delimited values (TDV) format not yet implemenented.");
122			return 1;
123			case DATAREADER_FORMAT_CSV:
124			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Comma-separated values (TDV) format not yet implemenented.");
125			return 1;
126	johnpye	809	default:
127	johnpye	812	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unknown file format '%s' specified",format);
128	johnpye	809	return 1;
129			}
130	johnpye	812
131	johnpye	809	return 0;
132			}
133
134	johnpye	819	typedef struct DataFileSearchData_struct{
135			struct FilePath fp; /< the relative path we're searching for /
136			ospath_stat_t buf; /*< saves memory allocation in the 'test' fn /
137			int error; /*< error code (in case stat or fopen failed) /
138			struct FilePath fp_found; /< the full path we found /
139			} DataFileSearchData;
140
141			FilePathTestFn datareader_searchpath_test;
142
143			/**
144			@return 1 on success
145			*/
146			int datareader_searchpath_test(struct FilePath path,void searchdata){
147			struct FilePath *fp1;
148			DataFileSearchData *sd;
149
150			sd = (DataFileSearchData *)searchdata;
151			assert(sd!=NULL);
152			assert(sd->fp!=NULL);
153
154			fp1 = ospath_concat(path,sd->fp);
155			if(fp1==NULL){
156			CONSOLE_DEBUG("Couldn't concatenate path");
157			return 0;
158			}
159
160			if(ospath_stat(fp1,&sd->buf)){
161			sd->error = errno;
162			ospath_free(fp1);
163			return 0;
164			}
165
166			sd->fp_found = fp1;
167			return 1;
168			};
169
170
171	johnpye	824	/**
172	johnpye	801	Initialise the datareader: open the file, check the number of columns, etc.
173			@return 0 on success
174	johnpye	812
175			@TODO search for the file in the ASCENDLIBRARY if not found immediately
176	johnpye	801	*/
177			int datareader_init(DataReader *d){
178	johnpye	812	ospath_stat_t s;
179	johnpye	819	char *tmp;
180			struct FilePath *sp1, fp2;
181			DataFileSearchData sd;
182	johnpye	812
183	johnpye	801	d->fp = ospath_new(d->fn);
184			if(d->fp==NULL){
185			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Invalid filepath");
186			return 1;
187			}
188	johnpye	809
189	johnpye	812	if(ospath_stat(d->fp,&s)){
190			if(errno==ENOENT){
191	johnpye	819	/* file doesn't exist: check the search path instead */
192	jpye	1524	tmp = Asc_GetEnv(ASC_ENV_LIBRARY);
193	johnpye	819	if(tmp==NULL){
194	jpye	1524	ERROR_REPORTER_HERE(ASC_PROG_ERROR,"No paths to search (is env var '%s' set?)",ASC_ENV_LIBRARY);
195	johnpye	819	return 1;
196			}
197
198			sp1 = ospath_searchpath_new(tmp);
199			if(sp1==NULL){
200	jpye	1524	ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Unable to process %s value '%s'",ASC_ENV_LIBRARY,tmp);
201	johnpye	819	/* memory error */
202			ascfree(tmp);
203			return -3;
204			}
205			ascfree(tmp);
206
207			sd.fp = d->fp;
208
209			fp2 = ospath_searchpath_iterate(sp1, &datareader_searchpath_test, &sd);
210
211			if(fp2==NULL){
212	johnpye	824	ERROR_REPORTER_HERE(ASC_USER_ERROR,"File '%s' not found in search path (error %d)",d->fn,sd.error);
213	johnpye	819	ospath_searchpath_free(sp1);
214			return -1;
215			}
216
217			ospath_searchpath_free(sp1);
218
219			/* replace our relative path with an absolute one */
220			ospath_free(d->fp);
221			d->fp = sd.fp_found;
222
223	johnpye	812	}else{
224	johnpye	819	ERROR_REPORTER_HERE(ASC_USER_ERROR,"The file '%s' cannot be accessed.\n"
225			"Check the file privileges, or try specifying an absolute path.",d->fn
226			);
227	johnpye	812	return 1;
228			}
229			}
230
231			CONSOLE_DEBUG("About to open the data file");
232
233	johnpye	811	d->f = ospath_fopen(d->fp,"r");
234			if(d->f == NULL){
235	johnpye	812	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Unable to open file '%s' for reading.",d->fn);
236	johnpye	809	return 1;
237			}
238
239	johnpye	812	CONSOLE_DEBUG("Data file open ok");
240
241			asc_assert(d->headerfn);
242			asc_assert(d->eoffn);
243			asc_assert(d->datafn);
244
245	johnpye	811	if((*d->headerfn)(d)){
246	johnpye	809	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error processing file header in '%s'",d->fn);
247	johnpye	811	fclose(d->f);
248	johnpye	809	return 1;
249			}
250	johnpye	824
251	johnpye	811	while(! (*d->eoffn)(d)){
252			if((*d->datafn)(d)){
253			ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error reading file data in '%s'",d->fn);
254			fclose(d->f);
255			return 1;
256			}
257	johnpye	809	}
258			fclose(d->f);
259
260	johnpye	812	d->i = 0; /* set current position to zero */
261
262	johnpye	809	return 0;
263	johnpye	801	}
264
265			/**
266			Shut down the data reader and deallocate any memory owned by it, then
267			free the memory at d.
268			*/
269			int datareader_delete(DataReader *d){
270			if(d->fp){
271			ospath_free(d->fp);
272			d->fp = NULL;
273			}
274			if(d->f){
275			fclose(d->f);
276			d->f = NULL;
277			}
278			ASC_FREE(d);
279			return 0;
280			}
281
282			/**
283			Return the number of inputs (independent variables) supplied in the
284			DataReader's current file. Can only be 1 at this stage.
285			*/
286			int datareader_num_inputs(const DataReader *d){
287	johnpye	812	return d->ninputs;
288	johnpye	801	}
289
290			/**
291			Return the number of outputs (dependent variables) found in the DataReader's
292			current file. Should be one or more.
293			*/
294			int datareader_num_outputs(const DataReader *d){
295			return d->noutputs;
296	johnpye	811	}
297	johnpye	801
298	johnpye	809
299	johnpye	812	int datareader_locate(DataReader d, double t, double t1, double *t2){
300			(*d->indepfn)(d,t1);
301			if(*t1 > t && d->i > 0){
302			/* start of current interval is too late */
303			do{
304			--(d->i);
305			(*d->indepfn)(d,t1);
306	johnpye	819	/CONSOLE_DEBUG("STEPPING BACK (d->i = %d currently), t1=%lf",d->i, t1); */
307	johnpye	812	}while(*t1 > t && d->i > 0);
308			}
309			/* now either d->i==0 or t1 < t*/
310	johnpye	819	/* CONSOLE_DEBUG("d->i==%d, t1=%lf",d->i,t1); /
311	johnpye	812	++d->i;
312			(*d->indepfn)(d,t2);
313			if(*t2 <= t){
314			/* end of current interface is too early */
315			do{
316	johnpye	819	/* CONSOLE_DEBUG("STEPPING FORWARD (d->i = %d currently), t1=%lf, t2=%lf",d->i,t1,t2); */
317			(*d->indepfn)(d,t1);
318	johnpye	812	++(d->i);
319			(*d->indepfn)(d,t2);
320			}while(*t2 < t && d->i < d->ndata);
321			}
322	johnpye	819	CONSOLE_DEBUG("d->i==%d, t1[0] = %lf, t2[0] = %lf",d->i,t1[0],t2[0]);
323	johnpye	812
324			if(d->i == d->ndata \|\| d->i == 0){
325			return 1;
326			}
327
328	johnpye	815	/* CONSOLE_DEBUG("INTERVAL OK"); */
329	johnpye	812	return 0;
330			}
331
332	johnpye	801	/**
333			Return an interpolated set of output values for the given input values.
334	johnpye	824	This should be computed such that the output values are smooth in their
335	johnpye	801	first derivatives.
336	johnpye	809
337			The required memory for the inputs and outputs must be allocated by the
338			caller, and indicated by the pointers 'inputs' and 'outputs'.
339
340	johnpye	801	@see datareader_deriv
341			@TODO implement this
342			*/
343			int datareader_func(DataReader d, double inputs, double *outputs){
344	johnpye	809	int i;
345	johnpye	812	double t1[1], t2[1];
346	johnpye	819	double v1[5], v2[5]; /** @TODO this is a fixed size... not good */
347	johnpye	812
348	johnpye	819	double t,g,dt;
349	johnpye	809	t = inputs[0];
350	johnpye	811
351	johnpye	819	CONSOLE_DEBUG("EVALUATING AT t = %lf",inputs[0]);
352	johnpye	812
353			asc_assert(d->indepfn);
354
355			if(datareader_locate(d,t,t1,t2)){
356			CONSOLE_DEBUG("LOCATION ERROR");
357			ERROR_REPORTER_HERE(ASC_USER_ERROR,"Time value t=%f is out of range",t);
358			return 1;
359			}
360	johnpye	819	CONSOLE_DEBUG("LOCATED AT t1 = %lf, t2 = %lf",t1[0], t2[0]);
361	johnpye	812
362			(*d->valfn)(d,v2);
363			--d->i;
364			(*d->valfn)(d,v1);
365	johnpye	819
366			CONSOLE_DEBUG("LOCATED OK, d->i = %d, t1 = %lf, t2 = %lf, v1=%lf, v2=%lf", d->i, t1[0], t2[0],v1[0],v2[0]);
367	johnpye	824
368	johnpye	812	for(i=0;i<d->noutputs;++i){
369	johnpye	819	dt = t2[0] - t1[0];
370			g = (v2[i]-v1[i])/dt;
371			outputs[i]=v1[i]+g(t-(t1));
372			CONSOLE_DEBUG("[%d]: DT = %lf, START = %lf, GRADIENT = %lf, VALUE=%lf",i,dt, v1[i],g,outputs[i]);
373	johnpye	812	}
374	johnpye	824
375	johnpye	812	return 0;
376	johnpye	801	}
377
378			/**
379			Return an interpolated set of output derivatives for the given input
380			values. These should be smooth.
381			@see datareader_func
382			@TODO implement this
383			*/
384			int datareader_deriv(DataReader d, double inputs, double *jacobian){
385	johnpye	812	int i;
386			double t1[1], t2[1];
387			double v1[2], v2[2];
388
389			double t;
390			t = inputs[0];
391
392	johnpye	815	/* CONSOLE_DEBUG("EVALUATING"); */
393	johnpye	812
394			asc_assert(d->indepfn);
395
396			if(datareader_locate(d,t,t1,t2)){
397	johnpye	815	CONSOLE_DEBUG("LOCATION ERROR");
398	johnpye	812	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Time value t=%f is out of range",t);
399			return 1;
400	johnpye	824	}
401	johnpye	812
402			(*d->valfn)(d,v2);
403			--d->i;
404			(*d->valfn)(d,v1);
405	johnpye	824
406	johnpye	812	for(i=0;i<d->noutputs;++i){
407			jacobian[i]=(v2[i]-v1[i])/(t2-t1);
408			}
409	johnpye	824
410	johnpye	812	return 0;
411	johnpye	801	}
412
413	johnpye	809