#!/usr/bin/env python
#	ASCEND modelling environment
#	Copyright (C) 2006, 2007 Carnegie Mellon University
#
#	This program is free software; you can redistribute it and/or modify
#	it under the terms of the GNU General Public License as published by
#	the Free Software Foundation; either version 2, or (at your option)
#	any later version.
#
#	This program is distributed in the hope that it will be useful,
#	but WITHOUT ANY WARRANTY; without even the implied warranty of
#	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#	GNU General Public License for more details.
#
#	You should have received a copy of the GNU General Public License
#	along with this program; if not, write to the Free Software
#	Foundation, Inc., 59 Temple Place - Suite 330,
#	Boston, MA 02111-1307, USA.

# This script gives a test suite for the high-level interface of ASCEND via
# Python. It is also planned to be a wrapper for the CUnit test suite, although
# this is still experimental.

import unittest
import os, sys
import math
import atexit

import platform
if platform.system() != "Windows":
	try:
		import dl
		_dlflags = dl.RTLD_GLOBAL|dl.RTLD_NOW
	except:
		# On platforms that unilaterally refuse to provide the 'dl' module
		# we'll just set the value and see if it works.
		print "Note: python 'dl' module not available on this system, guessing value of RTLD_* flags"
		_dlflags = 258

	sys.setdlopenflags(_dlflags)

class Ascend(unittest.TestCase):

	def setUp(self):
		import ascpy
		self.L = ascpy.Library()
	
	def tearDown(self):
		self.L.clear()
		del self.L

class AscendSelfTester(Ascend):

	def _run(self,modelname,solvername="QRSlv",filename=None,parameters={}):
		if filename==None:
			filename = 'johnpye/%s.a4c' % modelname
		self.L.load(filename)
		T = self.L.findType(modelname)
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver(solvername))
		for k,v in parameters.iteritems():
			M.setParameter(k,v)
		M.solve(ascpy.Solver(solvername),ascpy.SolverReporter())	
		M.run(T.getMethod('self_test'))
		return M

class TestCompiler(Ascend):

	def _run(self,filen,modeln=""):
		self.L.load('test/compiler/%s.a4c' % filen)
		T = self.L.findType('%s%s' % (filen,modeln))
		M = T.getSimulation('sim')
		M.build()

	def _runfail(self,filen,n,msg="failed"):
		try:
			self._run(filen,'fail%d' % n)
		except Exception,e:
			print "(EXPECTED) ERROR: %s" % e
			return
		self.fail(msg)


	def testloading(self):
		"""library startup"""
		pass

	def testsystema4l(self):
		"""loading system.a4l?"""
		self.L.load('system.a4l')

	def testatomsa4l(self):
		"""loading atoms.a4l?"""
		self.L.load('atoms.a4l')

	def testmissingreq(self):
		"""flagging a missing REQUIRE"""
		self._runfail('missingreq',1)

	def testmissingsubreq(self):
		"""flagging a subsidiary missing REQUIRE"""
		self._runfail('missingreq',1)

	def defaultmethodstest(self,modelname):
		self.L.load("test/defaultmethods.a4c")
		T = self.L.findType(modelname)
		M = T.getSimulation('sim')
		M.run(T.getMethod('on_load'))
		M.run(T.getMethod('self_test'))
		return M

	def testdefault1(self):
		self.defaultmethodstest('testdefault1')

	def testdefault2(self):
		self.defaultmethodstest('testdefault2')

	def testdefault3(self):
		self.defaultmethodstest('testdefault3')

	def testdefault4(self):
		self.defaultmethodstest('testdefault4')

	def testdefault5(self):
		self.defaultmethodstest('testdefault5')

	def testdefault6(self):
		self.defaultmethodstest('testdefault6')

	def testdefault7(self):
		self.defaultmethodstest('testdefault7')

	def testdefault8(self):
		self.defaultmethodstest('testdefault8')

	def testdefault9(self):
		self.defaultmethodstest('testdefault9')

	def testdefault10(self):
		self.defaultmethodstest('testdefault10')

	def testdefault11(self):
		self.defaultmethodstest('testdefault11')

	def testdefault12(self):
		self.defaultmethodstest('testdefault12')

	def testdefault13(self):
		self.defaultmethodstest('testdefault13')

	def testdefault14(self):
		self.defaultmethodstest('testdefault14')

	def testdefault15(self):
		self.defaultmethodstest('testdefault15')

	def testdefault16(self):
		self.defaultmethodstest('testdefault16')

	def testdefault17(self):
		self.defaultmethodstest('testdefault17')
	def testdefault18(self):
		self.defaultmethodstest('testdefault18')

	def testdefault19(self):
		self.defaultmethodstest('testdefault19')

	#def testdefault19fail(self):
	#	self.defaultmethodstest('testdefault19fail')

	def testdefault20(self):
		self.defaultmethodstest('testdefault20')

	#def testdefault20fail(self):
	#	self.defaultmethodstest('testdefault20fail')

