00001 00031 #include <itpp/base/specmat.h> 00032 #include <itpp/base/math/elem_math.h> 00033 #include <itpp/base/math/log_exp.h> 00034 #include <itpp/base/matfunc.h> 00035 00036 00037 namespace itpp 00038 { 00039 00040 ivec find(const bvec &invector) 00041 { 00042 it_assert(invector.size() > 0, "find(): vector cannot be empty"); 00043 ivec temp(invector.size()); 00044 int pos = 0; 00045 for (int i = 0;i < invector.size();i++) { 00046 if (invector(i) == bin(1)) { 00047 temp(pos) = i; 00048 pos++; 00049 } 00050 } 00051 temp.set_size(pos, true); 00052 return temp; 00053 } 00054 00056 00057 #define CREATE_SET_FUNS(typef,typem,name,value) \ 00058 typef name(int size) \ 00059 { \ 00060 typef t(size); \ 00061 t = value; \ 00062 return t; \ 00063 } \ 00064 \ 00065 typem name(int rows, int cols) \ 00066 { \ 00067 typem t(rows, cols); \ 00068 t = value; \ 00069 return t; \ 00070 } 00071 00072 #define CREATE_EYE_FUN(type,name,zero,one) \ 00073 type name(int size) { \ 00074 type t(size,size); \ 00075 t = zero; \ 00076 for (int i=0; i<size; i++) \ 00077 t(i,i) = one; \ 00078 return t; \ 00079 } 00080 00081 CREATE_SET_FUNS(vec, mat, ones, 1.0) 00082 CREATE_SET_FUNS(bvec, bmat, ones_b, bin(1)) 00083 CREATE_SET_FUNS(ivec, imat, ones_i, 1) 00084 CREATE_SET_FUNS(cvec, cmat, ones_c, std::complex<double>(1.0)) 00085 00086 CREATE_SET_FUNS(vec, mat, zeros, 0.0) 00087 CREATE_SET_FUNS(bvec, bmat, zeros_b, bin(0)) 00088 CREATE_SET_FUNS(ivec, imat, zeros_i, 0) 00089 CREATE_SET_FUNS(cvec, cmat, zeros_c, std::complex<double>(0.0)) 00090 00091 CREATE_EYE_FUN(mat, eye, 0.0, 1.0) 00092 CREATE_EYE_FUN(bmat, eye_b, bin(0), bin(1)) 00093 CREATE_EYE_FUN(imat, eye_i, 0, 1) 00094 CREATE_EYE_FUN(cmat, eye_c, std::complex<double>(0.0), std::complex<double>(1.0)) 00095 00096 template <class T> 00097 void eye(int size, Mat<T> &m) 00098 { 00099 m.set_size(size, size, false); 00100 m = T(0); 00101 for (int i = size - 1; i >= 0; i--) 00102 m(i, i) = T(1); 00103 } 00104 00106 00107 vec impulse(int size) 00108 { 00109 vec t(size); 00110 t.clear(); 00111 t[0] = 1.0; 00112 return t; 00113 } 00114 00115 vec linspace(double from, double to, int points) 00116 { 00117 if (points < 2) { 00118 // This is the "Matlab definition" of linspace 00119 vec output(1); 00120 output(0) = to; 00121 return output; 00122 } 00123 else { 00124 vec output(points); 00125 double step = (to - from) / double(points - 1); 00126 int i; 00127 for (i = 0; i < points; i++) 00128 output(i) = from + i * step; 00129 return output; 00130 } 00131 } 00132 00133 vec zigzag_space(double t0, double t1, int K) 00134 { 00135 it_assert(K > 0, "zigzag_space:() K must be positive"); 00136 ivec N = "0 1"; 00137 00138 int n = 2; 00139 for (int k = 0; k < K; k++) { 00140 ivec Nn = 2 * N; 00141 for (int i = 1; i < length(Nn); i += 2) { 00142 Nn = concat(Nn, i); 00143 n++; 00144 } 00145 N = Nn; 00146 } 00147 00148 vec T0 = linspace(t0, t1, n); 00149 vec Tt = zeros(n); 00150 for (int i = 0; i < n; i++) { 00151 Tt(i) = T0(N(i)); 00152 } 00153 return Tt; 00154 } 00155 00156 // Construct a Hadamard-imat of size "size" 00157 imat hadamard(int size) 00158 { 00159 it_assert(size > 0, "hadamard(): size is not a power of 2"); 00160 int logsize = ceil_i(::log2(static_cast<double>(size))); 00161 it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2"); 00162 00163 imat H(size, size); 00164 H(0, 0) = 1; 00165 00166 for (int i = 0; i < logsize; ++i) { 00167 int pow2 = 1 << i; 00168 for (int k = 0; k < pow2; ++k) { 00169 for (int l = 0; l < pow2; ++l) { 00170 H(k, l) = H(k, l); 00171 H(k + pow2, l) = H(k, l); 00172 H(k, l + pow2) = H(k, l); 00173 H(k + pow2, l + pow2) = (-1) * H(k, l); 00174 } 00175 } 00176 } 00177 return H; 00178 } 00179 00180 imat jacobsthal(int p) 00181 { 00182 int quadratic_residue; 00183 imat out(p, p); 00184 int i, j; 00185 00186 out = -1; // start with all elements equal to "-1" 00187 00188 // Generate a complete list of quadratic residues 00189 for (i = 0; i < (p - 1) / 2; i++) { 00190 quadratic_residue = ((i + 1) * (i + 1)) % p; 00191 // set this element in all rows (col-row) = quadratic_residue 00192 for (j = 0; j < p; j++) { 00193 out(j, (j + quadratic_residue) % p) = 1; 00194 } 00195 } 00196 00197 // set diagonal elements to zero 00198 for (i = 0; i < p; i++) { 00199 out(i, i) = 0; 00200 } 00201 return out; 00202 } 00203 00204 imat conference(int n) 00205 { 00206 it_assert_debug(n % 4 == 2, "conference(int n); wrong size"); 00207 int pm = n - 1; // p must be odd prime, not checked 00208 imat out(n, n); 00209 00210 out.set_submatrix(1, n - 1, 1, n - 1, jacobsthal(pm)); 00211 out.set_submatrix(0, 0, 1, n - 1, 1); 00212 out.set_submatrix(1, n - 1, 0, 0, 1); 00213 out(0, 0) = 0; 00214 00215 return out; 00216 } 00217 00218 00219 template <> 00220 const cmat toeplitz(const cvec &c) 00221 { 00222 int s = c.size(); 00223 cmat output(s, s); 00224 cvec c_conj = conj(c); 00225 for (int i = 1; i < s; ++i) { 00226 for (int j = 0; j < s - i; ++j) { 00227 output(i + j, j) = c_conj(i); 00228 } 00229 } 00230 // start from j = 0 here, because the main diagonal is not conjugated 00231 for (int j = 0; j < s; ++j) { 00232 for (int i = 0; i < s - j; ++i) { 00233 output(i, i + j) = c(j); 00234 } 00235 } 00236 return output; 00237 } 00238 00239 00240 mat rotation_matrix(int dim, int plane1, int plane2, double angle) 00241 { 00242 mat m; 00243 double c = std::cos(angle), s = std::sin(angle); 00244 00245 it_assert(plane1 >= 0 && plane2 >= 0 && 00246 plane1 < dim && plane2 < dim && plane1 != plane2, 00247 "Invalid arguments to rotation_matrix()"); 00248 00249 m.set_size(dim, dim, false); 00250 m = 0.0; 00251 for (int i = 0; i < dim; i++) 00252 m(i, i) = 1.0; 00253 00254 m(plane1, plane1) = c; 00255 m(plane1, plane2) = -s; 00256 m(plane2, plane1) = s; 00257 m(plane2, plane2) = c; 00258 00259 return m; 00260 } 00261 00262 void house(const vec &x, vec &v, double &beta) 00263 { 00264 double sigma, mu; 00265 int n = x.size(); 00266 00267 v = x; 00268 if (n == 1) { 00269 v(0) = 1.0; 00270 beta = 0.0; 00271 return; 00272 } 00273 sigma = sum(sqr(x(1, n - 1))); 00274 v(0) = 1.0; 00275 if (sigma == 0.0) 00276 beta = 0.0; 00277 else { 00278 mu = std::sqrt(sqr(x(0)) + sigma); 00279 if (x(0) <= 0.0) 00280 v(0) = x(0) - mu; 00281 else 00282 v(0) = -sigma / (x(0) + mu); 00283 beta = 2 * sqr(v(0)) / (sigma + sqr(v(0))); 00284 v /= v(0); 00285 } 00286 } 00287 00288 void givens(double a, double b, double &c, double &s) 00289 { 00290 double t; 00291 00292 if (b == 0) { 00293 c = 1.0; 00294 s = 0.0; 00295 } 00296 else { 00297 if (fabs(b) > fabs(a)) { 00298 t = -a / b; 00299 s = -1.0 / std::sqrt(1 + t * t); 00300 c = s * t; 00301 } 00302 else { 00303 t = -b / a; 00304 c = 1.0 / std::sqrt(1 + t * t); 00305 s = c * t; 00306 } 00307 } 00308 } 00309 00310 void givens(double a, double b, mat &m) 00311 { 00312 double t, c, s; 00313 00314 m.set_size(2, 2); 00315 00316 if (b == 0) { 00317 m(0, 0) = 1.0; 00318 m(1, 1) = 1.0; 00319 m(1, 0) = 0.0; 00320 m(0, 1) = 0.0; 00321 } 00322 else { 00323 if (fabs(b) > fabs(a)) { 00324 t = -a / b; 00325 s = -1.0 / std::sqrt(1 + t * t); 00326 c = s * t; 00327 } 00328 else { 00329 t = -b / a; 00330 c = 1.0 / std::sqrt(1 + t * t); 00331 s = c * t; 00332 } 00333 m(0, 0) = c; 00334 m(1, 1) = c; 00335 m(0, 1) = s; 00336 m(1, 0) = -s; 00337 } 00338 } 00339 00340 mat givens(double a, double b) 00341 { 00342 mat m(2, 2); 00343 givens(a, b, m); 00344 return m; 00345 } 00346 00347 00348 void givens_t(double a, double b, mat &m) 00349 { 00350 double t, c, s; 00351 00352 m.set_size(2, 2); 00353 00354 if (b == 0) { 00355 m(0, 0) = 1.0; 00356 m(1, 1) = 1.0; 00357 m(1, 0) = 0.0; 00358 m(0, 1) = 0.0; 00359 } 00360 else { 00361 if (fabs(b) > fabs(a)) { 00362 t = -a / b; 00363 s = -1.0 / std::sqrt(1 + t * t); 00364 c = s * t; 00365 } 00366 else { 00367 t = -b / a; 00368 c = 1.0 / std::sqrt(1 + t * t); 00369 s = c * t; 00370 } 00371 m(0, 0) = c; 00372 m(1, 1) = c; 00373 m(0, 1) = -s; 00374 m(1, 0) = s; 00375 } 00376 } 00377 00378 mat givens_t(double a, double b) 00379 { 00380 mat m(2, 2); 00381 givens_t(a, b, m); 00382 return m; 00383 } 00384 00386 template void eye(int, mat &); 00388 template void eye(int, bmat &); 00390 template void eye(int, imat &); 00392 template void eye(int, cmat &); 00393 00394 } // namespace itpp
Generated on Sun Jul 26 08:54:26 2009 for IT++ by Doxygen 1.5.9