Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 260 additions and 1299 deletions
doc/sphinx/simplifications.rst 0 → 100644
View file @ 5600b6b6
***************************************
Assignment Collection & Simplifications
***************************************
AssignmentCollection
====================
.. autoclass:: pystencils.AssignmentCollection
:members:
SimplificationStrategy
======================
.. autoclass:: pystencils.simp.SimplificationStrategy
:members:
Simplifications
===============
.. automodule:: pystencils.simp.simplifications
:members:
Subexpression insertion
=======================
The subexpression insertions have the goal to insert subexpressions which will not reduce the number of FLOPs.
For example a constant value kept as subexpression will lead to a new variable in the code which will occupy
a register slot. On the other side a single variable could just be inserted in all assignments.
.. automodule:: pystencils.simp.subexpression_insertion
:members:
doc/sphinx/stencil.rst 0 → 100644
View file @ 5600b6b6
*******
Stencil
*******
.. automodule:: pystencils.stencil
:members:
doc/sphinx/tutorials.rst 0 → 100644
View file @ 5600b6b6
Tutorials
=========
These tutorials are a good place to start if you are new to pystencils.
All tutorials and demos listed here are based on Jupyter notebooks that you can find in the pystencils repository.
It is a good idea to download them and run them directly to be able to play around with the code.
.. toctree::
:maxdepth: 1
/notebooks/01_tutorial_getting_started.ipynb
/notebooks/02_tutorial_basic_kernels.ipynb
/notebooks/03_tutorial_datahandling.ipynb
/notebooks/04_tutorial_advection_diffusion.ipynb
/notebooks/05_tutorial_phasefield_spinodal_decomposition.ipynb
/notebooks/06_tutorial_phasefield_dentritic_growth.ipynb
/notebooks/demo_assignment_collection.ipynb
/notebooks/demo_plotting_and_animation.ipynb
/notebooks/demo_derivatives.ipynb
/notebooks/demo_benchmark.ipynb
/notebooks/demo_wave_equation.ipynb
equationcollection/__init__.py deleted 100644 → 0
View file @ cfef0f9d
from pystencils.equationcollection.equationcollection import EquationCollection
from pystencils.equationcollection.simplificationstrategy import SimplificationStrategy
equationcollection/equationcollection.py deleted 100644 → 0
View file @ cfef0f9d
import sympy as sp
from pystencils.sympyextensions import fastSubs, countNumberOfOperations
class EquationCollection:
"""
A collection of equations with subexpression definitions, also represented as equations,
that are used in the main equations. EquationCollections can be passed to simplification methods.
These simplification methods can change the subexpressions, but the number and
left hand side of the main equations themselves is not altered.
Additionally a dictionary of simplification hints is stored, which are set by the functions that create
equation collections to transport information to the simplification system.
:ivar mainEquations: list of sympy equations
:ivar subexpressions: list of sympy equations defining subexpressions used in main equations
:ivar simplificationHints: dictionary that is used to annotate the equation collection with hints that are
used by the simplification system. See documentation of the simplification rules for
potentially required hints and their meaning.
"""
# ----------------------------------------- Creation ---------------------------------------------------------------
def __init__(self, equations, subExpressions, simplificationHints={}):
self.mainEquations = equations
self.subexpressions = subExpressions
self.simplificationHints = simplificationHints
def symbolGen():
"""Use this generator to create new unused symbols for subexpressions"""
counter = 0
while True:
counter += 1
newSymbol = sp.Symbol("xi_" + str(counter))
if newSymbol in self.boundSymbols:
continue
yield newSymbol
self.subexpressionSymbolNameGenerator = symbolGen()
def createNewWithAdditionalSubexpressions(self, newEquations, additionalSubExpressions):
assert len(self.mainEquations) == len(newEquations), "Number of update equations cannot be changed"
return EquationCollection(newEquations,
self.subexpressions + additionalSubExpressions,
self.simplificationHints)
def createNewWithSubstitutionsApplied(self, substitutionDict):
newSubexpressions = [fastSubs(eq, substitutionDict) for eq in self.subexpressions]
newEquations = [fastSubs(eq, substitutionDict) for eq in self.mainEquations]
return EquationCollection(newEquations, newSubexpressions, self.simplificationHints)
def addSimplificationHint(self, key, value):
assert key not in self.simplificationHints, "This hint already exists"
self.simplificationHints[key] = value
# ---------------------------------------------- Properties -------------------------------------------------------
@property
def allEquations(self):
return self.subexpressions + self.mainEquations
@property
def freeSymbols(self):
"""All symbols used in the equation collection, which have not been defined inside the equation system"""
freeSymbols = set()
for eq in self.allEquations:
freeSymbols.update(eq.rhs.atoms(sp.Symbol))
return freeSymbols - self.boundSymbols
@property
def boundSymbols(self):
"""Set of all symbols which occur on left-hand-sides i.e. all symbols which are defined."""
boundSymbolsSet = set([eq.lhs for eq in self.allEquations])
assert len(boundSymbolsSet) == len(self.subexpressions) + len(self.mainEquations), \
"Not in SSA form - same symbol assigned multiple times"
return boundSymbolsSet
@property
def definedSymbols(self):
"""All symbols that occur as left-hand-sides of the main equations"""
return set([eq.lhs for eq in self.mainEquations])
@property
def operationCount(self):
"""See :func:`countNumberOfOperations` """
return countNumberOfOperations(self.allEquations)
def get(self, symbols, fromMainEquationsOnly=False):
"""Return the equations which have symbols as left hand sides"""
if not hasattr(symbols, "__len__"):
symbols = list(symbols)
symbols = set(symbols)
if not fromMainEquationsOnly:
eqsToSearchIn = self.allEquations
else:
eqsToSearchIn = self.mainEquations
return [eq for eq in eqsToSearchIn if eq.lhs in symbols]
# ----------------------------------------- Display and Printing -------------------------------------------------
def _repr_html_(self):
def makeHtmlEquationTable(equations):
noBorder = 'style="border:none"'
htmlTable = '<table style="border:none; width: 100%; ">'
line = '<tr {nb}> <td {nb}>$${eq}$$</td> </tr> '
for eq in equations:
formatDict = {'eq': sp.latex(eq),
'nb': noBorder, }
htmlTable += line.format(**formatDict)
htmlTable += "</table>"
return htmlTable
result = ""
if len(self.subexpressions) > 0:
result += "<div>Subexpressions:</div>"
result += makeHtmlEquationTable(self.subexpressions)
result += "<div>Main Equations:</div>"
result += makeHtmlEquationTable(self.mainEquations)
return result
def __repr__(self):
return "Equation Collection for " + ",".join([str(eq.lhs) for eq in self.mainEquations])
def __str__(self):
result = "Subexpressions\n"
for eq in self.subexpressions:
result += str(eq) + "\n"
result += "Main Equations\n"
for eq in self.mainEquations:
result += str(eq) + "\n"
return result
# ------------------------------------- Manipulation ------------------------------------------------------------
def merge(self, other):
"""Returns a new collection which contains self and other. Subexpressions are renamed if they clash."""
ownDefs = set([e.lhs for e in self.mainEquations])
otherDefs = set([e.lhs for e in other.mainEquations])
assert len(ownDefs.intersection(otherDefs)) == 0, "Cannot merge, since both collection define the same symbols"
ownSubexpressionSymbols = {e.lhs: e.rhs for e in self.subexpressions}
substitutionDict = {}
processedOtherSubexpressionEquations = []
for otherSubexpressionEq in other.subexpressions:
if otherSubexpressionEq.lhs in ownSubexpressionSymbols:
if otherSubexpressionEq.rhs == ownSubexpressionSymbols[otherSubexpressionEq.lhs]:
continue # exact the same subexpression equation exists already
else:
# different definition - a new name has to be introduced
newLhs = self.subexpressionSymbolNameGenerator()
newEq = sp.Eq(newLhs, fastSubs(otherSubexpressionEq.rhs, substitutionDict))
processedOtherSubexpressionEquations.append(newEq)
substitutionDict[otherSubexpressionEq.lhs] = newLhs
else:
processedOtherSubexpressionEquations.append(fastSubs(otherSubexpressionEq, substitutionDict))
return EquationCollection(self.mainEquations + other.mainEquations,
self.subexpressions + processedOtherSubexpressionEquations)
def extract(self, symbolsToExtract):
"""
Creates a new equation collection with equations that have symbolsToExtract as left-hand-sides and
only the necessary subexpressions that are used in these equations
"""
symbolsToExtract = set(symbolsToExtract)
newEquations = []
subexprMap = {e.lhs: e.rhs for e in self.subexpressions}
handledSymbols = set()
queue = []
def addSymbolsFromExpr(expr):
dependentSymbols = expr.atoms(sp.Symbol)
for ds in dependentSymbols:
if ds not in handledSymbols:
queue.append(ds)
handledSymbols.add(ds)
for eq in self.mainEquations:
if eq.lhs in symbolsToExtract:
newEquations.append(eq)
addSymbolsFromExpr(eq.rhs)
while len(queue) > 0:
e = queue.pop(0)
if e not in subexprMap:
continue
else:
addSymbolsFromExpr(subexprMap[e])
newSubExpr = [eq for eq in self.subexpressions if eq.lhs in handledSymbols]
return EquationCollection(newEquations, newSubExpr)
def newWithoutUnusedSubexpressions(self):
"""Returns a new equation collection containing only the subexpressions that
are used/referenced in the equations"""
allLhs = [eq.lhs for eq in self.mainEquations]
return self.extract(allLhs)
def insertSubexpressions(self):
"""Returns a new equation collection by inserting all subexpressions into the main equations"""
if len(self.subexpressions) == 0:
return EquationCollection(self.mainEquations, self.subexpressions, self.simplificationHints)
subsDict = {self.subexpressions[0].lhs: self.subexpressions[0].rhs}
subExpr = [e for e in self.subexpressions]
for i in range(1, len(subExpr)):
subExpr[i] = fastSubs(subExpr[i], subsDict)
subsDict[subExpr[i].lhs] = subExpr[i].rhs
newEq = [fastSubs(eq, subsDict) for eq in self.mainEquations]
return EquationCollection(newEq, [], self.simplificationHints)
def lambdify(self, symbols, module=None, fixedSymbols={}):
"""
Returns a function to evaluate this equation collection
:param symbols: symbol(s) which are the parameter for the created function
:param module: same as sympy.lambdify paramter of same same, i.e. which module to use e.g. 'numpy'
:param fixedSymbols: dictionary with substitutions, that are applied before lambdification
"""
eqs = self.createNewWithSubstitutionsApplied(fixedSymbols).insertSubexpressions().mainEquations
lambdas = {eq.lhs: sp.lambdify(symbols, eq.rhs, module) for eq in eqs}
def f(*args, **kwargs):
return {s: f(*args, **kwargs) for s, f in lambdas.items()}
return f
equationcollection/simplifications.py deleted 100644 → 0
View file @ cfef0f9d
import sympy as sp
from pystencils.equationcollection import EquationCollection
from pystencils.sympyextensions import replaceAdditive
def sympyCSE(equationCollection):
"""
Searches for common subexpressions inside the equation collection, in both the existing subexpressions as well
as the equations themselves. It uses the sympy subexpression detection to do this. Return a new equation collection
with the additional subexpressions found
"""
symbolGen = equationCollection.subexpressionSymbolNameGenerator
replacements, newEq = sp.cse(equationCollection.subexpressions + equationCollection.mainEquations, symbols=symbolGen)
replacementEqs = [sp.Eq(*r) for r in replacements]
modifiedSubexpressions = newEq[:len(equationCollection.subexpressions)]
modifiedUpdateEquations = newEq[len(equationCollection.subexpressions):]
newSubexpressions = replacementEqs + modifiedSubexpressions
topologicallySortedPairs = sp.cse_main.reps_toposort([[e.lhs, e.rhs] for e in newSubexpressions])
newSubexpressions = [sp.Eq(a[0], a[1]) for a in topologicallySortedPairs]
return EquationCollection(modifiedUpdateEquations, newSubexpressions, equationCollection.simplificationHints)
def applyOnAllEquations(equationCollection, operation):
"""Applies sympy expand operation to all equations in collection"""
result = [operation(s) for s in equationCollection.mainEquations]
return equationCollection.createNewWithAdditionalSubexpressions(result, [])
def applyOnAllSubexpressions(equationCollection, operation):
return EquationCollection(equationCollection.mainEquations,
[operation(s) for s in equationCollection.subexpressions],
equationCollection.simplificationHints)
def subexpressionSubstitutionInExistingSubexpressions(equationCollection):
"""Goes through the subexpressions list and replaces the term in the following subexpressions"""
result = []
for outerCtr, s in enumerate(equationCollection.subexpressions):
newRhs = s.rhs
for innerCtr in range(outerCtr):
subExpr = equationCollection.subexpressions[innerCtr]
newRhs = replaceAdditive(newRhs, subExpr.lhs, subExpr.rhs, requiredMatchReplacement=1.0)
newRhs = newRhs.subs(subExpr.rhs, subExpr.lhs)
result.append(sp.Eq(s.lhs, newRhs))
return EquationCollection(equationCollection.mainEquations, result, equationCollection.simplificationHints)
def subexpressionSubstitutionInUpdateEquations(equationCollection):
"""Replaces already existing subexpressions in the equations of the equationCollection"""
result = []
for s in equationCollection.mainEquations:
newRhs = s.rhs
for subExpr in equationCollection.subexpressions:
newRhs = replaceAdditive(newRhs, subExpr.lhs, subExpr.rhs, requiredMatchReplacement=1.0)
result.append(sp.Eq(s.lhs, newRhs))
return equationCollection.createNewWithAdditionalSubexpressions(result, [])
equationcollection/simplificationstrategy.py deleted 100644 → 0
View file @ cfef0f9d
import sympy as sp
import textwrap
from collections import namedtuple
class SimplificationStrategy:
"""
A simplification strategy is an ordered collection of simplification rules.
Each simplification is a function taking an equation collection, and returning a new simplified
equation collection. The strategy can nicely print intermediate simplification stages and results
to Jupyter notebooks.
"""
def __init__(self):
self._rules = []
def addSimplificationRule(self, rule):
"""
Adds the given simplification rule to the end of the collection.
:param rule: function that taking one equation collection and returning a (simplified) equation collection
"""
self._rules.append(rule)
@property
def rules(self):
return self._rules
def apply(self, updateRule):
"""Applies all simplification rules to the given equation collection"""
for t in self._rules:
updateRule = t(updateRule)
return updateRule
def __call__(self, equationCollection):
"""Same as apply"""
return self.apply(equationCollection)
def createSimplificationReport(self, equationCollection):
"""
Returns a simplification report containing the number of operations at each simplification stage, together
with the run-time the simplification took.
"""
ReportElement = namedtuple('ReportElement', ['simplificationName', 'adds', 'muls', 'divs', 'runtime'])
class Report:
def __init__(self):
self.elements = []
def add(self, element):
self.elements.append(element)
def __str__(self):
try:
import tabulate
return tabulate(self.elements, headers=['Name', 'Adds', 'Muls', 'Divs', 'Runtime'])
except ImportError:
result = "Name, Adds, Muls, Divs, Runtime\n"
for e in self.elements:
result += ",".join(e) + "\n"
return result
def _repr_html_(self):
htmlTable = '<table style="border:none">'
htmlTable += "<tr> <th>Name</th> <th>Adds</th> <th>Muls</th> <th>Divs</th> <th>Runtime</th></tr>"
line = "<tr><td>{simplificationName}</td>" \
"<td>{adds}</td> <td>{muls}</td> <td>{divs}</td> <td>{runtime}</td> </tr>"
for e in self.elements:
htmlTable += line.format(**e._asdict())
htmlTable += "</table>"
return htmlTable
import time
report = Report()
op = equationCollection.operationCount
report.add(ReportElement("OriginalTerm", op['adds'], op['muls'], op['divs'], '-'))
for t in self._rules:
startTime = time.perf_counter()
equationCollection = t(equationCollection)
endTime = time.perf_counter()
op = equationCollection.operationCount
timeStr = "%.2f ms" % ((endTime - startTime) * 1000,)
report.add(ReportElement(t.__name__, op['adds'], op['muls'], op['divs'], timeStr))
return report
def showIntermediateResults(self, equationCollection, symbols=None):
class IntermediateResults:
def __init__(self, strategy, eqColl, resSyms):
self.strategy = strategy
self.equationCollection = eqColl
self.restrictSymbols = resSyms
def __str__(self):
def printEqCollection(title, eqColl):
text = title
if self.restrictSymbols:
text += "\n".join([str(e) for e in eqColl.get(self.restrictSymbols)])
else:
text += textwrap.indent(str(eqColl), " " * 3)
return text
result = printEqCollection("Initial Version", self.equationCollection)
eqColl = self.equationCollection
for rule in self.strategy.rules:
eqColl = rule(eqColl)
result += printEqCollection(rule.__name__, eqColl)
return result
def _repr_html_(self):
def printEqCollection(title, eqColl):
text = '<h5 style="padding-bottom:10px">%s</h5> <div style="padding-left:20px;">' % (title, )
if self.restrictSymbols:
text += "\n".join(["$$" + sp.latex(e) + '$$' for e in eqColl.get(self.restrictSymbols)])
else:
text += eqColl._repr_html_()
text += "</div>"
return text
result = printEqCollection("Initial Version", self.equationCollection)
eqColl = self.equationCollection
for rule in self.strategy.rules:
eqColl = rule(eqColl)
result += printEqCollection(rule.__name__, eqColl)
return result
return IntermediateResults(self, equationCollection, symbols)
def __repr__(self):
result = "Simplification Strategy:\n"
for t in self._rules:
result += " - %s\n" % (t.__name__,)
return result
field.py deleted 100644 → 0
View file @ cfef0f9d
from itertools import chain
import numpy as np
import sympy as sp
from sympy.core.cache import cacheit
from sympy.tensor import IndexedBase
from pystencils.types import TypedSymbol
class Field:
"""
With fields one can formulate stencil-like update rules on structured grids.
This Field class knows about the dimension, memory layout (strides) and optionally about the size of an array.
Creating Fields:
To create a field use one of the static create* members. There are two options:
1. create a kernel with fixed loop sizes i.e. the shape of the array is already known. This is usually the
case if just-in-time compilation directly from Python is done. (see :func:`Field.createFromNumpyArray`)
2. create a more general kernel that works for variable array sizes. This can be used to create kernels
beforehand for a library. (see :func:`Field.createGeneric`)
Dimensions:
A field has spatial and index dimensions, where the spatial dimensions come first.
The interpretation is that the field has multiple cells in (usually) two or three dimensional space which are
looped over. Additionally N values are stored per cell. In this case spatialDimensions is two or three,
and indexDimensions equals N. If you want to store a matrix on each point in a two dimensional grid, there
are four dimensions, two spatial and two index dimensions: ``len(arr.shape) == spatialDims + indexDims``
Indexing:
When accessing (indexing) a field the result is a FieldAccess which is derived from sympy Symbol.
First specify the spatial offsets in [], then in case indexDimension>0 the indices in ()
e.g. ``f[-1,0,0](7)``
Example without index dimensions:
>>> a = np.zeros([10, 10])
>>> f = Field.createFromNumpyArray("f", a, indexDimensions=0)
>>> jacobi = ( f[-1,0] + f[1,0] + f[0,-1] + f[0,1] ) / 4
Example with index dimensions: LBM D2Q9 stream pull
>>> stencil = np.array([[0,0], [0,1], [0,-1]])
>>> src = Field.createGeneric("src", spatialDimensions=2, indexDimensions=1)
>>> dst = Field.createGeneric("dst", spatialDimensions=2, indexDimensions=1)
>>> for i, offset in enumerate(stencil):
... sp.Eq(dst[0,0](i), src[-offset](i))
Eq(dst_C^0, src_C^0)
Eq(dst_C^1, src_S^1)
Eq(dst_C^2, src_N^2)
"""
@staticmethod
def createFromNumpyArray(fieldName, npArray, indexDimensions=0):
"""
Creates a field based on the layout, data type, and shape of a given numpy array.
Kernels created for these kind of fields can only be called with arrays of the same layout, shape and type.
:param fieldName: symbolic name for the field
:param npArray: numpy array
:param indexDimensions: see documentation of Field
"""
spatialDimensions = len(npArray.shape) - indexDimensions
if spatialDimensions < 1:
raise ValueError("Too many index dimensions. At least one spatial dimension required")
fullLayout = getLayoutFromNumpyArray(npArray)
spatialLayout = tuple([i for i in fullLayout if i < spatialDimensions])
assert len(spatialLayout) == spatialDimensions
strides = tuple([s // np.dtype(npArray.dtype).itemsize for s in npArray.strides])
shape = tuple([int(s) for s in npArray.shape])
return Field(fieldName, npArray.dtype, spatialLayout, shape, strides)
@staticmethod
def createGeneric(fieldName, spatialDimensions, dtype=np.float64, indexDimensions=0, layout='numpy'):
"""
Creates a generic field where the field size is not fixed i.e. can be called with arrays of different sizes
:param fieldName: symbolic name for the field
:param dtype: numpy data type of the array the kernel is called with later
:param spatialDimensions: see documentation of Field
:param indexDimensions: see documentation of Field
:param layout: tuple specifying the loop ordering of the spatial dimensions e.g. (2, 1, 0 ) means that
the outer loop loops over dimension 2, the second outer over dimension 1, and the inner loop
over dimension 0. Also allowed: the strings 'numpy' (0,1,..d) or 'reverseNumpy' (d, ..., 1, 0)
"""
if layout == 'numpy' or layout.lower() == 'c':
layout = tuple(range(spatialDimensions))
elif layout == 'reverseNumpy' or layout.lower() == 'f':
layout = tuple(reversed(range(spatialDimensions)))
if len(layout) != spatialDimensions:
raise ValueError("Layout")
shapeSymbol = IndexedBase(TypedSymbol(Field.SHAPE_PREFIX + fieldName, Field.SHAPE_DTYPE), shape=(1,))
strideSymbol = IndexedBase(TypedSymbol(Field.STRIDE_PREFIX + fieldName, Field.STRIDE_DTYPE), shape=(1,))
totalDimensions = spatialDimensions + indexDimensions
shape = tuple([shapeSymbol[i] for i in range(totalDimensions)])
strides = tuple([strideSymbol[i] for i in range(totalDimensions)])
return Field(fieldName, dtype, layout, shape, strides)
def __init__(self, fieldName, dtype, layout, shape, strides):
"""Do not use directly. Use static create* methods"""
self._fieldName = fieldName
self._dtype = numpyDataTypeToC(dtype)
self._layout = layout
self.shape = shape
self.strides = strides
@property
def spatialDimensions(self):
return len(self._layout)
@property
def indexDimensions(self):
return len(self.shape) - len(self._layout)
@property
def layout(self):
return self._layout
@property
def name(self):
return self._fieldName
@property
def spatialShape(self):
return self.shape[:self.spatialDimensions]
@property
def indexShape(self):
return self.shape[self.spatialDimensions:]
@property
def spatialStrides(self):
return self.strides[:self.spatialDimensions]
@property
def indexStrides(self):
return self.strides[self.spatialDimensions:]
@property
def dtype(self):
return self._dtype
def __repr__(self):
return self._fieldName
def neighbor(self, coordId, offset):
offsetList = [0] * self.spatialDimensions
offsetList[coordId] = offset
return Field.Access(self, tuple(offsetList))
def __getitem__(self, offset):
if type(offset) is np.ndarray:
offset = tuple(offset)
if type(offset) is str:
offset = tuple(directionStringToOffset(offset, self.spatialDimensions))
if type(offset) is not tuple:
offset = (offset,)
if len(offset) != self.spatialDimensions:
raise ValueError("Wrong number of spatial indices: "
"Got %d, expected %d" % (len(offset), self.spatialDimensions))
return Field.Access(self, offset)
def __call__(self, *args, **kwargs):
center = tuple([0]*self.spatialDimensions)
return Field.Access(self, center)(*args, **kwargs)
def __hash__(self):
return hash((self._layout, self.shape, self.strides, self._dtype, self._fieldName))
def __eq__(self, other):
selfTuple = (self.shape, self.strides, self.name, self.dtype)
otherTuple = (other.shape, other.strides, other.name, other.dtype)
return selfTuple == otherTuple
PREFIX = "f"
STRIDE_PREFIX = PREFIX + "stride_"
SHAPE_PREFIX = PREFIX + "shape_"
STRIDE_DTYPE = "const int *"
SHAPE_DTYPE = "const int *"
DATA_PREFIX = PREFIX + "d_"
class Access(sp.Symbol):
def __new__(cls, name, *args, **kwargs):
obj = Field.Access.__xnew_cached_(cls, name, *args, **kwargs)
return obj
def __new_stage2__(self, field, offsets=(0, 0, 0), idx=None):
fieldName = field.name
offsetsAndIndex = chain(offsets, idx) if idx is not None else offsets
constantOffsets = not any([isinstance(o, sp.Basic) for o in offsetsAndIndex])
if not idx:
idx = tuple([0] * field.indexDimensions)
if constantOffsets:
offsetName = offsetToDirectionString(offsets)
if field.indexDimensions == 0:
obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName)
elif field.indexDimensions == 1:
obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName + "^" + str(idx[0]))
else:
idxStr = ",".join([str(e) for e in idx])
obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName + "^" + idxStr)
else:
offsetName = "%0.10X" % (abs(hash(tuple(offsetsAndIndex))))
obj = super(Field.Access, self).__xnew__(self, fieldName + "_" + offsetName)
obj._field = field
obj._offsets = []
for o in offsets:
if isinstance(o, sp.Basic):
obj._offsets.append(o)
else:
obj._offsets.append(int(o))
obj._offsetName = offsetName
obj._index = idx
return obj
def __getnewargs__(self):
return self.name, self.offsets, self.index
__xnew__ = staticmethod(__new_stage2__)
__xnew_cached_ = staticmethod(cacheit(__new_stage2__))
def __call__(self, *idx):
if self._index != tuple([0]*self.field.indexDimensions):
print(self._index, tuple([0]*self.field.indexDimensions))
raise ValueError("Indexing an already indexed Field.Access")
idx = tuple(idx)
if len(idx) != self.field.indexDimensions and idx != (0,):
raise ValueError("Wrong number of indices: "
"Got %d, expected %d" % (len(idx), self.field.indexDimensions))
return Field.Access(self.field, self._offsets, idx)
@property
def field(self):
return self._field
@property
def offsets(self):
return self._offsets
@property
def requiredGhostLayers(self):
return int(np.max(np.abs(self._offsets)))
@property
def nrOfCoordinates(self):
return len(self._offsets)
@property
def offsetName(self):
return self._offsetName
@property
def index(self):
return self._index
def _hashable_content(self):
superClassContents = list(super(Field.Access, self)._hashable_content())
t = tuple(superClassContents + [hash(self._field), self._index] + self._offsets)
return t
def extractCommonSubexpressions(equations):
"""
Uses sympy to find common subexpressions in equations and returns
them in a topologically sorted order, ready for evaluation.
Usually called before list of equations is passed to :func:`createKernel`
"""
replacements, newEq = sp.cse(equations)
# Workaround for older sympy versions: here subexpressions (temporary = True) are extracted
# which leads to problems in Piecewise functions which have to a default case indicated by True
symbolsEqualToTrue = {r[0]: True for r in replacements if r[1] is sp.true}
replacementEqs = [sp.Eq(*r) for r in replacements if r[1] is not sp.true]
equations = replacementEqs + newEq
topologicallySortedPairs = sp.cse_main.reps_toposort([[e.lhs, e.rhs] for e in equations])
equations = [sp.Eq(a[0], a[1].subs(symbolsEqualToTrue)) for a in topologicallySortedPairs]
return equations
def getLayoutFromNumpyArray(arr):
"""
Returns a list indicating the memory layout (linearization order) of the numpy array.
Example:
>>> getLayoutFromNumpyArray(np.zeros([3,3,3]))
[0, 1, 2]
In this example the loop over the zeroth coordinate should be the outermost loop,
followed by the first and second. Elements arr[x,y,0] and arr[x,y,1] are adjacent in memory.
Normally constructed numpy arrays have this order, however by stride tricks or other frameworks, arrays
with different memory layout can be created.
"""
coordinates = list(range(len(arr.shape)))
return [x for (y, x) in sorted(zip(arr.strides, coordinates), key=lambda pair: pair[0], reverse=True)]
def numpyDataTypeToC(dtype):
"""Mapping numpy data types to C data types"""
if dtype == np.float64:
return "double"
elif dtype == np.float32:
return "float"
elif dtype == np.int32:
return "int"
raise NotImplementedError()
def offsetComponentToDirectionString(coordinateId, value):
"""
Translates numerical offset to string notation.
x offsets are labeled with east 'E' and 'W',
y offsets with north 'N' and 'S' and
z offsets with top 'T' and bottom 'B'
If the absolute value of the offset is bigger than 1, this number is prefixed.
:param coordinateId: integer 0, 1 or 2 standing for x,y and z
:param value: integer offset
Example:
>>> offsetComponentToDirectionString(0, 1)
'E'
>>> offsetComponentToDirectionString(1, 2)
'2N'
"""
nameComponents = (('W', 'E'), # west, east
('S', 'N'), # south, north
('B', 'T'), # bottom, top
)
if value == 0:
result = ""
elif value < 0:
result = nameComponents[coordinateId][0]
else:
result = nameComponents[coordinateId][1]
if abs(value) > 1:
result = "%d%s" % (abs(value), result)
return result
def offsetToDirectionString(offsetTuple):
"""
Translates numerical offset to string notation.
For details see :func:`offsetComponentToDirectionString`
:param offsetTuple: 3-tuple with x,y,z offset
Example:
>>> offsetToDirectionString([1, -1, 0])
'SE'
>>> offsetToDirectionString(([-3, 0, -2]))
'2B3W'
"""
names = ["", "", ""]
for i in range(len(offsetTuple)):
names[i] = offsetComponentToDirectionString(i, offsetTuple[i])
name = "".join(reversed(names))
if name == "":
name = "C"
return name
def directionStringToOffset(directionStr, dim=3):
"""
Reverse mapping of :func:`offsetToDirectionString`
:param directionStr: string representation of offset
:param dim: dimension of offset, i.e the length of the returned list
>>> directionStringToOffset('NW', dim=3)
array([-1, 1, 0])
>>> directionStringToOffset('NW', dim=2)
array([-1, 1])
>>> directionStringToOffset(offsetToDirectionString([3,-2,1]))
array([ 3, -2, 1])
"""
offsetMap = {
'C': np.array([0, 0, 0]),
'W': np.array([-1, 0, 0]),
'E': np.array([1, 0, 0]),
'S': np.array([0, -1, 0]),
'N': np.array([0, 1, 0]),
'B': np.array([0, 0, -1]),
'T': np.array([0, 0, 1]),
}
offset = np.array([0, 0, 0])
while len(directionStr) > 0:
factor = 1
firstNonDigit = 0
while directionStr[firstNonDigit].isdigit():
firstNonDigit += 1
if firstNonDigit > 0:
factor = int(directionStr[:firstNonDigit])
directionStr = directionStr[firstNonDigit:]
curOffset = offsetMap[directionStr[0]]
offset += factor * curOffset
directionStr = directionStr[1:]
return offset[:dim]
finitedifferences.py deleted 100644 → 0
View file @ cfef0f9d
import numpy as np
import sympy as sp
from pystencils.field import Field
from pystencils.transformations import fastSubs
def grad(var, dim=3):
r"""
Gradients are represented as a special symbol:
e.g. :math:`\nabla x = (x^{\Delta 0}, x^{\Delta 1}, x^{\Delta 2})`
This function takes a symbol and creates the gradient symbols according to convention above
:param var: symbol to take the gradient of
:param dim: dimension (length) of the gradient vector
"""
if hasattr(var, "__getitem__"):
return [[sp.Symbol("%s^Delta^%d" % (v.name, i)) for v in var] for i in range(dim)]
else:
return [sp.Symbol("%s^Delta^%d" % (var.name, i)) for i in range(dim)]
def discretizeCenter(term, symbolsToFieldDict, dx, dim=3):
"""
Expects term that contains given symbols and gradient components of these symbols and replaces them
by field accesses. Gradients are replaced centralized approximations:
``(upper neighbor - lower neighbor ) / ( 2*dx)``
:param term: term where symbols and gradient(symbol) should be replaced
:param symbolsToFieldDict: mapping of symbols to Field
:param dx: width and height of one cell
:param dim: dimension
Example:
>>> x = sp.Symbol("x")
>>> gradx = grad(x, dim=3)
>>> term = x * gradx[0]
>>> term
x*x^Delta^0
>>> f = Field.createGeneric('f', spatialDimensions=3)
>>> discretizeCenter(term, { x: f }, dx=1, dim=3)
f_C*(f_E/2 - f_W/2)
"""
substitutions = {}
for symbols, field in symbolsToFieldDict.items():
if not hasattr(symbols, "__getitem__"):
symbols = [symbols]
g = grad(symbols, dim)
substitutions.update({symbol: field(i) for i, symbol in enumerate(symbols)})
for d in range(dim):
up, down = __upDownOffsets(d, dim)
substitutions.update({g[d][i]: (field[up](i) - field[down](i)) / dx / 2 for i in range(len(symbols))})
return term.subs(substitutions)
def discretizeStaggered(term, symbolsToFieldDict, coordinate, coordinateOffset, dx, dim=3):
"""
Expects term that contains given symbols and gradient components of these symbols and replaces them
by field accesses. Gradients in coordinate direction are replaced by staggered version at cell boundary.
Symbols themselves and gradients in other directions are replaced by interpolated version at cell face.
:param term: input term where symbols and gradients are replaced
:param symbolsToFieldDict: mapping of symbols to Field
:param coordinate: id for coordinate (0 for x, 1 for y, ... ) defining cell boundary.
Only gradients in this direction are replaced e.g. if symbol^Delta^coordinate
:param coordinateOffset: either +1 or -1 for upper or lower face in coordinate direction
:param dx: width and height of one cell
:param dim: dimension
Example: Discretizing at right/east face of cell i.e. coordinate=0, offset=1)
>>> x, dx = sp.symbols("x dx")
>>> gradx = grad(x, dim=3)
>>> term = x * gradx[0]
>>> term
x*x^Delta^0
>>> f = Field.createGeneric('f', spatialDimensions=3)
>>> discretizeStaggered(term, symbolsToFieldDict={ x: f}, dx=dx, coordinate=0, coordinateOffset=1, dim=3)
(-f_C + f_E)*(f_C/2 + f_E/2)/dx
"""
assert coordinateOffset == 1 or coordinateOffset == -1
assert 0 <= coordinate < dim
substitutions = {}
for symbols, field in symbolsToFieldDict.items():
if not hasattr(symbols, "__getitem__"):
symbols = [symbols]
offset = [0] * dim
offset[coordinate] = coordinateOffset
offset = np.array(offset, dtype=np.int)
gradient = grad(symbols)[coordinate]
substitutions.update({s: (field[offset](i) + field(i)) / 2 for i, s in enumerate(symbols)})
substitutions.update({g: (field[offset](i) - field(i)) / dx * coordinateOffset for i, g in enumerate(gradient)})
for d in range(dim):
if d == coordinate:
continue
up, down = __upDownOffsets(d, dim)
for i, s in enumerate(symbols):
centerGrad = (field[up](i) - field[down](i)) / (2 * dx)
neighborGrad = (field[up+offset](i) - field[down+offset](i)) / (2 * dx)
substitutions[grad(s)[d]] = (centerGrad + neighborGrad) / 2
return fastSubs(term, substitutions)
def discretizeDivergence(vectorTerm, symbolsToFieldDict, dx):
"""
Computes discrete divergence of symbolic vector
:param vectorTerm: sequence of terms, interpreted as vector
:param symbolsToFieldDict: mapping of symbols to Field
:param dx: length of a cell
Example: Laplace stencil
>>> x, dx = sp.symbols("x dx")
>>> gradX = grad(x, dim=3)
>>> f = Field.createGeneric('f', spatialDimensions=3)
>>> sp.simplify(discretizeDivergence(gradX, {x : f}, dx))
(f_B - 6*f_C + f_E + f_N + f_S + f_T + f_W)/dx
"""
dim = len(vectorTerm)
result = 0
for d in range(dim):
for offset in [-1, 1]:
result += offset * discretizeStaggered(vectorTerm[d], symbolsToFieldDict, d, offset, dx, dim)
return result
def __upDownOffsets(d, dim):
coord = [0] * dim
coord[d] = 1
up = np.array(coord, dtype=np.int)
coord[d] = -1
down = np.array(coord, dtype=np.int)
return up, down
gpucuda/__init__.py deleted 100644 → 0
View file @ cfef0f9d
from pystencils.gpucuda.kernelcreation import createCUDAKernel
from pystencils.gpucuda.cudajit import makePythonFunction
from pystencils.backends.cbackend import generateCUDA
\ No newline at end of file
gpucuda/cudajit.py deleted 100644 → 0
View file @ cfef0f9d
import numpy as np
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
from pycuda.gpuarray import GPUArray
from pystencils.backends.cbackend import generateCUDA
from pystencils.transformations import symbolNameToVariableName
def numpyTypeFromString(typename, includePointers=True):
import ctypes as ct
typename = typename.replace("*", " * ")
typeComponents = typename.split()
basicTypeMap = {
'double': np.float64,
'float': np.float32,
'int': np.int32,
'long': np.int64,
}
resultType = None
for typeComponent in typeComponents:
typeComponent = typeComponent.strip()
if typeComponent == "const" or typeComponent == "restrict" or typeComponent == "volatile":
continue
if typeComponent in basicTypeMap:
resultType = basicTypeMap[typeComponent]
elif typeComponent == "*" and includePointers:
assert resultType is not None
resultType = ct.POINTER(resultType)
return resultType
def buildNumpyArgumentList(kernelFunctionNode, argumentDict):
argumentDict = {symbolNameToVariableName(k): v for k, v in argumentDict.items()}
result = []
for arg in kernelFunctionNode.parameters:
if arg.isFieldArgument:
field = argumentDict[arg.fieldName]
if arg.isFieldPtrArgument:
result.append(field.gpudata)
elif arg.isFieldShapeArgument:
strideArr = np.array(field.strides, dtype=np.int32) / field.dtype.itemsize
result.append(cuda.In(strideArr))
elif arg.isFieldStrideArgument:
shapeArr = np.array(field.shape, dtype=np.int32)
result.append(cuda.In(shapeArr))
else:
assert False
else:
param = argumentDict[arg.name]
expectedType = numpyTypeFromString(arg.dtype)
result.append(expectedType(param))
return result
def makePythonFunction(kernelFunctionNode, argumentDict={}):
mod = SourceModule(str(generateCUDA(kernelFunctionNode)))
func = mod.get_function(kernelFunctionNode.functionName)
def wrapper(**kwargs):
from copy import copy
fullArguments = copy(argumentDict)
fullArguments.update(kwargs)
shapes = set()
strides = set()
for argValue in fullArguments.values():
if isinstance(argValue, GPUArray):
shapes.add(argValue.shape)
strides.add(argValue.strides)
if len(strides) == 0:
raise ValueError("No GPU arrays passed as argument")
assert len(strides) < 2, "All passed arrays have to have the same strides"
assert len(shapes) < 2, "All passed arrays have to have the same size"
shape = list(shapes)[0]
dictWithBlockAndThreadNumbers = kernelFunctionNode.getCallParameters(shape)
args = buildNumpyArgumentList(kernelFunctionNode, fullArguments)
func(*args, **dictWithBlockAndThreadNumbers)
return wrapper
gpucuda/kernelcreation.py deleted 100644 → 0
View file @ cfef0f9d
import sympy as sp
from pystencils.transformations import resolveFieldAccesses, typeAllEquations, \
parseBasePointerInfo, typingFromSympyInspection
from pystencils.ast import Block, KernelFunction
from pystencils import Field
BLOCK_IDX = list(sp.symbols("blockIdx.x blockIdx.y blockIdx.z"))
THREAD_IDX = list(sp.symbols("threadIdx.x threadIdx.y threadIdx.z"))
def getLinewiseCoordinates(field, ghostLayers):
availableIndices = [THREAD_IDX[0]] + BLOCK_IDX
assert field.spatialDimensions <= 4, "This indexing scheme supports at most 4 spatial dimensions"
result = availableIndices[:field.spatialDimensions]
fastestCoordinate = field.layout[-1]
result[0], result[fastestCoordinate] = result[fastestCoordinate], result[0]
def getCallParameters(arrShape):
def getShapeOfCudaIdx(cudaIdx):
if cudaIdx not in result:
return 1
else:
return arrShape[result.index(cudaIdx)] - 2 * ghostLayers
return {'block': tuple([getShapeOfCudaIdx(idx) for idx in THREAD_IDX]),
'grid': tuple([getShapeOfCudaIdx(idx) for idx in BLOCK_IDX]) }
return [i + ghostLayers for i in result], getCallParameters
def createCUDAKernel(listOfEquations, functionName="kernel", typeForSymbol=None):
if not typeForSymbol or typeForSymbol == 'double':
typeForSymbol = typingFromSympyInspection(listOfEquations, "double")
elif typeForSymbol == 'float':
typeForSymbol = typingFromSympyInspection(listOfEquations, "float")
fieldsRead, fieldsWritten, assignments = typeAllEquations(listOfEquations, typeForSymbol)
readOnlyFields = set([f.name for f in fieldsRead - fieldsWritten])
allFields = fieldsRead.union(fieldsWritten)
code = KernelFunction(Block(assignments), fieldsRead.union(fieldsWritten), functionName)
code.globalVariables.update(BLOCK_IDX + THREAD_IDX)
fieldAccesses = code.atoms(Field.Access)
requiredGhostLayers = max([fa.requiredGhostLayers for fa in fieldAccesses])
coordMapping, getCallParameters = getLinewiseCoordinates(list(fieldsRead)[0], requiredGhostLayers)
allFields = fieldsRead.union(fieldsWritten)
basePointerInfo = [['spatialInner0']]
basePointerInfos = {f.name: parseBasePointerInfo(basePointerInfo, [2, 1, 0], f) for f in allFields}
coordMapping = {f.name: coordMapping for f in allFields}
resolveFieldAccesses(code, readOnlyFields, fieldToFixedCoordinates=coordMapping,
fieldToBasePointerInfo=basePointerInfos)
# add the function which determines #blocks and #threads as additional member to KernelFunction node
# this is used by the jit
code.getCallParameters = getCallParameters
return code
llvm/__init__.py deleted 100644 → 0
View file @ cfef0f9d
from .kernelcreation import createKernel
llvm/control_flow.py deleted 100644 → 0
View file @ cfef0f9d
import llvmlite.ir as ir
class Loop:
def __init__(self, builder, start_val, stop_val, step_val=1, loop_name='loop', phi_name="_phi"):
self.builder = builder
self.start_val = start_val
self.stop_val = stop_val
self.step_val = step_val
self.loop_name = loop_name
self.phi_name = phi_name
def __enter__(self):
self.loop_end, self.after, phi = self._for_loop(self.start_val, self.stop_val, self.step_val, self.loop_name,
self.phi_name)
return phi
def _for_loop(self, start_val, stop_val, step_val, loop_name, phi_name):
# TODO size of int??? unisgned???
integer = ir.IntType(64)
# Loop block
pre_loop_bb = self.builder.block
loop_bb = self.builder.append_basic_block(name='loop_' + loop_name)
self.builder.branch(loop_bb)
# Insert an explicit fall through from the current block to loop_bb
self.builder.position_at_start(loop_bb)
# Add phi
phi = self.builder.phi(integer, name=phi_name)
phi.add_incoming(start_val, pre_loop_bb)
loop_end_bb = self.builder.append_basic_block(name=loop_name + "_end_bb")
self.builder.position_at_start(loop_end_bb)
next_var = self.builder.add(phi, step_val, name=loop_name + '_next_it')
cond = self.builder.icmp_unsigned('<', next_var, stop_val, name=loop_name + "_cond")
after_bb = self.builder.append_basic_block(name=loop_name + "_after_bb")
self.builder.cbranch(cond, loop_bb, after_bb)
phi.add_incoming(next_var, loop_end_bb)
self.builder.position_at_end(loop_bb)
return loop_end_bb, after_bb, phi
def __exit__(self, exc_type, exc, exc_tb):
self.builder.branch(self.loop_end)
self.builder.position_at_end(self.after)
llvm/jit.py deleted 100644 → 0
View file @ cfef0f9d
import llvmlite.binding as llvm
import logging.config
class Eval(object):
def __init__(self):
llvm.initialize()
llvm.initialize_all_targets()
llvm.initialize_native_target()
llvm.initialize_native_asmprinter()
self.target = llvm.Target.from_default_triple()
def compile(self, module):
logger.debug('=============Preparse')
logger.debug(str(module))
llvmmod = llvm.parse_assembly(str(module))
llvmmod.verify()
logger.debug('=============Function in IR')
logger.debug(str(llvmmod))
# TODO cpu, features, opt
cpu = llvm.get_host_cpu_name()
features = llvm.get_host_cpu_features()
logger.debug('=======Things')
logger.debug(cpu)
logger.debug(features.flatten())
target_machine = self.target.create_target_machine(cpu=cpu, features=features.flatten(), opt=2)
logger.debug('Machine = ' + str(target_machine.target_data))
with open('gen.ll', 'w') as f:
f.write(str(llvmmod))
optimize = True
if optimize:
pmb = llvm.create_pass_manager_builder()
pmb.opt_level = 2
pmb.disable_unit_at_a_time = False
pmb.loop_vectorize = True
pmb.slp_vectorize = True
# TODO possible to pass for functions
pm = llvm.create_module_pass_manager()
pm.add_instruction_combining_pass()
pm.add_function_attrs_pass()
pm.add_constant_merge_pass()
pm.add_licm_pass()
pmb.populate(pm)
pm.run(llvmmod)
logger.debug("==========Opt")
logger.debug(str(llvmmod))
with open('gen_opt.ll', 'w') as f:
f.write(str(llvmmod))
with llvm.create_mcjit_compiler(llvmmod, target_machine) as ee:
ee.finalize_object()
logger.debug('==========Machine code')
logger.debug(target_machine.emit_assembly(llvmmod))
with open('gen.S', 'w') as f:
f.write(target_machine.emit_assembly(llvmmod))
with open('gen.o', 'wb') as f:
f.write(target_machine.emit_object(llvmmod))
# fptr = CFUNCTYPE(c_double, c_double, c_double)(ee.get_function_address('add2'))
# result = fptr(2, 3)
# print(result)
return 0
if __name__ == "__main__":
logger = logging.getLogger(__name__)
else:
logger = logging.getLogger(__name__)
llvm/kernelcreation.py deleted 100644 → 0
View file @ cfef0f9d
import sympy as sp
from pystencils.transformations import resolveFieldAccesses, makeLoopOverDomain, typingFromSympyInspection, \
typeAllEquations, getOptimalLoopOrdering, parseBasePointerInfo, moveConstantsBeforeLoop, splitInnerLoop
from pystencils.types import TypedSymbol, DataType
from pystencils.field import Field
import pystencils.ast as ast
def createKernel(listOfEquations, functionName="kernel", typeForSymbol=None, splitGroups=(),
iterationSlice=None, ghostLayers=None):
"""
Creates an abstract syntax tree for a kernel function, by taking a list of update rules.
Loops are created according to the field accesses in the equations.
:param listOfEquations: list of sympy equations, containing accesses to :class:`pystencils.field.Field`.
Defining the update rules of the kernel
:param functionName: name of the generated function - only important if generated code is written out
:param typeForSymbol: a map from symbol name to a C type specifier. If not specified all symbols are assumed to
be of type 'double' except symbols which occur on the left hand side of equations where the
right hand side is a sympy Boolean which are assumed to be 'bool' .
:param splitGroups: Specification on how to split up inner loop into multiple loops. For details see
transformation :func:`pystencils.transformation.splitInnerLoop`
:param iterationSlice: if not None, iteration is done only over this slice of the field
:param ghostLayers: a sequence of pairs for each coordinate with lower and upper nr of ghost layers
if None, the number of ghost layers is determined automatically and assumed to be equal for a
all dimensions
:return: :class:`pystencils.ast.KernelFunction` node
"""
if not typeForSymbol:
typeForSymbol = typingFromSympyInspection(listOfEquations, "double")
def typeSymbol(term):
if isinstance(term, Field.Access) or isinstance(term, TypedSymbol):
return term
elif isinstance(term, sp.Symbol):
return TypedSymbol(term.name, DataType(typeForSymbol[term.name]))
else:
raise ValueError("Term has to be field access or symbol")
fieldsRead, fieldsWritten, assignments = typeAllEquations(listOfEquations, typeForSymbol)
allFields = fieldsRead.union(fieldsWritten)
readOnlyFields = set([f.name for f in fieldsRead - fieldsWritten])
body = ast.Block(assignments)
loopOrder = getOptimalLoopOrdering(allFields)
code = makeLoopOverDomain(body, functionName, iterationSlice=iterationSlice,
ghostLayers=ghostLayers, loopOrder=loopOrder)
if splitGroups:
typedSplitGroups = [[typeSymbol(s) for s in splitGroup] for splitGroup in splitGroups]
splitInnerLoop(code, typedSplitGroups)
basePointerInfo = [['spatialInner0'], ['spatialInner1']]
basePointerInfos = {field.name: parseBasePointerInfo(basePointerInfo, loopOrder, field) for field in allFields}
resolveFieldAccesses(code, readOnlyFields, fieldToBasePointerInfo=basePointerInfos)
moveConstantsBeforeLoop(code)
return code
\ No newline at end of file
pyproject.toml 0 → 100644
View file @ 5600b6b6
[project]
name = "pystencils"
description = "Speeding up stencil computations on CPUs and GPUs"
dynamic = ["version"]
readme = "README.md"
authors = [
{ name = "Martin Bauer" },
{ name = "Jan Hönig " },
{ name = "Markus Holzer" },
{ name = "Frederik Hennig" },
{ email = "cs10-codegen@fau.de" },
]
license = { file = "COPYING.txt" }
requires-python = ">=3.10"
dependencies = ["sympy>=1.9,<=1.12.1", "numpy>=1.8.0", "appdirs", "joblib", "pyyaml", "fasteners"]
classifiers = [
"Development Status :: 4 - Beta",
"Framework :: Jupyter",
"Topic :: Software Development :: Code Generators",
"Topic :: Scientific/Engineering :: Physics",
"Intended Audience :: Developers",
"Intended Audience :: Science/Research",
"License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)",
]
[project.urls]
"Bug Tracker" = "https://i10git.cs.fau.de/pycodegen/pystencils/-/issues"
"Documentation" = "https://pycodegen.pages.i10git.cs.fau.de/pystencils/"
"Source Code" = "https://i10git.cs.fau.de/pycodegen/pystencils"
[project.optional-dependencies]
gpu = ['cupy']
alltrafos = ['islpy', 'py-cpuinfo']
bench_db = ['blitzdb', 'pymongo', 'pandas']
interactive = [
'matplotlib',
'ipy_table',
'imageio',
'jupyter',
'pyevtk',
'rich',
'graphviz',
]
use_cython = [
'Cython'
]
doc = [
'sphinx',
'sphinx_rtd_theme',
'nbsphinx',
'sphinxcontrib-bibtex',
'sphinx_autodoc_typehints',
'pandoc',
]
tests = [
'pytest',
'pytest-cov',
'pytest-html',
'ansi2html',
'pytest-xdist',
'flake8',
'nbformat',
'nbconvert',
'ipython',
'matplotlib',
'py-cpuinfo',
'randomgen>=1.18',
]
[build-system]
requires = [
"setuptools>=61",
"versioneer[toml]>=0.29",
# 'Cython'
]
build-backend = "setuptools.build_meta"
[tool.setuptools.package-data]
pystencils = [
"include/*.h",
"boundaries/createindexlistcython.pyx"
]
[tool.setuptools.packages.find]
where = ["src"]
include = ["pystencils", "pystencils.*"]
namespaces = false
[tool.versioneer]
# See the docstring in versioneer.py for instructions. Note that you must
# re-run 'versioneer.py setup' after changing this section, and commit the
# resulting files.
VCS = "git"
style = "pep440"
versionfile_source = "src/pystencils/_version.py"
versionfile_build = "pystencils/_version.py"
tag_prefix = "release/"
parentdir_prefix = "pystencils-"
pytest.ini 0 → 100644
View file @ 5600b6b6
[pytest]
testpaths = src tests doc/notebooks
pythonpath = src
python_files = test_*.py *_test.py scenario_*.py
norecursedirs = *.egg-info .git .cache .ipynb_checkpoints htmlcov
addopts = --doctest-modules --durations=20 --cov-config pytest.ini
markers =
longrun: tests only run at night since they have large execution time
notebook: mark for notebooks
# these warnings all come from third party libraries.
filterwarnings =
ignore:an integer is required:DeprecationWarning
ignore:\s*load will be removed, use:PendingDeprecationWarning
ignore:the imp module is deprecated in favour of importlib:DeprecationWarning
ignore:.*is a deprecated alias for the builtin `bool`:DeprecationWarning
ignore:'contextfilter' is renamed to 'pass_context':DeprecationWarning
ignore:Using or importing the ABCs from 'collections' instead of from 'collections.abc':DeprecationWarning
ignore:Animation was deleted without rendering anything:UserWarning
[run]
branch = True
source = src/pystencils
tests
omit = doc/*
tests/*
setup.py
quicktest.py
conftest.py
versioneer.py
src/pystencils/jupytersetup.py
src/pystencils/cpu/msvc_detection.py
src/pystencils/sympy_gmpy_bug_workaround.py
src/pystencils/cache.py
src/pystencils/pacxx/benchmark.py
src/pystencils/_version.py
venv/
[report]
exclude_lines =
# Have to re-enable the standard pragma
pragma: no cover
def __repr__
def _repr_html_
# Don't complain if tests don't hit defensive assertion code:
raise AssertionError
raise NotImplementedError
NotImplementedError()
#raise ValueError
# Don't complain if non-runnable code isn't run:
if 0:
if False:
if __name__ == .__main__.:
skip_covered = True
fail_under = 85
[html]
directory = coverage_report
quicktest.py 0 → 100644
View file @ 5600b6b6
#!/usr/bin/env python3
from contextlib import redirect_stdout
import io
from tests.test_quicktests import (
test_basic_kernel,
test_basic_blocking_staggered,
test_basic_vectorization,
)
quick_tests = [
test_basic_kernel,
test_basic_blocking_staggered,
test_basic_vectorization,
]
if __name__ == "__main__":
print("Running pystencils quicktests")
for qt in quick_tests:
print(f" -> {qt.__name__}")
with redirect_stdout(io.StringIO()):
qt()
release.sh 0 → 100644
View file @ 5600b6b6
#!/bin/bash
echo "Existing versions"
git tag -l | grep release
echo "Enter the next version"
read new_version
git tag -s release/${new_version}
git push origin master release/${new_version}
rm -rf dist
python setup.py sdist
twine upload dist/*