class TestSolver(AscendSelfTester):
	
	def testlog10(self):
		M = self._run('testlog10')

	def testrootsofpoly(self):
		self._run('roots_of_poly',filename="roots_of_poly.a4c")

	def testcollapsingcan(self):
		self._run('collapsingcan',filename="collapsingcan.a4c")

	def testdistancecalc(self):
		self._run('distance_calc',filename="distance_calc.a4c")

	def testconopt(self):
		self._run('conopttest',"CONOPT",filename="conopttest.a4c")				

	def testcmslv2(self):
		self._run('testcmslv2',"CMSlv")	

	def testsunpos1(self):
		self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")

	def testsunpos2(self):
		self._run('example_1_6_2',"QRSlv","johnpye/sunpos.a4c")

	def testsunpos3(self):
		self._run('example_1_7_1',"QRSlv","johnpye/sunpos.a4c")

	def testsunpos4(self):
		self._run('example_1_7_2',"QRSlv","johnpye/sunpos.a4c")

	def testsunpos5(self):
		self._run('example_1_7_3',"QRSlv","johnpye/sunpos.a4c")

	def testsunpos6(self):
		self._run('example_1_8_1',"QRSlv","johnpye/sunpos.a4c")

	def testinstanceas(self):
		M = self._run('example_1_6_1',"QRSlv","johnpye/sunpos.a4c")
		self.assertAlmostEqual( float(M.t_solar), M.t_solar.as("s"))
		self.assertAlmostEqual( float(M.t_solar)/3600, M.t_solar.as("h"))

	def testwritegraph(self):
		M = self._run('testlog10')
		if platform.system!="Windows":
			M.write(sys.stderr,"dot")		
		else:
			self.fail("not implemented on windows")

	def testrelinclude(self):
		self.L.load('test/relinclude.a4c')
		T = self.L.findType('relinclude')
		M = T.getSimulation('sim')
		M.eq1.setIncluded(True)
		M.eq2.setIncluded(False)
		M.eq3.setIncluded(False)
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		self.assertAlmostEqual( float(M.z), 2.0)
		M.eq1.setIncluded(False)
		M.eq2.setIncluded(True)
		M.eq3.setIncluded(False)
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		self.assertAlmostEqual( float(M.z), 4.0)
		M.eq1.setIncluded(False)
		M.eq2.setIncluded(False)
		M.eq3.setIncluded(True)
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		self.assertAlmostEqual( float(M.z), 4.61043629206)


class TestLRSlv(AscendSelfTester):
	def testonerel(self):
		self._run('onerel',"LRSlv","test/lrslv/onerel.a4c")

# need to migrate to 'FIX boolvar', currently not supported...
#	def testonerel(self):
#		self._run('onerel',"LRSlv","test/lrslv/onerel.a4c")

	def testsequencecrash(self):
		try:
			self._run('onerel',"LRSlv","test/lrslv/sequencecrash.a4c")
		except:
			pass
			# it just has to not crash, that's all

	def testsequence(self):
		self._run('onerel',"LRSlv","test/lrslv/sequence.a4c")


class TestCMSlv(AscendSelfTester):
	def testsonic(self):
		M = self._run('sonic',"CMSlv","sonic.a4c")
		assert(M.sonic_flow.getBoolValue())

		# other side of boundary...
		M.D.setRealValueWithUnits(4.,"cm")	
		T = self.L.findType('sonic')
		M.solve(ascpy.Solver('CMSlv'),ascpy.SolverReporter())
		M.run(T.getMethod('self_test'))
		assert(not M.sonic_flow.getBoolValue())

	def testheatex(self):
		self._run('heatex',"CMSlv","heatex.a4c")
	def testphaseeq(self):
		self._run('phaseq',"CMSlv","phaseq.a4c")
	def testpipeline(self):
		self._run('pipeline',"CMSlv","pipeline.a4c"
			,{'infinity':3.2e9}
		)
	def testrachford(self):
		self._run('rachford',"CMSlv","rachford.a4c")
	def testlinmassbal(self):
		self._run('linmassbal',"CMSlv","linmassbal.a4c")


class TestMatrix(AscendSelfTester):
	def testlog10(self):
		M = self._run('testlog10')
		print "FETCHING MATRIX................."
		X = M.getMatrix()
# this stuff crashes Windows because the FILE* structure used by Python is not the same
# as used by MinGW...
		#print "GOT MATRIX"
		#sys.stderr.flush()
		#sys.stdout.flush()
		#F = os.tmpfile()
		#X.write(F.fileno,"mmio")
		#F.seek(0)
		#print F.read()
		
class TestIntegrator(Ascend):

	def testListIntegrators(self):
		I = ascpy.Integrator.getEngines()
		s1 = sorted([str(i) for i in I])
		assert 'IDA' in s1
		assert 'LSODE' in s1

	# this routine is reused by both testIDA and testLSODE
	def _testIntegrator(self,integratorname):
		self.L.load('johnpye/shm.a4c')
		M = self.L.findType('shm').getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		P = M.getParameters()
		M.setParameter('feastol',1e-12)
		print M.getChildren()
		assert float(M.x) == 10.0
		assert float(M.v) == 0.0
		t_end = math.pi

		I = ascpy.Integrator(M)
		I.setReporter(ascpy.IntegratorReporterNull(I))
		I.setEngine(integratorname);
		I.setLinearTimesteps(ascpy.Units("s"), 0.0, t_end, 100);
		I.setMinSubStep(0.0001); # these limits are required by IDA at present (numeric diff)
		I.setMaxSubStep(0.1);
		I.setInitialSubStep(0.001);
		I.setMaxSubSteps(200);
		if(integratorname=='IDA'):
			I.setParameter('autodiff',False)
		for p in M.getParameters():
			print p.getName(),"=",p.getValue()
		I.analyse();
		I.solve();
		print "At end of simulation,"
		print "x = %f" % M.x
		print "v = %f" % M.v
		assert abs(float(M.x) + 10) < 1e-2
		assert abs(float(M.v)) < 1e-2
		assert I.getNumObservedVars() == 3

	def testInvalidIntegrator(self):
		self.L.load('johnpye/shm.a4c') 
		M = self.L.findType('shm').getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		try:
			I.setEngine('___NONEXISTENT____')
		except IndexError:
			return
		self.fail("setEngine did not raise error!")

	def testLSODE(self):
		self._testIntegrator('LSODE')

	def testIDA(self):
		self._testIntegrator('IDA')

	def testparameters(self):
		self.L.load('johnpye/shm.a4c')
		M = self.L.findType('shm').getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		P = I.getParameters()
		for p in P:
			print p.getName(),"=",p.getValue()
		assert len(P)==12
		assert P[0].isStr()
		assert P[0].getName()=="linsolver"
		assert P[0].getValue()=='DENSE'
		assert P[2].getName()=="maxord"
		assert P[3].getName()=="autodiff"
		assert P[3].getValue()==True
		assert P[8].getName()=="atolvect"
		assert P[8].getBoolValue() == True
		P[3].setBoolValue(False)
		assert P[3].getBoolValue()==False
		I.setParameters(P)
		assert I.getParameterValue('autodiff')==False
		I.setParameter('autodiff',True)
		try:
			v = I.getParameterValue('nonexist')
		except KeyError:
			pass
		else:
			self.fail('Failed to trip invalid Integrator parameter')

