generic/solver/diffvars.c

/*  ASCEND modelling environment
    Copyright (C) 2007 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.
*//** @file
    Routines for creating 'derivative chains' during system build.
*//*
    By John Pye, Jan 2007
*/
#include "diffvars.h"

#include <utilities/ascMalloc.h>
#include <utilities/ascPanic.h>
#include <utilities/error.h>

#include <compiler/instance_io.h>

#include "analyse_impl.h"
#include "system_impl.h"

/* #define DIFFVARS_DEBUG */

/*------------------------------------------------------------------------------
    DERIVATIVES & DIFFERENTIAL VARIABLES
*/

SolverDiffVarCollection *system_get_diffvars(slv_system_t sys){
    return sys->diffvars;
}

static 
int CmpDiffVars(const struct solver_ipdata *a, const struct solver_ipdata *b){
    if(a->u.v.odeid < b->u.v.odeid)return -1;
    if(a->u.v.odeid > b->u.v.odeid)return 1;
    /* now, the ode_id is equal: next level of sorting... */
    if(a->u.v.deriv < b->u.v.deriv)return -1;
    if(a->u.v.deriv > b->u.v.deriv)return 1;
    return 0;
}

/**
    This function steals a lot of what was in integrator.c, but we try to make
    it more general and capable of extracting info about high-order derivative
    chains. One eye open to the possibility of doing Pantelides (or other?) 
    style index reduction.

    Also, note that we try to only use data from the struct problem_t, NOT
    to go back to the instance hierarchy. The idea is that we have already got
    what we need in classify_instance, etc.

    @return 0 on success
*/
int system_generate_diffvars(slv_system_t sys, struct problem_t *prob){
    SolverDiffVarCollection *diffvars = NULL;
    struct solver_ipdata *vip, *vipnext;
    SolverDiffVarSequence *seq;
    struct gl_list_t *seqs;
    long i, seqstart, nalg, ndiff;
    short j;
    char cont;
    short maxorder = 0;

    asc_assert(prob);
    
    if(gl_length(prob->diffvars)==0){
        CONSOLE_DEBUG("No differential variables were seen. Skipping generation of diffvars struct.");
        return 0;
    }

    seqs = gl_create(prob->nr);
    
    /* add the list of algebraic variables to the structure too */
    asc_assert(prob->algebvars);
    nalg = gl_length(prob->algebvars);
    for(i=0; i<nalg; ++i){
        vip = (struct solver_ipdata *)gl_fetch(prob->algebvars,i+1);    
        seq = ASC_NEW(SolverDiffVarSequence);
        seq->n = 1;
        seq->ode_id = 0;
        seq->vars = ASC_NEW_ARRAY(struct var_variable*,1);
        seq->vars[0] = vip->u.v.data;
        gl_append_ptr(seqs,(void *)seq);
    }
#ifdef DIFFVARS_DEBUG
    CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);
    CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
