pygtk/interface/simulation.cpp

#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <sstream>
using namespace std;

extern "C"{
#include <utilities/ascConfig.h>
#include <utilities/ascSignal.h>
#include <utilities/ascMalloc.h>
#include <general/dstring.h>
#include <general/tm_time.h>
#include <compiler/instance_enum.h>
#include <compiler/fractions.h>
#include <compiler/compiler.h>
#include <compiler/dimen.h>
#include <compiler/symtab.h>
#include <compiler/instance_io.h>
#include <compiler/instantiate.h>
#include <compiler/bintoken.h>
#include <compiler/instance_enum.h>
#include <compiler/instquery.h>
#include <compiler/check.h>
#include <compiler/name.h>
#include <compiler/pending.h>

#include <utilities/readln.h>
#include <solver/mtx.h>
#include <solver/slv_types.h>
#include <solver/var.h>
#include <solver/rel.h>
#include <solver/discrete.h>
#include <solver/conditional.h>
#include <solver/logrel.h>
#include <solver/bnd.h>
#include <solver/calc.h>
#include <solver/relman.h>
#include <solver/slv_common.h>
#include <solver/linsol.h>
#include <solver/linsolqr.h>
#include <solver/slv_client.h>
#include <solver/system.h>
#include <solver/slv_interface.h>
#include <solver/slvDOF.h>
#include <solver/slv3.h>
#include <solver/slv_stdcalls.h>
}

#include "simulation.h"
#include "solver.h"
#include "solverparameters.h"
#include "name.h"
#include "incidencematrix.h"
#include "variable.h"

/**
    Create an instance of a type (call compiler etc)

    @TODO fix mutex on compile command filenames
*/
Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
    is_built = false;
    // Create an Instance object for the 'simulation root' (we'll call
    // it the 'simulation model') and it can be fetched using 'getModel()'
    // any time later.
    //simroot = Instanc(GetSimulationRoot(i),name);
}

Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
    is_built = old.is_built;
    sys = old.sys;
    bin_srcname = old.bin_srcname;
    bin_objname = old.bin_objname;
    bin_libname = old.bin_libname;
    bin_cmd = old.bin_cmd;
    bin_rm = old.bin_rm;
}

Simulation::~Simulation(){
    //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
}

Instanc &
Simulation::getModel(){
    if(!simroot.getInternalType()){
        throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
    }
    return simroot;
}

void
Simulation::checkDoF() const{
        cerr << "CHECKING DOF..." << endl;
        int dof, status;
        if(!sys){
                throw runtime_error("System not yet built");
        }
        slvDOF_status(sys, &status, &dof);
        switch(status){
                case 1: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Underspecified; %d degrees of freedom",dof); break;
                case 2: ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Square"); break;
                case 3: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally singular"); break;
                case 4: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Overspecified"); break;
                case 5:
                        throw runtime_error("Unable to resolve degrees of freedom"); break;
                default:
                        throw runtime_error("Invalid return status from slvDOF_status");
        }
}