class TestLSODE(Ascend):

	def testzill(self):
		self.L.load('johnpye/zill.a4c')
		T = self.L.findType('zill')
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('LSODE')
		I.setMinSubStep(1e-7)
		I.setMaxSubStep(0.001)
		I.setMaxSubSteps(10000)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
		I.analyse()
		I.solve()
		M.run(T.getMethod('self_test'))

	def testnewton(self):
		sys.stderr.write("STARTING TESTNEWTON\n")
		self.L.load('johnpye/newton.a4c')
		T = self.L.findType('newton')
		M = T.getSimulation('sim')
		M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())	
		I = ascpy.Integrator(M)
		I.setEngine('LSODE')
		I.setParameter('rtolvect',False)
		I.setParameter('rtol',1e-7)
		I.setParameter('atolvect',False)
		I.setParameter('atol',1e-7)
		I.setMinSubStep(1e-7)
		I.setMaxSubStep(0.001)
		I.setMaxSubSteps(10000)
		
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
		I.analyse()
		I.solve()
		print "At end of simulation,"
		print "x = %f" % M.x
		print "v = %f" % M.v
		M.run(T.getMethod('self_test'))

	def testlotka(self):
		self.L.load('johnpye/lotka.a4c')
		M = self.L.findType('lotka').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		I = ascpy.Integrator(M)
		I.setEngine('LSODE')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
		I.analyse()
		print "Number of vars = %d" % I.getNumVars()
		assert I.getNumVars()==2
		I.solve()
		assert I.getNumObservedVars() == 3;
		assert abs(M.R - 832) < 1.0
		assert abs(M.F - 21.36) < 0.1

	def testwritegraph(self):
		self.L.load('johnpye/lotka.a4c')
		M = self.L.findType('lotka').getSimulation('sim')
		F = file('lotka.dot','w')
		M.build()
		M.write(F,"dot")		


#-------------------------------------------------------------------------------
# Testing of a external blackbox functions

class TestBlackBox(AscendSelfTester):
	def testparsefail0(self):
		try:
			self.L.load('test/blackbox/parsefail0.a4c')
			self.fail("parsefail0 should not have loaded without errors")
		except:
			pass

	def testparsefail1(self):
		try:
			self.L.load('test/blackbox/parsefail1.a4c')
			self.fail("parsefail1 should not have loaded without errors")
		except:
			pass

	def testparsefail2(self):
		try:
			self.L.load('test/blackbox/parsefail2.a4c')
			self.fail("parsefail2 should not have loaded without errors")
		except:
			pass

	def testparsefail3(self):
		try:
			self.L.load('test/blackbox/parsefail3.a4c')
			self.fail("parsefail3 should not have loaded without errors")
		except:
			pass

	def testparsefail4(self):
		try:
			self.L.load('test/blackbox/parsefail4.a4c')
			self.fail("parsefail4 should not have loaded")
		except:
			pass

	def testfail1(self):
		"""Mismatched arg counts check-- tests bbox, not ascend."""
		self.L.load('test/blackbox/fail1.a4c')
		try:
			M = self.L.findType('fail1').getSimulation('sim')
			self.fail("expected exception was not raised")
		except RuntimeError,e:
			print "Caught exception '%s', assumed ok" % e

	def testfail2(self):
		"""Incorrect data arg check -- tests bbox, not ascend"""
		self.L.load('test/blackbox/fail2.a4c')
		try:
			M = self.L.findType('fail2').getSimulation('sim')
			self.fail("expected exception was not raised")
		except RuntimeError,e:
			print "Caught exception '%s', assumed ok (should mention errors during instantiation)" % e

	def testpass1(self):
		"""simple single bbox forward solve"""
		M = self._run('pass1',filename='test/blackbox/pass.a4c')

	def testpass2(self):
		"""simple single bbox reverse solve"""
		M = self._run('pass2',filename='test/blackbox/pass.a4c')

	def testpass3(self):
		"""simple double bbox solve"""
		M = self._run('pass3',filename='test/blackbox/pass3.a4c')

	def testpass4(self):
		"""simple double bbox reverse solve"""
		M = self._run('pass4',filename='test/blackbox/pass3.a4c')

	def testpass5(self):
		M = self._run('pass5',filename='test/blackbox/pass5.a4c')

	def testpass6(self):
		M = self._run('pass6',filename='test/blackbox/pass5.a4c')

	def testpass7(self):
		M = self._run('pass7',filename='test/blackbox/passmerge.a4c')

	def testpass8(self):
		M = self._run('pass8',filename='test/blackbox/passmerge.a4c')

	def testpass9(self):
		M = self._run('pass9',filename='test/blackbox/passmerge.a4c')

	def testpass10(self):
		M = self._run('pass10',filename='test/blackbox/passmerge.a4c')

	def testpass11(self):
		M = self._run('pass11',filename='test/blackbox/passmerge.a4c')

	def testpass12(self):
		M = self._run('pass12',filename='test/blackbox/passmerge.a4c')