#endif
    
    /* first sort the list of diffvars */
    gl_sort(prob->diffvars, (CmpFunc)CmpDiffVars);
    
    /* scan the list for derivs that are missing vars */
    i = 1; cont = TRUE; seqstart =1; ndiff = 0;
    vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,i);
    if(vip->u.v.deriv > 1){
        ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
        FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
        WriteInstanceName(ASCERR,vip->i,prob->root);
        FPRINTF(ASCERR,"')");
        error_reporter_end_flush();
        return 2;
    }
    while(cont){
        /* CONSOLE_DEBUG("Working, seqstart=%ld",seqstart); */
        if(i >= gl_length(prob->diffvars))cont = FALSE;
        else vipnext = (struct solver_ipdata *)gl_fetch(prob->diffvars,i+1);

        if(cont && vipnext->u.v.odeid == vip->u.v.odeid){
            /* same sequence, check that it's the next derivative */
            asc_assert(vip->u.v.deriv >= vip->u.v.deriv);

            if(vipnext->u.v.deriv == vip->u.v.deriv){
                ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
                FPRINTF(ASCERR,"Repeated ode_type %d for ode_id %d (var '%s')",vip->u.v.deriv,vip->u.v.odeid);
                WriteInstanceName(ASCERR,vip->i,prob->root);
                FPRINTF(ASCERR,"')");
                error_reporter_end_flush();
                return 1;
            }else if(vipnext->u.v.deriv > vip->u.v.deriv + 1){
                ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
                FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
                WriteInstanceName(ASCERR,vip->i,prob->root);
                FPRINTF(ASCERR,"')");
                error_reporter_end_flush();
                return 2;
            }

            /* allow this var into the sequence, loop again */
            vip = vipnext; ++i; continue;
        }

        /* close the current sequence */
        if(i == seqstart){
            ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
            FPRINTF(ASCERR,"Variable '");
            WriteInstanceName(ASCWAR,vip->i,prob->root);
            FPRINTF(ASCERR,"' declared differential without derivative being identified (ode_id=%d)",vip->u.v.odeid);
            error_reporter_end_flush();
            return 3;
        }
        seq = ASC_NEW(SolverDiffVarSequence);
        seq->n = i - seqstart + 1;
        seq->ode_id = vip->u.v.odeid;
        seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
#ifdef DIFFVARS_DEBUG
        CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n);
#endif
        for(j=0;j<seq->n;++j){
            vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
            seq->vars[j]=vip->u.v.data;
            /* set the VAR_ALGEB flag as req */
            var_set_diff(seq->vars[j],j==0 && seq->n > 1);
#ifdef DIFFVARS_DEBUG
            CONSOLE_DEBUG("seq, j=%d: is %s",j
                ,(var_diff(seq->vars[j]) ? "DIFF" : 
                    (var_deriv(seq->vars[j]) ? "deriv" : "alg")
                ) 
            );
#endif
        }
        gl_append_ptr(seqs,(void *)seq);
        ndiff++;
#ifdef DIFFVARS_DEBUG
        CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs));
#endif

        if(vip->u.v.deriv > maxorder){
            maxorder = vip->u.v.deriv;
            /* CONSOLE_DEBUG("maxorder increased to %d",maxorder); */
        }

        ++i; seqstart=i; vip=vipnext;
        if(cont && vip->u.v.deriv != 1){
            ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
            FPRINTF(ASCERR,"Missing undifferentiated form of var '");
            WriteInstanceName(ASCWAR,vip->i,prob->root);
            FPRINTF(ASCERR,"' (ode_id = %d, ode_type = %d)",vip->u.v.odeid,vip->u.v.deriv);
            error_reporter_end_flush();
            return 4;
        }
        continue;
    }
    
#ifdef DIFFVARS_DEBUG
    CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);
#endif

    diffvars = ASC_NEW(SolverDiffVarCollection);
    diffvars->maxorder = maxorder;
    diffvars->nseqs = gl_length(seqs);
    diffvars->seqs = ASC_NEW_ARRAY(SolverDiffVarSequence,diffvars->nseqs);
    for(i=0;i<diffvars->nseqs;++i){
        diffvars->seqs[i] = *(SolverDiffVarSequence *)gl_fetch(seqs,i+1);
    }
    gl_destroy(seqs);
    diffvars->nalg = nalg;
    diffvars->ndiff = ndiff;

    diffvars->nindep = gl_length(prob->indepvars);
    diffvars->indep = ASC_NEW_ARRAY(struct var_variable *,diffvars->nindep);
    for(i=0;i<diffvars->nindep;++i){
        vip = (struct solver_ipdata *)gl_fetch(prob->indepvars,i+1);
        diffvars->indep[i] = vip->u.v.data;
    }
#ifdef DIFFVARS_DEBUG
    CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);
