generic/integrator/idaanalyse.c

#include "idaanalyse.h"

#include <utilities/ascPanic.h>

#ifdef ASC_IDA_NEW_ANALYSE
# include <solver/diffvars.h>
# include <solver/slv_stdcalls.h>
# include <solver/block.h>
#include <solver/slvDOF.h>
#endif

#define ANALYSE_DEBUG

static int integrator_ida_check_partitioning(IntegratorSystem *sys);
static int integrator_ida_check_diffindex(IntegratorSystem *sys);
static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys);

const var_filter_t integrator_ida_nonderiv = {
    VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
    VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0         | 0
};

const var_filter_t integrator_ida_deriv = {
    VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
    VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0         | VAR_DERIV
};

const rel_filter_t integrator_ida_rel = {
    REL_INCLUDED | REL_EQUALITY | REL_ACTIVE,
    REL_INCLUDED | REL_EQUALITY | REL_ACTIVE 
};

/**
    Check derivative chains: if a var or derivative is inelligible,
    follow down the chain fixing & setting zero any derivatives.
*/
static int integrator_ida_check_vars(IntegratorSystem *sys){
    const SolverDiffVarCollection *diffvars;
    char *varname, *derivname;
    int n_y = 0;
    int i, j;
    struct var_variable *v;

    SolverDiffVarSequence seq;

    /* we shouldn't have allocated these yet: just be sure */
    asc_assert(sys->y==NULL);
    asc_assert(sys->ydot==NULL);
    asc_assert(sys->n_y==0);

    /* get the the dervative chains from the system */
    diffvars = slv_get_diffvars(sys->system);
    if(diffvars==NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
        return 1;
    }
    
    /* add the variables from the derivative chains */
    for(i=0; i<diffvars->nseqs; ++i){
        seq = diffvars->seqs[i];
        asc_assert(seq.n >= 1);
        v = seq.vars[0];
        if(!var_apply_filter(v,&integrator_ida_nonderiv)){
            varname = var_make_name(sys->system,v);
            CONSOLE_DEBUG("'%s' fails non-deriv filter",varname);
            ASC_FREE(varname);          
            for(j=1;j<seq.n;++j){
                v = seq.vars[j];
                var_set_active(v,FALSE);
                var_set_value(v,0);
                varname = var_make_name(sys->system,v);             
                CONSOLE_DEBUG("Derivative '%s' SET ZERO AND INACTIVE",varname);
                ASC_FREE(v);
            }
            continue;
        }

        varname = var_make_name(sys->system,v);
        CONSOLE_DEBUG("We will use var '%s'",varname);
        ASC_FREE(varname);
        if(seq.n > 2){
            varname = var_make_name(sys->system,v);             
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'");
            ASC_FREE(varname);
            return 2;
        }else if(seq.n==2){
            /* diff var */
            if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
                /* add the diff & deriv vars to the lists */
                n_y++;
                varname = var_make_name(sys->system,v);
                derivname = var_make_name(sys->system,seq.vars[1]);
                CONSOLE_DEBUG("Added '%s' and '%s'",varname,derivname);
                ASC_FREE(varname);
                ASC_FREE(derivname);
                continue;
            }
            varname = var_make_name(sys->system,v);
            derivname = var_make_name(sys->system,seq.vars[1]);
            CONSOLE_DEBUG("Derivative '%s' of '%s' fails filter; convert '%s' to algebraic",derivname,varname,varname);
            ASC_FREE(varname);
            ASC_FREE(derivname);
            /* fall through */
        }

        varname = var_make_name(sys->system,v);
        CONSOLE_DEBUG("Adding '%s' to algebraic",varname);
        ASC_FREE(varname);
        n_y++;
    }

    /* we assert that all vars in y meet the integrator_ida_nonderiv filter */
    /* we assert that all vars in ydot meet the integrator_ida_deriv filter */

    CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
    sys->n_y = n_y;
    return 0;
}