void
Simulation::run(const Method &method){
    cerr << "RUNNING PROCEDURE " << method.getName() << endl;
    Nam name = Nam(method.getSym());
    //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
    Proc_enum pe;
    pe = Initialize(
        &*(getModel().getInternalType()) ,name.getInternalType(), "__not_named__"
        ,ASCERR
        ,0, NULL, NULL
    );

    if(pe == Proc_all_ok){
        ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
        //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
    }else{
        stringstream ss;
        ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
        switch(pe){
            case Proc_CallOK: ss << "Call OK"; break;
            case Proc_CallError: ss << "Error occurred in call"; break;
            case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
            case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
            case Proc_CallContinue: ss << "Skip to next iteration"; break;

            case Proc_break: ss << "Break"; break;
            case Proc_continue: ss << "Continue"; break;
            case Proc_fallthru: ss << "Fall-through"; break;
            case Proc_return: ss << "Return"; break;
            case Proc_stop: ss << "Stop"; break;
            case Proc_stack_exceeded: ss << "Stack exceeded"; break;
            case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
            case Proc_case_matched: ss << "Case matched"; break;
            case Proc_case_unmatched: ss << "Case unmatched"; break;

            case Proc_case_undefined_value: ss << "Undefined value in case"; break;
            case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
            case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
            case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
            case Proc_case_wrong_index: ss << "Wrong index in case"; break;
            case Proc_case_wrong_value: ss << "Wrong value in case"; break;
            case Proc_case_extra_values: ss << "Extra values in case"; break;
            case Proc_bad_statement: ss << "Bad statement"; break;
            case Proc_bad_name: ss << "Bad name"; break;
            case Proc_for_duplicate_index: ss << "Duplicate index"; break;
            case Proc_for_set_err: ss << "For set error"; break;
            case Proc_for_not_set: ss << "For not set"; break;
            case Proc_illegal_name_use: ss << "Illegal name use"; break;
            case Proc_name_not_found: ss << "Name not found"; break;
            case Proc_instance_not_found: ss << "Instance not found"; break;
            case Proc_type_not_found: ss << "Type not found"; break;
            case Proc_illegal_type_use: ss << "Illegal use"; break;
            case Proc_proc_not_found: ss << "Method not found"; break;
            case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
            case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
            case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
            case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
            case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
            case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
            case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
            case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
            case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
            case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
            case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
            case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
            case Proc_if_not_logical: ss << "If expression is not logical"; break;
            case Proc_user_interrupt: ss << "User interrupt"; break;
            case Proc_infinite_loop: ss << "Infinite loop"; break;
            case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
            case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
            case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
            case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
            case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
            case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
            case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
            case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
            case Proc_lhs_error: ss << "Left-hand-side error"; break;
            case Proc_rhs_error: ss << "Right-hand-side error"; break;
            case Proc_unknown_error: ss << "Unknown error"; break;
            default:
                ss << "Invalid error code";
        }


        ss << " (" << int(pe) << ")";
        throw runtime_error(ss.str());
    }
}

const bool
Simulation::check(){
    cerr << "CHECKING SIMULATION" << endl;
    Instance *i1 = getModel().getInternalType();
    CheckInstance(stderr, &*i1);
    cerr << "...DONE CHECKING" << endl;
}

void
Simulation::build(){
    cerr << "BUILDING SIMULATION..." << endl;
    Instance *i1 = getModel().getInternalType();
    sys = system_build(&*i1);
    if(!sys){
        throw runtime_error("Unable to build system");
    }
    is_built = true;
    cerr << "...DONE BUILDING" << endl;
}

vector<Variable>
Simulation::getFixableVariables(){
    cerr << "GETTING FIXABLE VARIABLES..." << endl;
    vector<Variable> vars;
    vars.reserve(100);

    if(!sys){
        throw runtime_error("Simulation system not yet built");
    }

    int32 **vip; /** TODO ensure 32 bit integers are used */

    // Get IDs of elegible variables in array at vip...
    if(!slvDOF_eligible(sys,vip)){
        ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
    }else{
        //cerr << "FIXABLE VARS FOUND" << endl;
        struct var_variable **vp = slv_get_solvers_var_list(sys);

        /*struct var_variable *first_var = vp[0];
        char *first_var_name = var_make_name(sys,first_var);
        cerr << "FIRST SYS VAR IS NAMED " << var_make_name(s,first_var) << endl;
        ascfree(first_var_name);*/

        if(vp==NULL){
            throw runtime_error("Simulation variable list is null");
        }

        // iterate through this list until we find a -1:
        int i=0;
        int var_index = (*vip)[i];
        while(var_index >= 0){
            //cerr << "FOUND VARIABLE var_index = " << var_index << endl;
            struct var_variable *var = vp[var_index];
            //cerr << "VARIABLE " << var_index << " IS ELIGIBLE" << endl;
            char *var_name = var_make_name(sys,var);
            //cerr << "ELIGIBLE VAR: " << var_name << endl;
            ascfree(var_name);
            vars.push_back( Variable(this, var) );
            ++i;
            var_index = (*vip)[i];
        }
        ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
        //cerr << "END ELEGIBLE VARS LIST" << endl;
        ascfree(*vip);
        //cerr << "FREED VIP LIST" << endl;
    }

    //cerr << "FINISHED WITH FINDING ELEGIBLE VARIABLES" << endl;
    return vars;
}

    
void
Simulation::solve(Solver solver){
    if(!is_built){
        throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
    }

    cerr << "SIMULATION::SOLVE STARTING..." << endl;
    enum inst_t k = getModel().getKind();
    if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");

    Instance *i1 = getInternalType();
    int npend = NumberPendingInstances(&*i1);
    if(npend)throw runtime_error("Can't solve: There are still %d pending instances");

    if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");

    cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
    setSolver(solver);


    cerr << "PRESOLVING SYSTEM...";
    slv_presolve(sys);
    cerr << "DONE" << endl;

    cerr << "SOLVING SYSTEM..." << endl;
    // Add some stuff here for cleverer iteration....
    unsigned niter = 1000;
    double updateinterval = 0.02;

    double starttime = tm_cpu_time();
    double lastupdate = starttime;
    slv_status_t status;
    int solved_vars=0;
    bool stop=false;

    for(int iter = 1; iter <= niter && !stop; ++iter){
        slv_get_status(sys,&status);
        if(status.ready_to_solve){
            slv_iterate(sys);
        }
        if(tm_cpu_time() - lastupdate > updateinterval && iter > 0){
            if(solved_vars < status.block.previous_total_size){
                solved_vars = status.block.previous_total_size;
                CONSOLE_DEBUG("Solved %5d vars. Now block %3d/%-3d...", solved_vars, status.block.current_block, status.block.number_of);
            }else{
                CONSOLE_DEBUG("Total %5d iterations (%5d in current block of %d vars)", iter, status.block.iteration, status.block.current_size);
            }
            lastupdate = tm_cpu_time();
        }
    }
    double elapsed = tm_cpu_time() - starttime;

    if(status.ok){
        cerr << "... DONE SOLVING SYSTEM" << endl;
    }else{
        ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver failed");
    }

    cerr << "SOLVER PERFORMED " << status.iteration << " ITERATIONS IN " << elapsed << "s" << endl;

    if(status.iteration_limit_exceeded){
        ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Exceeded interation limit");
    }

    if(status.converged){
        ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Solver converged: %d iterations (%.2f s)"
            ,status.iteration,elapsed);
    }else{
        ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver not converged after %d iterations (%.2f s).",status.iteration,elapsed);
    }

}

