Non-rectangular arrays and data repacking

Download 410,71 Kb. Pdf ko'rish
bet	8/12
Sana	16.11.2023
Hajmi	410,71 Kb.
	#100091

1 ... 4 5 6 7 8 9 10 11 12

Bog'liq
cython cise

Non-rectangular arrays and data repacking
Sometimes data does not fit naturally in rectangu-
lar arrays, and Cython is especially well-suited to
this situation. One such example arises in cosmol-
ogy. Satellite experiments such as the Wilkinson
Microwave Anisotropy Probe have produced high-
resolution images of the Cosmic Microwave Back-
ground, a primary source of information about the
early universe. The resulting images are spheri-
cal, as they contain values for all directions on the
sky.
The spherical harmonic transform of these maps,
a “fourier transform on the sphere”, is especially
important. It has complex coefficients a
`m
where
the indices run over 0 ≤ ` ≤ `
max
, −` ≤ m ≤ `.
An average of the entire map is stored in a
0,0
,
followed by three elements to describe the dipole
component, a
1,−1
, a
1,0
, a
1,1
, and so on. Data like
this can be stored in a one-dimensional array and
elements looked up at position `
2
+ ` + m.
It is possible, but not trivial, to operate on such
data using NumPy whole-array operations. The
problem is that NumPy functions, such as find-
ing the variance, are primarily geared towards
rectangular arrays. If the data was rectangular,
one could estimate the variance per `, averaging
over m, by calling np.var(data, axis=1). This
doesn’t work for non-rectangular data.
While
there are workarounds, such as the reduceat
method and masked arrays, we have found it
much more straightforward to write the obvious
loops over ` and m using Cython. For compar-
ison, with Python and NumPy one could loop
over ` and call repeatedly call np.var for sub-
slices of the data, which was 27 times slower in
our case (`
max
= 1500). Using a naive double
loop over both ` and m was more than a 1000
times slower in Python than in Cython. (Inciden-
tally, the variance per `, or power spectrum, is
the primary quantity of interest to observational
cosmologists.)
The spherical harmonic transform mentioned
above is computed using the Fortran library
HEALPix
3
, which can readily be called from
Cython with the help of Fwrap.
However,
HEALPix spits out the result as a 2D array, with
roughly half of the elements unoccupied.
The
waste of storage aside, 2D arrays are often incon-
venient – with 1D arrays one can treat each set of
coefficients as a vector, and perform linear alge-
bra, estimate covariance matrices and so on, the
usual way. Again, it is possible to quickly reorder
the data the way we want it with a Cython loop.
With all the existing code out there wanting data
in slightly different order and formats, for loops
are not about to disappear.
Fwrap
Whereas C and C++ integrate closely with
Cython, Fortran wrappers in Cython are gener-
ated with Fwrap, a separate utility that is dis-
tributed separately from Cython. Fwrap
[fwrap]
is a tool that automates wrapping Fortran source
in C, Cython and Python, allowing Fortran code
to benefit from the dynamism and flexibility of
Python. Fwrapped code can be seamlessly inte-
grated into a C, Cython or Python project. The
utility transparently supports most of the features
introduced in Fortran 90/95/2003, and will han-
dle nearly all Fortran 77 source as well. Fwrap
does not currently support derived types or func-
tion callbacks, but support for these features is
scheduled in an upcoming release.
Thanks to the C interoperability features supplied
in the Fortran 2003 standard – and supported in
recent versions of all widely-used Fortran 90/95
compilers – Fwrap generates wrappers that are
portable across platforms and compilers. Fwrap is
intended to be as friendly as possible, and handles
the Fortran parsing and generation automatically.
It also generates a build script for the project that
will portably build a Python extension module
from the wrapper files.
Fwrap is similar in intent to other Fortran-Python
3
Hierarchical Equal Area isoLatitude Pixelization,
G´
orski et al,
http://healpix.jpl.nasa.gov/
6
IEEE Computing in science and Engineering

tools such as F2PY, PyFort and Forthon. F2PY is
distributed with NumPy and is a capable tool for
wrapping Fortran 77 codes. Fwrap’s approach dif-
fers in that it leverages Cython to create Python
bindings. Manual tuning of the wrapper can be
easily accomplished by simply modifying the gen-
erated Cython code, rather than using a restricted
domain-specific language. Another benefit is re-
duced overhead when calling Fortran code from
Cython.
Consider a real world example: wrapping a sub-
routine from netlib’s LAPACK Fortran 90 source.
We will use the Fortran 90 subroutine interface
for dgesdd, used to compute the singular value de-
composition arrays U, S, and VT of a real array A,
such that A = U * DIAG(S) * VT. This routine
is typical of Fortran 90 source code – it has scalar
and array arguments with different intents and
different datatypes. We have augmented the ar-
gument declarations with INTENT attributes and
removed extraneous work array arguments for il-
lustration purposes:
SUBROUTINE DGESDD(JOBZ, M, N, A, LDA, S, &
& U, LDU, VT, LDVT, INFO)
! .. Scalar Arguments ..
CHARACTER, INTENT(IN) :: JOBZ
INTEGER, INTENT(OUT)
:: INFO
INTEGER, INTENT(IN)
:: LDA, LDU, LDVT &
& M, N
! .. Array Arguments ..
DOUBLE PRECISION, INTENT(INOUT) :: &
& A(LDA, *)
DOUBLE PRECISION, INTENT(OUT)
:: &
& S(*), U(LDU, *), VT(LDVT, *)
! DGESDD subroutine body
END SUBROUTINE DGESDD
When invoked on the above Fortran code, Fwrap
parses the code and makes it available to C,
Cython and Python. If desired, we can generate
a deployable package for use on computers
that don’t have Fwrap or Cython installed.
To use the wrapped code from Python, we
must first set up the subroutine arguments—in

Download 410,71 Kb.

1 ... 4 5 6 7 8 9 10 11 12

Download 410,71 Kb.

Pdf ko'rish