# this test doesn't work: 'system is inconsistent' -- and structurally singular
#	def testpass13(self):
#		"""cross-merged input/output solve"""
#		M = self._run('pass13',filename='test/blackbox/passmerge.a4c')

	def testpass14(self):
		"""cross-merged input/output reverse solve"""
		M = self._run('pass14',filename='test/blackbox/passmerge.a4c')

	def testpass20(self):
		M = self._run('pass20',filename='test/blackbox/passarray.a4c')

	def testparsefail21(self):
		"""dense array of black boxes wrong syntax"""
		try:
			self.L.load('test/blackbox/parsefail21.a4c')
			self.fail("parsefail21 should not have loaded without errors")
		except:
			pass

	def testpass22(self):
		M = self._run('pass22',filename='test/blackbox/passarray.a4c')

	def testpass23(self):
		M = self._run('pass23',filename='test/blackbox/passarray.a4c')

	def testpass61(self):
		M = self._run('pass61',filename='test/blackbox/reinstantiate.a4c')

	def testpass62(self):
		M = self._run('pass62',filename='test/blackbox/reinstantiate.a4c')

	def testpass64(self):
		M = self._run('pass64',filename='test/blackbox/reinstantiate.a4c')

	def testpass65(self):
		M = self._run('pass65',filename='test/blackbox/reinstantiate.a4c')

	def testpass66(self):
		M = self._run('pass66',filename='test/blackbox/reinstantiate.a4c')

	def testpass67(self):
		M = self._run('pass67',filename='test/blackbox/reinstantiate.a4c')

class TestExtFn(AscendSelfTester):
	def testextfntest(self):
		M = self._run('extfntest',filename='johnpye/extfn/extfntest.a4c')
		self.assertAlmostEqual(M.y, 2);
		self.assertAlmostEqual(M.x, 1);
		self.assertAlmostEqual(M.y, M.x + 1);

	def testextrelfor(self):
		M = self._run('extrelfor',filename='johnpye/extfn/extrelfor.a4c')

## @TODO fix bug with badly-named bbox rel in a loop (Ben, maybe)
#	def testextrelforbadnaming(self):
#		self.L.load('johnpye/extfn/extrelforbadnaming.a4c')
#		T = self.L.findType('extrelfor')
#		M = T.getSimulation('sim')
#		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
#		print "x[1] = %f" % M.x[1]
#		print "x[2] = %f" % M.x[2]
#		print "x[3] = %f" % M.x[3]
#		print "x[4] = %f" % M.x[4]
#		print "x[5] = %f" % M.x[5]
#		M.run(T.getMethod('self_test'))

	def testextrelrepeat(self):
		M = self._run('extrelrepeat',filename='johnpye/extfn/extrelrepeat.a4c')

#-------------------------------------------------------------------------------
# Testing of Sensitivity module

class TestSensitivity(AscendSelfTester):
	def test1(self):
		self.L.load('sensitivity_test.a4c')
		T = self.L.findType('sensitivity_test')
		M = T.getSimulation('sim',False)
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		M.run(T.getMethod('analyse'))
		M.run(T.getMethod('self_test'))

#	def testall(self):
#		self.L.load('sensitivity_test.a4c')
#		T = self.L.findType('sensitivity_test_all')
#		M = T.getSimulation('sim',False)
#		M.run(T.getMethod('on_load'))
#		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
#		M.run(T.getMethod('analyse'))
#		M.run(T.getMethod('self_test'))
# CAUSES CRASH
				
#-------------------------------------------------------------------------------
# Testing of a ExtPy - external python methods

class TestExtPy(AscendSelfTester):
	def test1(self):
		self.L.load('johnpye/extpy/extpytest.a4c')
		T = self.L.findType('extpytest')
		M = T.getSimulation('sim')
		M.run(T.getMethod('self_test'))
		
	def test2(self):
		self.L.load('johnpye/extpy/extpytest.a4c')
		T = self.L.findType('extpytest')
		M = T.getSimulation('sim')
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		M.run(T.getMethod('pythonthing'))
		# causes crash!

#-------------------------------------------------------------------------------
# Testing of saturated steam properties library (iapwssatprops.a4c)

class TestSteam(AscendSelfTester):

	def testiapwssatprops1(self):
		M = self._run('testiapwssatprops1',filename='steam/iapwssatprops.a4c')
	def testiapwssatprops2(self):
		M = self._run('testiapwssatprops2',filename='steam/iapwssatprops.a4c')
	def testiapwssatprops3(self):
		M = self._run('testiapwssatprops3',filename='steam/iapwssatprops.a4c')

	# test the stream model basically works
	def testsatsteamstream(self):
		M = self._run('satsteamstream',filename='steam/satsteamstream.a4c')

	# test that we can solve in terms of various (rho,u)
	def testsatuv(self):
		self.L.load('steam/iapwssat.a4c')
		T = self.L.findType('testiapwssatuv')
		M = T.getSimulation('sim',False)
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		print "p = %f bar" % M.p.as('bar');
		print "T = %f C" % (M.T.as('K') - 273.15);
		print "x = %f" % M.x;
		M.run(T.getMethod('self_test'))
		M.run(T.getMethod('values2'))
#		M.v.setRealValueWithUnits(1.0/450,"m^3/kg");
#		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		print "p = %f bar" % M.p.as('bar');
		print "T = %f C" % (M.T.as('K') - 273.15);
		print "x = %f" % M.x;
		M.run(T.getMethod('self_test2'))
		

## @TODO fix error capture from bounds checking during initialisation
#	def testiapwssat1(self):
#		M = self._run('testiapwssat1',filename='steam/iapwssat.a4c')

	def testdsgsat(self):
		self.L.load('steam/dsgsat3.a4c')
		T = self.L.findType('dsgsat3')
		M = T.getSimulation('sim',False)
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		self.assertAlmostEqual(M.dTw_dt[2],0.0);
		M.run(T.getMethod('configure_dynamic'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		return M

	def testdsgsatrepeat(self):
		self.L.load('steam/dsgsat3.a4c')
		T = self.L.findType('dsgsat3')
		M = T.getSimulation('sim',False)
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())

	def testvary(self):
		self.L.load('steam/dsgsat3.a4c')
		T = self.L.findType('dsgsat3')
		M = T.getSimulation('sim',False)
		M.run(T.getMethod('on_load'))
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		print "----- setting qdot_s -----"
		M.qdot_s.setRealValueWithUnits(1000,"W/m")
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		print "----- setting qdot_s -----"
		M.qdot_s.setRealValueWithUnits(2000,"W/m")
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())

	def teststeadylsode(self):
		"test that steady conditions are stable with LSODE"
		M = self.testdsgsat()
		#M.qdot_s.setRealValueWithUnits(1000,"W/m")
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		#M.setParameter('
	 	I = ascpy.Integrator(M)
		I.setEngine('LSODE')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
		I.analyse()	
		I.solve()

