Enable linting with ruff

.editorconfig

@@ -7,6 +7,7 @@ root = true
 [*]
 end_of_line = lf
 insert_final_newline = true
+trim_trailing_whitespace = false
 
 # Matches multiple files with brace expansion notation
 # Set default charset

.github/workflows/ci.yml

@@ -14,6 +14,20 @@ concurrency:
   cancel-in-progress: true
 
 jobs:
+  typos:
+    name: Spelling (typos)
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+    - uses: crate-ci/typos@master
+
+  ruff:
+    name: Linting (ruff)
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+    - uses: chartboost/ruff-action@v1
+
   test:
     name: Unittest (${{ matrix.os }}-${{ matrix.compiler }}-py${{ matrix.python-version }})
     runs-on: ${{ matrix.os }}

docs/conf.py

@@ -1,65 +1,71 @@
-import sys
 import os
+import sys
 from unittest.mock import MagicMock as Mock
+
 from setuptools_scm import get_version
 
+
 MOCK_MODULES = [
-	'numpy', 'numpy.testing', 'numpy.random',
-	'scipy', 'scipy.integrate', 'scipy.integrate._ode', 'scipy.stats', 'scipy.integrate._ivp.ivp',
-	'symengine', 'symengine.printing', 'symengine.lib.symengine_wrapper',
-	'jitcxde_common.helpers','jitcxde_common.numerical','jitcxde_common.symbolic','jitcxde_common.transversal',
+	"numpy", "numpy.testing", "numpy.random",
+	"scipy", "scipy.integrate", "scipy.integrate._ode", "scipy.stats", "scipy.integrate._ivp.ivp",
+	"symengine", "symengine.printing", "symengine.lib.symengine_wrapper",
+	"jitcxde_common.helpers","jitcxde_common.numerical","jitcxde_common.symbolic","jitcxde_common.transversal",
 	]
 sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
 
-class ode_mock(object): pass
-sys.modules['scipy.integrate'] = Mock(ode=ode_mock)
-class GroupHandler_mock(object): pass
-sys.modules['jitcxde_common.transversal'] = Mock(GroupHandler=GroupHandler_mock)
+class ode_mock:
+	pass
+
+class GroupHandler_mock:
+	pass
+
+sys.modules["scipy.integrate"] = Mock(ode=ode_mock)
+sys.modules["jitcxde_common.transversal"] = Mock(GroupHandler=GroupHandler_mock)
 
 sys.path.insert(0,os.path.abspath("../examples"))
 sys.path.insert(0,os.path.abspath("../jitcode"))
 
-needs_sphinx = '1.3'
+needs_sphinx = "1.3"
 
 extensions = [
-    'sphinx.ext.autodoc',
-    'sphinx.ext.autosummary',
-    'sphinx.ext.mathjax',
-    'numpydoc',
-    'sphinx.ext.graphviz'
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx.ext.mathjax",
+    "sphinx.ext.graphviz"
+    "numpydoc",
 ]
 
-source_suffix = '.rst'
+source_suffix = ".rst"
 
-master_doc = 'index'
+master_doc = "index"
 
-project = u'JiTCODE'
-copyright = u'2016, Gerrit Ansmann'
+project = "JiTCODE"
+copyright = "2016, Gerrit Ansmann"
 
-release = version = get_version(root='..', relative_to=__file__)
+release = version = get_version(root="..", relative_to=__file__)
 
 default_role = "any"
 
 add_function_parentheses = True
 
 add_module_names = False
 
-html_theme = 'nature'
-pygments_style = 'colorful'
-htmlhelp_basename = 'JiTCODEdoc'
+html_theme = "nature"
+pygments_style = "colorful"
+htmlhelp_basename = "JiTCODEdoc"
 
 numpydoc_show_class_members = False
-autodoc_member_order = 'bysource'
+autodoc_member_order = "bysource"
 
 graphviz_output_format = "svg"
 
-toc_object_entries_show_parents = 'hide'
+toc_object_entries_show_parents = "hide"
 
 def on_missing_reference(app, env, node, contnode):
-	if node['reftype'] == 'any':
+	if node["reftype"] == "any":
 		return contnode
 	else:
 		return None
 
 def setup(app):
-	app.connect('missing-reference', on_missing_reference)
+	app.connect("missing-reference", on_missing_reference)

examples/SW_of_Roesslers.py

@@ -1,5 +1,4 @@
 #!/usr/bin/python3
