 |
MAOS
Multithreaded Adaptive Optics Simulator
|
|
MAOS
Multithreaded Adaptive Optics Simulator
|
Macros | |
| #define | AOS_LIB_TYPE |
| #define | AOS_DMAT(A) d##A |
Typedefs | |
| typedef real(* | dminsearch_fun) (const real *x, void *info) |
Calls X(new) with a filename to be saved to.
| dmat * dmat_cast | ( | const_anyarray | A | ) |
check and cast an object to matrix
check the size of matrix if exist. Otherwise create it. content is not zeroed if already exists.
creat a reference to an existing matrix. header is duplicated if exists.
create an new reference another with different shape.
creat a new matrix referencing columns in existing matrix. reference counted.
Override a stack matrix struct with pointers to columns of another matrix. Does not add reference count the original data.
Resize a matrix by adding or removing columns or rows. Data is kept whenever possible.
Create a new sub matrix of nx*ny starting from(sx,sy) and copy the data.
Concatenate two matrixes into 1 along "dim" (1 for x, 2 for y)
| int dzerocol | ( | dmat * | A, |
| int | icol | ||
| ) |
Set column elements to zero.
Permute the vector so that out(:)=in(p);
Reverse permute the vector so that out(p)=in(:);
| int dflip | ( | dmat * | A, |
| int | axis | ||
| ) |
Flip the matrix along the set axis. 0 to flip both, 1 along x, 2 along y.
| real dvecsum | ( | const double *__restrict | p, |
| long | nx | ||
| ) |
sum of a vector using pointer and length
Sort all columns of A, in ascending order if ascend is non zero, otherwise in descending order.
| void dvecmaxmin | ( | const double *__restrict | p, |
| long | N, | ||
| double * | max, | ||
| double * | min | ||
| ) |
Compute max, min and sum. Has to handle NAN nicely. Complex values are converted into magnitude during comparison.
| real dvecdot | ( | const double *__restrict | p1, |
| const double *__restrict | p2, | ||
| const double *__restrict | p3, | ||
| long | n | ||
| ) |
compute sum(p1.*p2.*p3)
| void dvecnormalize_sum | ( | double *__restrict | p, |
| long | nx, | ||
| double | norm | ||
| ) |
normalize vector so that the sum is norm;
| void dvecnormalize_sumabs | ( | double *__restrict | p, |
| long | nx, | ||
| double | norm | ||
| ) |
normalize vector so that the sum of abs is norm;
| real dvecmaxabs | ( | const double *__restrict | p, |
| long | n | ||
| ) |
find the maximum of abs of all elements using pointer and count.
Find unique values and sort from small to large.
| input | input array |
| ascend | sorting direction |
reorder B and embed/crop into center of A .
4 * * 3 * * * * * * * * 2 * * 1
to
1 2 3 4
cyclic shift A by nx and ny to B.
4 3 1 2 2 1 to 3 4
create an new cell similar to A in shape. When a cell is empty, it is created with an emtry matrix and cannot be overriden.
Create an new cell with matrix specified. Each block is stored continuously in memory.
Create an new cell with matrix specified. Each block is stored continuously in memory.
| dcell * dcellnew_file | ( | long | nx, |
| long | ny, | ||
| long * | nnx, | ||
| long * | nny, | ||
| const char * | keywords, | ||
| const char * | format, | ||
| ... | |||
| ) |
Calls cellnew3 with a filename to be saved to.
| dcell * dcellnewsame_file | ( | long | nx, |
| long | ny, | ||
| long | mx, | ||
| long | my, | ||
| const char * | keywords, | ||
| const char * | format, | ||
| ... | |||
| ) |
Calls cellnew_same with a filename to be saved to.
creat a cell reference an existing cell by referencing the elements.
duplicate a X(cell) object. Not implemented for lmat or imat as X(cellcp) is defined in spmm.c
Copy B to A.
Takes parameters of matrix, sparse matrix, or cell array of them.
reduce nx*ny cell matrix to 1*ny if dim=1 and nx*1 if dim=2 by merging the cells.
concatenate two cell matrices along dimenstion 'dim'.
concatenate two cell matrices along dimenstion 'dim'.
concatenate coresponding elements of each X(cell). They must have the same shape.
convert a vector to cell using dimensions specified in dims. Reference the vector
make A a cell array using shape information from ref if *B is NULL.
Notice that empty cells may be replaced by zeros if it is within valid column or row.
input is nsa*ncol cell. each cell has npix=nx*ny elements. Extract icol of cell as npix*nsa array.
Fill A with random uniform numbers between [0, 1]*max
Fill A with random normal distribution numbers with standard deviation of sigma.
compute weighted dot product a'*(w*b)
| void dcwmcol | ( | dmat *__restrict | A, |
| const dmat *__restrict | B | ||
| ) |
Compute component wise multiply A=A.*(B1*wt1+B2*wt2) component-wise multiply of three matrices. A=A.*W.*(B1*wt1+B2*wt2) Component-wise multiply each column of A with B A(:,i)=A(:,i).*B;
| void dcwmrow | ( | dmat *__restrict | A, |
| const dmat *__restrict | B | ||
| ) |
component-wise multiply of columns of A with combination of B1 and B2: A(:,i)=A(:,i)*(B1*wt1+B2*wt2); component wise multiply of 2d complex matrix A,W and 1d vector B. A(:,i)=A(:,i).*W(:,i).*B; Component wise multiply each row of A with B. A(i,:)=A(i,:)*B
component-wise multiply of rows of A with combination of B1 and B2: A(i,:)=A(i,:).*(B1*wt1+B2*wt2); Component wise division B=B./A. 0/0 is replace by 'value';
| void dmulvec | ( | double *__restrict | y, |
| const dmat *__restrict | A, | ||
| const double *__restrict | x, | ||
| const double | alpha | ||
| ) |
multiply a X(mat) matrix with a vector and accumulate to y: y+=A*x*alpha
| void dmm | ( | dmat ** | C0, |
| const double | beta, | ||
| const dmat * | A, | ||
| const dmat * | B, | ||
| const char | trans[2], | ||
| const double | alpha | ||
| ) |
compute matrix product using blas dgemm.; C=beta*C+ alpha *trans(A)*trans(B); if C exist.
inplace invert a small square SPD matrix using lapack dposv_, usually (A'*w*A). by solving Ax=I; copy x to A. dposv_ modifies A also. be careful
Compute the cholesky decomposition of a symmetric semi-definit dense matrix.
compute the pseudo inverse of matrix A with weigthing of full matrix W or sparse matrix weighting Wsp. For full matrix, wt can be either W or diag (W) for diagonal weighting. B=inv(A'*W*A+tikcr)*A'*W; thres is the threshold to truncate eigenvalues.
compute the relative difference betwee two vectors. sqrt(||A-B||/||A||) using norm2. for debugging purpose.
| int dcircle | ( | dmat * | A, |
| double | cx, | ||
| double | cy, | ||
| double | dx, | ||
| double | dy, | ||
| double | r, | ||
| double | val | ||
| ) |
Generate a new gray pixel map based on bilinear influence functions used in mkw. It creates slightly larger map than an filled circle to account for the sampling effect. The Center and Radius cx,cy,r are in unit of meter, The sampling dx,dy specify spacing of the points in meter along x or y dimension. Each full point received value of 'val'
| int drectangle | ( | dmat * | A, |
| double | cx, | ||
| double | cy, | ||
| double | rx, | ||
| double | ry, | ||
| double | theta, | ||
| double | val | ||
| ) |
Generate a new gray pixel map based on bilinear influence functions used in mkw. It creates slightly larger map than an filled rectangle to account for the sampling effect. The Center and half width cx,cy,rx,ry are in unit of meter, The sampling dx,dy specify spacing of the points in meter along x or y dimension. Each full point received value of 'val'
rotate the column vectors CCW, equivalent as rotate coordinate theta CW. A is nx2 while A(:,1) is x, A(:,2) is y.
rotate the row vectors CCW. same as rotate coordinate theta CW. A is 2xn A(:,1) is x, A(:,2) is y.
rotate a 2x2 covariance matrix A by theta CCW (coordinate rotate -theta CCW) or from ra to xy coordinate. R*A*R';
T matrix vector multiply optimized for just three values. y=A*x;
Compute the parabolic fit center using 3x3 points around the peak.
| void dcog | ( | double * | grad, |
| const dmat * | im, | ||
| double | offsetx, | ||
| double | offsety, | ||
| double | thres, | ||
| double | bkgrnd, | ||
| double | flux | ||
| ) |
Compute thresholded center of gravity. The threshold is absolute value. bkgrnd is removed from i0 when computing cog. offset is the offset of the reference point (cog=0) from the physical center. all length are given in terms of pixel.
Shift the image in A to center on physical center+[offsetx,offsety] using cog and fft.
Compute SVD of a general matrix A. A=U*diag(S)*V'; diag(S) is returned.
computes pow(A,power) in place using svd. positive thres: Drop eigenvalues that are smaller than thres * max eigen value negative thres: Drop eigenvalues that are smaller than thres * previous eigen value (sorted descendantly).
| A | [in/out] The matrix | |
| [in] | power | power of eigen values. usually -1 for inverse |
| [in] | thres1 | SVD inverse absolute threshold |
| [in] | thres2 | SVD inverse relative threshold |
computes out=out*alpha+exp(A*beta) using scaling and squaring method. Larger scaling is more accurate but slower. Tested against matlab expm
First, scaling the matrix by 2^-n times so it is smaller than 1/threshold. Then compute exp(A*beta*2^-n) using Tylor expansion. Then compute the final answer by taking exponential of 2^n.
compute B=bc*B+ac*A behavior changed on 2009-11-02. if A is NULL, don't do anything.
compute B=bc*B+ac*A. Use smallest dimension of both. behavior changed on 2009-11-02. if A is NULL, don't do anything.
Compute B[:,icol]=bc*B[:,icol]+ac*A where A is vector.
Interpolation of 1d array
Interpolate using linear interp. xin is the coordinate of yin. xnew is the coordinate of the output.
Interpolate using log(xin) and log(xnew) xin is the coordinate of yin. xnew is the coordinate of the output.
| void dblend | ( | dmat *__restrict | A, |
| dmat *__restrict | B, | ||
| int | overlap | ||
| ) |
blend B into center of A with width of overlap. The center (size is B->nx-overlap, B->ny-overlap) of A is replaced by center of B . The overlapping area is blended
For each entry in A, call repeatly to collect its histogram, centered at center, spaced by spacing, for n bins in total. center if at bin n/2.
1D Cubic spline interpolation preparation. if x has only 1 column: x is the coordinate. y is the function value. if x has two columns: first column is the coordinate, y is null.
It is upto the user to make sure that the coordinate is increasingly ordered and evenly spaced .
If the values of a function \(f(x)\) and its derivative are know at x=0, and x=1 (normalized coordinate), then the function can be interpolated on the interval [0,1] using a third degree polynomial. This is called cubic interpolation. The formula of this polynomial can be easily derived.
A third degree polynomial and its derivative:
\[ f(x)=ax^3+bx^2+cx+d \]
\[ f(x)=3ax^3+2bx+c \]
The coefficients can be derived from the value and derivatives:
\begin{eqnarray*} a&=&2f(0)-2f(1)+f^\prime (0)+f^\prime(1)\\ b&=&-3f(0)+3f(1)-2f^\prime(0)-f^\prime(0)\\ c&=&f^\prime(0)\\ d&=&f(0)\\ \end{eqnarray*}
the derivatives can be computed as
\begin{eqnarray*} f^\prime(0)&=&\frac{f(1)-f(-1)}{2}\\ f^\prime(1)&=&\frac{f(2)-f(0)}{2}\\ \end{eqnarray*}
so we have the formula
\begin{eqnarray*} a&=&-0.5 f(-1) + 1.5 f(0) - 1.5 f(1) + 0.5 f(2)\\ b&=& f(-1) - 2.5 f(0) + 2 f(1) - 0.5 f(2)\\ c&=&-0.5 f(-1) + 0.5 f(1) \\ d&=& f(0) \\ \end{eqnarray*}
for the boundary pints, replace
\[f^\prime(0)=(f(1)-f(-1))/2\]
by
\[f^\prime(0)=(f(1)-f(0))\]
Otehr type of boundaries are handled in the same way.
see www.paulinternet.nl/?page=bicubicx
1D Cubic spline interpolation preparation. if x has only 1 column: x is the coordinate. y is the function value. if x has two columns: first column is the coordinate, y is null.
It is upto the user to make sure that the coordinate is increasingly ordered and evenly spaced .
If the values of a function \(f(x)\) and its derivative are know at x=0, and x=1 (normalized coordinate), then the function can be interpolated on the interval [0,1] using a third degree polynomial. This is called cubic interpolation. The formula of this polynomial can be easily derived.
A third degree polynomial and its derivative:
\[ f(x)=ax^3+bx^2+cx+d \]
\[ f(x)=3ax^3+2bx+c \]
The coefficients can be derived from the value and derivatives:
\begin{eqnarray*} a&=&2f(0)-2f(1)+f^\prime (0)+f^\prime(1)\\ b&=&-3f(0)+3f(1)-2f^\prime(0)-f^\prime(0)\\ c&=&f^\prime(0)\\ d&=&f(0)\\ \end{eqnarray*}
the derivatives can be computed as
\begin{eqnarray*} f^\prime(0)&=&\frac{f(1)-f(-1)}{2}\\ f^\prime(1)&=&\frac{f(2)-f(0)}{2}\\ \end{eqnarray*}
so we have the formula
\begin{eqnarray*} a&=&-0.5 f(-1) + 1.5 f(0) - 1.5 f(1) + 0.5 f(2)\\ b&=& f(-1) - 2.5 f(0) + 2 f(1) - 0.5 f(2)\\ c&=&-0.5 f(-1) + 0.5 f(1) \\ d&=& f(0) \\ \end{eqnarray*}
for the boundary pints, replace
\[f^\prime(0)=(f(1)-f(-1))/2\]
by
\[f^\prime(0)=(f(1)-f(0))\]
Otehr type of boundaries are handled in the same way.
see www.paulinternet.nl/?page=bicubicx
Evluate the cubic spline represented by nx5 matrix coeff, at location array xnew.
Do 1D cubic spline all at once by calling X(spline_prep) and X(spline_evald)
| void dembed | ( | dmat *__restrict | out, |
| const dmat *__restrict | in, | ||
| const double | theta | ||
| ) |
Embed a ninx*niny matrix in into A with optional rotation by theta CCW (output points are rotated by theta -CCW) around the fft center. Used to rotate the PSF from x-y to radial-azimuthal coordinate in radial format CCD. A may be bigger or smaller than B.
| void dembedd | ( | dmat *__restrict | out, |
| const dmat *__restrict | in, | ||
| const double | theta | ||
| ) |
embed a ninx*niny matrix in into A with optional rotation by -theta CCW (coordinate rotate theta CCW) around the fft center. Used to rotate the PSF from x-y to radial-azimuthal coordinate in radial format CCD. A may be bigger or smaller than B.
| real dfwhm | ( | dmat * | A | ) |
Calculate number of pixels having values larger than or equal to half of maximum and convert to diameter.
For more accurate calculation of Gaussian beams, use the sqrt(2*log(2))*X(gauss_width).
Calculate D4s width of Gaussian 2D function
| mr | Equivalent radius of the same area |
| ma | major axis |
| mb | minor axi |
| angle | angle |
| A | The irradiance (intensity) |
| thres | The threshold relative to peak. |
| real dgauss_width | ( | dmat * | A, |
| double | thres | ||
| ) |
A convenient wrapper
| real dfwhm_gauss | ( | dmat * | A | ) |
Use gauss fit to compute fwhm
Compute the enclosed energy or azimuthal average of a.
| psf | The input array |
| dvec | The diameter for enclosed energy, or radius for azimuthal average |
| type | The type. -1: azimuthal average, 0: within a square, 1: within a circle, 2: within a slit |
| nthread | Number of CPU threads to use. |
| int dminsearch | ( | double * | x, |
| int | nmod, | ||
| double | ftol, | ||
| int | nmax, | ||
| dminsearch_fun | fun, | ||
| void * | info | ||
| ) |
Search minimum along multiple dimenstions. scale contains the scale of each dimension. x contains initial warm restart values.
| void dbessik | ( | double | x, |
| double | xnu, | ||
| double * | ri, | ||
| double * | rk, | ||
| double * | rip, | ||
| double * | rkp | ||
| ) |
Returns the modified Bessel functions ri = Iν , rk = Kν and their derivatives rip = Iν , rkp = Kν, for positive x and for xnu = ν ≥ 0.
| void dcellscale | ( | anyarray | A_, |
| double | w | ||
| ) |
scale each element of A.
invert each component of the dcell. Each cell is treated as an individual matrix.
compute the pseudo inverse of block matrix A. A is n*p cell, wt n*n cell or sparse cell. \(B=inv(A'*W*A)*A'*W\)
| [in] | A | The matrix to pseudo invert |
| [in] | W_ | The weighting matrix. dense or sparse |
| [in] | thres1 | SVD absolute inverse threshold |
| [in] | thres2 | SVD relative inverse threshold |
compute the power of a block matrix using svd method. First convert it do X(mat), do the power, and convert back to block matrix.
| int dcellclip | ( | dcell * | Ac, |
| double | min, | ||
| double | max | ||
| ) |
clip a X(cell) array to max at 'max', min at 'min'
Create a new sub cell matrix of nx*ny starting from(sx,sy)
2D cubic spline interpolation preparation. x is the x coordinate vector of the 2-d grid. y is the y coordinate vector of the 2-d grid. z is defined on the 2-d grid. It is upto the user to make sure that the coordinate is increasingly ordered and evenly spaced .
The boundaries are handled in the same way is X(spline). i.e. replace
\[f^\prime(0)=(f(1)-f(-1))/2\]
by
\[f^\prime(0)=(f(1)-f(0))\]
Otehr type of boundaries are handled in the same way.
Evaluate 2D cubic spline at location defined 2-d arrays by xnew, ynew
| void dmaprot | ( | anyarray | A_, |
| real | theta | ||
| ) |
Rotate a or multiple 2d array ccw by angle
Routines for input/output to bin/fits file. Function to write dense matrix data into a file pointer. Generally used by library developer
Function to read dense matrix into memory from file pointer. Generally used by library developer.
Create a new X(mat) matrix object, mmapped from file. The file is truncated if already exists in rw mode.
| dcell * dcellnew_mmap | ( | long | nx, |
| long | ny, | ||
| long * | nnx, | ||
| long * | nny, | ||
| const char * | keywords, | ||
| const char * | format, | ||
| ... | |||
| ) |
Create a new dense matrix cell object, mmapped from file. The file is truncated if already exists.
| dcell * dcellnewsame_mmap | ( | long | nx, |
| long | ny, | ||
| long | mx, | ||
| long | my, | ||
| const char * | keywords, | ||
| const char * | format, | ||
| ... | |||
| ) |
Create a new dense matrix cell object, with identical blocks, mmapped from file. be aware that the data is not 8-byte aligned. The file is truncated if already exists.
Map the file to memory in read only, shared mode.
Create a nx*ny sparse matrix with memory for nmax max elements allocated.
Create a new sparse matrix and fill in uniform random numbers with filling factor of 'fill'
Check a sparse array for wrong orders. Return 1 if is lower triangle, 2 if upper triangle, and 3 if diagonal.
Create a new sparse matrix with diagonal elements set to vec*alpha
| void dspmuldiag | ( | dsp *__restrict | A, |
| const double * | w, | ||
| double | alpha | ||
| ) |
Multiply a sparse matrix inplace with a diagonal weighting matrix whose diagonal values are stored in w. W_ii=w_i; W_ij=0 if i!=j A=A*W*alpha; W is a diagonal sparse matrix. diag(W) is w multiply w[i] to all numbers in column[i]
sparse matrix multiply with a vector: y=y+op(A)*x*alpha op(A)=A if trans=='n' op(A)=A' if trans=='t' op(A)=conj(A') if trans=='c'
| void dspmulcreal | ( | double *__restrict | y, |
| const dsp * | A, | ||
| const dcomplex *__restrict | x, | ||
| double | alpha | ||
| ) |
Multiply a sparse matrix with the real part of a complex vector: y=y+A*creal(x)*alpha
sparse matrix multiply with dense matrix
| void dmulsp | ( | dmat ** | yout, |
| const dmat * | x, | ||
| const dsp * | A, | ||
| const char | trans[2], | ||
| const double | alpha | ||
| ) |
dense matrix multiply with sparse matrix, implemented by transpose all: y=x*A-->y'=A'*x'; conjugation is handled carefully.
Multiply two vectors with weighting by sparse matrix. return y'*(A*x)
Multiply two cell arrays with weighting by sparse matrix
| void dcellmm | ( | panyarray | C0_, |
| const_anyarray | A_, | ||
| const_anyarray | B_, | ||
| const char | trans[2], | ||
| const double | alpha | ||
| ) |
Compute A*B and add to C0.
C0=C0+op(A)*op(B)*alpha; op(A)=A if trans[0]=='n' op(A)=A' if trans[0]=='t' op(A)=conj(A') if trans[0]=='c'
op(B)=B if trans[1]=='n' op(B)=B' if trans[1]=='t' op(B)=conj(B') if trans[1]=='c'
A may be dense or sparse matrix.
2009-11-09: There was initially a beta parameter It was implemented wrongly for beta!=1 because for every call to dmm, the already accumulated ones are scaled. removed beta.
a different interface for multiplying cells.
Convert sparse matrix into dense matrix and add to output: out0=out0+full(A)*alpha
Expand sparse matrix cell to dense matrix cell.
Added two sparse matrices: return A*a+B*b
Add a sparse matrix to another: A0=A0+B
| void dcelladd | ( | panyarray | pA_, |
| double | ac, | ||
| const_anyarray | B_, | ||
| double | bc | ||
| ) |
Calculate A_=A_*ac+B_*bc;
Takes parameters of matrix, sparse matrix, or cell array of them.
Add alpha times identity to a sparse matrix. If X(sp) is not symmetric, only add diagonal to first nm*nm block for nm=min(nx,ny) Assume the the sparse array is sported correctly.
| void dcelladdI | ( | anyarray | A_, |
| double | alpha | ||
| ) |
Add alpha to diagnonal elements of A_. A_ can be matrix or sparse matrix.
Multiply two sparse arrays and return the result op(A)*op(B)
Multiply two sparse arrays and add to the third: C0=C0+op(A)*op(B)*scale
Concatenate two sparse array along dim dimension
Sum abs of elements of sparse array along dimension dim
Drop elements that are EPS times the largest value.
Make sure the elements are sorted correctly. Does not change the location of data. can be done without harm.
Create a sparse convolution operator C with C(i,j)=A(i-j); A must be very sparse with only a view non-zero value otherwise C will be too full.
Permute rows and columns of sparse matrix A;
Invert a SPD sparse matrix that is block diagonal with block sizes of bs.
Extrat the diagonal blocks of size bs into cell arrays.
| dmat * dcell2m | ( | const_anyarray | A_ | ) |
Convert dense or sparse matrix cell to matrix.
Function to write sparse matrix data into file pointed using a file pointer. Generally used by library developer. We do not convert data during saving, but rather do the conversion during reading.