#	def testpeturblsode(self):
#		"test that steady conditions are stable with LSODE"
#		M = self.testdsgsat()
#		# here is the peturbation...
#		M.qdot_s.setRealValueWithUnits(1000,"W/m")
#		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
#	 	I = ascpy.Integrator(M)
#		I.setEngine('LSODE')
#		I.setReporter(ascpy.IntegratorReporterConsole(I))
#		I.setLinearTimesteps(ascpy.Units("s"), 0, 5, 1)
#		I.analyse()	
#		I.solve()

	def teststeadyida(self):	
		M = self.testdsgsat()
		self.assertAlmostEqual(M.dTw_dt[2],0.0)
		Tw1 = float(M.T_w[2])
		T = self.L.findType('dsgsat3')
		M.run(T.getMethod('free_states'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setParameter('safeeval',True)
		I.setParameter('rtol',1e-4)
		I.setParameter('atolvect',False)
		I.setParameter('atol',1e-4)
		I.setParameter('maxord',3)		
		I.setInitialSubStep(0.001)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 3600, 10)
		I.analyse()
		I.solve()
		self.assertAlmostEqual(float(M.T_w[2]),Tw1)
		M.qdot_s.setRealValueWithUnits(1000,"W/m")
		self.assertAlmostEqual(M.qdot_s.as("W/m"),1000)
		M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
		print "dTw/dt = %f" % M.dTw_dt[2]
		self.assertNotAlmostEqual(M.dTw_dt[2],0.0)
		F=file('dsgsat.dot','w')
		M.write(F,'dot')

	def testpeturbida(self):	
		M = self.testdsgsat()
		self.assertAlmostEqual(M.dTw_dt[2],0.0)
		T = self.L.findType('dsgsat3')
		M.run(T.getMethod('free_states'))
		# here is the peturbation...
		qdot_s = float(M.qdot_s)
		print "OLD QDOT_S = %f" % qdot_s
		M.qdot_s.setRealValueWithUnits(6000,"W/m")
		# IDA has its own initial conditions solver, so no need to call QRSlv here
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		#I.setLinearTimesteps(ascpy.Units("s"), 0,300,300)
		I.setLogTimesteps(ascpy.Units("s"), 0.009, 1200, 150)
		I.analyse()
		F = file('ga.mm','w')
		I.writeMatrix(F,'dg/dz')
		F = file('gd.mm','w')
		I.writeMatrix(F,'dg/dx')
		F = file('fa.mm','w')
		I.writeMatrix(F,'df/dz')
		F = file('fd.mm','w')
		I.writeMatrix(F,'df/dx')
		F = file('fdp.mm','w')
		I.writeMatrix(F,"df/dx'")
		I.solve()
		
#-------------------------------------------------------------------------------
# Testing of freesteam external steam properties functions

with_freesteam = True
try:
	# we assume that if the freesteam python module is installed, the ASCEND
	# external library will also be.
	import freesteam
	have_freesteam = True
except ImportError,e:
	have_freesteam = False

if with_freesteam and have_freesteam:
	class TestFreesteam(AscendSelfTester):
#		def testfreesteamtest(self):
#			"""run the self-test cases bundled with freesteam"""
#			self._run('testfreesteam',filename='testfreesteam.a4c')

		def testload(self):
			"""check that we can load 'thermalequilibrium2' (IMPORT "freesteam", etc)"""
			self.L.load('johnpye/thermalequilibrium2.a4c')

		def testinstantiate(self):
			"""load an instantiate 'thermalequilibrium2'"""
			self.testload()
			M = self.L.findType('thermalequilibrium2').getSimulation('sim')
			return M

		def testintegrate(self):
			"""integrate transfer of heat from one mass of water/steam to another
			according to Newton's law of cooling"""
			M = self.testinstantiate()
			M.setSolver(ascpy.Solver("QRSlv"))
			I = ascpy.Integrator(M)
			I.setEngine('LSODE')
			I.setReporter(ascpy.IntegratorReporterConsole(I))
			I.setLinearTimesteps(ascpy.Units("s"), 0, 3000, 30)
			I.setMinSubStep(0.01)
			I.setInitialSubStep(1)
			I.analyse()
			print "Number of vars = %d" % I.getNumVars()
			assert I.getNumVars()==2
			I.solve()
			assert I.getNumObservedVars() == 3;
			print "S[1].T = %f K" % M.S[1].T
			print "S[2].T = %f K" % M.S[2].T
			print "Q = %f W" % M.Q		
			self.assertAlmostEqual(float(M.S[1].T),506.77225109,4);
			self.assertAlmostEqual(float(M.S[2].T),511.605173967,5);
			self.assertAlmostEqual(float(M.Q),-48.32922877329,3);
			self.assertAlmostEqual(float(M.t),3000);
			print "Note that the above values have not been verified analytically"

		def testcollapsingcan2(self):
			""" solve the collapsing can model using IAPWS-IF97 steam props """
			M = self._run("collapsingcan2",filename="collapsingcan2.a4c");

#-------------------------------------------------------------------------------
# Testing of the brent-solver EXTERNAL method

class TestBrent(AscendSelfTester):
	def testbrent(self):
		M = self._run('brent1',filename='test/brent.a4c')

#-------------------------------------------------------------------------------
# Testing of IDA's analysis module

class TestIDA(Ascend):
	def _run(self,filen,modeln=""):
		self.L.load('test/ida/%s.a4c' % filen)
		T = self.L.findType('%s%s' % (filen,modeln))
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.analyse()
		return M;

	def _runfail(self,filen,n,msg="failed"):
		try:
			self._run(filen,'fail%d' % n)
		except Exception,e:
			print "(EXPECTED) ERROR: %s" % e
			return
		self.fail(msg)

	def testsinglederiv(self):
		self._run('singlederiv')

	def testsinglederivfail1(self):
		self._runfail('singlederiv',1
			,"t.ode_id=-1 did not trigger error")

	def testsinglederivfail2(self):
		self._runfail('singlederiv',2
			,"dy_dt.ode_id=2 did not trigger error")

	def testsinglederivfail3(self):
		self._runfail('singlederiv',3
			,"dy_dt.ode_type=3 did not trigger error")

	def testsinglederivfail4(self):
		self._runfail('singlederiv',4
			,"duplicate ode_type=1 did not trigger error")

	def testsinglederivfail5(self):
		self._runfail('singlederiv',5
			,"duplicate ode_type=1 did not trigger error")

	def testsinglederivfail6(self):
		self._runfail('singlederiv',6
			,"duplicate ode_type=1 did not trigger error")

	def testtwoderiv(self):
		self._run('twoderiv')

	def testtwoderivfail1(self):
		self._runfail('twoderiv',1)

	def testtwoderivfail2(self):
		self._runfail('twoderiv',2)

	def testtwoderivfail3(self):
		self._runfail('twoderiv',3)
	def testtwoderivfail4(self):
		self._runfail('twoderiv',4)
	def testtwoderivfail5(self):
		self._runfail('twoderiv',5)

	def testnoderivs(self):
		self._runfail('noderivs',1)

	def testnoindeps(self):
		self._runfail('indeps',1)

	def testtwoindeps(self):
		self._runfail('indeps',2)

	def testfixedvars(self):
		self._run('fixedvars')

	def testfixedvars1(self):
		self._run('fixedvars',1)

# fails the index check
#	def testfixedvars2(self):
#		self._run('fixedvars',2)

# fails the index check
#	def testfixedvars3(self):
#		self._run('fixedvars',3)

	def testincidence(self):
		self._run('incidence')

	def testincidence1(self):
		self._run('incidence',1)
	def testincidence2(self):
		self._run('incidence',2)
	def testincidence3(self):
		M = self._run('incidence',3)

	def testincidence4(self):
		self._run('incidence',4)
	def testincidencefail5(self):
		self._runfail('incidence',5)

	def testwritematrix(self):
		self.L.load('test/ida/writematrix.a4c')
		T = self.L.findType('writematrix')
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.analyse()
# this stuff fails on Windows because FILE structure is different python vs mingw
#		F = os.tmpfile()
#		I.writeMatrix(F,"dF/dy")
#		F.seek(0)
#		print F.read()
#		F1 = os.tmpfile()
#		I.writeMatrix(F1,"dF/dy'")
#		F1.seek(0)
#		print F1.read()
#		F1 = os.tmpfile()
#		I.writeMatrix(F1,"dg/dx")
#		F1.seek(0)
#		print F1.read()
#		# for the moment you'll have to check these results manually.

	def testwritematrix2(self):
		self.L.load('test/ida/writematrix.a4c')
		T = self.L.findType('writematrix2')
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.analyse()
# this stuff fails on Windows because FILE structure is different python vs mingw
#		F = os.tmpfile()
#		I.writeMatrix(F,"dF/dy")
#		F.seek(0)
#		print F.read()
#		F1 = os.tmpfile()
#		I.writeMatrix(F1,"dF/dy'")
#		F1.seek(0)
#		print F1.read()
#		F1 = os.tmpfile()
#		I.writeMatrix(F1,"dg/dx")
#		F1.seek(0)
#		print F1.read()
		#F1 = os.tmpfile()
		#I.writeMatrix(F1,"dydp/dyd")
		#F1.seek(0)
		#print F1.read()
		# for the moment you'll have to check these results manually.

	def testindexproblem(self):
		self.L.load('test/ida/indexproblem.a4c')
		T = self.L.findType('indexproblem')
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.analyse()
		pass

	def testindexproblem2(self):
		self.L.load('test/ida/indexproblem.a4c')
		T = self.L.findType('indexproblem2')
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		try:
			I.analyse()
		except Exception,e:
			return
		self.fail("Index problem not detected")

	def testboundaries(self):
		self.L.load('test/ida/boundaries.a4c')
		T = self.L.findType('boundaries')
		M = T.getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.analyse()
		I.setLogTimesteps(ascpy.Units("s"), 0.1, 20, 20)
		I.setParameter('linsolver','DENSE')
		I.setParameter('calcic','Y')
		I.setParameter('linsolver','DENSE')
		I.setParameter('safeeval',False)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.solve()

# doesn't work yet:
#	def testincidence5(self):
#		self._run('incidence',5)


#-------------------------------------------------------------------------------
# Testing of IDA models using DENSE linear solver

class TestIDADENSE(Ascend):
	"""IDA DAE integrator, DENSE linear solver"""

	def testlotka(self):
		self.L.load('johnpye/lotka.a4c')
		M = self.L.findType('lotka').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5);
		I.setParameter('linsolver','DENSE')
		I.setParameter('rtol',1e-8);
		I.analyse()
		assert I.getNumVars()==2
		assert abs(M.R - 1000) < 1e-300
		I.solve()
		assert I.getNumObservedVars() == 3
		assert abs(M.R - 832) < 1.0
		assert abs(M.F - 21.36) < 0.1
		
	def testdenx(self):
		print "-----------------------------====="
		self.L.load('johnpye/idadenx.a4c')
		M = self.L.findType('idadenx').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('calcic','YA_YDP')
		I.setParameter('linsolver','DENSE')
		I.setParameter('safeeval',True)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
		I.setMaxSubStep(0);
		I.setInitialSubStep(0)
		I.setMaxSubSteps(0)
		I.setParameter('autodiff',True)
		I.analyse()
		I.solve()
		assert abs(float(M.y1) - 5.1091e-08) < 2e-9
		assert abs(float(M.y2) - 2.0437e-13) < 2e-14
		assert abs(float(M.y3) - 1.0) < 1e-5

	def testhires(self):
		self.L.load('test/hires.a4c')
		T = self.L.findType('hires')
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setParameter('rtol',1.1e-15)
		I.setParameter('atolvect',0)
		I.setParameter('atol',1.1e-15)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLogTimesteps(ascpy.Units(""), 1, 321.8122, 5)
		I.setInitialSubStep(1e-5)
		I.setMaxSubSteps(10000)
		I.analyse()
		I.solve()
		for i in range(8):
			print "y[%d] = %.20g" % (i+1, M.y[i+1])
		M.run(T.getMethod('self_test'))

	def testchemakzo(self):
		self.L.load('test/chemakzo.a4c')
		T = self.L.findType('chemakzo')
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setParameter('rtol',1e-15)
		I.setParameter('atolvect',0)
		I.setParameter('atol',1e-15)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 1, 180, 5)
		I.setInitialSubStep(1e-13)
		I.setMaxSubSteps(10000)
		I.analyse()
		I.solve()
		for i in range(6):
			print "y[%d] = %.20g" % (i+1, M.y[i+1])
		M.run(T.getMethod('self_test'))

	def testtransamp(self):
		self.L.load('test/transamp.a4c')
		T = self.L.findType('transamp')
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setParameter('rtol',1e-7)
		I.setParameter('atolvect',0)
		I.setParameter('atol',1e-7)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0.05, 0.2, 20)
		I.setInitialSubStep(0.00001)
		I.setMaxSubSteps(10000)
		I.analyse()
		I.solve()
		for i in range(6):
			print "y[%d] = %.20g" % (i+1, M.y[i+1])
		M.run(T.getMethod('self_test'))