#endif

    /* I know, I know, this might not be the best place for this... */
    diffvars->nobs = gl_length(prob->obsvars);
    diffvars->obs = ASC_NEW_ARRAY(struct var_variable *,diffvars->nobs);
    for(i=0;i<diffvars->nobs;++i){
        vip = (struct solver_ipdata *)gl_fetch(prob->obsvars,i+1);
        diffvars->obs[i] = vip->u.v.data;
    }
#ifdef DIFFVARS_DEBUG
    CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);
    CONSOLE_DEBUG("There were %ld rels",prob->nr);
#endif

    slv_set_diffvars(sys,(void *)diffvars);

    return 0;
}


int system_diffvars_debug(slv_system_t sys,FILE *fp){
    int i, j;
    char *varname;
    const SolverDiffVarCollection *diffvars;
    SolverDiffVarSequence seq;
    diffvars = sys->diffvars;
    if(diffvars==NULL){
        fprintf(fp,"NO DIFFVARS (NULL)");
        return 0;
    }
    for(i=0; i<diffvars->nseqs;++i){
        seq = diffvars->seqs[i];
        fprintf(fp,"%d: ",i);
        for(j=0; j<seq.n; ++j){
            if(j)fprintf(fp," <-- ");
            fprintf(fp,"%d: (%p)",var_sindex(seq.vars[j]),seq.vars[j]);
            varname = var_make_name(sys,seq.vars[j]);
            fprintf(fp,"'%s'",varname);
            ASC_FREE(varname);
        }
        fprintf(fp,"\n");
    }
    return 0;
}