void
Simulation::write(){
    simroot.write();
}

//------------------------------------------
// ASSIGNING SOLVER TO SIMULATION

void
Simulation::setSolver(Solver &solver){
    cerr << "SETTING SOLVER ON SIMULATION TO " << solver.getName() << endl;

    if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
    // Update the solver object because sometimes an alternative solver can be returned, apparently.

    int selected = slv_select_solver(sys, solver.getIndex());
    //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;

    if(selected<0){
        ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
        throw runtime_error("Failed to select solver");
    }

    if(selected!=solver.getIndex()){
        solver = Solver(slv_solver_name(selected));
        ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
    }

    if( slv_eligible_solver(sys) <= 0){
        ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
        throw runtime_error("Inelegible solver");
    }
}

const Solver
Simulation::getSolver() const{
    int index = slv_get_selected_solver(sys);
    //cerr << "Simulation::getSolver: index = " << index << endl;
    if(index<0)throw runtime_error("No solver selected");

    return Solver(slv_solver_name(index));
}


/**
    Get solver parameters struct wrapped up as a SolverParameters class.
*/
SolverParameters
Simulation::getSolverParameters() const{
    if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");

    slv_parameters_t p;
    slv_get_parameters(sys,&p);
    return SolverParameters(p);
}

/**
    Update the solver parameters by passing a new set back
*/
void
Simulation::setSolverParameters(SolverParameters &P){
    if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
    slv_set_parameters(sys, &(P.getInternalType()));
}

slv_system_structure *
Simulation::getSystem(){
    if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
    return sys;
}