/**
    Sort the lists. First we will put the non-derivative vars, then we will
    put the derivative vars. Then we will put all the others.
*/
static int integrator_ida_sort_rels_and_vars(IntegratorSystem *sys){
    int ny1, nydot, nr;

    /* we should not have allocated y or ydot yet */
    asc_assert(sys->y==NULL && sys->ydot==NULL);

    /* but we should have found some variables (and know how many) */
    asc_assert(sys->n_y);

    if(system_cut_vars(sys->system, 0, &integrator_ida_nonderiv, &ny1)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
        return 1;
    }

    CONSOLE_DEBUG("cut_vars: ny1 = %d, sys->n_y = %d",ny1,sys->n_y);
    asc_assert(ny1 == sys->n_y);

    if(system_cut_vars(sys->system, ny1, &integrator_ida_deriv, &nydot)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
        return 1;
    }

    if(system_cut_rels(sys->system, 0, &integrator_ida_rel, &nr)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
        return 1;
    }

    if(ny1 != nr){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
        return 2;
    }

    return 0;
}

/**
    Build up the lists 'ydot' and 'y_id.
    'ydot' allows us to find the derivative of a variable given its sindex.
    'y_id' allows us to locate the variable of a derivative given its sindex.

    Hence
    y_id[var_sindex(derivvar)-sys->n_y] = var_sindex(diffvar)
    ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)

    Note that there is 'shifting' required when looking up y_id.

    Note that we want to get rid of 'y' but for the moment we are also assigning
    that. We want to be rid of it because if the vars are ordered correctly
    it shouldn't be needed.
*/
static int integrator_ida_create_lists(IntegratorSystem *sys){
    const SolverDiffVarCollection *diffvars;
    int i, j;
    struct var_variable *v;

    SolverDiffVarSequence seq;

    asc_assert(sys->y==NULL);
    asc_assert(sys->ydot==NULL);
    asc_assert(sys->y_id== NULL);

    /* get the the dervative chains from the system */
    diffvars = slv_get_diffvars(sys->system);
    if(diffvars==NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
        return 1;
    }
    
    /* allocate space (more than enough) */
    sys->y = ASC_NEW_ARRAY(struct var_variable *,sys->n_y);
    sys->y_id = ASC_NEW_ARRAY_CLEAR(int,sys->n_y);
    sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,sys->n_y);
    j = 0;

    /* add the variables from the derivative chains */
    for(i=0; i<diffvars->nseqs; ++i){
        seq = diffvars->seqs[i];
        asc_assert(seq.n >= 1);
        v = seq.vars[0];
        if(!var_apply_filter(v,&integrator_ida_nonderiv)){
            continue;
        }

        if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
            asc_assert(var_sindex(seq.vars[1])-sys->n_y >= 0);
            asc_assert(var_sindex(seq.vars[1])-sys->n_y < sys->n_y);
            sys->y_id[var_sindex(seq.vars[1])-sys->n_y] = j;
            sys->ydot[j] = seq.vars[1];
        }else{
            asc_assert(sys->ydot[j]==NULL);
        }

        sys->y[j] = v;
        j++;
    }

    return 0;
}

