1 | /*
|
---|
2 | * vector_ops.cpp
|
---|
3 | *
|
---|
4 | * Created on: Apr 1, 2010
|
---|
5 | * Author: crueger
|
---|
6 | */
|
---|
7 |
|
---|
8 | #include "vector.hpp"
|
---|
9 | #include "Plane.hpp"
|
---|
10 | #include "log.hpp"
|
---|
11 | #include "verbose.hpp"
|
---|
12 | #include "gslmatrix.hpp"
|
---|
13 | #include "leastsquaremin.hpp"
|
---|
14 | #include "info.hpp"
|
---|
15 | #include "Helpers/fast_functions.hpp"
|
---|
16 | #include "Exceptions/LinearDependenceException.hpp"
|
---|
17 | #include "Exceptions/SkewException.hpp"
|
---|
18 |
|
---|
19 | #include <gsl/gsl_linalg.h>
|
---|
20 | #include <gsl/gsl_matrix.h>
|
---|
21 | #include <gsl/gsl_permutation.h>
|
---|
22 | #include <gsl/gsl_vector.h>
|
---|
23 |
|
---|
24 | /**
|
---|
25 | * !@file
|
---|
26 | * These files defines several common operation on vectors that should not
|
---|
27 | * become part of the main vector class, because they are either to complex
|
---|
28 | * or need methods from other subsystems that should not be moved to
|
---|
29 | * the LinAlg-Subsystem
|
---|
30 | */
|
---|
31 |
|
---|
32 | /** Creates a new vector as the one with least square distance to a given set of \a vectors.
|
---|
33 | * \param *vectors set of vectors
|
---|
34 | * \param num number of vectors
|
---|
35 | * \return true if success, false if failed due to linear dependency
|
---|
36 | */
|
---|
37 | bool LSQdistance(Vector &res,const Vector **vectors, int num)
|
---|
38 | {
|
---|
39 | int j;
|
---|
40 |
|
---|
41 | for (j=0;j<num;j++) {
|
---|
42 | Log() << Verbose(1) << j << "th atom's vector: " << vectors[j] << endl;
|
---|
43 | }
|
---|
44 |
|
---|
45 | int np = 3;
|
---|
46 | struct LSQ_params par;
|
---|
47 |
|
---|
48 | const gsl_multimin_fminimizer_type *T =
|
---|
49 | gsl_multimin_fminimizer_nmsimplex;
|
---|
50 | gsl_multimin_fminimizer *s = NULL;
|
---|
51 | gsl_vector *ss, *y;
|
---|
52 | gsl_multimin_function minex_func;
|
---|
53 |
|
---|
54 | size_t iter = 0, i;
|
---|
55 | int status;
|
---|
56 | double size;
|
---|
57 |
|
---|
58 | /* Initial vertex size vector */
|
---|
59 | ss = gsl_vector_alloc (np);
|
---|
60 | y = gsl_vector_alloc (np);
|
---|
61 |
|
---|
62 | /* Set all step sizes to 1 */
|
---|
63 | gsl_vector_set_all (ss, 1.0);
|
---|
64 |
|
---|
65 | /* Starting point */
|
---|
66 | par.vectors = vectors;
|
---|
67 | par.num = num;
|
---|
68 |
|
---|
69 | for (i=NDIM;i--;)
|
---|
70 | gsl_vector_set(y, i, (vectors[0]->at(i) - vectors[1]->at(i))/2.);
|
---|
71 |
|
---|
72 | /* Initialize method and iterate */
|
---|
73 | minex_func.f = &LSQ;
|
---|
74 | minex_func.n = np;
|
---|
75 | minex_func.params = (void *)∥
|
---|
76 |
|
---|
77 | s = gsl_multimin_fminimizer_alloc (T, np);
|
---|
78 | gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
|
---|
79 |
|
---|
80 | do
|
---|
81 | {
|
---|
82 | iter++;
|
---|
83 | status = gsl_multimin_fminimizer_iterate(s);
|
---|
84 |
|
---|
85 | if (status)
|
---|
86 | break;
|
---|
87 |
|
---|
88 | size = gsl_multimin_fminimizer_size (s);
|
---|
89 | status = gsl_multimin_test_size (size, 1e-2);
|
---|
90 |
|
---|
91 | if (status == GSL_SUCCESS)
|
---|
92 | {
|
---|
93 | printf ("converged to minimum at\n");
|
---|
94 | }
|
---|
95 |
|
---|
96 | printf ("%5d ", (int)iter);
|
---|
97 | for (i = 0; i < (size_t)np; i++)
|
---|
98 | {
|
---|
99 | printf ("%10.3e ", gsl_vector_get (s->x, i));
|
---|
100 | }
|
---|
101 | printf ("f() = %7.3f size = %.3f\n", s->fval, size);
|
---|
102 | }
|
---|
103 | while (status == GSL_CONTINUE && iter < 100);
|
---|
104 |
|
---|
105 | for (i=(size_t)np;i--;)
|
---|
106 | res[i] = gsl_vector_get(s->x, i);
|
---|
107 | gsl_vector_free(y);
|
---|
108 | gsl_vector_free(ss);
|
---|
109 | gsl_multimin_fminimizer_free (s);
|
---|
110 |
|
---|
111 | return true;
|
---|
112 | };
|
---|
113 |
|
---|
114 | /** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
|
---|
115 | * \param *axis rotation axis
|
---|
116 | * \param alpha rotation angle in radian
|
---|
117 | */
|
---|
118 | Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha)
|
---|
119 | {
|
---|
120 | Vector a,y;
|
---|
121 | Vector res;
|
---|
122 | // normalise this vector with respect to axis
|
---|
123 | a = vec;
|
---|
124 | a.ProjectOntoPlane(axis);
|
---|
125 | // construct normal vector
|
---|
126 | try {
|
---|
127 | y = Plane(axis,a,0).getNormal();
|
---|
128 | }
|
---|
129 | catch (MathException &excp) {
|
---|
130 | // The normal vector cannot be created if there is linar dependency.
|
---|
131 | // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
|
---|
132 | return vec;
|
---|
133 | }
|
---|
134 | y.Scale(vec.Norm());
|
---|
135 | // scale normal vector by sine and this vector by cosine
|
---|
136 | y.Scale(sin(alpha));
|
---|
137 | a.Scale(cos(alpha));
|
---|
138 | res = vec.Projection(axis);
|
---|
139 | // add scaled normal vector onto this vector
|
---|
140 | res += y;
|
---|
141 | // add part in axis direction
|
---|
142 | res += a;
|
---|
143 | return res;
|
---|
144 | };
|
---|