-# -*- coding: utf-8 -*-
 
 """
 This example showcases several advanced features of JiTCODE that are relevant for an efficient integration of more complex systems as well as how to deal with some special situations. Therefore it is pretty bizarre from a dynamical perspective.
@@ -13,7 +12,7 @@
 	\\dot{z}_i &= b + z_i (x_i -c) + k \\sum_{j=0}^N (z_j-z_i)
 	\\end{alignedat}
 
-The control parameters shall be :math:`a = 0.165`, :math:`b = 0.2`, :math:`c = 10.0`, and :math:`k = 0.01`. The (frequency) parameter :math:`ω_i` shall be picked randomly from the uniform distribution on :math:`[0.8,1.0]` for each :math:`i`. :math:`A∈ℝ^{N×N}` shall be the adjacency matrix of a one-dimensional small-world network (which shall be provided by a function `small_world_network` in the following example code). So, the :math:`x` compenents are coupled diffusively with a small-world coupling topology, while the :math:`z` components are coupled diffusively to their mean field, with the coupling term being modulated with :math:`\\sin(t)`.
+The control parameters shall be :math:`a = 0.165`, :math:`b = 0.2`, :math:`c = 10.0`, and :math:`k = 0.01`. The (frequency) parameter :math:`ω_i` shall be picked randomly from the uniform distribution on :math:`[0.8,1.0]` for each :math:`i`. :math:`A∈ℝ^{N×N}` shall be the adjacency matrix of a one-dimensional small-world network (which shall be provided by a function `small_world_network` in the following example code). So, the :math:`x` components are coupled diffusively with a small-world coupling topology, while the :math:`z` components are coupled diffusively to their mean field, with the coupling term being modulated with :math:`\\sin(t)`.
 
 Without further ado, here is the example code (`complete running example <https://raw.githubusercontent.com/neurophysik/jitcode/master/examples/SW_of_Roesslers.py>`_); highlighted lines will be commented upon below:
 
@@ -47,10 +46,10 @@ def small_world_network(number_of_nodes, nearest_neighbours, rewiring_probabilit
 		# rewiring
 		for i in range(n):
 			for j in range(i):
-				if A[i,j] and (np.random.random() < rewiring_probability):
+				if A[i,j] and (rng.random() < rewiring_probability):
 					A[j,i] = A[i,j] = False
 					while True:
-						i_new,j_new = np.random.randint(0,n,2)
+						i_new,j_new = rng.integers(0,n,2)
 						if A[i_new,j_new] or i_new==j_new:
 							continue
 						else:
@@ -60,15 +59,17 @@ def small_world_network(number_of_nodes, nearest_neighbours, rewiring_probabilit
 		return A
 
 	# example-start
-	from jitcode import jitcode, y, t
 	import numpy as np
 	import symengine
+
+	from jitcode import jitcode, t, y
 
 	# parameters
 	# ----------
 
+	rng = np.random.default_rng()
 	N = 500
-	ω = np.random.uniform(0.8,1.0,N)
+	ω = rng.uniform(0.8,1.0,N)
 	a = 0.165
 	b = 0.2
 	c = 10.0
@@ -99,11 +100,11 @@ def f():
 	# integrate
 	# ---------
 
-	initial_state = np.random.random(3*N)
+	initial_state = rng.random(3*N)
 
 	ODE = jitcode(f, helpers=helpers, n=3*N)
 	ODE.generate_f_C(simplify=False, do_cse=False, chunk_size=150)
-	ODE.set_integrator('dopri5')
+	ODE.set_integrator("dopri5")
 	ODE.set_initial_value(initial_state,0.0)
 
 	# data structure: x[0], v[0], z[0], x[1], …, x[N], v[N], z[N]

examples/benchmark.py

@@ -1,16 +1,18 @@
-import numpy as np
-from numpy.random import choice, uniform
 from time import process_time
-from scipy.integrate import ode, solve_ivp, odeint
+
+import numpy as np
+from scipy.integrate import ode, odeint, solve_ivp
 from scipy.integrate._ivp.ivp import METHODS
-from jitcode import jitcode, y
 from symengine import sin
 
+from jitcode import jitcode, y
+
+
 solver_ode = "dopri5"
 solver_ivp = "RK45"
 
 # Context manager for timing
-class timer(object):
+class timer:
 	def __init__(self,name):
 		self.name = name
 
@@ -20,105 +22,105 @@ def __enter__(self):
 	def __exit__(self,*args):
 		end = process_time()
 		duration = end-self.start
-		print("%s took %.5f s" % (self.name,duration))
+		print(f"{self.name} took {duration:.5f}s")
 
 # The actual test
 def test_scenario(name,fun,initial,times,rtol,atol):
 	print(40*"-",name,40*"-",sep="\n")
 
-	with timer("ode (%s)"%solver_ode):
-		I = ode(fun)
-		I.set_integrator(solver_ode,rtol=rtol,atol=atol,nsteps=10**8)
-		I.set_initial_value(initial,0.0)
-		result = np.vstack([I.integrate(time) for time in times])
-	assert I.successful()
+	with timer(f"ode ({solver_ode})"):
+		solver = ode(fun)
+		solver.set_integrator(solver_ode,rtol=rtol,atol=atol,nsteps=10**8)
+		solver.set_initial_value(initial,0.0)
+		_result = np.vstack([solver.integrate(time) for time in times])
+	assert solver.successful()
 
 	inv_fun = lambda y,t: fun(t,y)
 	with timer("odeint with suboptimal function (LSODA)"):
-		result = odeint(
+		_result = odeint(
 				func=inv_fun,
-				y0=initial, t=[0.0]+list(times),
+				y0=initial, t=[0.0, *times],
 				rtol=rtol, atol=atol,
 				mxstep=10**8
 				)
 
-	with timer("solve_ivp (%s) without result"%solver_ivp):
-		I = solve_ivp(
+	with timer(f"solve_ivp ({solver_ivp}) without result"):
+		solver = solve_ivp(
 				fun,
 				t_span=(0,times[-1]),
 				y0=initial,
 				method=solver_ivp, rtol=rtol, atol=atol
 			)
-	assert I.status != -1
+	assert solver.status != -1
 
-	with timer("solve_ivp (%s)"%solver_ivp):
-		I = solve_ivp(
+	with timer(f"solve_ivp ({solver_ivp})"):
+		solver = solve_ivp(
 				fun,
 				t_span=(0,times[-1]), t_eval=times,
 				y0=initial,
 				method=solver_ivp, rtol=rtol, atol=atol
 			)
-		result = I.y
-	assert I.status != -1
+		_result = solver.y
+	assert solver.status != -1
 
-	with timer("solve_ivp (%s) with dense_output"%solver_ivp):
-		I = solve_ivp(
+	with timer(f"solve_ivp ({solver_ivp}) with dense_output"):
+		solver = solve_ivp(
 				fun,
 				t_span=(0,times[-1]),
 				y0=initial,
 				method=solver_ivp, rtol=rtol, atol=atol,
 				dense_output=True
 			)
-		result = np.vstack([I.sol(time) for time in times])
-	assert I.status != -1
+		_result = np.vstack([solver.sol(time) for time in times])
+	assert solver.status != -1
 
-	with timer("%s with dense output"%solver_ivp):
-		I = METHODS[solver_ivp](
+	with timer(f"{solver_ivp} with dense output"):
+		solver = METHODS[solver_ivp](
 				fun=fun,
 				y0=initial, t0=0.0, t_bound=times[-1],
 				rtol=rtol, atol=atol
 			)
 		def solutions():
 			for time in times:
-				while I.t < time:
-					I.step()
-				yield I.dense_output()(time)
-		result = np.vstack(list(solutions()))
-	assert I.status != "failed"
+				while solver.t < time:
+					solver.step()
+				yield solver.dense_output()(time)
+		_result = np.vstack(list(solutions()))
+	assert solver.status != "failed"
 
-	with timer("%s with manual resetting"%solver_ivp):
-		I = METHODS[solver_ivp](
+	with timer(f"{solver_ivp} with manual resetting"):
+		solver = METHODS[solver_ivp](
 				fun=fun,
 				y0=initial, t0=0.0, t_bound=times[-1],
 				rtol=rtol, atol=atol
 			)
 		def solutions():
 			for time in times:
-				I.t_bound = time
-				I.status = "running"
-				while I.status == "running":
-					I.step()
-				yield I.y
-		result = np.vstack(list(solutions()))
-	assert I.status != "failed"
+				solver.t_bound = time
+				solver.status = "running"
+				while solver.status == "running":
+					solver.step()
+				yield solver.y
+		_result = np.vstack(list(solutions()))
+	assert solver.status != "failed"
 
-	with timer("%s with reinitialising"%solver_ivp):
+	with timer(f"{solver_ivp} with reinitialising"):
 		def solutions():
 			current_time = 0.0
 			state = initial
 			for time in times:
-				I = METHODS[solver_ivp](
+				solver = METHODS[solver_ivp](
 						fun=fun,
 						y0=state, t0=current_time, t_bound=time,
 						rtol=rtol, atol=atol
 					)
-				while I.status == "running":
-					I.step()
-				assert I.status != "failed"
+				while solver.status == "running":
+					solver.step()
+				assert solver.status != "failed"
 				current_time = time
-				state = I.y
+				state = solver.y
 				yield state
-		result = np.vstack(list(solutions()))
+		_result = np.vstack(list(solutions()))
 
 # Using compiled functions to make things faster
 def get_compiled_function(f):
@@ -142,9 +144,12 @@ def get_compiled_function(f):
 )
 
 
+rng = np.random.default_rng()
+
 n, c, q = 100, 3.0, 0.2
-A = choice( [1,0], size=(n,n), p=[q,1-q] )
-omega = uniform(-0.5,0.5,n)
+A = rng.choice( [1,0], size=(n,n), p=[q,1-q] )
+omega = rng.uniform(-0.5,0.5,n)
+
 def kuramotos_f():
 	for i in range(n):
 		coupling_sum = sum(
@@ -157,10 +162,8 @@ def kuramotos_f():
 test_scenario(
 	name = "random network of Kuramoto oscillators",
 	fun = get_compiled_function(kuramotos_f),
-	initial = uniform(0,2*np.pi,n),
+	initial = rng.uniform(0,2*np.pi,n),
 	times = range(1,10001,10),
 	rtol = 1e-13,
 	atol = 1e-6,
 )
-
-