# MODEL FAILS ANALYSIS: we need to add support for non-incident differential vars
#	def testpollution(self):
#		self.L.load('test/pollution.a4c')
#		T = self.L.findType('pollution')
#		M = T.getSimulation('sim')
#		M.setSolver(ascpy.Solver('QRSlv'))
#		I = ascpy.Integrator(M)
#		I.setEngine('IDA')
#		I.setParameter('linsolver','DENSE')
#		I.setParameter('rtol',1.1e-15)
#		I.setParameter('atolvect',0)
#		I.setParameter('atol',1.1e-15)
#		I.setReporter(ascpy.IntegratorReporterConsole(I))
#		I.setLogTimesteps(ascpy.Units("s"), 1, 60, 5)
#		I.setInitialSubStep(1e-5)
#		I.setMaxSubSteps(10000)
#		I.analyse()
#		I.solve()
#		for i in range(20):
#			print "y[%d] = %.20g" % (i+1, M.y[i+1])
#		M.run(T.getMethod('self_test'))

## @TODO fails during IDACalcIC (model too big?)
#	def testkryx(self):
#		self.L.load('johnpye/idakryx.a4c')
#		ascpy.getCompiler().setUseRelationSharing(False)
#		M = self.L.findType('idakryx').getSimulation('sim')
#		M.setSolver(ascpy.Solver('QRSlv'))
#		M.build()
#		I = ascpy.Integrator(M)
#		I.setEngine('IDA')
#		I.setReporter(ascpy.IntegratorReporterConsole(I))
#		I.setParameter('linsolver','DENSE')
#		I.setParameter('maxl',8)
#		I.setParameter('gsmodified',False)
#		I.setParameter('autodiff',True)
#		I.setParameter('rtol',0)
#		I.setParameter('atol',1e-3);
#		I.setParameter('atolvect',False)
#		I.setParameter('calcic','YA_YDP')
#		I.analyse()
#		I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 11)
#		I.solve()
#		assert abs(M.u[2][2].getValue()) < 1e-5
	
