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/system_impl.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);

/**
    A var can be non-incident. If it *is* non incident and we're going to
    keep it, it will have to have derivative that *is* incident, and that
    meets the following filter.

    If it doesn't have a valid derivative (eg the derivative is fixed, or
    the variable doesn't HAVE a derivative), we will mark the non-deriv
    var non-ACTIVE, so anyway it will end up meeting this filter after we've
    run integrator_ida_check_vars.
*/
const var_filter_t integrator_ida_nonderiv = {
    VAR_SVAR | VAR_ACTIVE | VAR_FIXED | VAR_DERIV,
    VAR_SVAR | 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;
    int n_y = 0;
    int i, j;
    struct var_variable *v;
    SolverDiffVarSequence seq;
    int vok;

#define VARMSG(MSG) \
    varname = var_make_name(sys->system,v); \
    CONSOLE_DEBUG(MSG,varname); \
    ASC_FREE(varname)

    CONSOLE_DEBUG("BEFORE CHECKING VARS");
    integrator_ida_analyse_debug(sys,stderr);

    /* 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 = sys->system->diffvars;
    if(diffvars==NULL){
        ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
        return 1;
    }

    system_var_list_debug(sys->system);
    
    /* 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];

        vok = var_apply_filter(v,&integrator_ida_nonderiv);

        if(vok && !var_incident(v)){
            VARMSG("good var '%s' is not incident");
            /* var meets our filter, but perhaps it's not incident? */
            if(seq.n == 1 || var_apply_filter(seq.vars[1],&integrator_ida_nonderiv)){
                VARMSG("DEACTIVATING NON-INCIDENT VAR '%s' (NO DERIVATIVE)");
                var_set_active(v,0);
                vok = 0;
            }else{
                VARMSG("'%s' has a derivative that's OK");
                ERROR_REPORTER_HERE(ASC_USER_ERROR,"Non-incident var with an incident derivative. ASCEND can't handle this case at the moment, but we hope to fix it.");
                return 1;
            }
        }       

        if(!vok){
            VARMSG("'%s' fails non-deriv filter");
            for(j=1;j<seq.n;++j){
                v = seq.vars[j];
                var_set_active(v,FALSE);
                var_set_value(v,0);
                VARMSG("Derivative '%s' SET ZERO AND INACTIVE");
            }
            continue;
        }       

        VARMSG("We will use var '%s'");
        if(seq.n > 2){
            varname = var_make_name(sys->system,v);             
            ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'",varname);
            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++;
                VARMSG("Added diff var '%s'");
                v = seq.vars[1]; VARMSG("and its derivative '%s'");
                continue;
            }
            VARMSG("Diff var '%s' being converted to alg var...");
            v = seq.vars[1]; VARMSG("...because deriv var '%s' fails filter");
            /* fall through */
        }

        VARMSG("Adding '%s' to algebraic");
        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;


    CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS");
    integrator_ida_analyse_debug(sys,stderr);

    /* 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;
    char *varname;

    SolverDiffVarSequence seq;

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

    CONSOLE_DEBUG("BEFORE MAKING LISTS");
    integrator_ida_analyse_debug(sys,stderr);

    /* get the the dervative chains from the system */
    diffvars = sys->system->diffvars;
    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;

    CONSOLE_DEBUG("Passing through chains...");

    /* add the variables from the derivative chains */
    for(i=0; i<diffvars->nseqs; ++i){

        CONSOLE_DEBUG("i = %d",i);
            
        seq = diffvars->seqs[i];
        asc_assert(seq.n >= 1);
        v = seq.vars[0];

        varname = var_make_name(sys->system, v);
        CONSOLE_DEBUG("alg '%s'",varname);
        ASC_FREE(varname);


        if(!var_apply_filter(v,&integrator_ida_nonderiv)){
            continue;
        }

        varname = var_make_name(sys->system, v);
        CONSOLE_DEBUG("alg '%s' is GOOD",varname);
        ASC_FREE(varname);

        if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
            asc_assert(var_sindex(seq.vars[1]) >= sys->n_y);
            asc_assert(var_sindex(seq.vars[1])-sys->n_y < sys->n_y);

            varname = var_make_name(sys->system, seq.vars[1]);
            CONSOLE_DEBUG("diff '%s' IS GOOD",varname);
            ASC_FREE(varname);

            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++;
    }

    CONSOLE_DEBUG("AFTER MAKING LISTS");
    integrator_ida_analyse_debug(sys,stderr);

    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 = sys->system->diffvars;

    /* 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
    int *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 = %d -------------",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 %d 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;
}


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