00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_ANALYTIC_MATRIX_HPP
00014 #define __MMX_ANALYTIC_MATRIX_HPP
00015 #include <continewz/analytic_vector.hpp>
00016 #include <algebramix/series_matrix.hpp>
00017 namespace mmx {
00018 #define TMPL_DEF template<typename C, typename V>
00019 #define TMPL template<typename C, typename V>
00020 #define Analytic analytic<C,V>
00021 #define Analytic_rep analytic_rep<C,V>
00022 #define Assume assume_bound<typename unvectorize<R >::val>
00023 #define Scalar typename unvectorize<C>::val
00024 #define Radius Abs_type(Scalar)
00025 #define Series series<C>
00026 #define R Abs_type(C)
00027 #define MR Abs_type(matrix<C> )
00028
00029
00030
00031
00032
00033 TMPL Analytic access (const analytic<matrix<C>,V>& f, nat i, nat j);
00034 TMPL matrix<Analytic > as_matrix (const analytic<matrix<C>,V>& f);
00035 TMPL analytic<matrix<C>,V> as_analytic (const matrix<Analytic >& m);
00036
00037 TMPL analytic<matrix<C>,V>
00038 solve_matrix_lde_init (const analytic<matrix<C>,V>& f, const matrix<C>& c);
00039 TMPL matrix<Analytic >
00040 solve_lde (const matrix<Analytic >& f);
00041 TMPL matrix<Analytic >
00042 solve_lde_init (const matrix<Analytic >& f, const matrix<C>& c);
00043
00044
00045
00046
00047
00048 #define Analytic_matrix analytic<matrix<C>,V>
00049 #define Matrix_analytic matrix<analytic<C,V> >
00050 STYPE_TO_TYPE(TMPL,as_matrix_type,Analytic_matrix,Matrix_analytic);
00051 #undef Matrix_analytic
00052 #undef Analytic_matrix
00053
00054 template<typename C>
00055 struct unvectorize<matrix<C> > {
00056
00057 typedef C val;
00058 static inline vector<val> encode (const matrix<C>& x) {
00059 nat nr= rows (x), nc= cols (x);
00060 vector<val> r (C(0), nr*nc);
00061 for (nat i=0; i<nr; i++)
00062 for (nat j=0; j<nc; j++)
00063 r[i*nc+j]= x(i,j);
00064 return r; }
00065 static inline matrix<C>
00066 decode (const vector<val>& v, const matrix<C>& gauge) {
00067 nat nr= rows (gauge), nc= cols (gauge);
00068 matrix<C> r (C(0), nr, nc);
00069 for (nat i=0; i<nr; i++)
00070 for (nat j=0; j<nc; j++)
00071 r (i, j)= v[i*nc+j];
00072 return r; }
00073 };
00074
00075 TMPL
00076 class matrix_access_analytic_rep: public Analytic_rep {
00077 protected:
00078 const analytic<matrix<C>,V> f;
00079 nat i, j;
00080 public:
00081 matrix_access_analytic_rep
00082 (const analytic<matrix<C>,V>& f2, nat i2, nat j2):
00083 Analytic_rep (get_format1 (CF(f2))), f (f2), i (i2), j (j2) {}
00084 void Clear_cache (nat which) const {
00085 Analytic_rep::Clear_cache (which);
00086 clear_cache (f, which); }
00087 Series Expand () const {
00088 return access (expand (f), i, j); }
00089 Radius Radius_bound (nat order) const {
00090 return radius_bound (f, order); }
00091 R Tail_bound (const Radius& r, nat order, Assume& a) const {
00092 return read (tail_bound (this->f, r, order, a), i, j); }
00093 Analytic Move (const Scalar& z) const {
00094 return access (move (f, z), i, j); }
00095 C Eval (const Scalar& z) const {
00096 return eval (f, z) (i, j); }
00097 Analytic Derive () const {
00098 return access (derive (f), i, j); }
00099 };
00100
00101 TMPL Analytic
00102 access (const analytic<matrix<C>,V>& f, nat i, nat j) {
00103 return (Analytic_rep*) new matrix_access_analytic_rep<C,V> (f, i, j);
00104 }
00105
00106 TMPL matrix<Analytic >
00107 as_matrix (const analytic<matrix<C>,V>& f) {
00108 nat nr= rows (f[0]), nc= cols (f[0]);
00109 matrix<Analytic > m (Analytic (0, get_format1 (CF(f))), nr, nc);
00110 for (nat i=0; i<nr; i++)
00111 for (nat j=0; j<nc; j++)
00112 m (i, j)= access (f, i, j);
00113 return m;
00114 }
00115
00116 TMPL matrix<R >
00117 head_extremum (const analytic<matrix<C>,V>& f,
00118 const Radius& r, nat order, int type) {
00119 nat i, j, nr= rows (f[0]), nc= cols (f[0]);
00120 matrix<R > b (R (0), nr, nc);
00121 for (i=0; i<nr; i++)
00122 for (j=0; j<nc; j++)
00123 b (i, j)= head_extremum (access (f, i, j), r, order, type);
00124 return b;
00125 }
00126
00127 TMPL
00128 class matrix_analytic_rep: public analytic_rep<matrix<C>,V> {
00129 const matrix<Analytic > m;
00130 public:
00131 matrix_analytic_rep (const matrix<Analytic >& m2):
00132 analytic_rep<matrix<C>,V> (format<matrix<C> > (get_format1 (CF(m2)))),
00133 m (m2) {}
00134 void Clear_cache (nat which) const {
00135 analytic_rep<vector<C>,V>::Clear_cache (which);
00136 for (nat i=0; i<rows (m); i++)
00137 for (nat j=0; j<cols (m); j++)
00138 clear_cache (m (i, j), which); }
00139 series<matrix<C> > Expand () const {
00140 nat i, j, nr= rows (m), nc= cols (m);
00141 matrix<Series > r (Series (0), nr, nc);
00142 for (i=0; i<nr; i++)
00143 for (j=0; j<nc; j++)
00144 r (i, j)= expand (m (i, j));
00145 return as_series (r); }
00146 Radius Radius_bound (nat order) const {
00147 Radius r= Maximal (Radius);
00148 for (nat i=0; i<rows (m); i++)
00149 for (nat j=0; j<cols (m); j++)
00150 r= min (r, radius_bound (m (i, j), order));
00151 return r; }
00152 MR Tail_bound (const Radius& r, nat order, Assume& a) const {
00153 nat i, j, nr= rows (m), nc= cols (m);
00154 MR b (R (0), nr, nc);
00155 for (i=0; i<nr; i++)
00156 for (j=0; j<nc; j++)
00157 b (i, j)= tail_bound (m (i, j), r, order, a);
00158 return b; }
00159 analytic<matrix<C>,V> Move (const Scalar& z) const {
00160 nat i, j, nr= rows (m), nc= cols (m);
00161 matrix<Analytic > r (Analytic (0, get_format1 (this->tfm())), nr, nc);
00162 for (i=0; i<nr; i++)
00163 for (j=0; j<nc; j++)
00164 r (i,j)= move (m (i, j), z);
00165 return as_analytic (r); }
00166 matrix<C> Eval (const Scalar& z) const {
00167 nat i, j, nr= rows (m), nc= cols (m);
00168 matrix<C> r (C (0), nr, nc);
00169 for (i=0; i<nr; i++)
00170 for (j=0; j<nc; j++)
00171 r (i, j)= eval (m (i, j), z);
00172 return r; }
00173 analytic<matrix<C>,V> Derive () const {
00174 return as_analytic (derive (m)); }
00175 };
00176
00177 TMPL analytic<matrix<C>,V>
00178 as_analytic (const matrix<Analytic >& m) {
00179 return (analytic_rep<matrix<C>,V>*) new matrix_analytic_rep<C,V> (m);
00180 }
00181
00182 TMPL analytic<matrix<C>,V>
00183 from_matrix (const matrix<Analytic >& m) {
00184 return (analytic_rep<matrix<C>,V>*) new matrix_analytic_rep<C,V> (m);
00185 }
00186
00187 TMPL matrix<Analytic >
00188 move (const matrix<Analytic >& m, const C& z) {
00189 nat nr= rows (m), nc= cols (m);
00190 matrix<Analytic > r (Analytic (0, get_format (z)), nr, nc);
00191 for (nat i=0; i<nr; i++)
00192 for (nat j=0; j<nc; j++)
00193 r (i, j)= move (m (i, j), z);
00194 return r;
00195 }
00196
00197 TMPL matrix<C>
00198 eval (const matrix<Analytic >& m, const C& z) {
00199 nat nr= rows (m), nc= cols (m);
00200 matrix<C> r (C (0), nr, nc);
00201 for (nat i=0; i<nr; i++)
00202 for (nat j=0; j<nc; j++)
00203 r (i, j)= eval (m (i, j), z);
00204 return r;
00205 }
00206
00207
00208
00209
00210
00211 TMPL analytic<matrix<C>,V>
00212 solve_matrix_lde_init (const analytic<matrix<C>,V>& f, const matrix<C>& c) {
00213 return unary_recursive_analytic<solve_matrix_lde_op> (f, c);
00214 }
00215
00216 TMPL matrix<Analytic >
00217 solve_lde (const matrix<Analytic >& f) {
00218 ASSERT (is_square_matrix (f), "square matrix expected");
00219 matrix<C> c (C(1), rows (f), cols (f));
00220 return as_matrix (solve_matrix_lde_init (as_analytic (f), c));
00221 }
00222
00223 TMPL matrix<Analytic >
00224 solve_lde_init (const matrix<Analytic >& f, const matrix<C>& c) {
00225 ASSERT (is_square_matrix (f), "square matrix expected");
00226 ASSERT (rows (f) == rows (c), "dimensions don't match");
00227 return as_matrix (solve_matrix_lde_init (as_analytic (f), c));
00228 }
00229
00230
00231
00232
00233
00234 template<typename C> list<C>
00235 straight_path (const C& start, const C& end, nat steps= 5) {
00236 list<C> r;
00237 for (nat i=steps; i>0; i--)
00238 r= cons (start + (i*(end-start)) / steps, r);
00239 return r;
00240 }
00241
00242 template<typename C> list<C>
00243 circle_path (const C& start, const C& center, int orient= 1, nat steps= 17) {
00244 list<C> r;
00245 for (int i=0; ((nat) i)<steps; i++) {
00246 C z= exp ((C(-orient*2*i) * Pi(C) * Imaginary(C)) / C(steps));
00247 r= cons (center + z * (start - center), r);
00248 }
00249 return r;
00250 }
00251
00252 template<typename C> list<C>
00253 delta_to_absolute (const list<C>& l, const C& pos= C(0)) {
00254 if (is_nil (l)) return l;
00255 else return cons (pos + car (l), delta_to_absolute (cdr (l), pos + car (l)));
00256 }
00257
00258 template<typename C> list<C>
00259 absolute_to_delta (const list<C>& l, const C& pos= C(0)) {
00260 if (is_nil (l)) return l;
00261 else return cons (car (l) - pos, absolute_to_delta (cdr (l), car (l)));
00262 }
00263
00264 template<typename T, typename C> T
00265 move (const T& x, const list<C>& p) {
00266 if (is_nil (p)) return x;
00267 else return move (move (x, car (p)), cdr (p));
00268 }
00269
00270 TMPL matrix<C>
00271 connection_matrix (const matrix<Analytic >& m, const list<C>& p) {
00272 nat n= N(p);
00273 if (n <= 1) {
00274 ASSERT (is_square_matrix (m), "square matrix expected");
00275 if (n == 0) return matrix<C> (C(1), rows (m), cols (m));
00276 else return eval (solve_lde (m), car (p));
00277 }
00278 else {
00279 list<C> p1= range (p, 0, n>>1);
00280 list<C> p2= range (p, n>>1, n);
00281 matrix<C> r1= connection_matrix (m, p1);
00282 matrix<C> r2= connection_matrix (move (m, p1), p2);
00283 return r2 * r1;
00284 }
00285 }
00286
00287 #undef TMPL_DEF
00288 #undef TMPL
00289 #undef Analytic
00290 #undef Analytic_rep
00291 #undef Assume
00292 #undef Scalar
00293 #undef Radius
00294 #undef Series
00295 #undef R
00296 #undef MR
00297 }
00298 #endif // __MMX_ANALYTIC_MATRIX_HPP
