%-*- mode: Noweb; noweb-code-mode: python-mode -*- \section{Typedef} \label{sec:typedef} <>= cdef extern from "utilities.h": ctypedef double rtd ctypedef struct vector: rtd x rtd y rtd z void set_device(int) cdef extern from "cuda.h": ctypedef struct float2: float x float y @ <>= def setDevice(int id): set_device(id) @ \section{JSONAbstract} \label{sec:jsonabstract} \index{utilities!python!JSONAbstract} <>= import numpy as np import math from libc cimport math cimport numpy as np cimport utilities import os import json import boto3 import re import zipfile class JSONAbstract: def __init__(self, jprms = None, jsonfile = None, zipdata = None, cache = '', verbose = True): if jsonfile is not None: if (not os.path.isfile(jsonfile)) and (zipdata is not None): matches = re.match(r'(.*)://(.*)/(.*)/(.*.zip)',zipdata) if matches.group(1)=='s3': print("Downloading "+matches.group(4)+" from s3://"+matches.group(2)+"/"+matches.group(3)+" to "+cache+matches.group(4)+"...") s3 = boto3.resource('s3') s3.Bucket(matches.group(2)).download_file(matches.group(3)+"/"+matches.group(4),cache+matches.group(4)) try: with zipfile.ZipFile(cache+matches.group(4), "r") as z: print("Unzipping "+cache+matches.group(4)+"...") z.extractall(cache) except IOError as e: print("Unable to open %s!"%(cache+matches.group(4))) try: with open(jsonfile) as f: self.jprms = json.loads(f.read()) except IOError as e: print("Unable to open %s!"%(jsonfile)) else: self.jprms = jprms if verbose: print(self) def __str__(self): return json.dumps(self.jprms,indent=2, separators=(',', ': ')) def __getitem__(self, key): return self.jprms[key] @ \section{Sensor} \label{sec:sensor} This is an abstract class for all CEO sensors <>= cdef class Sensor: pass <>= cdef class Sensor: def calibrate(self,*args,**kwargs): pass def reset(self,**kwargs): pass def analyze(self,*args,**kwargs): pass def propagate(self,*args,**kwargs): pass def process(self,**kwargs): pass def get_measurement(self,**kwargs): pass def __pos__(self): pass def __invert__(self): pass property data: def __get__(self): return np.zeros(1) property Data: def __get__(self): return np.zeros((1,1)) @ \section{Stopwatch} \label{sec:stopwatch-1} The cython wrapper for stopwatch (Sec.~\ref{sec:stopwatch}) is: \index{utilities!python!StopWatch} \subsection{PXD file} \label{sec:pxd-file} <>= ## stopwatch structure cdef cppclass stopwatch: void tic() void toc() void toc(float *) @ \subsection{PYX file} \label{sec:pyx-file} <>= ## stopwatch structure cdef class StopWatch: """ Creates a stopwatch object to time cuda based routines Attributes ---------- elapsedTime : float The time elapsed between a tic() and a toc() calls """ cdef stopwatch *_c_stopwatch cdef public float elapsedTime def __cinit__(self): self._c_stopwatch = new stopwatch() def tic(self): """ Starts the stopwatch """ self._c_stopwatch.tic() def toc(self): """ Stops the stopwatch """ self._c_stopwatch.toc(&(self.elapsedTime)) def ping(self): """ Stops the stopwatch and display the elapsed time in ms """ self._c_stopwatch.toc(&(self.elapsedTime)) print "Elapsed time: %gms"%(self.elapsedTime) @ \section{Mask} \label{sec:mask} The cython wrapper for the mask structure (Sec.~\ref{sec:mask-structure}) is \subsection{PXD file} \label{sec:pxd-file-1} <>= ## mask structure cdef cppclass mask: int nel, nnz char *m float *f float area int *d__piston_mask void setup(int) void setup(int, float) void setup(int, float, int, int, int) void setup(float , float , float , float , float , float , float , int) void setup_circular(int, float) void setup_circular(int, float, float) void setup_GMT(int, float) void set_gmt_piston(float *, float *) void reset() void alter(float *) void add(char *, int) void cleanup() @ \subsubsection{Class definition} \label{sec:class-definition} <>= # MaskAbstract cdef class MaskAbstract: cdef mask *_c_mask cdef readonly int n cdef readonly cuFloatArray f cdef cuIntArray _piston_mask_ cdef class Mask(MaskAbstract): pass @ \subsection{PYX file} \label{sec:pyx-file-1} \subsubsection{MaskAbstract} \label{sec:maskabstract} \index{utilities!python!MaskAbstract} <>= ## mask structure cdef class MaskAbstract: """ The abstract class for all other mask classes Attributes ---------- f : cuFloatArray the mask as a float array nnz : int, readonly the number of non zeros values in the mask nel : int, readonly the total number of elements in the mask area : float, readonly the area covered by the mask See also -------- cuFloatArray : an interface for GPU host and device float data """ def __cinit__(self, *args, **kwargs): self._c_mask = new mask() def __init__(self, int n): self.n = n self.f = cuFloatArray(shape=(self.n,1)) self.f._c_gpu.dev_data = self._c_mask.f def alter(self, cuFloatArray filter): """ Modifies the mask pattern Parameters ---------- filter : cuFloatArray The filter to apply to the mask """ self._c_mask.alter(filter._c_gpu.dev_data) def alter(self, MaskAbstract M): """ Modifies the mask pattern Parameters ---------- M : MaskAbstract The filter to apply to the mask """ self._c_mask.alter(M._c_mask.f) property nnz: def __get__(self): return self._c_mask.nnz property nel: def __get__(self): return self._c_mask.nel property area: def __get__(self): return self._c_mask.area @ \subsubsection{Mask} \label{sec:mask-1} \index{utilities!python!Mask} <>= cdef class Mask(MaskAbstract): """ Creates a square mask object. Parameters ---------- n : an int the linear sampling in pixel of the pupil L : float the linear size of the mask in meter i_s : float, optional the row pixel coordinate of the mask center, default: 0 j_s : float, optional the column pixel coordinate of the mask center, default: 0 n_s : float, optional the mask size in pixel, default: 0 means n_s=n theta : float the rotation angle [rd] of the mask with respect to i_s and j_s, default: 0 See also -------- MaskAbstract """ def __init__(self,int n, float L=0.0, float i_0=0.0, float j_0=0.0, float n_s=0.0, float theta=0.0, float i_s=0.0, float j_s=0.0): if n_s==0: if (L>0): super(Mask,self).__init__(n*n) self._c_mask.setup(n,L) else: super(Mask,self).__init__(n) self._c_mask.setup(n) else: super(Mask,self).__init__(n_s*n_s) self._c_mask.setup(n_s,L,i_0,j_0,theta,i_s,j_s,n) self.f = cuFloatArray(shape=(self.n,1)) self.f._c_gpu.dev_data = self._c_mask.f def __dealloc__(self): self._c_mask.cleanup() @ \subsubsection{Telescope} \label{sec:telescope} \index{utilities!python!Telescope} <>= cdef class Telescope(MaskAbstract): """ Creates a circular mask object. Parameters ---------- n : int the sampling in pixel of the pupil. D : float the diamter in meter of the pupil See also -------- MaskAbstract """ def __init__(self, int n, float D, float scale=1.0): super(Telescope,self).__init__(n*n) self._c_mask.setup_circular(n,D,scale) self.f = cuFloatArray(shape=(self.n,1)) self.f._c_gpu.dev_data = self._c_mask.f def __dealloc__(self): self._c_mask.cleanup() @ \subsubsection{GMT} \label{sec:gmt} \index{utilities!python!GMT} <>= cdef class GMT(MaskAbstract): """ Creates a mask object for the Giant Magellan Telescope pupil. Parameters ---------- n : int The sampling in pixel of the pupil. S : float The size of the pupil in meter. Attributes ---------- piston_mask : list A list of piston mask, one per segment See also -------- MaskAbstract """ def __init__(self, int n, float S=25.5): super(GMT,self).__init__(n*n) self._c_mask.setup_GMT(n,S) self.f = cuFloatArray(shape=(self.n,1)) self.f._c_gpu.dev_data = self._c_mask.f self._piston_mask_ = cuIntArray(shape=(1,self.n)) self._piston_mask_._c_gpu.dev_data = self._c_mask.d__piston_mask; def __dealloc__(self): self._c_mask.cleanup() def set_gmt_piston(self, cuFloatArray phase, float[::1] value): """ Adds segment piston to the phase Parameters ---------- phase : cuFloatArray The source wavefront phase. value : float[::1] The 1D array of segment piston in meter """ cdef cuFloatArray p p = cuFloatArray(host_data=value) self._c_mask.set_gmt_piston(phase._c_gpu.dev_data, p._c_gpu.dev_data) def get_gmt_piston(self, cuFloatArray phase): """ Retrieves the segment piston from the phase Parameters ---------- phase : cuFloatArray The source wavefront phase. Returns ------- numpy array : The 1D array of segment piston in meter """ Q = self.piston_mask[0] return np.dot(phase.host().ravel(),Q.T)/np.sum(Q,axis=1) property piston_mask: def __get__(self): P = self._piston_mask_.host() F = self.f.host() return [ np.array( [(P[kk,:]*F.ravel())==k for k in range(1,8)] ) for kk in range(self._piston_mask_.shape[0])] @ \section{Device to host} \label{sec:device-host} The cython interface between GPU and CPU array (Sec.~\ref{sec:dev-to-host}). \subsection{PXD file} \label{sec:pxd-file-2} <>= ## device to host void dev2host( float *host_data, float *dev_data, int N) void host2dev( float **dev_data, float *host_data, int N) void freedev( float **dev_data ) @ \section{cu$<$Type$>$Array} \label{sec:cuarray} The cu$<$Type$>$Array class is a cython class that is a container for CUDA variables: <>= cdef extern from "utilities.h": cdef cppclass gpu_t[T]: T *dev_data T *host_data int N void setup() void setup(int) void dev_malloc() void free_dev() void dev2host() void host2dev() void reset() # void double2float(gpu_t[float] * ) @ \subsection{cuFloatArray} \label{sec:cufloatarray} \index{utilities!python!cuFloatArray} <>= # cuFloatArray from operator import mul from numpy cimport ndarray, dtype from libc.stdint cimport intptr_t cdef class cuFloatArray: cdef gpu_t[float] *_c_gpu <> @ <>= cdef public dtype type cdef public int size cdef public shape cdef public order cdef dev_malloc cdef public ndarray host_data cdef units @ <>= from functools import reduce cdef class cuFloatArray: """ A class to interact with a cuda float array <> """ def __cinit__(self,shape=None,order='C',host_data=None,dev_malloc=False): self._c_gpu = new gpu_t[float]() self.type = dtype(np.float32) <> <> def assign_ptr(self): self._c_gpu.host_data = self.host_data.data; @ <>= Parameters ---------- shape : tuple of int, optional Number of rows and colums; default: None order : string, optional Set the memory layout, "C" or "F"; default: "C" host_data : list or tuple or numpy array, optional Host data to copy to the device; default: None dev_malloc : bool, optional Flag to allocate memory on device; default: False Attributes ---------- shape : tuple The shape of the data array order : string, optional The memory layout, "C" or "F" size : float The number of elements in the data array host_data : ndarray The data on the host dev_malloc : bool True is memory has been allocated on the device @ <>= self._c_gpu.setup(); self.dev_malloc = dev_malloc self.order = order if shape is not None: self.shape = shape self.size = np.prod(shape) self._c_gpu.N = self.size if self.dev_malloc: self._c_gpu.N = self.size self._c_gpu.dev_malloc() if host_data is not None: self.host_data = np.array( host_data , dtype=self.type , order=self.order) self.shape = host_data.shape self.size = np.prod(self.shape) self._c_gpu.N = self.size if not self.dev_malloc: self._c_gpu.dev_malloc() self.assign_ptr() self._c_gpu.host2dev() self.dev_malloc = True self.units = {'micron': 1e6, 'nm': 1e9, 'arcsec': 180*3600/math.pi, 'mas': 1e3*180*3600/math.pi} @ <>= def __dealloc__(self): if self.dev_malloc: self._c_gpu.free_dev() def host(self,units=None,zm=False,mask=None,shape=None, order=None): <> def reset(self): self._c_gpu.reset() @property def dev_ptr(self): cdef intptr_t a a = self._c_gpu.dev_data return a @property def nbytes(self): cdef size_t nb nb = sizeof(self.type)*self._c_gpu.N return nb def __repr__(self): repr(self.host()) def __str__(self): str(self.host()) @ <>= <> if self.shape==(0,0) or self.shape is None: raise ValueError("Data must have a size greater than zero!!") if shape is not None: prod = lambda x,y : x*y assert reduce(prod,shape,)==reduce(prod,self.shape), "Total size of array must be unchanged!" self.shape = shape if order is not None: self.order = order self.host_data = np.zeros(self.shape,order=self.order,dtype=self.type) self.assign_ptr() self._c_gpu.dev2host() if zm: if mask is None: self.host_data -= self.host_data.mean() else: mask = np.reshape(mask,self.shape) mean_host_data = np.sum(self.host_data*mask)/np.sum(mask) self.host_data = mask*(self.host_data-mean_host_data) if units is not None: self.host_data = self.host_data*self.units[units] return self.host_data @ <>= """ Transfers data on the GPU device to the memory of the host Parameters ---------- shape : tuple, optional The shape of the data array, default to None order : string, optional Set the memory layout, "C" or "F"; default: None units : string, optional The data converstion units: ""micron"", ""nm"", ""arcsec" and ""mas"", default to None zm : boolean, optional If True, the data mean is removed; default to False mask : boolean If zm is True, the mean is computed for data values where the mask is True Returns ------- numpy array : a copy of the device data on the host """ @ \subsection{cuIntArray} \label{sec:cuintarray} \index{utilities!python!cuIntArray} <>= # cuIntArray cdef class cuIntArray: cdef gpu_t[int] *_c_gpu <> <>= cdef class cuIntArray: """ A class to interact with a cuda integer array <> """ def __cinit__(self,shape=None,host_data=None,dev_malloc=False,order='C'): self._c_gpu = new gpu_t[int]() self.type = dtype(np.int32) <> <> def assign_ptr(self): self._c_gpu.host_data = self.host_data.data; @ \subsection{cuDoubleArray} \label{sec:cudoublearray} \index{utilities!python!cuDoubleArray} <>= # cuDoubleArray cdef class cuDoubleArray: cdef gpu_t[double] *_c_gpu <> <>= cdef class cuDoubleArray: """ A class to interact with a cuda double array <> """ def __cinit__(self,shape=None,host_data=None,dev_malloc=False,order='C'): self._c_gpu = new gpu_t[double]() self.type = dtype(np.double) <> <> def assign_ptr(self): self._c_gpu.host_data = self.host_data.data; """ @property def __single__(self): cdef cuFloatArray data data = cuFloatArray(shape=self.shape, dev_malloc=True) self._c_gpu.double2float(data._c_gpu) return data """ @ \subsection{cuFloatComplexArray} \label{sec:cufloatarray} \index{utilities!python!cuFloatComplexArray} <>= # cuFloatComplexArray from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free cdef class cuFloatComplexArray: cdef gpu_t[float2] *_c_gpu <> @ <>= cdef: dtype type public int size public shape public order dev_malloc float2 *_host_data_ public ndarray host_data units @ <>= cdef class cuFloatComplexArray: """ A class to interact with a cuda float complex array <> """ def __cinit__(self,shape=None,host_data=None,dev_malloc=False,order='C'): self._c_gpu = new gpu_t[float2]() self.type = dtype(np.complex64) <> <> def assign_ptr(self): self._host_data_ = PyMem_Malloc(self.size*sizeof(float2)) self._c_gpu.host_data = self._host_data_ @ <>= Parameters ---------- shape : tuple of int, optional Number of rows and colums; default: None order : string, optional Set the memory layout, "C" or "F"; default: "C" host_data : list or tuple or numpy array, optional Host data to copy to the device; default: None dev_malloc : bool, optional Flag to allocate memory on device; default: False Attributes ---------- shape : tuple The shape of the data array order : string, optional The memory layout, "C" or "F" size : float The number of elements in the data array host_data : ndarray The data on the host dev_malloc : bool True is memory has been allocated on the device @ <>= cdef int k cdef tuple subs self._c_gpu.setup(); self.dev_malloc = dev_malloc self.order = order if shape is not None: self.shape = shape self.size = np.prod(shape) self._c_gpu.N = self.size if self.dev_malloc: self._c_gpu.N = self.size self._c_gpu.dev_malloc() self.assign_ptr() if host_data is not None: self.host_data = np.array( host_data , dtype=self.type , order=self.order) self.shape = host_data.shape self.size = np.prod(self.shape) self._c_gpu.N = self.size if not self.dev_malloc: self._c_gpu.dev_malloc() self.assign_ptr() for k in range(self.size): subs = np.unravel_index(k,self.shape) self._host_data_[k].x = self.host_data.real[subs] self._host_data_[k].y = self.host_data.imag[subs] self._c_gpu.host2dev() self.dev_malloc = True self.units = {'micron': 1e6, 'nm': 1e9, 'arcsec': 180*3600/math.pi, 'mas': 1e3*180*3600/math.pi} @ <>= def __dealloc__(self): if self.dev_malloc: self._c_gpu.free_dev() PyMem_Free(self._host_data_) def host(self,units=None,zm=False,mask=None,shape=None,order=None): <> def reset(self): self._c_gpu.reset() @ <>= <> cdef int k cdef tuple subs if self.shape==(0,0) or self.shape is None: raise ValueError("Data must have a size greater than zero!!") if shape is not None: prod = lambda x,y : x*y assert reduce(prod,shape,)==reduce(prod,self.shape), "Total size of array must be unchanged!" self.shape = shape if order is not None: self.order = order self.host_data = np.zeros(self.shape,order=self.order,dtype=self.type) self.assign_ptr() self._c_gpu.dev2host() for k in range(self.size): subs = np.unravel_index(k,self.shape) self.host_data.real[subs] = self._host_data_[k].x self.host_data.imag[subs] = self._host_data_[k].y if units is not None: self.host_data = self.host_data*self.units[units] return self.host_data @ <>= """ Transfers data on the GPU device to the memory of the host Parameters ---------- shape : tuple, optional The shape of the data array, default to None order : string, optional Set the memory layout, "C" or "F"; default: None units : string, optional The data converstion units: ""micron"", ""nm"", ""arcsec" and ""mas"", default to None Returns ------- numpy array : a copy of the device data on the host """ @ \section{Atmospheric dispersion} \label{sec:atmosph-disp-1} The cython wrapper for the atmospheric dispersion routine (Sec.~\ref{sec:atmosph-disp}) is \index{utilities!python!atmosphereRefractiveIndex} \index{utilities!python!atmosphericDispersion} <>= ## atmosphere dispersion double atmosphere_refractive_index(float wavelength, float altitude, float temperature, float humidity) double atmospheric_dispersion(float wavelength, float delta_wavelength, float zenith, float altitude, float temperature, float humidity) @ <>= ## atmosphere dispersion cpdef double atmosphereRefractiveIndex(float wavelength, float altitude, float temperature, float humidity ): return atmosphere_refractive_index(wavelength, altitude, temperature, humidity) cpdef double atmosphericDispersion(float wavelength, float delta_wavelength, float zenith, float altitude, float temperature, float humidity): return atmospheric_dispersion(wavelength, delta_wavelength, zenith, altitude, temperature, humidity) @ \section{Wavefront differentiation} \label{sec:wavefr-diff-1} \index{utilities!python!wavefrontFiniteDifference} The cython wrapper for the wavefront differentiation routines (Sec.~\ref{sec:wavefr-diff}) is <>= ## wavefront differentiation void wavefront_finite_difference(float *sx, float *sy, int NL, float *phi, int n, float d, int N_GS) void wavefront_finite_difference(float *sx, float *sy, int NL, float *phi, int n, float d, mask *pupil, int N_GS) @ <>= ## wavefront differentiation cpdef wavefrontFiniteDifference(cuFloatArray sx, cuFloatArray sy, int NL, cuFloatArray phi, int n, float d, MaskAbstract M=None, int N_GS=1): """ Computes the average finite difference of the wavefront Parameters ---------- sx : cuFloatArray The x-axis slopes, memory must be allocated with dev_malloc=True sy : cuFloatArray The y-axis slopes, memory must be allocated with dev_malloc=True NL : int The linear size of the lenslet array phi : cuFloatArray The wavefront n : int The number of pixel per lenslet d : float The lenslet pitch M : MaskAbstract, optional The valid lenslet mask N_GS : int, optional The number of WFS guide stars """ if M is None: wavefront_finite_difference(sx._c_gpu.dev_data, sy._c_gpu.dev_data, NL, phi._c_gpu.dev_data, n, d, N_GS) else: wavefront_finite_difference(sx._c_gpu.dev_data, sy._c_gpu.dev_data, NL, phi._c_gpu.dev_data, n, d, M._c_mask, N_GS) @ \subsection{Special functions} \label{sec:spec_fun} \subsubsection{Irregular Modified Bessel Functions} \label{sec:irreg_modif_bessel} \index{utilities!python!Knu} <>= double _K_nu_(double, double) @ <>= cpdef double Knu(double nu, double x): return _K_nu_(nu,x) @ \section{Geometry} \label{sec:geometry-1} \index{utilities!python!polywind} The cython wrapper for the routine that computes the winding number of a polygon (Sec.~\ref{sec:geometry}) is <>= ## polygon winding number int polywind(double Px, double Py, double *Vx, double *Vy, int NV) void polywinds(int *W, double *Px, double *Py, int NP, double *Vx, double *Vy, int NV) @ <>= ## polygon winding number cpdef polyWind(double Px, double Py, double[:] Vx, double[:] Vy): """ Computes the winding number of a polygon around a point P Parameters ---------- Px : double The x coordinates of the point Py : double The y coordinates of the point Vx : ndarray The x coordinate array of the polygon vertices Vy : ndarray The y coordinate array of the polygon vertices """ return polywind(Px, Py, &Vx[0], &Vy[0], Vx.size) cpdef polyWinds(cuDoubleArray Px, cuDoubleArray Py, cuDoubleArray Vx, cuDoubleArray Vy): """ Computes the winding numbers of a polygon around a array of points P Parameters ---------- Px : cuDoubleArray The x coordinate array of the points Py : cuDoubleArray The y coordinate array of the points Vx : cuDoubleArray The x coordinate array of the polygon vertices Vy : cuDoubleArray The y coordinate array of the polygon vertices """ cdef cuIntArray W cdef int NP, NV NP = Px.size NV = Vx.size W = cuIntArray(shape=Px.shape, dev_malloc=True) polywinds(W._c_gpu.dev_data, Px._c_gpu.dev_data, Py._c_gpu.dev_data, NP, Vx._c_gpu.dev_data, Vy._c_gpu.dev_data, NV) return W @ \section{Sparse matrix} \label{sec:sparse-matrix-1} \index{utilities!python!SparseMatrix} <>= ## sparsematrix structure cdef cppclass sparseMatrix: int nnz, n_row, n_col int *csrColIndH int *csrRowPtrH float *csrValH void gradient(int , int, float ) void interpolation(int, int, mask *,float i0, float) void MVM(float *, float *) void add(sparseMatrix *C, sparseMatrix *A) void cleanup() @ <>= cdef class SparseMatrix: cdef: sparseMatrix *_c_sparseMatrix readonly int nnz, n_row, n_col readonly cuFloatArray csrValH readonly cuIntArray csrColIndH, csrRowPtrH cdef class SparseGradient(SparseMatrix): cdef: readonly int NL, N_LENSLET_PX readonly float d #from source cimport Complex_amplitude #cdef class SparseInterpolation(SparseMatrix): # pass @ <>= ## SparseMatrix cdef class SparseMatrix: def __cinit__(self, *args, **kwargs): print('SparseMatrix.cinit', args, kwargs) self._c_sparseMatrix = new sparseMatrix(); def __init__(self, int nnz, int n_row, int n_col, *args, **kwargs): print('SparseMatrix.init', nnz, n_row, n_col, args, kwargs) self.nnz = nnz self.n_row = n_row self.n_col = n_col self.csrValH = cuFloatArray(shape=(self.nnz,1)) self.csrValH._c_gpu.dev_data \ = self._c_sparseMatrix.csrValH self.csrColIndH = cuIntArray(shape=(self.nnz,1)) self.csrColIndH._c_gpu.dev_data \ = self._c_sparseMatrix.csrColIndH self.csrRowPtrH = cuIntArray(shape=(self.n_row+1,1)) self.csrRowPtrH._c_gpu.dev_data \ = self._c_sparseMatrix.csrRowPtrH def __dealloc__(self): self._c_sparseMatrix.cleanup() def __mul__(_S_,_X_): if isinstance(_S_,SparseMatrix) and isinstance(_X_,cuFloatArray): S = _S_ X = _X_ Y = cuFloatArray(shape=(S.n_row,1),dev_malloc=True) S._c_sparseMatrix.MVM(Y._c_gpu.dev_data,X._c_gpu.dev_data) return Y else: return NotImplemented def __add__(_A_,_B_): if isinstance(_A_,SparseMatrix) and isinstance(_B_,SparseMatrix): A = _A_ B = _B_ C = SparseMatrix(A.nnz+B.nnz, A.n_row, A.n_col) A._c_sparseMatrix.add(C._c_sparseMatrix,B._c_sparseMatrix) return C else: return NotImplemented property H: def __get__(self): from scipy.sparse import csr_matrix return csr_matrix((self.csrValH.host().ravel(), self.csrColIndH.host().ravel(), self.csrRowPtrH.host().ravel()), shape=(self.n_row,self.n_col)) @ \subsection{SparseGradient} \label{sec:sparsegradient} \index{utilities!python!SparseGradient} <>= ## SparseGradient cdef class SparseGradient(SparseMatrix): def __init__(self, int NL, int N_LENSLET_PX, float d): self.NL = NL self.N_LENSLET_PX = N_LENSLET_PX self.d = d self._c_sparseMatrix.gradient(NL, N_LENSLET_PX, d) super().__init__(self._c_sparseMatrix.nnz, self._c_sparseMatrix.n_row, self._c_sparseMatrix.n_col) ## SparseInterpolation #cdef class SparseInterpolation(SparseMatrix): # def __init__(self,int NI, int NP, Complex_amplitude W, # float i0, float j0): # self._c_sparseMatrix.interpolation(NI,NP,W._c_complex_amplitude.M,i0,j0) # super().__init__(self._c_sparseMatrix.nnz, # self._c_sparseMatrix.n_row, # self._c_sparseMatrix.n_col) @ Adding the class definitions to the definition file <>= <>