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 <compiler/instance_io.h>

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

/*------------------------------------------------------------------------------
    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/*, *varname*/;
    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);
    }
    CONSOLE_DEBUG("Added %ld algebraic vars to chains",nalg);

    CONSOLE_DEBUG("Sorting %ld differential & derivative vars...",gl_length(prob->diffvars));
    
    /* 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(ASCWAR,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(vipnext->u.v.odeid == vip->u.v.odeid);
            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 '",vip->u.v.deriv,vip->u.v.odeid);
                WriteInstanceName(ASCWAR,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(ASCWAR,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;
        /* CONSOLE_DEBUG("Saving sequence ode_id = %ld, n = %d",seq->ode_id,seq->n); */
        seq->vars = ASC_NEW_ARRAY(struct var_variable*,seq->n);
        for(j=0;j<seq->n;++j){
            vip = (struct solver_ipdata *)gl_fetch(prob->diffvars,seqstart+j);
            seq->vars[j]=vip->u.v.data;
            /* varname = var_make_name(sys,seq->vars[j]);
            CONSOLE_DEBUG("seq, j=%d: '%s'",j,varname);
            ASC_FREE(varname); */
        }
        gl_append_ptr(seqs,(void *)seq);
        ndiff++;
        /* CONSOLE_DEBUG("Completed seq %ld",gl_length(seqs)); */

        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;
    }
    
    CONSOLE_DEBUG("Identified %ld derivative chains, maximum length %d...",gl_length(seqs),maxorder);

    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;
    }
    CONSOLE_DEBUG("Identified %ld indep vars",diffvars->nindep);

    /* 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;
    }
    CONSOLE_DEBUG("Identified %ld obs vers",diffvars->nobs);

    CONSOLE_DEBUG("There were %ld rels",prob->nr);

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