Skip to content

Commit

Permalink
Merge pull request #413 from MethodicalAcceleratorDesign/dev
Browse files Browse the repository at this point in the history 
Dev to Master for 0.9.7-1
  • Loading branch information
@ldeniau
ldeniau authored Dec 4, 2023
2 parents 3cf6908 + bdc0daf commit a06df86
Show file tree
Hide file tree
Showing 13 changed files with 272 additions and 95 deletions.
120 changes: 114 additions & 6 deletions src/mad_mat.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1144,6 +1144,114 @@ int
mad_cmat_invc_r (const cpx_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)
{ return mad_cmat_invc(y, CPX(x_re,x_im), r, m, n, rcond); }

// -- pseudo-inverse ----------------------------------------------------------o

// SVD decomposition A = U * S * V.t()
// A:[m x n], U:[m x m], S:[min(m,n)], V:[n x n]

int
mad_mat_pinvn (const num_t y[], num_t x, num_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{
CHKYR; // compute x/Y[m x n]
ssz_t mn = MIN(m,n);
mad_alloc_tmp(num_t, s, mn );
mad_alloc_tmp(num_t, U, m*m);
mad_alloc_tmp(num_t, V, n*n);
mad_alloc_tmp(num_t, S, n*m); mad_vec_fill(0, S, n*m);

int rnk = 0;
int info = mad_mat_svd(y, U, s, V, m, n);
if (info != 0) goto finalize;

// Remove ncond (largest) singular values
idx_t k = 0;
FOR(i,MIN(mn,-ncond)) s[k++] = 0;

// Tolerance on keeping singular values.
rcond = MAX(fabs(rcond), DBL_EPSILON);

// Keep relevant singular values and reject ncond (smallest) singular values
FOR(i,k,mn)
if (mn-i >= ncond && s[i] >= rcond*s[k]) S[i*m+i] = (++rnk, 1/s[i]);
else break;

mad_mat_muld(V, S, r, n, m, n);
mad_mat_mult(r, U, r, n, m, m);

if (x != 1) mad_vec_muln(r, x, r, m*n);

finalize:
mad_free_tmp(s);
mad_free_tmp(U);
mad_free_tmp(V);
mad_free_tmp(S);

return rnk;
}

int // without complex-by-value version
mad_mat_pinvc_r (const num_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{ return mad_mat_pinvc(y, CPX(x_re,x_im), r, m, n, rcond, ncond); }

int
mad_mat_pinvc (const num_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{
CHKYR; // compute x/Y[m x n]
mad_alloc_tmp(num_t, rr, m*n);
int info = mad_mat_pinvn(y, 1, rr, m, n, rcond, ncond);
if (!info) mad_vec_mulc(rr, x, r, m*n);
mad_free_tmp(rr);
return info;
}

int
mad_cmat_pinvc_r (const cpx_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{ return mad_cmat_pinvc(y, CPX(x_re,x_im), r, m, n, rcond, ncond); }

int
mad_cmat_pinvn (const cpx_t y[], num_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{ return mad_cmat_pinvc(y, CPX(x,0), r, m, n, rcond, ncond); }

int
mad_cmat_pinvc (const cpx_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond)
{
CHKYR; // compute x/Y[m x n]
ssz_t mn = MIN(m,n);
mad_alloc_tmp(num_t, s, mn );
mad_alloc_tmp(cpx_t, U, m*m);
mad_alloc_tmp(cpx_t, V, n*n);
mad_alloc_tmp(num_t, S, n*m); mad_vec_fill(0, S, n*m);

int rnk = 0;
int info = mad_cmat_svd(y, U, s, V, m, n);
if (info != 0) goto finalize;

// Remove ncond (largest) singular values
idx_t k = 0;
FOR(i,MIN(mn,-ncond)) s[k++] = 0;

// Tolerance on keeping singular values.
rcond = MAX(fabs(rcond), DBL_EPSILON);

// Keep relevant singular values and reject ncond (smallest) singular values
FOR(i,k,mn)
if (mn-i >= ncond && s[i] >= rcond*s[k]) S[i*m+i] = (++rnk, 1/s[i]);
else break;

mad_cmat_muldm(V, S, r, n, m, n);
mad_cmat_mult (r, U, r, n, m, m);

if (x != 1) mad_cvec_mulc(r, x, r, m*n);

finalize:
mad_free_tmp(s);
mad_free_tmp(U);
mad_free_tmp(V);
mad_free_tmp(S);

return rnk;
}

// -- divide ------------------------------------------------------------------o

// note:
Expand Down Expand Up @@ -2200,10 +2308,10 @@ mad_mat_svdcnd(const num_t a[], idx_t c[], ssz_t m, ssz_t n,
if (N > mn || N <= 0) N = mn;

// Tolerance on components similarity in V columns.
tol = MAX(tol, DBL_EPSILON);
tol = MAX(fabs(tol), DBL_EPSILON);

// Tolerance on keeping singular values.
rcond = MAX(rcond, DBL_EPSILON);
rcond = MAX(fabs(rcond), DBL_EPSILON);

// Number of columns to remove.
idx_t nc = 0;
Expand All @@ -2212,7 +2320,7 @@ mad_mat_svdcnd(const num_t a[], idx_t c[], ssz_t m, ssz_t n,

// Loop over increasing singular values.
RFOR(i,mn,mn-N) {
// Singular value is large, stop checking.
// Singular value is large enough, stop checking.
if (S[i] > rcond*S[0]) break;

// Loop over rows of V (i.e. columns of V^T)
Expand Down Expand Up @@ -2266,10 +2374,10 @@ mad_cmat_svdcnd(const cpx_t a[], idx_t c[], ssz_t m, ssz_t n,
if (N > mn || N <= 0) N = mn;

// Tolerance on components similarity in V columns.
tol = MAX(tol, DBL_EPSILON);
tol = MAX(fabs(tol), DBL_EPSILON);

// Tolerance on keeping singular values.
rcond = MAX(rcond, DBL_EPSILON);
rcond = MAX(fabs(rcond), DBL_EPSILON);

// Number of columns to remove.
idx_t nc = 0;
Expand All @@ -2278,7 +2386,7 @@ mad_cmat_svdcnd(const cpx_t a[], idx_t c[], ssz_t m, ssz_t n,

// Loop over increasing singular values.
RFOR(i,mn,mn-N) {
// Singular value is large, stop checking.
// Singular value is large enough, stop checking.
if (S[i] > rcond*S[0]) break;

// Loop over rows of V (i.e. columns of V^T)
Expand Down
6 changes: 6 additions & 0 deletions src/mad_mat.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,9 @@ int mad_mat_det (const num_t x[], num_t *r ,
int mad_mat_invn (const num_t y[], num_t x, num_t r[], ssz_t m, ssz_t n, num_t rcond); // num / mat
int mad_mat_invc (const num_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / mat
int mad_mat_invc_r (const num_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / mat
int mad_mat_pinvn (const num_t y[], num_t x, num_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // num / mat (pseudo-inverse)
int mad_mat_pinvc (const num_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / mat
int mad_mat_pinvc_r (const num_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / mat
int mad_mat_div (const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // mat / mat
int mad_mat_divm (const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // mat / cmat
int mad_mat_solve (const num_t a[], const num_t b[], num_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // min|b-ax| (QR)
Expand Down Expand Up @@ -90,6 +93,9 @@ int mad_cmat_det (const cpx_t x[], cpx_t *r ,
int mad_cmat_invn (const cpx_t y[], num_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // num / cmat
int mad_cmat_invc (const cpx_t y[], cpx_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / cmat
int mad_cmat_invc_r (const cpx_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / cmat
int mad_cmat_pinvn (const cpx_t y[], num_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // num / cmat (pseudo-inverse)
int mad_cmat_pinvc (const cpx_t y[], cpx_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / cmat
int mad_cmat_pinvc_r (const cpx_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / cmat
int mad_cmat_div (const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // cmat / cmat
int mad_cmat_divm (const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // cmat / mat
int mad_cmat_solve (const cpx_t a[], const cpx_t b[], cpx_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // min|b-ax| (QR)
Expand Down
4 changes: 4 additions & 0 deletions src/madl_cmad.mad
View file
Original file line number Diff line number Diff line change
Expand Up @@ -367,6 +367,8 @@ void mad_mat_muldm (const num_t x[], const cpx_t y[], cpx_t r[], ssz_t
int mad_mat_det (const num_t x[], num_t *r , ssz_t n); // det(mat)
int mad_mat_invn (const num_t y[], num_t x , num_t r[], ssz_t m, ssz_t n, num_t rcond); // num / mat
int mad_mat_invc_r (const num_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / mat
int mad_mat_pinvn (const num_t y[], num_t x, num_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // num / mat (pseudo-inverse)
int mad_mat_pinvc_r (const num_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / mat
int mad_mat_div (const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // mat / mat
int mad_mat_divm (const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // mat / cmat
int mad_mat_solve (const num_t a[], const num_t b[], num_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // min|b-ax| (QR)
Expand Down Expand Up @@ -407,6 +409,8 @@ void mad_cmat_muldm (const cpx_t x[], const num_t y[], cpx_t r[], ssz_t
int mad_cmat_det (const cpx_t x[], cpx_t *r , ssz_t n); // det(cmat)
int mad_cmat_invn (const cpx_t y[], num_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // num / cmat
int mad_cmat_invc_r (const cpx_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond); // cpx / cmat
int mad_cmat_pinvn (const cpx_t y[], num_t x , cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // num / cmat (pseudo-inverse)
int mad_cmat_pinvc_r (const cpx_t y[], num_t x_re,num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond, int ncond); // cpx / cmat
int mad_cmat_div (const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // cmat / cmat
int mad_cmat_divm (const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // cmat / mat
int mad_cmat_solve (const cpx_t a[], const cpx_t b[], cpx_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond); // min|b-ax| (QR)
Expand Down
62 changes: 28 additions & 34 deletions src/madl_cofind.mad
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,16 +114,10 @@ end
-- cofind using jacobian ------------------------------------------------------o

local function cofind_jac (self, mflw)
local cotol, coiter, codiff, totalpath in self

-- sanity check and adjustment
if cotol < codiff*10 then
warn("cotol < codiff*10, adjusting codiff to ctol/100")
codiff = cotol/100
end
local coitr, cotol, costp, totalpath in self

-- save current orbits, extend mflw of n particles to n*(1+6) particles
local n, X0 = mflw.npar, table.new(mflw.npar,0)
local n, X0, dH = mflw.npar, table.new(mflw.npar,0), table.new(mflw.npar,0)
for i=n,1,-1 do
local m, ii = mflw[i], 7*(i-1)+1
local mt = {__index=m, __newindex=
Expand All @@ -134,18 +128,20 @@ local function cofind_jac (self, mflw)
setmetatable(mc,mt) -- connect secondary particles to primary particle
end
mflw[ii], m.id, m.coid = m, ii, m.id
X0[i] = par2vec(m) -- save current orbit
X0[i] = par2vec(m) -- save current orbit
dH[i] = costp * X0[i]:norm() -- save current diff step
if dH[i] == 0 then dH[i] = max(costp, cotol) end
end
mflw.npar, mflw.tpar, n = 7*n, 7*n, 7*n

-- Note: coid = primary particle original id (i.e. from cofind)
-- id = primary and secondary particle id (i.e. in track)

-- save finite differences, final translation, temporaries
local dH, O1, X, R = par2vec(codiff), par2vec(self.O1), vector(6), matrix(6)
local O1, X, R = par2vec(self.O1), vector(6), matrix(6)

-- search for fix points
for itr=1,coiter do
for itr=1,coitr do
-- sanity check
assert(n%7 == 0, "unexpected corrupted set of particle blocks")
-- for i=1,n do print("itr="..itr, mflw[i]) end
Expand All @@ -159,7 +155,7 @@ local function cofind_jac (self, mflw)
m = mflw[i+j]
assert(m.coid == coid, "unexpected corrupted set of particle blocks")
vec2par(X0[coid], m)
m[vn[j]] = X0[coid][j] + dH[j]
m[vn[j]] = X0[coid][j] + dH[coid]
end
end

Expand All @@ -170,7 +166,7 @@ local function cofind_jac (self, mflw)
if n ~= mflw.npar then
warn("cofind: lost %d particles at iteration %d", n-mflw.npar, itr)
-- save information (spos, turn and status set by lostpar)
for i=mflw.npar+1,n do mflw[i].coiter = itr end
for i=mflw.npar+1,n do mflw[i].coitr = itr end
n = mflw.npar

-- filter out particles with at least one lost particle in their block
Expand Down Expand Up @@ -198,9 +194,9 @@ local function cofind_jac (self, mflw)
local m = mflw[i]

-- retrieve previous orbit, current orbit and jacobian
local X0 = X0[m.coid] ; par2vec(m, X)
local X0, dh = X0[m.coid], dH[m.coid] ; par2vec(m, X)
for j=1,6 do ; for k=1,6 do -- compute jacobian R_jk = df(x_j)/dx_k
R:set(j,k, (mflw[i+k][vn[j]] - X[j])/dH[k])
R:set(j,k, (mflw[i+k][vn[j]] - X[j])/dh)
end end

-- update X0 = X0-dx if |dx| > cotol, where dx solves (R-I)dx = (X-O1)-X0
Expand All @@ -212,12 +208,12 @@ local function cofind_jac (self, mflw)

if typ then -- "stable/singular"
if typ == "stable"
then vec2par(X0, m)
then vec2par(X0, m) ; dH[m.coid] = costp * X0:norm()
else warn("cofind: singular matrix (rnk=%d) at iteration %d \z
for particle %d.", rnk, itr, m.coid)
end
-- save information in stable/singular primary particle
m.rank, m.status, m.coiter = rnk, typ, itr
m.rank, m.status, m.coitr = rnk, typ, itr
-- swap with last tracked block
for j=0,6 do mflw[i+j], mflw[n-6+j] = mflw[n-6+j], mflw[i+j] end
n = n - 7
Expand All @@ -238,9 +234,9 @@ local function cofind_jac (self, mflw)

-- 5. mark remaining particles as unstable
if n ~= 0 then
warn("cofind: closed orbit(s) did not converge in %d iterations", coiter)
warn("cofind: closed orbit(s) did not converge in %d iterations", coitr)
for i=1,n do
mflw[i].coiter, mflw[i].status = coiter, "unstable"
mflw[i].coitr, mflw[i].status = coitr, "unstable"
end
end

Expand All @@ -266,7 +262,7 @@ end
-- cofind using map -----------------------------------------------------------o

local function cofind_map (self, mflw)
local cotol, coiter, totalpath in self
local coitr, cotol, totalpath in self

-- save current orbits and truncate computation at order 1
local X0, n = table.new(mflw.npar,0), mflw.npar
Expand All @@ -280,7 +276,7 @@ local function cofind_map (self, mflw)
local O1, X, R = par2vec(self.O1), vector(6), matrix(6)

-- search for fix points
for itr=1,coiter do
for itr=1,coitr do

-- 1. set order 0 to orbit, order 1 to I, higher orders to 0 (if any)
for i=1,n do mflw[i]:setvar(X0[mflw[i].id]) end
Expand All @@ -292,7 +288,7 @@ local function cofind_map (self, mflw)
if n ~= mflw.npar then
warn("cofind: lost %d particles at iteration %d", n-mflw.npar, itr)
-- save information (spos, turn and status set by lostpar)
for i=mflw.npar+1,n do mflw[i].coiter = itr end
for i=mflw.npar+1,n do mflw[i].coitr = itr end
n = mflw.npar
end

Expand All @@ -319,7 +315,7 @@ local function cofind_map (self, mflw)
for damap %d.", rnk, itr, m.id)
end
-- save information in stable/singular damap
m.rank, m.status, m.coiter = rnk, typ, itr
m.rank, m.status, m.coitr = rnk, typ, itr
-- swap with last tracked damap
mflw[i], mflw[n] = mflw[n], mflw[i]
n = n - 1
Expand All @@ -341,9 +337,9 @@ local function cofind_map (self, mflw)

-- 5. mark remaining damap as unstable
if n ~= 0 then
warn("cofind: closed orbit(s) did not converge in %d iterations", coiter)
warn("cofind: closed orbit(s) did not converge in %d iterations", coitr)
for i=1,n do
mflw[i].coiter, mflw[i].status = coiter, "unstable"
mflw[i].coitr, mflw[i].status = coitr, "unstable"
end
end

Expand All @@ -363,10 +359,10 @@ end
-- output status: stable, unstable, singular, lost.

local function exec (self)
local mapdef, radiate, cotol, coiter, codiff in self
assertf(is_positive(cotol ), "invalid cotol %.15g (positive number expected)" , cotol )
assertf(is_positive(coiter), "invalid coiter %d (positive number expected)" , coiter)
assertf(is_positive(codiff), "invalid codiff %.15g (positive number expected)", codiff)
local mapdef, radiate, coitr, cotol, costp in self
assertf(is_positive(coitr), "invalid coitr %d (positive number expected)" , coitr)
assertf(is_positive(cotol), "invalid cotol %.15g (positive number expected)", cotol)
assertf(is_positive(costp), "invalid costp %.15g (positive number expected)", costp)

-- prepare template for tracking, block quantum radiation and photon tracking
local _, mflw = track { exec=false } :copy_variables(self)
Expand Down Expand Up @@ -425,9 +421,9 @@ local cofind = command 'cofind' {
savesel=nil, -- save selector (predicate) (trck)
apersel=nil, -- aper selector (predicate, default atentry) (trck)

coiter=20, -- maximum number of iterations (cofn)
coitr=20, -- maximum number of iterations (cofn)
cotol=1e-8, -- closed orbit tolerance (i.e. min |dx|) (cofn)
codiff=1e-10, -- finite differences step for jacobian (cofn)
costp=1e-8, -- relative finite differences step for jacobian (cofn)
O1=0, -- optional final coordinates translation (cofn)

info=nil, -- information level (output on terminal) (trck)
Expand All @@ -438,9 +434,7 @@ local cofind = command 'cofind' {
exec=exec, -- command to execute upon children creation

__attr = tblcat( -- list of all setup attributes
track.__attr,
{'coiter', 'cotol', 'codiff', 'O1'},
{noeval=track.__attr.noeval}
track.__attr, {}, {noeval=track.__attr.noeval}
)
} :set_readonly() -- reference cofind command is readonly

Expand Down
Loading

0 comments on commit a06df86

Please sign in to comment.