00001 00033 #include <itpp/base/specmat.h> 00034 #include <itpp/base/elmatfunc.h> 00035 #include <itpp/base/matfunc.h> 00036 #include <itpp/base/stat.h> 00037 00038 00039 namespace itpp { 00040 00041 ivec find(const bvec &invector) 00042 { 00043 it_assert(invector.size()>0,"find(): vector cannot be empty"); 00044 ivec temp(invector.size()); 00045 int pos=0; 00046 for (int i=0;i<invector.size();i++) { 00047 if (invector(i)==bin(1)) { 00048 temp(pos)=i;pos++; 00049 } 00050 } 00051 temp.set_size(pos, true); 00052 return temp; 00053 } 00054 00055 #ifndef DOXYGEN_SHOULD_SKIP_THIS 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 00105 #endif //DOXYGEN_SHOULD_SKIP_THIS 00106 00107 vec impulse(int size) { 00108 vec t(size); 00109 t.clear(); 00110 t[0]=1.0; 00111 return t; 00112 } 00113 00114 vec linspace(double from, double to, int points) { 00115 if (points<2){ 00116 // This is the "Matlab definition" of linspace 00117 vec output(1); 00118 output(0)=to; 00119 return output; 00120 } 00121 else{ 00122 vec output(points); 00123 double step = (to - from) / double(points-1); 00124 int i; 00125 for (i=0; i<points; i++) 00126 output(i) = from + i*step; 00127 return output; 00128 } 00129 } 00130 00131 00132 // Construct a Hadamard-imat of size "size" 00133 imat hadamard(int size) { 00134 it_assert(size > 0, "hadamard(): size is not a power of 2"); 00135 int logsize = ceil_i(log2(static_cast<double>(size))); 00136 it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2"); 00137 00138 imat H(size, size); H(0,0) = 1; 00139 00140 for (int i = 0; i < logsize; ++i) { 00141 int pow2 = 1 << i; 00142 for (int k = 0; k < pow2; ++k) { 00143 for (int l = 0; l < pow2; ++l) { 00144 H(k, l) = H(k, l); 00145 H(k+pow2, l) = H(k, l); 00146 H(k, l+pow2) = H(k, l); 00147 H(k+pow2, l+pow2) = (-1) * H(k, l); 00148 } 00149 } 00150 } 00151 return H; 00152 } 00153 00154 imat jacobsthal(int p) 00155 { 00156 int quadratic_residue; 00157 imat out(p,p); 00158 int i, j; 00159 00160 out = -1; // start with all elements equal to "-1" 00161 00162 // Generate a complete list of quadratic residues 00163 for (i=0; i<(p-1)/2; i++) { 00164 quadratic_residue=((i+1)*(i+1))%p; 00165 // set this element in all rows (col-row) = quadratic_residue 00166 for (j=0; j<p; j++) { 00167 out(j, (j+quadratic_residue)%p)=1; 00168 } 00169 } 00170 00171 // set diagonal elements to zero 00172 for (i=0; i<p; i++) { 00173 out(i,i)=0; 00174 } 00175 return out; 00176 } 00177 00178 imat conference(int n) 00179 { 00180 it_assert1(n%4 == 2, "conference(int n); wrong size"); 00181 int pm=n-1; // p must be odd prime, not checked 00182 imat out(n,n); 00183 00184 out.set_submatrix(1,n-1,1,n-1, jacobsthal(pm)); 00185 out.set_submatrix(0,0,1,n-1, 1); 00186 out.set_submatrix(1,n-1,0,0, 1); 00187 out(0,0)=0; 00188 00189 return out; 00190 } 00191 00192 cmat toeplitz(const cvec &c, const cvec &r) { 00193 int size = c.size(); 00194 it_assert(size == r.size(), 00195 "toeplitz(): Incorrect sizes of input vectors."); 00196 cmat output(size, size); 00197 cvec c_conj = conj(c); 00198 c_conj[0] = conj(c_conj[0]); 00199 00200 for (int i = 0; i < size; i++) { 00201 cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1); 00202 output.set_submatrix(i, size - 1, i, i, tmp); 00203 } 00204 for (int i = 0; i < size - 1; i++) { 00205 cmat tmp = reshape(r(1, size - 1 - i), 1, size - 1 - i); 00206 output.set_submatrix(i, i, i + 1, size - 1, tmp); 00207 } 00208 00209 return output; 00210 } 00211 00212 cmat toeplitz(const cvec &c) { 00213 int size = c.size(); 00214 cmat output(size, size); 00215 cvec c_conj = conj(c); 00216 c_conj[0] = conj(c_conj[0]); 00217 00218 for (int i = 0; i < size; i++) { 00219 cmat tmp = reshape(c_conj(0, size - 1 - i), size - i, 1); 00220 output.set_submatrix(i, size - 1, i, i, tmp); 00221 } 00222 for (int i = 0; i < size - 1; i++) { 00223 cmat tmp = reshape(c(1, size - 1 - i), 1, size - 1 - i); 00224 output.set_submatrix(i, i, i + 1, size - 1, tmp); 00225 } 00226 00227 return output; 00228 } 00229 00230 mat toeplitz(const vec &c, const vec &r) { 00231 00232 mat output(c.size(), r.size()); 00233 00234 for (int i=0; i<c.size(); i++) { 00235 for(int j=0; j<std::min(r.size(), c.size()-i); j++) 00236 output(i+j,j) = c(i); 00237 } 00238 00239 for (int j=1; j<r.size(); j++) { 00240 for(int i=0; i<std::min(c.size(), r.size()-j); i++) 00241 output(i,i+j) = r(j); 00242 } 00243 00244 return output; 00245 } 00246 00247 mat toeplitz(const vec &c) { 00248 mat output(c.size(), c.size()); 00249 00250 for (int i=0; i<c.size(); i++) { 00251 for(int j=0; j<c.size()-i; j++) 00252 output(i+j,j) = c(i); 00253 } 00254 00255 for (int j=1; j<c.size(); j++) { 00256 for(int i=0; i<c.size()-j; i++) 00257 output(i,i+j) = c(j); 00258 } 00259 00260 return output; 00261 } 00262 00263 mat rotation_matrix(int dim, int plane1, int plane2, double angle) 00264 { 00265 mat m; 00266 double c = std::cos(angle), s = std::sin(angle); 00267 00268 it_assert(plane1>=0 && plane2>=0 && 00269 plane1<dim && plane2<dim && plane1!=plane2, 00270 "Invalid arguments to rotation_matrix()"); 00271 00272 m.set_size(dim, dim, false); 00273 m = 0.0; 00274 for (int i=0; i<dim; i++) 00275 m(i,i) = 1.0; 00276 00277 m(plane1, plane1) = c; 00278 m(plane1, plane2) = -s; 00279 m(plane2, plane1) = s; 00280 m(plane2, plane2) = c; 00281 00282 return m; 00283 } 00284 00285 void house(const vec &x, vec &v, double &beta) 00286 { 00287 double sigma, mu; 00288 int n = x.size(); 00289 00290 v = x; 00291 if (n == 1) { 00292 v(0) = 1.0; 00293 beta = 0.0; 00294 return; 00295 } 00296 sigma = energy(x(1, n-1)); 00297 v(0) = 1.0; 00298 if (sigma == 0.0) 00299 beta = 0.0; 00300 else { 00301 mu = std::sqrt(sqr(x(0)) + sigma); 00302 if (x(0) <= 0.0) 00303 v(0) = x(0) - mu; 00304 else 00305 v(0) = -sigma / (x(0) + mu); 00306 beta = 2 * sqr(v(0)) / (sigma + sqr(v(0))); 00307 v /= v(0); 00308 } 00309 } 00310 00311 void givens(double a, double b, double &c, double &s) 00312 { 00313 double t; 00314 00315 if (b == 0) { 00316 c = 1.0; 00317 s = 0.0; 00318 } 00319 else { 00320 if (fabs(b) > fabs(a)) { 00321 t = -a/b; 00322 s = -1.0 / std::sqrt(1 + t*t); 00323 c = s * t; 00324 } 00325 else { 00326 t = -b/a; 00327 c = 1.0 / std::sqrt(1 + t*t); 00328 s = c * t; 00329 } 00330 } 00331 } 00332 00333 void givens(double a, double b, mat &m) 00334 { 00335 double t, c, s; 00336 00337 m.set_size(2,2); 00338 00339 if (b == 0) { 00340 m(0,0) = 1.0; 00341 m(1,1) = 1.0; 00342 m(1,0) = 0.0; 00343 m(0,1) = 0.0; 00344 } 00345 else { 00346 if (fabs(b) > fabs(a)) { 00347 t = -a/b; 00348 s = -1.0 / std::sqrt(1 + t*t); 00349 c = s * t; 00350 } 00351 else { 00352 t = -b/a; 00353 c = 1.0 / std::sqrt(1 + t*t); 00354 s = c * t; 00355 } 00356 m(0,0) = c; 00357 m(1,1) = c; 00358 m(0,1) = s; 00359 m(1,0) = -s; 00360 } 00361 } 00362 00363 mat givens(double a, double b) 00364 { 00365 mat m(2,2); 00366 givens(a, b, m); 00367 return m; 00368 } 00369 00370 00371 void givens_t(double a, double b, mat &m) 00372 { 00373 double t, c, s; 00374 00375 m.set_size(2,2); 00376 00377 if (b == 0) { 00378 m(0,0) = 1.0; 00379 m(1,1) = 1.0; 00380 m(1,0) = 0.0; 00381 m(0,1) = 0.0; 00382 } 00383 else { 00384 if (fabs(b) > fabs(a)) { 00385 t = -a/b; 00386 s = -1.0 / std::sqrt(1 + t*t); 00387 c = s * t; 00388 } 00389 else { 00390 t = -b/a; 00391 c = 1.0 / std::sqrt(1 + t*t); 00392 s = c * t; 00393 } 00394 m(0,0) = c; 00395 m(1,1) = c; 00396 m(0,1) = -s; 00397 m(1,0) = s; 00398 } 00399 } 00400 00401 mat givens_t(double a, double b) 00402 { 00403 mat m(2,2); 00404 givens_t(a, b, m); 00405 return m; 00406 } 00407 00409 template void eye(int, mat &); 00411 template void eye(int, bmat &); 00413 template void eye(int, imat &); 00415 template void eye(int, cmat &); 00416 00417 } // namespace itpp
Generated on Fri Jan 11 08:51:37 2008 for IT++ by Doxygen 1.3.9.1