/** 
    Perform additional problem analysis to prepare problem for integration with 
    IDA.

    We assume that the analyse_generate_diffvars routine has been called, so 
    that we just need to call slv_get_diffvars to access the derivative
    chains.

    We can also assume that the independent variable has been found.

    See mailing list, ~Jan 2007.

    Note, the stuff for identifying the static and output sub-problems should
    be part of integrator.c, not this file. We will assume this is handled 

    @return 0 on success 
    @see integrator_analyse
*/
int integrator_ida_analyse(IntegratorSystem *sys){
    int res;
    const SolverDiffVarCollection *diffvars;
    char *varname;
    int i;

    asc_assert(sys->engine==INTEG_IDA);

#ifdef ANALYSE_DEBUG
    CONSOLE_DEBUG("Starting IDA analysis");
#endif

    res = integrator_ida_check_vars(sys);
    if(res){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
        return 1;
    }

#ifdef ANALYSE_DEBUG
    CONSOLE_DEBUG("Sorting rels and vars");
#endif

    res = integrator_ida_sort_rels_and_vars(sys);
    if(res){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
        return 1;
    }

#ifdef ANALYSE_DEBUG
    CONSOLE_DEBUG("Creating lists");
#endif

    res = integrator_ida_create_lists(sys);
    if(res){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
        return 1;
    }

#ifdef ANALYSE_DEBUG
    CONSOLE_DEBUG("Checking");
#endif

    asc_assert(sys->y);
    asc_assert(sys->ydot);
    asc_assert(sys->y_id);

    integrator_ida_analyse_debug(sys,stderr);

    if(integrator_ida_check_diffindex(sys)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
        return 360;
    }

    /* the following causes TestIDA.testincidence3 to fail:
    if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
        return 2;
    }*/

#if 0
    /* check structural singularity for the two IDACalcIC scenarios */

    /* ...(1) FIX the derivatives */
    CONSOLE_DEBUG("Checking system with derivatives fixed...");
    for(i=0;i<n_y;++i){
        if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
    }

    slv_block_partition(sys->system);
    res = integrator_ida_block_check(sys);
    if(res)return 100 + res;

    /* ...(2) FREE the derivatives, FIX the diffvars */
    CONSOLE_DEBUG("Checking system with differential variables fixed...");
    for(i=0;i<n_y;++i){
        if(sys->ydot[i]){
            var_set_fixed(sys->ydot[i],0);
            var_set_fixed(sys->y[i],1);
        }
    }
    slv_block_partition(sys->system);
    res = integrator_ida_block_check(sys);
    if(res)return 200 + res;
    
    /* ...(3) restore as it was, FREE the diffvars */
    for(i=0;i<n_y;++i){
        if(sys->ydot[i]){
            var_set_fixed(sys->y[i],0);
        }
    }   
#endif

    /* get the the dervative chains from the system */
    diffvars = slv_get_diffvars(sys->system);

    /* check the indep var is same as was located elsewhere */
    if(diffvars->nindep>1){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
        return 3;
    }else if(diffvars->nindep<1){
        ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
        return 4;
    }else if(diffvars->indep[0]!=sys->x){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
        return 5;
    }

    /* get the observations */
    /** @TODO this should use a 'problem provider' API instead of static data
    from the system object */
    sys->n_obs = diffvars->nobs;
    sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
    for(i=0;i<sys->n_obs;++i){
        sys->obs[i] = diffvars->obs[i];
        varname = var_make_name(sys->system,sys->obs[i]);
        CONSOLE_DEBUG("'%s' is observation",varname);
        ASC_FREE(varname);
    }

    /*   - 'y' list as [ya|yd] */
    /*   - sparsity pattern for dF/dy and dF/dy' */
    /*   - sparsity pattern for union of above */
    /*   - block decomposition based on above */
    /*   - block decomposition results in reordering of y and y' */
    /*   - boundaries (optional) */
    /* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
    return -1; */

    return 0;
}

/*------------------------------------------------------------------------------
  ANALYSIS ROUTINE (new implementation)
*/

static int integrator_ida_check_partitioning(IntegratorSystem *sys){
    long i, nv;
    int err=0;
    char *varname;
    struct var_variable **vlist, *v;
    nv = slv_get_num_solvers_vars(sys->system);
    vlist = slv_get_solvers_var_list(sys->system);
    var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
                      ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
    for(i=0;i<nv;++i){
        v=vlist[i];
        if(!var_apply_filter(v,&vf))continue;
        varname = var_make_name(sys->system,v);
        if(!var_deriv(v)){
            fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
            if(i>=sys->n_y){
                ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
                err++;
            }
        }else{
            fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
            if(i<sys->n_y){
                ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
                    ,varname, sys->n_y, i
                );
                err++;
            }
        }
        ASC_FREE(varname);
    }
    if(err){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
        return 1;
    }
    return 0;
}

int integrator_ida_block_check(IntegratorSystem *sys){
    int res;
    int dof;
#ifdef ANALYSE_DEBUG
    long *vlist, *vp, i, nv, nv_ok;
    char *varname;
    struct var_variable **solversvars;

    var_filter_t vfilt = {
        VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
        VAR_ACTIVE | VAR_INCIDENT | 0
    };
            
    nv = slv_get_num_solvers_vars(sys->system);
    solversvars = slv_get_solvers_var_list(sys->system);
    CONSOLE_DEBUG("-------------- nv = %ld -------------",nv);
    for(nv_ok=0, i=0; i < nv; ++i){
        if(var_apply_filter(solversvars[i],&vfilt)){
            varname = var_make_name(sys->system,solversvars[i]);
            fprintf(stderr,"%s\n",varname);
            ASC_FREE(varname);
            nv_ok++;
        }
    }
    CONSOLE_DEBUG("----------- got %ld ok -------------",nv_ok);
#endif

    if(!slvDOF_status(sys->system, &res, &dof)){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
        return -1;
    }
    switch(res){
        case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
        case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
        case 3: CONSOLE_DEBUG("System is structurally singular"); break;
        case 4: CONSOLE_DEBUG("System is overspecified"); break;
        default: 
            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
            return -2;
    }

#ifdef ANALYSE_DEBUG
    /* if it was underspecified, what vars could be fixed? */
    if(res==1){
        CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
        solversvars = slv_get_solvers_var_list(sys->system);
        if(!slvDOF_eligible(sys->system, &vlist)){
            ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
            return -3;
        }
        for(vp=vlist;*vp!=-1;++vp){
            varname = var_make_name(sys->system, solversvars[*vp]);
            CONSOLE_DEBUG("Fixable var: %s",varname);
            ASC_FREE(varname);
        }
        CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
        return 1;
    }
#endif

    return res;
}

