IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
specmat.cpp
Go to the documentation of this file.
1 
30 #include <itpp/base/specmat.h>
32 #include <itpp/base/math/log_exp.h>
33 #include <itpp/base/itcompat.h>
34 #include <itpp/base/matfunc.h>
35 
36 
37 namespace itpp
38 {
39 
40 ivec find(const bvec &invector)
41 {
42  it_assert(invector.size() > 0, "find(): vector cannot be empty");
43  ivec temp(invector.size());
44  int pos = 0;
45  for (int i = 0;i < invector.size();i++) {
46  if (invector(i) == bin(1)) {
47  temp(pos) = i;
48  pos++;
49  }
50  }
51  temp.set_size(pos, true);
52  return temp;
53 }
54 
56 
57 #define CREATE_SET_FUNS(typef,typem,name,value) \
58  typef name(int size) \
59  { \
60  typef t(size); \
61  t = value; \
62  return t; \
63  } \
64  \
65  typem name(int rows, int cols) \
66  { \
67  typem t(rows, cols); \
68  t = value; \
69  return t; \
70  }
71 
72 #define CREATE_EYE_FUN(type,name,zero,one) \
73  type name(int size) { \
74  type t(size,size); \
75  t = zero; \
76  for (int i=0; i<size; i++) \
77  t(i,i) = one; \
78  return t; \
79  }
80 
81 CREATE_SET_FUNS(vec, mat, ones, 1.0)
82 CREATE_SET_FUNS(bvec, bmat, ones_b, bin(1))
83 CREATE_SET_FUNS(ivec, imat, ones_i, 1)
84 CREATE_SET_FUNS(cvec, cmat, ones_c, std::complex<double>(1.0))
85 
86 CREATE_SET_FUNS(vec, mat, zeros, 0.0)
87 CREATE_SET_FUNS(bvec, bmat, zeros_b, bin(0))
88 CREATE_SET_FUNS(ivec, imat, zeros_i, 0)
89 CREATE_SET_FUNS(cvec, cmat, zeros_c, std::complex<double>(0.0))
90 
91 CREATE_EYE_FUN(mat, eye, 0.0, 1.0)
92 CREATE_EYE_FUN(bmat, eye_b, bin(0), bin(1))
93 CREATE_EYE_FUN(imat, eye_i, 0, 1)
94 CREATE_EYE_FUN(cmat, eye_c, std::complex<double>(0.0), std::complex<double>(1.0))
95 
96 template <class T>
97 void eye(int size, Mat<T> &m)
98 {
99  m.set_size(size, size, false);
100  m = T(0);
101  for (int i = size - 1; i >= 0; i--)
102  m(i, i) = T(1);
103 }
104 
106 
107 vec impulse(int size)
108 {
109  vec t(size);
110  t.clear();
111  t[0] = 1.0;
112  return t;
113 }
114 
115 vec linspace(double from, double to, int points)
116 {
117  if (points < 2) {
118  // This is the "Matlab definition" of linspace
119  vec output(1);
120  output(0) = to;
121  return output;
122  }
123  else {
124  vec output(points);
125  double step = (to - from) / double(points - 1);
126  int i;
127  for (i = 0; i < points; i++)
128  output(i) = from + i * step;
129  return output;
130  }
131 }
132 
133 vec zigzag_space(double t0, double t1, int K)
134 {
135  it_assert(K > 0, "zigzag_space:() K must be positive");
136  ivec N = "0 1";
137 
138  int n = 2;
139  for (int k = 0; k < K; k++) {
140  ivec Nn = 2 * N;
141  for (int i = 1; i < length(Nn); i += 2) {
142  Nn = concat(Nn, i);
143  n++;
144  }
145  N = Nn;
146  }
147 
148  vec T0 = linspace(t0, t1, n);
149  vec Tt = zeros(n);
150  for (int i = 0; i < n; i++) {
151  Tt(i) = T0(N(i));
152  }
153  return Tt;
154 }
155 
156 // Construct a Hadamard-imat of size "size"
157 imat hadamard(int size)
158 {
159  it_assert(size > 0, "hadamard(): size is not a power of 2");
160  int logsize = ceil_i(::log2(static_cast<double>(size)));
161  it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2");
162 
163  imat H(size, size);
164  H(0, 0) = 1;
165 
166  for (int i = 0; i < logsize; ++i) {
167  int pow2 = 1 << i;
168  for (int k = 0; k < pow2; ++k) {
169  for (int l = 0; l < pow2; ++l) {
170  H(k, l) = H(k, l);
171  H(k + pow2, l) = H(k, l);
172  H(k, l + pow2) = H(k, l);
173  H(k + pow2, l + pow2) = (-1) * H(k, l);
174  }
175  }
176  }
177  return H;
178 }
179 
180 imat jacobsthal(int p)
181 {
182  int quadratic_residue;
183  imat out(p, p);
184  int i, j;
185 
186  out = -1; // start with all elements equal to "-1"
187 
188  // Generate a complete list of quadratic residues
189  for (i = 0; i < (p - 1) / 2; i++) {
190  quadratic_residue = ((i + 1) * (i + 1)) % p;
191  // set this element in all rows (col-row) = quadratic_residue
192  for (j = 0; j < p; j++) {
193  out(j, (j + quadratic_residue) % p) = 1;
194  }
195  }
196 
197  // set diagonal elements to zero
198  for (i = 0; i < p; i++) {
199  out(i, i) = 0;
200  }
201  return out;
202 }
203 
204 imat conference(int n)
205 {
206  it_assert_debug(n % 4 == 2, "conference(int n); wrong size");
207  int pm = n - 1; // p must be odd prime, not checked
208  imat out(n, n);
209 
210  out.set_submatrix(1, 1, jacobsthal(pm));
211  out.set_submatrix(0, 0, 1, n - 1, 1);
212  out.set_submatrix(1, n - 1, 0, 0, 1);
213  out(0, 0) = 0;
214 
215  return out;
216 }
217 
218 
219 template <>
220 const cmat toeplitz(const cvec &c)
221 {
222  int s = c.size();
223  cmat output(s, s);
224  cvec c_conj = conj(c);
225  for (int i = 1; i < s; ++i) {
226  for (int j = 0; j < s - i; ++j) {
227  output(i + j, j) = c_conj(i);
228  }
229  }
230  // start from j = 0 here, because the main diagonal is not conjugated
231  for (int j = 0; j < s; ++j) {
232  for (int i = 0; i < s - j; ++i) {
233  output(i, i + j) = c(j);
234  }
235  }
236  return output;
237 }
238 
239 
240 mat rotation_matrix(int dim, int plane1, int plane2, double angle)
241 {
242  mat m;
243  double c = std::cos(angle), s = std::sin(angle);
244 
245  it_assert(plane1 >= 0 && plane2 >= 0 &&
246  plane1 < dim && plane2 < dim && plane1 != plane2,
247  "Invalid arguments to rotation_matrix()");
248 
249  m.set_size(dim, dim, false);
250  m = 0.0;
251  for (int i = 0; i < dim; i++)
252  m(i, i) = 1.0;
253 
254  m(plane1, plane1) = c;
255  m(plane1, plane2) = -s;
256  m(plane2, plane1) = s;
257  m(plane2, plane2) = c;
258 
259  return m;
260 }
261 
262 void house(const vec &x, vec &v, double &beta)
263 {
264  double sigma, mu;
265  int n = x.size();
266 
267  v = x;
268  if (n == 1) {
269  v(0) = 1.0;
270  beta = 0.0;
271  return;
272  }
273  sigma = sum(sqr(x(1, n - 1)));
274  v(0) = 1.0;
275  if (sigma == 0.0)
276  beta = 0.0;
277  else {
278  mu = std::sqrt(sqr(x(0)) + sigma);
279  if (x(0) <= 0.0)
280  v(0) = x(0) - mu;
281  else
282  v(0) = -sigma / (x(0) + mu);
283  beta = 2 * sqr(v(0)) / (sigma + sqr(v(0)));
284  v /= v(0);
285  }
286 }
287 
288 void givens(double a, double b, double &c, double &s)
289 {
290  double t;
291 
292  if (b == 0) {
293  c = 1.0;
294  s = 0.0;
295  }
296  else {
297  if (fabs(b) > fabs(a)) {
298  t = -a / b;
299  s = -1.0 / std::sqrt(1 + t * t);
300  c = s * t;
301  }
302  else {
303  t = -b / a;
304  c = 1.0 / std::sqrt(1 + t * t);
305  s = c * t;
306  }
307  }
308 }
309 
310 void givens(double a, double b, mat &m)
311 {
312  double t, c, s;
313 
314  m.set_size(2, 2);
315 
316  if (b == 0) {
317  m(0, 0) = 1.0;
318  m(1, 1) = 1.0;
319  m(1, 0) = 0.0;
320  m(0, 1) = 0.0;
321  }
322  else {
323  if (fabs(b) > fabs(a)) {
324  t = -a / b;
325  s = -1.0 / std::sqrt(1 + t * t);
326  c = s * t;
327  }
328  else {
329  t = -b / a;
330  c = 1.0 / std::sqrt(1 + t * t);
331  s = c * t;
332  }
333  m(0, 0) = c;
334  m(1, 1) = c;
335  m(0, 1) = s;
336  m(1, 0) = -s;
337  }
338 }
339 
340 mat givens(double a, double b)
341 {
342  mat m(2, 2);
343  givens(a, b, m);
344  return m;
345 }
346 
347 
348 void givens_t(double a, double b, mat &m)
349 {
350  double t, c, s;
351 
352  m.set_size(2, 2);
353 
354  if (b == 0) {
355  m(0, 0) = 1.0;
356  m(1, 1) = 1.0;
357  m(1, 0) = 0.0;
358  m(0, 1) = 0.0;
359  }
360  else {
361  if (fabs(b) > fabs(a)) {
362  t = -a / b;
363  s = -1.0 / std::sqrt(1 + t * t);
364  c = s * t;
365  }
366  else {
367  t = -b / a;
368  c = 1.0 / std::sqrt(1 + t * t);
369  s = c * t;
370  }
371  m(0, 0) = c;
372  m(1, 1) = c;
373  m(0, 1) = -s;
374  m(1, 0) = s;
375  }
376 }
377 
378 mat givens_t(double a, double b)
379 {
380  mat m(2, 2);
381  givens_t(a, b, m);
382  return m;
383 }
384 
386 template void eye(int, mat &);
388 template void eye(int, bmat &);
390 template void eye(int, imat &);
392 template void eye(int, cmat &);
393 
394 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2