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;

    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 = 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){


        CONSOLE_DEBUG("BEFORE CHAIN %d",i);
        integrator_ida_analyse_debug(sys,stderr);

        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(varname);
            }
            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'",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++;
                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;

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

    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 = 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;

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


/**
    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 analyse_diffvars_debug(sys->system,fp);
}

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