IncidenceMatrix
Simulation::getIncidenceMatrix(){
    return IncidenceMatrix(*this);
}
1	#include <iostream>
2	#include <iomanip>
3	#include <stdexcept>
4	#include <sstream>
5	using namespace std;
6
7	extern "C"{
8	#include <utilities/ascConfig.h>
9	#include <utilities/ascSignal.h>
10	#include <utilities/ascMalloc.h>
11	#include <general/dstring.h>
12	#include <general/tm_time.h>
13	#include <compiler/instance_enum.h>
14	#include <compiler/fractions.h>
15	#include <compiler/compiler.h>
16	#include <compiler/dimen.h>
17	#include <compiler/symtab.h>
18	#include <compiler/instance_io.h>
19	#include <compiler/instantiate.h>
20	#include <compiler/bintoken.h>
21	#include <compiler/instance_enum.h>
22	#include <compiler/instquery.h>
23	#include <compiler/check.h>
24	#include <compiler/name.h>
25	#include <compiler/pending.h>
26
27	#include <utilities/readln.h>
28	#include <solver/mtx.h>
29	#include <solver/slv_types.h>
30	#include <solver/var.h>
31	#include <solver/rel.h>
32	#include <solver/discrete.h>
33	#include <solver/conditional.h>
34	#include <solver/logrel.h>
35	#include <solver/bnd.h>
36	#include <solver/calc.h>
37	#include <solver/relman.h>
38	#include <solver/slv_common.h>
39	#include <solver/linsol.h>
40	#include <solver/linsolqr.h>
41	#include <solver/slv_client.h>
42	#include <solver/system.h>
43	#include <solver/slv_interface.h>
44	#include <solver/slvDOF.h>
45	#include <solver/slv3.h>
46	#include <solver/slv_stdcalls.h>
47	}
48
49	#include "simulation.h"
50	#include "solver.h"
51	#include "solverparameters.h"
52	#include "name.h"
53	#include "incidencematrix.h"
54	#include "variable.h"
55
56	/**
57	Create an instance of a type (call compiler etc)
58
59	@TODO fix mutex on compile command filenames
60	*/
61	Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
62	is_built = false;
63	// Create an Instance object for the 'simulation root' (we'll call
64	// it the 'simulation model') and it can be fetched using 'getModel()'
65	// any time later.
66	//simroot = Instanc(GetSimulationRoot(i),name);
67	}
68
69	Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
70	is_built = old.is_built;
71	sys = old.sys;
72	bin_srcname = old.bin_srcname;
73	bin_objname = old.bin_objname;
74	bin_libname = old.bin_libname;
75	bin_cmd = old.bin_cmd;
76	bin_rm = old.bin_rm;
77	}
78
79	Simulation::~Simulation(){
80	//CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
81	}
82
83	Instanc &
84	Simulation::getModel(){
85	if(!simroot.getInternalType()){
86	throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
87	}
88	return simroot;
89	}
90
91	void
92	Simulation::checkDoF() const{
93	cerr << "CHECKING DOF..." << endl;
94	int dof, status;
95	if(!sys){
96	throw runtime_error("System not yet built");
97	}
98	slvDOF_status(sys, &status, &dof);
99	switch(status){
100	case 1: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Underspecified; %d degrees of freedom",dof); break;
101	case 2: ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Square"); break;
102	case 3: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally singular"); break;
103	case 4: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Overspecified"); break;
104	case 5:
105	throw runtime_error("Unable to resolve degrees of freedom"); break;
106	default:
107	throw runtime_error("Invalid return status from slvDOF_status");
108	}
109	}
110
111	void
112	Simulation::run(const Method &method){
113	cerr << "RUNNING PROCEDURE " << method.getName() << endl;
114	Nam name = Nam(method.getSym());
115	//cerr << "CREATED NAME '" << name.getName() << "'" << endl;
116	Proc_enum pe;
117	pe = Initialize(
118	&*(getModel().getInternalType()) ,name.getInternalType(), "__not_named__"
119	,ASCERR
120	,0, NULL, NULL
121	);
122
123	if(pe == Proc_all_ok){
124	ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
125	//cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
126	}else{
127	stringstream ss;
128	ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
129	switch(pe){
130	case Proc_CallOK: ss << "Call OK"; break;
131	case Proc_CallError: ss << "Error occurred in call"; break;
132	case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
133	case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
134	case Proc_CallContinue: ss << "Skip to next iteration"; break;
135
136	case Proc_break: ss << "Break"; break;
137	case Proc_continue: ss << "Continue"; break;
138	case Proc_fallthru: ss << "Fall-through"; break;
139	case Proc_return: ss << "Return"; break;
140	case Proc_stop: ss << "Stop"; break;
141	case Proc_stack_exceeded: ss << "Stack exceeded"; break;
142	case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
143	case Proc_case_matched: ss << "Case matched"; break;
144	case Proc_case_unmatched: ss << "Case unmatched"; break;
145
146	case Proc_case_undefined_value: ss << "Undefined value in case"; break;
147	case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
148	case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
149	case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
150	case Proc_case_wrong_index: ss << "Wrong index in case"; break;
151	case Proc_case_wrong_value: ss << "Wrong value in case"; break;
152	case Proc_case_extra_values: ss << "Extra values in case"; break;
153	case Proc_bad_statement: ss << "Bad statement"; break;
154	case Proc_bad_name: ss << "Bad name"; break;
155	case Proc_for_duplicate_index: ss << "Duplicate index"; break;
156	case Proc_for_set_err: ss << "For set error"; break;
157	case Proc_for_not_set: ss << "For not set"; break;
158	case Proc_illegal_name_use: ss << "Illegal name use"; break;
159	case Proc_name_not_found: ss << "Name not found"; break;
160	case Proc_instance_not_found: ss << "Instance not found"; break;
161	case Proc_type_not_found: ss << "Type not found"; break;
162	case Proc_illegal_type_use: ss << "Illegal use"; break;
163	case Proc_proc_not_found: ss << "Method not found"; break;
164	case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
165	case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
166	case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
167	case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
168	case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
169	case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
170	case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
171	case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
172	case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
173	case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
174	case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
175	case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
176	case Proc_if_not_logical: ss << "If expression is not logical"; break;
177	case Proc_user_interrupt: ss << "User interrupt"; break;
178	case Proc_infinite_loop: ss << "Infinite loop"; break;
179	case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
180	case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
181	case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
182	case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
183	case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
184	case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
185	case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
186	case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
187	case Proc_lhs_error: ss << "Left-hand-side error"; break;
188	case Proc_rhs_error: ss << "Right-hand-side error"; break;
189	case Proc_unknown_error: ss << "Unknown error"; break;
190	default:
191	ss << "Invalid error code";
192	}
193
194
195	ss << " (" << int(pe) << ")";
196	throw runtime_error(ss.str());
197	}
198	}
199
200	const bool
201	Simulation::check(){
202	cerr << "CHECKING SIMULATION" << endl;
203	Instance *i1 = getModel().getInternalType();
204	CheckInstance(stderr, &*i1);
205	cerr << "...DONE CHECKING" << endl;
206	}
207
208	void
209	Simulation::build(){
210	cerr << "BUILDING SIMULATION..." << endl;
211	Instance *i1 = getModel().getInternalType();
212	sys = system_build(&*i1);
213	if(!sys){
214	throw runtime_error("Unable to build system");
215	}
216	is_built = true;
217	cerr << "...DONE BUILDING" << endl;
218	}
219
220	vector<Variable>
221	Simulation::getFixableVariables(){
222	cerr << "GETTING FIXABLE VARIABLES..." << endl;
223	vector<Variable> vars;
224	vars.reserve(100);
225
226	if(!sys){
227	throw runtime_error("Simulation system not yet built");
228	}
229
230	int32 vip; / TODO ensure 32 bit integers are used */
231
232	// Get IDs of elegible variables in array at vip...
233	if(!slvDOF_eligible(sys,vip)){
234	ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
235	}else{
236	//cerr << "FIXABLE VARS FOUND" << endl;
237	struct var_variable **vp = slv_get_solvers_var_list(sys);
238
239	/struct var_variable first_var = vp[0];
240	char *first_var_name = var_make_name(sys,first_var);
241	cerr << "FIRST SYS VAR IS NAMED " << var_make_name(s,first_var) << endl;
242	ascfree(first_var_name);*/
243
244	if(vp==NULL){
245	throw runtime_error("Simulation variable list is null");
246	}
247
248	// iterate through this list until we find a -1:
249	int i=0;
250	int var_index = (*vip)[i];
251	while(var_index >= 0){
252	//cerr << "FOUND VARIABLE var_index = " << var_index << endl;
253	struct var_variable *var = vp[var_index];
254	//cerr << "VARIABLE " << var_index << " IS ELIGIBLE" << endl;
255	char *var_name = var_make_name(sys,var);
256	//cerr << "ELIGIBLE VAR: " << var_name << endl;
257	ascfree(var_name);
258	vars.push_back( Variable(this, var) );
259	++i;
260	var_index = (*vip)[i];
261	}
262	ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
263	//cerr << "END ELEGIBLE VARS LIST" << endl;
264	ascfree(*vip);
265	//cerr << "FREED VIP LIST" << endl;
266	}
267
268	//cerr << "FINISHED WITH FINDING ELEGIBLE VARIABLES" << endl;
269	return vars;
270	}
271
272
273	void
274	Simulation::solve(Solver solver){
275	if(!is_built){
276	throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
277	}
278
279	cerr << "SIMULATION::SOLVE STARTING..." << endl;
280	enum inst_t k = getModel().getKind();
281	if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");
282
283	Instance *i1 = getInternalType();
284	int npend = NumberPendingInstances(&*i1);
285	if(npend)throw runtime_error("Can't solve: There are still %d pending instances");
286
287	if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
288
289	cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
290	setSolver(solver);
291
292
293	cerr << "PRESOLVING SYSTEM...";
294	slv_presolve(sys);
295	cerr << "DONE" << endl;
296
297	cerr << "SOLVING SYSTEM..." << endl;
298	// Add some stuff here for cleverer iteration....
299	unsigned niter = 1000;
300	double updateinterval = 0.02;
301
302	double starttime = tm_cpu_time();
303	double lastupdate = starttime;
304	slv_status_t status;
305	int solved_vars=0;
306	bool stop=false;
307
308	for(int iter = 1; iter <= niter && !stop; ++iter){
309	slv_get_status(sys,&status);
310	if(status.ready_to_solve){
311	slv_iterate(sys);
312	}
313	if(tm_cpu_time() - lastupdate > updateinterval && iter > 0){
314	if(solved_vars < status.block.previous_total_size){
315	solved_vars = status.block.previous_total_size;
316	CONSOLE_DEBUG("Solved %5d vars. Now block %3d/%-3d...", solved_vars, status.block.current_block, status.block.number_of);
317	}else{
318	CONSOLE_DEBUG("Total %5d iterations (%5d in current block of %d vars)", iter, status.block.iteration, status.block.current_size);
319	}
320	lastupdate = tm_cpu_time();
321	}
322	}
323	double elapsed = tm_cpu_time() - starttime;
324
325	if(status.ok){
326	cerr << "... DONE SOLVING SYSTEM" << endl;
327	}else{
328	ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver failed");
329	}
330
331	cerr << "SOLVER PERFORMED " << status.iteration << " ITERATIONS IN " << elapsed << "s" << endl;
332
333	if(status.iteration_limit_exceeded){
334	ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Exceeded interation limit");
335	}
336
337	if(status.converged){
338	ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Solver converged: %d iterations (%.2f s)"
339	,status.iteration,elapsed);
340	}else{
341	ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver not converged after %d iterations (%.2f s).",status.iteration,elapsed);
342	}
343
344	}
345
346	void
347	Simulation::write(){
348	simroot.write();
349	}
350
351	//------------------------------------------
352	// ASSIGNING SOLVER TO SIMULATION
353
354	void
355	Simulation::setSolver(Solver &solver){
356	cerr << "SETTING SOLVER ON SIMULATION TO " << solver.getName() << endl;
357
358	if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
359	// Update the solver object because sometimes an alternative solver can be returned, apparently.
360
361	int selected = slv_select_solver(sys, solver.getIndex());
362	//cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
363
364	if(selected<0){
365	ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
366	throw runtime_error("Failed to select solver");
367	}
368
369	if(selected!=solver.getIndex()){
370	solver = Solver(slv_solver_name(selected));
371	ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
372	}
373
374	if( slv_eligible_solver(sys) <= 0){
375	ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
376	throw runtime_error("Inelegible solver");
377	}
378	}
379
380	const Solver
381	Simulation::getSolver() const{
382	int index = slv_get_selected_solver(sys);
383	//cerr << "Simulation::getSolver: index = " << index << endl;
384	if(index<0)throw runtime_error("No solver selected");
385
386	return Solver(slv_solver_name(index));
387	}
388
389
390	/**
391	Get solver parameters struct wrapped up as a SolverParameters class.
392	*/
393	SolverParameters
394	Simulation::getSolverParameters() const{
395	if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
396
397	slv_parameters_t p;
398	slv_get_parameters(sys,&p);
399	return SolverParameters(p);
400	}
401
402	/**
403	Update the solver parameters by passing a new set back
404	*/
405	void
406	Simulation::setSolverParameters(SolverParameters &P){
407	if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
408	slv_set_parameters(sys, &(P.getInternalType()));
409	}
410
411	slv_system_structure *
412	Simulation::getSystem(){
413	if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
414	return sys;
415	}
416
417	IncidenceMatrix
418	Simulation::getIncidenceMatrix(){
419	return IncidenceMatrix(*this);
420	}