#-------------------------------------------------------------------------------
# Testing of IDA models using SPGMR linear solver (Krylov)
	
# these tests are disabled until SPGMR preconditioning has been implemented
class TestIDASPGMR:#(Ascend):
	def testlotka(self):
		self.L.load('johnpye/lotka.a4c')
		M = self.L.findType('lotka').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 5)
		I.setParameter('rtol',1e-8)
		I.analyse()
		assert I.getNumVars()==2
		assert abs(M.R - 1000) < 1e-300
		I.solve()
		assert I.getNumObservedVars() == 3
		assert abs(M.R - 832) < 1.0
		assert abs(M.F - 21.36) < 0.1


	def testkryx(self):
		self.L.load('johnpye/idakryx.a4c')
		M = self.L.findType('idakryx').getSimulation('sim')
		M.build()
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setParameter('linsolver','SPGMR')
		I.setParameter('prec','JACOBI')
		I.setParameter('maxl',8)
		I.setParameter('gsmodified',False)
		I.setParameter('autodiff',True)
		I.setParameter('gsmodified',True)
		I.setParameter('rtol',0)
		I.setParameter('atol',1e-3);
		I.setParameter('atolvect',False)
		I.setParameter('calcic','Y')
		I.analyse()
		I.setLogTimesteps(ascpy.Units("s"), 0.01, 10.24, 10);
		print M.udot[1][3]
		I.solve()
		assert 0

	def testzill(self):
		self.L.load('johnpye/zill.a4c')
		T = self.L.findType('zill')
		M = T.getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('safeeval',False)
		I.setMinSubStep(1e-7)
		I.setMaxSubStep(0.001)
		I.setMaxSubSteps(10000)
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units(), 1.0, 1.5, 5)
		I.analyse()
		I.solve()
		M.run(T.getMethod('self_test'))

	def testdenxSPGMR(self):
		self.L.load('johnpye/idadenx.a4c')
		M = self.L.findType('idadenx').getSimulation('sim')
		M.setSolver(ascpy.Solver('QRSlv'))
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLogTimesteps(ascpy.Units("s"), 0.4, 4e10, 11)
		I.setMaxSubStep(0);
		I.setInitialSubStep(0);
		I.setMaxSubSteps(0);
		I.setParameter('autodiff',True)
		I.setParameter('linsolver','SPGMR')
		I.setParameter('gsmodified',False)
		I.setParameter('maxncf',10)
		I.analyse()
		I.solve()
		assert abs(float(M.y1) - 5.1091e-08) < 1e-10
		assert abs(float(M.y2) - 2.0437e-13) < 1e-15
		assert abs(float(M.y3) - 1.0) < 1e-5