static int integrator_ida_check_diffindex(IntegratorSystem *sys){
    int i, nv, err = 0;
    struct var_variable **vlist;
    int diffindex;

    CONSOLE_DEBUG("Checking diffindex vector");

    if(sys->y_id == NULL){
        CONSOLE_DEBUG("y_id not allocated");
        return 1;
    }

    vlist = slv_get_solvers_var_list(sys->system);
    nv = slv_get_num_solvers_vars(sys->system);

    for(i=0; i<nv; ++i){
        if(var_deriv(vlist[i])){
            if(i < sys->n_y){
                ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv in wrong position");
                err++;
            }else{
                if(var_sindex(vlist[i])!=i){
                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv var with incorrect sindex");
                    err++;
                }
                diffindex = sys->y_id[i - sys->n_y];
                if(diffindex >= sys->n_y){
                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff y_id[%d] value too big",i-sys->n_y);
                    err++;
                }
                if(var_sindex(vlist[diffindex])!=diffindex){
                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff var with incorrect sindex");
                    err++;
                }
            }
        }else{
            if(i!=var_sindex(vlist[i])){
                ERROR_REPORTER_HERE(ASC_PROG_ERR,"var_sindex incorrect");
                err++;
            }
        }
    }
    if(err){
        CONSOLE_DEBUG("Errors found in diffindex");
    }
    return err;
}

/**
    @TODO provide a macro implementation of this
*/
int integrator_ida_diffindex(const IntegratorSystem *sys, const struct var_variable *deriv){
    int diffindex;
#ifdef DIFFINDEX_DEBUG
    asc_assert( var_deriv    (deriv));
    asc_assert( var_active   (deriv));
    asc_assert( var_incident (deriv));
    asc_assert( var_svar     (deriv));

    asc_assert(!var_fixed    (deriv));

    asc_assert(var_sindex(deriv) >= sys->n_y);
    asc_assert(diffindex == var_sindex(sys->y[diffindex]));
#endif
    asc_assert(var_sindex(deriv) >= sys->n_y);
    diffindex = sys->y_id[var_sindex(deriv) - sys->n_y];

    return diffindex;
}


int integrator_ida_analyse_debug(const IntegratorSystem *sys,FILE *fp){
    return analyse_diffvars_debug(sys->system,fp);
}