1	/* ASCEND modelling environment
2	Copyright (C) 2007 Carnegie Mellon University
3
4	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
9	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
14	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	//* @file
19	Routines for creating 'derivative chains' during system build.
20	//
21	By John Pye, Jan 2007
22	*/
23	#include "diffvars.h"
24
25	#include <utilities/ascMalloc.h>
26	#include <utilities/ascPanic.h>
27	#include <utilities/error.h>
28
29	#include <compiler/instance_io.h>
30
31	#include "analyse_impl.h"
32	#include "system_impl.h"
33
34	/* #define DIFFVARS_DEBUG */
35
36	/*------------------------------------------------------------------------------
37	DERIVATIVES & DIFFERENTIAL VARIABLES
38	*/
39
40	SolverDiffVarCollection *system_get_diffvars(slv_system_t sys){
41	return sys->diffvars;
42	}
43
44	static
45	int CmpDiffVars(const struct solver_ipdata a, const struct solver_ipdata b){
46	if(a->u.v.odeid < b->u.v.odeid)return -1;
47	if(a->u.v.odeid > b->u.v.odeid)return 1;
48	/* now, the ode_id is equal: next level of sorting... */
49	if(a->u.v.deriv < b->u.v.deriv)return -1;
50	if(a->u.v.deriv > b->u.v.deriv)return 1;
51	return 0;
52	}
53
54	/**
55	This function steals a lot of what was in integrator.c, but we try to make
56	it more general and capable of extracting info about high-order derivative
57	chains. One eye open to the possibility of doing Pantelides (or other?)
58	style index reduction.
59
60	Also, note that we try to only use data from the struct problem_t, NOT
61	to go back to the instance hierarchy. The idea is that we have already got
62	what we need in classify_instance, etc.
63
64	@return 0 on success
65	*/
66	int system_generate_diffvars(slv_system_t sys, struct problem_t *prob){
67	SolverDiffVarCollection *diffvars = NULL;
68	struct solver_ipdata vip, vipnext;
69	SolverDiffVarSequence *seq;
70	struct gl_list_t *seqs;
71	long i, seqstart, nalg, ndiff;
72	short j;
73	char cont;
74	short maxorder = 0;
75
76	asc_assert(prob);
77
78	if(gl_length(prob->diffvars)==0){
79	CONSOLE_DEBUG("No differential variables were seen. Skipping generation of diffvars struct.");
80	return 0;
81	}
82
83	seqs = gl_create(prob->nr);
84
85	/* add the list of algebraic variables to the structure too */
86	asc_assert(prob->algebvars);
87	nalg = gl_length(prob->algebvars);
88	for(i=0; i<nalg; ++i){
89	vip = (struct solver_ipdata *)gl_fetch(prob->algebvars,i+1);
90	seq = ASC_NEW(SolverDiffVarSequence);
91	seq->n = 1;
92	seq->ode_id = 0;
93	seq->vars = ASC_NEW_ARRAY(struct var_variable*,1);
94	seq->vars[0] = vip->u.v.data;
95	gl_append_ptr(seqs,(void *)seq);
96	}
97	#ifdef DIFFVARS_DEBUG
98	CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);
99	CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
100	#endif
101
102	/* first sort the list of diffvars */
103	gl_sort(prob->diffvars, (CmpFunc)CmpDiffVars);
104
105	/* scan the list for derivs that are missing vars */
106	i = 1; cont = TRUE; seqstart =1; ndiff = 0;
107	vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,i);
108	if(vip->u.v.deriv > 1){
109	ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
110	FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
111	WriteInstanceName(ASCERR,vip->i,prob->root);
112	FPRINTF(ASCERR,"')");
113	error_reporter_end_flush();
114	return 2;
115	}
116	while(cont){
117	/* CONSOLE_DEBUG("Working, seqstart=%ld",seqstart); */
118	if(i >= gl_length(prob->diffvars))cont = FALSE;
119	else vipnext = (struct solver_ipdata *)gl_fetch(prob->diffvars,i+1);
120
121	if(cont && vipnext->u.v.odeid == vip->u.v.odeid){
122	/* same sequence, check that it's the next derivative */
123	asc_assert(vip->u.v.deriv >= vip->u.v.deriv);
124
125	if(vipnext->u.v.deriv == vip->u.v.deriv){
126	ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
127	FPRINTF(ASCERR,"Repeated ode_type %d for ode_id %d (var '%s')",vip->u.v.deriv,vip->u.v.odeid);
128	WriteInstanceName(ASCERR,vip->i,prob->root);
129	FPRINTF(ASCERR,"')");
130	error_reporter_end_flush();
131	return 1;
132	}else if(vipnext->u.v.deriv > vip->u.v.deriv + 1){
133	ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
134	FPRINTF(ASCERR,"Missing ode_type %d for ode_id %d (check var '",vip->u.v.deriv+1,vip->u.v.odeid);
135	WriteInstanceName(ASCERR,vip->i,prob->root);
136	FPRINTF(ASCERR,"')");
137	error_reporter_end_flush();
138	return 2;
139	}
140
141	/* allow this var into the sequence, loop again */
142	vip = vipnext; ++i; continue;
143	}
144
145	/* close the current sequence */
146	if(i == seqstart){
147	ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
148	FPRINTF(ASCERR,"Variable '");
149	WriteInstanceName(ASCWAR,vip->i,prob->root);
150	FPRINTF(ASCERR,"' declared differential without derivative being identified (ode_id=%d)",vip->u.v.odeid);
151	error_reporter_end_flush();
152	return 3;
153	}
154	seq = ASC_NEW(SolverDiffVarSequence);
155	seq->n = i - seqstart + 1;
156	seq->ode_id = vip->u.v.odeid;
157	seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
158	#ifdef DIFFVARS_DEBUG
159	CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n);
160	#endif
161	for(j=0;j<seq->n;++j){
162	vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
163	seq->vars[j]=vip->u.v.data;
164	/* set the VAR_ALGEB flag as req */
165	var_set_diff(seq->vars[j],j==0 && seq->n > 1);
166	#ifdef DIFFVARS_DEBUG
167	CONSOLE_DEBUG("seq, j=%d: is %s",j
168	,(var_diff(seq->vars[j]) ? "DIFF" :
169	(var_deriv(seq->vars[j]) ? "deriv" : "alg")
170	)
171	);
172	#endif
173	}
174	gl_append_ptr(seqs,(void *)seq);
175	ndiff++;
176	#ifdef DIFFVARS_DEBUG
177	CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs));
178	#endif
179
180	if(vip->u.v.deriv > maxorder){
181	maxorder = vip->u.v.deriv;
182	/* CONSOLE_DEBUG("maxorder increased to %d",maxorder); */
183	}
184
185	++i; seqstart=i; vip=vipnext;
186	if(cont && vip->u.v.deriv != 1){
187	ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
188	FPRINTF(ASCERR,"Missing undifferentiated form of var '");
189	WriteInstanceName(ASCWAR,vip->i,prob->root);
190	FPRINTF(ASCERR,"' (ode_id = %d, ode_type = %d)",vip->u.v.odeid,vip->u.v.deriv);
191	error_reporter_end_flush();
192	return 4;
193	}
194	continue;
195	}
196
197	#ifdef DIFFVARS_DEBUG
198	CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);
199	#endif
200
201	diffvars = ASC_NEW(SolverDiffVarCollection);
202	diffvars->maxorder = maxorder;
203	diffvars->nseqs = gl_length(seqs);
204	diffvars->seqs = ASC_NEW_ARRAY(SolverDiffVarSequence,diffvars->nseqs);
205	for(i=0;i<diffvars->nseqs;++i){
206	diffvars->seqs[i] = (SolverDiffVarSequence )gl_fetch(seqs,i+1);
207	}
208	gl_destroy(seqs);
209	diffvars->nalg = nalg;
210	diffvars->ndiff = ndiff;
211
212	diffvars->nindep = gl_length(prob->indepvars);
213	diffvars->indep = ASC_NEW_ARRAY(struct var_variable *,diffvars->nindep);
214	for(i=0;i<diffvars->nindep;++i){
215	vip = (struct solver_ipdata *)gl_fetch(prob->indepvars,i+1);
216	diffvars->indep[i] = vip->u.v.data;
217	}
218	#ifdef DIFFVARS_DEBUG
219	CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);
220	#endif
221
222	/* I know, I know, this might not be the best place for this... */
223	diffvars->nobs = gl_length(prob->obsvars);
224	diffvars->obs = ASC_NEW_ARRAY(struct var_variable *,diffvars->nobs);
225	for(i=0;i<diffvars->nobs;++i){
226	vip = (struct solver_ipdata *)gl_fetch(prob->obsvars,i+1);
227	diffvars->obs[i] = vip->u.v.data;
228	}
229	#ifdef DIFFVARS_DEBUG
230	CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);
231	CONSOLE_DEBUG("There were %ld rels",prob->nr);
232	#endif
233
234	slv_set_diffvars(sys,(void *)diffvars);
235
236	return 0;
237	}
238
239
240	int system_diffvars_debug(slv_system_t sys,FILE *fp){
241	int i, j;
242	char *varname;
243	const SolverDiffVarCollection *diffvars;
244	SolverDiffVarSequence seq;
245	diffvars = sys->diffvars;
246	if(diffvars==NULL){
247	fprintf(fp,"NO DIFFVARS (NULL)");
248	return 0;
249	}
250	for(i=0; i<diffvars->nseqs;++i){
251	seq = diffvars->seqs[i];
252	fprintf(fp,"%d: ",i);
253	for(j=0; j<seq.n; ++j){
254	if(j)fprintf(fp," <-- ");
255	fprintf(fp,"%d: (%p)",var_sindex(seq.vars[j]),seq.vars[j]);
256	varname = var_make_name(sys,seq.vars[j]);
257	fprintf(fp,"'%s'",varname);
258	ASC_FREE(varname);
259	}
260	fprintf(fp,"\n");
261	}
262	return 0;
263	}
264