class TestDOPRI5(Ascend):
	def testlotka(self):
		self.L.load('johnpye/dopri5/dopri5test.a4c')
		M = self.L.findType('dopri5test').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())	
		I = ascpy.Integrator(M)
		I.setEngine('DOPRI5')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 200, 20)
		I.setParameter('rtol',1e-8)
		I.analyse()
		assert I.getNumVars()==1
		I.solve()
	def testaren(self):
		self.L.load('johnpye/dopri5/aren.a4c')
		M = self.L.findType('aren').getSimulation('sim')
		M.setSolver(ascpy.Solver("QRSlv"))
		M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())
		I = ascpy.Integrator(M)
		I.setEngine('DOPRI5')
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		#xend = 17.0652165601579625588917206249
		I.setLinearTimesteps(ascpy.Units("s"), 0, 17.0652165601579625588917206249, 10)
		I.setParameter('rtol',1e-7)
		I.setParameter('atol',1e-7)
		I.setParameter('tolvect',False)
		I.setMinSubStep(0);
		I.setMaxSubStep(0);
		I.setInitialSubStep(0);
		I.analyse()
		I.solve()
		assert abs(float(M.y[0]) - 0.994) < 1e-5
		assert abs(float(M.y[1]) - 0.0) < 1e-5

# move code above down here if you want to temporarily avoid testing it
class NotToBeTested:
	def nothing(self):
		pass

	def testnewton(self):
		sys.stderr.write("STARTING TESTNEWTON\n")
		self.L.load('johnpye/newton.a4c')
		T = self.L.findType('newton')
		M = T.getSimulation('sim')
		M.solve(ascpy.Solver("QRSlv"),ascpy.SolverReporter())	
		I = ascpy.Integrator(M)
		I.setEngine('IDA')
		I.setParameter('linsolver','DENSE')
		I.setParameter('safeeval',True)
		I.setParameter('rtol',1e-8)
		I.setMaxSubStep(0.001)
		I.setMaxSubSteps(10000)
		
		I.setReporter(ascpy.IntegratorReporterConsole(I))
		I.setLinearTimesteps(ascpy.Units("s"), 0, 2*float(M.v)/float(M.g), 2)
		I.analyse()
		I.solve()
		print "At end of simulation,"
		print "x = %f" % M.x
		print "v = %f" % M.v
		M.run(T.getMethod('self_test'))

if __name__=='__main__':
	# a whole bag of tricks to make sure we get the necessary dirs in our ascend, python and ld path vars
	restart = 0

	if platform.system()=="Windows":
		LD_LIBRARY_PATH="PATH"
		SEP = ";"
	else:
		LD_LIBRARY_PATH="LD_LIBRARY_PATH"
		SEP = ":"

	freesteamdir = os.path.expanduser("~/freesteam/ascend")
	modeldirs = [os.path.abspath(os.path.join(sys.path[0],"models")),os.path.abspath(freesteamdir)]
	solverdir = os.path.abspath(os.path.join(sys.path[0],"solvers"))
	solverdirs = [os.path.join(solverdir,s) for s in "qrslv","cmslv","lrslv","conopt"]
	if not os.environ.get('ASCENDLIBRARY'):
		os.environ['ASCENDLIBRARY'] = SEP.join(modeldirs+solverdirs)
		restart = 1
	else:
		envmodelsdir = [os.path.abspath(i) for i in os.environ['ASCENDLIBRARY'].split(SEP)]
		for l in modeldirs:
			if l in envmodelsdir[len(modeldirs):]:
				envmodelsdir.remove(l)
				restart = 1
		for l in modeldirs:
			if l not in envmodelsdir:
				envmodelsdir.insert(0,l)
				restart = 1
		os.environ['ASCENDLIBRARY'] = SEP.join(envmodelsdir)	

	libdirs = ["pygtk","."]
	libdirs = [os.path.normpath(os.path.join(sys.path[0],l)) for l in libdirs]
	if not os.environ.get(LD_LIBRARY_PATH):
		os.environ[LD_LIBRARY_PATH]=SEP.join(libdirs)
		restart = 1
	else:
		envlibdirs = [os.path.normpath(i) for i in os.environ[LD_LIBRARY_PATH].split(SEP)]
		for l in libdirs:
			if l in envlibdirs[len(libdirs):]:
				envlibdirs.remove(l)
				restart = 1
		for l in libdirs:
			if l not in envlibdirs:
				envlibdirs.insert(0,l)
				restart = 1		
		os.environ[LD_LIBRARY_PATH] = SEP.join(envlibdirs)

	pypath = os.path.normpath(os.path.join(sys.path[0],"pygtk"))
	if not os.environ.get('PYTHONPATH'):
		os.environ['PYTHONPATH']=pypath
	else:
		envpypath = os.environ['PYTHONPATH'].split(SEP)
		if pypath not in envpypath:
			envpypath.insert(0,pypath)
			os.environ['PYTHONPATH']=SEP.join(envpypath)
			restart = 1

	if restart:
		if platform.system()=="Windows":
			pass
		else:
			script = os.path.join(sys.path[0],"test.py")					
			sys.stderr.write("Restarting with...\n")
			sys.stderr.write("  export LD_LIBRARY_PATH=%s\n" % os.environ.get(LD_LIBRARY_PATH))
			sys.stderr.write("  export PYTHONPATH=%s\n" % os.environ.get('PYTHONPATH'))
			sys.stderr.write("  export ASCENDLIBRARY=%s\n" % os.environ.get('ASCENDLIBRARY'))
			sys.stderr.flush()
			os.execvp("python",[script] + sys.argv)

	import ascpy

	try:
		import cunit
	except:
		pass

	atexit.register(ascpy.shutdown)
	#suite = unittest.TestSuite()
	#suite = unittest.defaultTestLoader.loadTestsFromName('__main__')
	#unittest.TextTestRunner(verbosity=2).run(suite)
	unittest.main()	