1	#include "idaanalyse.h"
2
3	#include <utilities/ascPanic.h>
4
5	#ifdef ASC_IDA_NEW_ANALYSE
6	# include <solver/diffvars.h>
7	# include <solver/slv_stdcalls.h>
8	# include <solver/block.h>
9	#include <solver/slvDOF.h>
10	#endif
11
12	#define ANALYSE_DEBUG
13
14	static int integrator_ida_check_partitioning(IntegratorSystem *sys);
15	static int integrator_ida_check_diffindex(IntegratorSystem *sys);
16	static int integrator_ida_rebuild_diffindex(IntegratorSystem *sys);
17
18	const var_filter_t integrator_ida_nonderiv = {
19	VAR_SVAR \| VAR_INCIDENT \| VAR_ACTIVE \| VAR_FIXED \| VAR_DERIV,
20	VAR_SVAR \| VAR_INCIDENT \| VAR_ACTIVE \| 0 \| 0
21	};
22
23	const var_filter_t integrator_ida_deriv = {
24	VAR_SVAR \| VAR_INCIDENT \| VAR_ACTIVE \| VAR_FIXED \| VAR_DERIV,
25	VAR_SVAR \| VAR_INCIDENT \| VAR_ACTIVE \| 0 \| VAR_DERIV
26	};
27
28	const rel_filter_t integrator_ida_rel = {
29	REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE,
30	REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE
31	};
32
33	/**
34	Check derivative chains: if a var or derivative is inelligible,
35	follow down the chain fixing & setting zero any derivatives.
36	*/
37	static int integrator_ida_check_vars(IntegratorSystem *sys){
38	const SolverDiffVarCollection *diffvars;
39	char varname, derivname;
40	int n_y = 0;
41	int i, j;
42	struct var_variable *v;
43
44	SolverDiffVarSequence seq;
45
46	/* we shouldn't have allocated these yet: just be sure */
47	asc_assert(sys->y==NULL);
48	asc_assert(sys->ydot==NULL);
49	asc_assert(sys->n_y==0);
50
51	/* get the the dervative chains from the system */
52	diffvars = slv_get_diffvars(sys->system);
53	if(diffvars==NULL){
54	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
55	return 1;
56	}
57
58	/* add the variables from the derivative chains */
59	for(i=0; i<diffvars->nseqs; ++i){
60	seq = diffvars->seqs[i];
61	asc_assert(seq.n >= 1);
62	v = seq.vars[0];
63	if(!var_apply_filter(v,&integrator_ida_nonderiv)){
64	varname = var_make_name(sys->system,v);
65	CONSOLE_DEBUG("'%s' fails non-deriv filter",varname);
66	ASC_FREE(varname);
67	for(j=1;j<seq.n;++j){
68	v = seq.vars[j];
69	var_set_active(v,FALSE);
70	var_set_value(v,0);
71	varname = var_make_name(sys->system,v);
72	CONSOLE_DEBUG("Derivative '%s' SET ZERO AND INACTIVE",varname);
73	ASC_FREE(v);
74	}
75	continue;
76	}
77
78	varname = var_make_name(sys->system,v);
79	CONSOLE_DEBUG("We will use var '%s'",varname);
80	ASC_FREE(varname);
81	if(seq.n > 2){
82	varname = var_make_name(sys->system,v);
83	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'");
84	ASC_FREE(varname);
85	return 2;
86	}else if(seq.n==2){
87	/* diff var */
88	if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
89	/* add the diff & deriv vars to the lists */
90	n_y++;
91	varname = var_make_name(sys->system,v);
92	derivname = var_make_name(sys->system,seq.vars[1]);
93	CONSOLE_DEBUG("Added '%s' and '%s'",varname,derivname);
94	ASC_FREE(varname);
95	ASC_FREE(derivname);
96	continue;
97	}
98	varname = var_make_name(sys->system,v);
99	derivname = var_make_name(sys->system,seq.vars[1]);
100	CONSOLE_DEBUG("Derivative '%s' of '%s' fails filter; convert '%s' to algebraic",derivname,varname,varname);
101	ASC_FREE(varname);
102	ASC_FREE(derivname);
103	/* fall through */
104	}
105
106	varname = var_make_name(sys->system,v);
107	CONSOLE_DEBUG("Adding '%s' to algebraic",varname);
108	ASC_FREE(varname);
109	n_y++;
110	}
111
112	/* we assert that all vars in y meet the integrator_ida_nonderiv filter */
113	/* we assert that all vars in ydot meet the integrator_ida_deriv filter */
114
115	CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
116	sys->n_y = n_y;
117	return 0;
118	}
119
120
121	/**
122	Sort the lists. First we will put the non-derivative vars, then we will
123	put the derivative vars. Then we will put all the others.
124	*/
125	static int integrator_ida_sort_rels_and_vars(IntegratorSystem *sys){
126	int ny1, nydot, nr;
127
128	/* we should not have allocated y or ydot yet */
129	asc_assert(sys->y==NULL && sys->ydot==NULL);
130
131	/* but we should have found some variables (and know how many) */
132	asc_assert(sys->n_y);
133
134	if(system_cut_vars(sys->system, 0, &integrator_ida_nonderiv, &ny1)){
135	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
136	return 1;
137	}
138
139	CONSOLE_DEBUG("cut_vars: ny1 = %d, sys->n_y = %d",ny1,sys->n_y);
140	asc_assert(ny1 == sys->n_y);
141
142	if(system_cut_vars(sys->system, ny1, &integrator_ida_deriv, &nydot)){
143	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
144	return 1;
145	}
146
147	if(system_cut_rels(sys->system, 0, &integrator_ida_rel, &nr)){
148	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
149	return 1;
150	}
151
152	if(ny1 != nr){
153	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
154	return 2;
155	}
156
157	return 0;
158	}
159
160	/**
161	Build up the lists 'ydot' and 'y_id.
162	'ydot' allows us to find the derivative of a variable given its sindex.
163	'y_id' allows us to locate the variable of a derivative given its sindex.
164
165	Hence
166	y_id[var_sindex(derivvar)-sys->n_y] = var_sindex(diffvar)
167	ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)
168
169	Note that there is 'shifting' required when looking up y_id.
170
171	Note that we want to get rid of 'y' but for the moment we are also assigning
172	that. We want to be rid of it because if the vars are ordered correctly
173	it shouldn't be needed.
174	*/
175	static int integrator_ida_create_lists(IntegratorSystem *sys){
176	const SolverDiffVarCollection *diffvars;
177	int i, j;
178	struct var_variable *v;
179
180	SolverDiffVarSequence seq;
181
182	asc_assert(sys->y==NULL);
183	asc_assert(sys->ydot==NULL);
184	asc_assert(sys->y_id== NULL);
185
186	/* get the the dervative chains from the system */
187	diffvars = slv_get_diffvars(sys->system);
188	if(diffvars==NULL){
189	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
190	return 1;
191	}
192
193	/* allocate space (more than enough) */
194	sys->y = ASC_NEW_ARRAY(struct var_variable *,sys->n_y);
195	sys->y_id = ASC_NEW_ARRAY_CLEAR(int,sys->n_y);
196	sys->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,sys->n_y);
197	j = 0;
198
199	/* add the variables from the derivative chains */
200	for(i=0; i<diffvars->nseqs; ++i){
201	seq = diffvars->seqs[i];
202	asc_assert(seq.n >= 1);
203	v = seq.vars[0];
204	if(!var_apply_filter(v,&integrator_ida_nonderiv)){
205	continue;
206	}
207
208	if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
209	asc_assert(var_sindex(seq.vars[1])-sys->n_y >= 0);
210	asc_assert(var_sindex(seq.vars[1])-sys->n_y < sys->n_y);
211	sys->y_id[var_sindex(seq.vars[1])-sys->n_y] = j;
212	sys->ydot[j] = seq.vars[1];
213	}else{
214	asc_assert(sys->ydot[j]==NULL);
215	}
216
217	sys->y[j] = v;
218	j++;
219	}
220
221	return 0;
222	}
223
224	/**
225	Perform additional problem analysis to prepare problem for integration with
226	IDA.
227
228	We assume that the analyse_generate_diffvars routine has been called, so
229	that we just need to call slv_get_diffvars to access the derivative
230	chains.
231
232	We can also assume that the independent variable has been found.
233
234	See mailing list, ~Jan 2007.
235
236	Note, the stuff for identifying the static and output sub-problems should
237	be part of integrator.c, not this file. We will assume this is handled
238
239	@return 0 on success
240	@see integrator_analyse
241	*/
242	int integrator_ida_analyse(IntegratorSystem *sys){
243	int res;
244	const SolverDiffVarCollection *diffvars;
245	char *varname;
246	int i;
247
248	asc_assert(sys->engine==INTEG_IDA);
249
250	#ifdef ANALYSE_DEBUG
251	CONSOLE_DEBUG("Starting IDA analysis");
252	#endif
253
254	res = integrator_ida_check_vars(sys);
255	if(res){
256	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
257	return 1;
258	}
259
260	#ifdef ANALYSE_DEBUG
261	CONSOLE_DEBUG("Sorting rels and vars");
262	#endif
263
264	res = integrator_ida_sort_rels_and_vars(sys);
265	if(res){
266	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
267	return 1;
268	}
269
270	#ifdef ANALYSE_DEBUG
271	CONSOLE_DEBUG("Creating lists");
272	#endif
273
274	res = integrator_ida_create_lists(sys);
275	if(res){
276	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
277	return 1;
278	}
279
280	#ifdef ANALYSE_DEBUG
281	CONSOLE_DEBUG("Checking");
282	#endif
283
284	asc_assert(sys->y);
285	asc_assert(sys->ydot);
286	asc_assert(sys->y_id);
287
288	integrator_ida_analyse_debug(sys,stderr);
289
290	if(integrator_ida_check_diffindex(sys)){
291	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
292	return 360;
293	}
294
295	/* the following causes TestIDA.testincidence3 to fail:
296	if(sys->n_y != slv_get_num_solvers_rels(sys->system)){
297	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
298	return 2;
299	}*/
300
301	#if 0
302	/* check structural singularity for the two IDACalcIC scenarios */
303
304	/* ...(1) FIX the derivatives */
305	CONSOLE_DEBUG("Checking system with derivatives fixed...");
306	for(i=0;i<n_y;++i){
307	if(sys->ydot[i])var_set_fixed(sys->ydot[i],1);
308	}
309
310	slv_block_partition(sys->system);
311	res = integrator_ida_block_check(sys);
312	if(res)return 100 + res;
313
314	/* ...(2) FREE the derivatives, FIX the diffvars */
315	CONSOLE_DEBUG("Checking system with differential variables fixed...");
316	for(i=0;i<n_y;++i){
317	if(sys->ydot[i]){
318	var_set_fixed(sys->ydot[i],0);
319	var_set_fixed(sys->y[i],1);
320	}
321	}
322	slv_block_partition(sys->system);
323	res = integrator_ida_block_check(sys);
324	if(res)return 200 + res;
325
326	/* ...(3) restore as it was, FREE the diffvars */
327	for(i=0;i<n_y;++i){
328	if(sys->ydot[i]){
329	var_set_fixed(sys->y[i],0);
330	}
331	}
332	#endif
333
334	/* get the the dervative chains from the system */
335	diffvars = slv_get_diffvars(sys->system);
336
337	/* check the indep var is same as was located elsewhere */
338	if(diffvars->nindep>1){
339	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
340	return 3;
341	}else if(diffvars->nindep<1){
342	ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
343	return 4;
344	}else if(diffvars->indep[0]!=sys->x){
345	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
346	return 5;
347	}
348
349	/* get the observations */
350	/** @TODO this should use a 'problem provider' API instead of static data
351	from the system object */
352	sys->n_obs = diffvars->nobs;
353	sys->obs = ASC_NEW_ARRAY(struct var_variable *,sys->n_obs);
354	for(i=0;i<sys->n_obs;++i){
355	sys->obs[i] = diffvars->obs[i];
356	varname = var_make_name(sys->system,sys->obs[i]);
357	CONSOLE_DEBUG("'%s' is observation",varname);
358	ASC_FREE(varname);
359	}
360
361	/* - 'y' list as [ya\|yd] */
362	/* - sparsity pattern for dF/dy and dF/dy' */
363	/* - sparsity pattern for union of above */
364	/* - block decomposition based on above */
365	/* - block decomposition results in reordering of y and y' */
366	/* - boundaries (optional) */
367	/* ERROR_REPORTER_HERE(ASC_PROG_ERR,"Implementation incomplete");
368	return -1; */
369
370	return 0;
371	}
372
373	/*------------------------------------------------------------------------------
374	ANALYSIS ROUTINE (new implementation)
375	*/
376
377	static int integrator_ida_check_partitioning(IntegratorSystem *sys){
378	long i, nv;
379	int err=0;
380	char *varname;
381	struct var_variable *vlist, v;
382	nv = slv_get_num_solvers_vars(sys->system);
383	vlist = slv_get_solvers_var_list(sys->system);
384	var_filter_t vf = {VAR_SVAR\|VAR_ACTIVE\|VAR_INCIDENT\|VAR_FIXED
385	,VAR_SVAR\|VAR_ACTIVE\|VAR_INCIDENT\|0 };
386	for(i=0;i<nv;++i){
387	v=vlist[i];
388	if(!var_apply_filter(v,&vf))continue;
389	varname = var_make_name(sys->system,v);
390	if(!var_deriv(v)){
391	fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
392	if(i>=sys->n_y){
393	ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
394	err++;
395	}
396	}else{
397	fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
398	if(i<sys->n_y){
399	ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
400	,varname, sys->n_y, i
401	);
402	err++;
403	}
404	}
405	ASC_FREE(varname);
406	}
407	if(err){
408	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
409	return 1;
410	}
411	return 0;
412	}
413
414	int integrator_ida_block_check(IntegratorSystem *sys){
415	int res;
416	int dof;
417	#ifdef ANALYSE_DEBUG
418	long vlist, vp, i, nv, nv_ok;
419	char *varname;
420	struct var_variable **solversvars;
421
422	var_filter_t vfilt = {
423	VAR_ACTIVE \| VAR_INCIDENT \| VAR_FIXED,
424	VAR_ACTIVE \| VAR_INCIDENT \| 0
425	};
426
427	nv = slv_get_num_solvers_vars(sys->system);
428	solversvars = slv_get_solvers_var_list(sys->system);
429	CONSOLE_DEBUG("-------------- nv = %ld -------------",nv);
430	for(nv_ok=0, i=0; i < nv; ++i){
431	if(var_apply_filter(solversvars[i],&vfilt)){
432	varname = var_make_name(sys->system,solversvars[i]);
433	fprintf(stderr,"%s\n",varname);
434	ASC_FREE(varname);
435	nv_ok++;
436	}
437	}
438	CONSOLE_DEBUG("----------- got %ld ok -------------",nv_ok);
439	#endif
440
441	if(!slvDOF_status(sys->system, &res, &dof)){
442	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
443	return -1;
444	}
445	switch(res){
446	case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
447	case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
448	case 3: CONSOLE_DEBUG("System is structurally singular"); break;
449	case 4: CONSOLE_DEBUG("System is overspecified"); break;
450	default:
451	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
452	return -2;
453	}
454
455	#ifdef ANALYSE_DEBUG
456	/* if it was underspecified, what vars could be fixed? */
457	if(res==1){
458	CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
459	solversvars = slv_get_solvers_var_list(sys->system);
460	if(!slvDOF_eligible(sys->system, &vlist)){
461	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
462	return -3;
463	}
464	for(vp=vlist;*vp!=-1;++vp){
465	varname = var_make_name(sys->system, solversvars[*vp]);
466	CONSOLE_DEBUG("Fixable var: %s",varname);
467	ASC_FREE(varname);
468	}
469	CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
470	return 1;
471	}
472	#endif
473
474	return res;
475	}
476
477	static int integrator_ida_check_diffindex(IntegratorSystem *sys){
478	int i, nv, err = 0;
479	struct var_variable **vlist;
480	int diffindex;
481
482	CONSOLE_DEBUG("Checking diffindex vector");
483
484	if(sys->y_id == NULL){
485	CONSOLE_DEBUG("y_id not allocated");
486	return 1;
487	}
488
489	vlist = slv_get_solvers_var_list(sys->system);
490	nv = slv_get_num_solvers_vars(sys->system);
491
492	for(i=0; i<nv; ++i){
493	if(var_deriv(vlist[i])){
494	if(i < sys->n_y){
495	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv in wrong position");
496	err++;
497	}else{
498	if(var_sindex(vlist[i])!=i){
499	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Deriv var with incorrect sindex");
500	err++;
501	}
502	diffindex = sys->y_id[i - sys->n_y];
503	if(diffindex >= sys->n_y){
504	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff y_id[%d] value too big",i-sys->n_y);
505	err++;
506	}
507	if(var_sindex(vlist[diffindex])!=diffindex){
508	ERROR_REPORTER_HERE(ASC_PROG_ERR,"Diff var with incorrect sindex");
509	err++;
510	}
511	}
512	}else{
513	if(i!=var_sindex(vlist[i])){
514	ERROR_REPORTER_HERE(ASC_PROG_ERR,"var_sindex incorrect");
515	err++;
516	}
517	}
518	}
519	if(err){
520	CONSOLE_DEBUG("Errors found in diffindex");
521	}
522	return err;
523	}
524
525	/**
526	@TODO provide a macro implementation of this
527	*/
528	int integrator_ida_diffindex(const IntegratorSystem sys, const struct var_variable deriv){
529	int diffindex;
530	#ifdef DIFFINDEX_DEBUG
531	asc_assert( var_deriv (deriv));
532	asc_assert( var_active (deriv));
533	asc_assert( var_incident (deriv));
534	asc_assert( var_svar (deriv));
535
536	asc_assert(!var_fixed (deriv));
537
538	asc_assert(var_sindex(deriv) >= sys->n_y);
539	asc_assert(diffindex == var_sindex(sys->y[diffindex]));
540	#endif
541	asc_assert(var_sindex(deriv) >= sys->n_y);
542	diffindex = sys->y_id[var_sindex(deriv) - sys->n_y];
543
544	return diffindex;
545	}
546
547
548	int integrator_ida_analyse_debug(const IntegratorSystem sys,FILE fp){
549	return analyse_diffvars_debug(sys->system,fp);
550	}
551