← Back to branch summary

~max-mcginley/third-order-lie-map/generator

« back to all changes in this revision

Viewing changes to tests/LieOperator_unittest.cc

  • Committer: Max McGinley
  • Date: 2015-09-10 17:13:14 UTC
  • Revision ID: max.mcginley@stfc.ac.uk-20150910171314-dsy6iqs2u3qhm3u4
More tests added

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/LieOperator_unittest.cc
13
13
 
14
14
namespace MAUS {
15
15
 
 
16
 
 
17
TEST(LieOperatorTest, PoissonBracket) {
 
18
        Matrix<double> poissonBase(4,4,0.0);
 
19
        for (int i = 0; i < 4; i++) {
 
20
                for (int j = 0; j < 4; j++) {
 
21
                        lie_polynomial u1;
 
22
                        addCoefficient(std::vector<int>(1,i),1.0,u1);
 
23
                        lie_polynomial u2;
 
24
                        addCoefficient(std::vector<int>(1,j),1.0,u2);
 
25
                        lie_polynomial p = LieOperator::poissonBracket(u1,u2);
 
26
                        poissonBase(i+1,j+1) = p[0L];
 
27
                }
 
28
        }
 
29
        EXPECT_EQ(poissonBase, LieMap::get_symplectic_matrix(4));
 
30
        
 
31
        //Testing PB of {p_x - x, sin(x)} = -cos(x)
 
32
        lie_polynomial sinp;
 
33
        lie_polynomial sinq;
 
34
        lie_polynomial cosp;
 
35
        lie_polynomial cosq;
 
36
        lie_polynomial expq;
 
37
        lie_polynomial expmp;
 
38
        lie_polynomial sin2q;
 
39
        int fact = 1;
 
40
        for (int i = 0; i < 6; i++) {
 
41
                addCoefficient(std::vector<int>(2*i,0),std::pow(-1,i)/fact,cosq);
 
42
                addCoefficient(std::vector<int>(2*i,2),std::pow(-1,i)/fact,cosp);
 
43
                addCoefficient(std::vector<int>(2*i,0),1./fact,expq);
 
44
                addCoefficient(std::vector<int>(2*i,2),1./fact,expmp);
 
45
                fact *= 2*i + 1;
 
46
                addCoefficient(std::vector<int>(2*i + 1,0),std::pow(-1,i)/fact,sinq);
 
47
                addCoefficient(std::vector<int>(2*i + 1,2),std::pow(-1,i)/fact,sinp);
 
48
                addCoefficient(std::vector<int>(2*i + 1,0),2.*std::pow(-4,i)/fact,sin2q);
 
49
                addCoefficient(std::vector<int>(2*i+1,0),1./fact,expq);
 
50
                addCoefficient(std::vector<int>(2*i+1,2),-1./fact,expmp);
 
51
                fact *= 2*i + 2;
 
52
        }
 
53
        
 
54
        lie_polynomial sinsqq = LieOperator::multiplyPolynomials(sinq,sinq);
 
55
        lie_polynomial pbCos = LieOperator::poissonBracket(LieOperator::multiplyPolynomials(expq,cosp),LieOperator::multiplyPolynomials(expmp,sinsqq));
 
56
        lie_polynomial expect = LieOperator::multiplyPolynomials(expq,LieOperator::multiplyPolynomials(sinp,LieOperator::multiplyPolynomials(sin2q,expmp)));
 
57
        addToPolynomial(LieOperator::multiplyPolynomials(expq,LieOperator::multiplyPolynomials(cosp,LieOperator::multiplyPolynomials(expmp,sinsqq))),expect,-1.);
 
58
        lie_polynomial expectTrunc;
 
59
        for (lie_polynomial::const_iterator it = expect.cbegin(); it != expect.cend(); it++) {
 
60
                if (getIndexFromLongCode(it->first).size() <= 11) addCoefficient(it->first, it->second,expectTrunc);
 
61
        }
 
62
        lie_polynomial pbTrunc;
 
63
        for (lie_polynomial::const_iterator it = pbCos.cbegin(); it != pbCos.cend(); it++) {
 
64
                if (getIndexFromLongCode(it->first).size() <= 11) addCoefficient(it->first, it->second,pbTrunc);
 
65
        }
 
66
        for (int i = 0; i < 12; i++) {
 
67
                for (int j = 0; j < 12 - i; j++) {
 
68
                        std::vector<int> index(i,0);
 
69
                        index.insert(index.end(),j,2);
 
70
                        EXPECT_NEAR(getCoefficient(index,expectTrunc).second,getCoefficient(index,pbTrunc).second,1.e-15) << "Poisson bracket up to index " << i+j << " has failed\nexpected:\n" << expectTrunc << "\nactual:\n" << pbTrunc << "\n";
 
71
                }
 
72
        }
 
73
}
 
74
 
 
75
TEST(LieOperatorTest, MultiplyPolynomial) {
 
76
        lie_polynomial sinp;
 
77
        lie_polynomial cosp;
 
78
        lie_polynomial sin2p;
 
79
        int fact = 1;
 
80
        for (int i = 0; i < 8; i++) {
 
81
                addCoefficient(std::vector<int>(2*i,2),std::pow(-1,i)/fact,cosp);
 
82
                fact *= 2*i + 1;
 
83
                addCoefficient(std::vector<int>(2*i + 1,2),std::pow(-1,i)/fact,sinp);
 
84
                addCoefficient(std::vector<int>(2*i + 1,2),std::pow(-4,i)/(fact),sin2p);
 
85
                fact *= 2*i + 2;
 
86
        }
 
87
        lie_polynomial sincos = LieOperator::multiplyPolynomials(sinp,cosp);
 
88
        for (int i = 0; i < 13; i++) {
 
89
                for (int j = 0; j < 13 - i; j++) {
 
90
                        std::vector<int> index(i,2);
 
91
                        index.insert(index.end(),j,2);
 
92
                        EXPECT_NEAR(getCoefficient(index,sin2p).second,getCoefficient(index,sincos).second,1.e-16) << "Multiply polynomial up to index " << i+j << " has failed";
 
93
                }
 
94
        }
 
95
}
 
96
 
16
97
TEST(LieOperatorTest, Symplecticity) {
17
98
        lie_polynomial pol;
18
99
        for (int j = 0; j < 4; j++) {
30
111
        LieOperator l(pol);
31
112
        Matrix<double> S(4,4,0.0);
32
113
        for (int i = 1; i <= 4; i++) {
33
 
                for (int j = 1; j < i; j++) {
 
114
                for (int j = 1; j <= i; j++) {
34
115
                        double f = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*2.0;
35
116
                        S(i,j) = f;
36
117
                        S(j,i) = f;
37
118
                }
38
119
                S(i,i) = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*2.0;
39
120
        }
40
 
        Matrix<double> B = exponentMatrix(Matrix<double>(4,4,LieMap::get_symplectic_matrix(4).data()) * S, 20,1.e-10);
 
121
        Matrix<double> B = exponentMatrix(LieMap::get_symplectic_matrix(4) * S, 20,1.e-10);
41
122
        
42
123
        
43
124
        for (int i = 0; i < 3; i++) {
53
134
                
54
135
                free(v1); free(v2);
55
136
                
56
 
                ASSERT_NEAR(prodIn, prodOut,(10.*varScale*varScale*varScale*varScale)) << "Randomly generates symplectic matrix exp(JS) was not symplectic";
 
137
                ASSERT_NEAR(prodIn, prodOut,(10.*varScale*varScale*varScale*varScale)) << "Randomly generated symplectic matrix exp(JS) was not symplectic";
57
138
        }
58
139
        
59
140
        PolynomialMap map = l.applyTransformToMatrix(B);
60
 
        Matrix<double> sym(4,4,LieMap::get_symplectic_matrix(4).data());
 
141
        Matrix<double> sym(LieMap::get_symplectic_matrix(4));
61
142
        Matrix<double> id(4,4,0.0);
62
143
        for (int i =1; i<=4; i++) id(i,i) = 1.0;
63
144
        
66
147
                
67
148
                for (int i = 0; i < 4; i++) coordsIn[i] = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*varScale*2.;
68
149
                
69
 
                Matrix<double> thirdJacobian(4,4,0.0);
70
 
                for (int i = 0; i < 4; i++) {
71
 
                        for (int j = 0; j < 4; j++) {
72
 
                                lie_polynomial mapPol = getLieFromPolynomialMap(map,i,3);
73
 
                                lie_polynomial der = LieOperator::differentiate(mapPol, j);
74
 
                                PolynomialMap derivative = getPolynomialMapFromLie(der,0);
75
 
                                double* out = new double[1];
76
 
                                derivative.F(coordsIn,out);
77
 
                                thirdJacobian(i+1,j+1) = out[0];
78
 
                                delete[] out;
79
 
                        }
80
 
                }
81
 
 
82
 
                Matrix<double> jacobian = B + thirdJacobian;
 
150
                Matrix<double> jacobian = LieOperator::getJacobianMatrix(map,coordsIn,4,4);
83
151
                
84
152
                Matrix<double> mjm = transpose(jacobian) * sym * jacobian;
85
 
                EXPECT_TRUE(AreMatricesNear(sym,mjm,0.01*std::pow(varScale,2))) << "Tangent matrix for third " << thirdJacobian << "\nPLUS first order\n" << B << "\nis not symplectic...\nM^T J M = " << mjm;
 
153
                EXPECT_TRUE(AreMatricesNear(sym,mjm,std::pow(20.,4)*std::pow(varScale,9))) << "B =\n" << B << "\nJacobian matrix\n" << jacobian << "\nis not symplectic. mjm =\n" << mjm << "\n";
86
154
        
87
155
                delete[] coordsIn;
88
156
        }
89
157
}
90
158
 
91
 
static const std::vector<std::vector<int> > indices = {{0,0,0,0}, {0,1,2,3}, {1,3,1,3}, {2,3,0,0}, {1,1,1,3}, {1,1,2}, {1,1,3,0}};
92
 
static const std::vector<double> coeffs = {0.325123461, 1.21464257, -2.56326413, 0.5215215, -0.6754357, 0.1848755, -0.9365346};
93
 
static const std::vector<int> vars = {0,2,1,0,1,1,2};
94
 
static const std::vector<int> pows = {4,1,2,2,3,2,0};
95
 
static const std::vector<std::vector<int> > indices_diff = {{0,0,0}, {0,1,3}, {1,3,3}, {3,2,0}, {1,1,3}, {1,2}, {-1,-1,-1}};
 
159
TEST(LieOperatorTest, SimpleHarmonicMotion) {
 
160
        
 
161
}
 
162
 
 
163
//If lie polynomial is independent of x, p_x should be conserved.
 
164
TEST(LieOperatorTest, Conservation) {
 
165
        Matrix<double> S(4,4,0.0);
 
166
        for (int i = 2; i <= 4; i++) {
 
167
                for (int j = 2; j <= i; j++) {
 
168
                        double f = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*2.0;
 
169
                        S(i,j) = f;
 
170
                        S(j,i) = f;
 
171
                }
 
172
                S(i,i) = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*2.0;
 
173
        }
 
174
        
 
175
        lie_polynomial pol;
 
176
        for (int j = 1; j < 4; j++) {
 
177
                for (int k = 1; k <= j; k++) {
 
178
                        for (int m = 1; m <= k; m++) {
 
179
                                for (int n = 1; n <= m; n++) {
 
180
                                        std::vector<int> perm(4,0); perm[0] = j; perm[1] = k; perm[2] = m; perm[3] = n;
 
181
                                        double co = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*20.0;
 
182
                                        addCoefficient(makeCoefficient(perm,co),pol);
 
183
                                }
 
184
                        }
 
185
                }
 
186
        }
 
187
        LieOperator l(pol);
 
188
        
 
189
        Matrix<double> B = exponentMatrix(LieMap::get_symplectic_matrix(4) * S, 20,1.e-10);
 
190
 
 
191
        for (int i = 0; i < 3; i++) {
 
192
                double* v1 = BASEparticle::generateRandomCoordinates();
 
193
                double* v2 = BASEparticle::generateRandomCoordinates();
 
194
                
 
195
                Vector<double> vec1 = Vector<double>(v1,4);
 
196
                Vector<double> vec2 = Vector<double>(v2,4);
 
197
                
 
198
                double prodIn = LieOperator::symplecticInnerProduct(v1,v2);
 
199
 
 
200
                double prodOut = LieOperator::symplecticInnerProduct(B*vec1,B*vec2);
 
201
                
 
202
                free(v1); free(v2);
 
203
                
 
204
                ASSERT_NEAR(prodIn, prodOut,(10.*varScale*varScale*varScale*varScale)) << "Randomly generated symplectic matrix exp(JS) was not symplectic";
 
205
        }
 
206
        
 
207
        PolynomialMap map = l.applyTransformToMatrix(B);
 
208
        
 
209
        for (int i = 0; i < 10; i++) {
 
210
                double* coordsIn = new double[4];
 
211
                double* coordsOut = new double[4];
 
212
                for (int i = 0; i < 4; i++) coordsIn[i] = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*varScale*2.;
 
213
                map.F(coordsIn,coordsOut);
 
214
                EXPECT_DOUBLE_EQ(coordsIn[2],coordsOut[2]) << "Hamiltonian independent of x does not preserve p_x";
 
215
        }
 
216
}
 
217
 
 
218
static const std::vector<std::vector<int> > indices = {{0,0,0,0}, {0,1,2,3}, {1,3,1,3}, {2,3,0,0}, {1,1,1,3}, {1,1,2}, {1,1,3,0}, {0,0,0,0,0,0,0,1,1,2,3}};
 
219
static const std::vector<double> coeffs = {0.325123461, 1.21464257, -2.56326413, 0.5215215, -0.6754357, 0.1848755, -0.9365346,1.376238759};
 
220
static const std::vector<int> vars = {0,2,1,0,1,1,2,0};
 
221
static const std::vector<int> pows = {4,1,2,2,3,2,0,7};
 
222
static const std::vector<std::vector<int> > indices_diff = {{0,0,0}, {0,1,3}, {1,3,3}, {3,2,0}, {1,1,3}, {1,2}, {-1,-1,-1},{0,0,0,0,0,0,1,1,2,3}};
96
223
 
97
224
TEST(LieOperatorTest, Differentiate) {
98
225
        lie_polynomial pol;
108
235
        }
109
236
}
110
237
 
111
 
TEST(LieOperatorTest, ApplyMatrix) {
112
 
        MAUS::lie_polynomial pol;
113
 
        addCoefficient(indices[1],coeffs[1],pol); //should be a {0,1,2,3} mononomial
114
 
        LieOperator op(pol);
115
 
        Matrix<double> mat(4,4,0.0);
116
 
        for (int i = 0; i < 4; i++) {
117
 
                mat(1,i+1) = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*2.0;
 
238
TEST(LieOperatorTest,SextupoleKick) {
 
239
        Matrix<double> id(4,4,0.0);
 
240
        for (int i = 0; i < 4; i++) id(i+1,i+1) = 1.0;
 
241
        lie_polynomial pol;
 
242
        double s = (static_cast <double> (rand()) / (static_cast <float> (RAND_MAX)) - 0.5)*varScale*200.;
 
243
        addCoefficient(std::vector<int>(3,0),s,pol);
 
244
        LieOperator l(pol);
 
245
        PolynomialMap m = l.applyTransformToMatrix(id);
 
246
        std::vector<PolynomialMap::PolynomialCoefficient> coeffs = m.GetCoefficientsAsVector();
 
247
        long expCode1 = getLongCodeFromIndex(std::vector<int>(2,0));
 
248
        long expCode2 = getLongCodeFromIndex(std::vector<int>(1,0));
 
249
        long expCode3 = getLongCodeFromIndex(std::vector<int>(1,1));
 
250
        long expCode4 = getLongCodeFromIndex(std::vector<int>(1,2));
 
251
        long expCode5 = getLongCodeFromIndex(std::vector<int>(1,3));
 
252
        for (std::vector<PolynomialMap::PolynomialCoefficient>::const_iterator it = coeffs.cbegin(); it != coeffs.cend(); it++) {
 
253
                long code = getLongCodeFromIndex(it->InVariables());
 
254
                if (code == expCode1 && it->OutVariable() == 2) EXPECT_DOUBLE_EQ(3.*s,it->Coefficient());
 
255
                else if (code == expCode2 && it->OutVariable() == 0) EXPECT_DOUBLE_EQ(1.,it->Coefficient());
 
256
                else if (code == expCode3 && it->OutVariable() == 1) EXPECT_DOUBLE_EQ(1.,it->Coefficient());
 
257
                else if (code == expCode4 && it->OutVariable() == 2) EXPECT_DOUBLE_EQ(1.,it->Coefficient());
 
258
                else if (code == expCode5 && it->OutVariable() == 3) EXPECT_DOUBLE_EQ(1.,it->Coefficient());
 
259
                else EXPECT_DOUBLE_EQ(0.0,it->Coefficient());
118
260
        }
119
 
        PolynomialMap map = op.applyTransformToMatrix(mat);
120
 
        Matrix<double> mapMat = map.GetCoefficientsAsMatrix();
121
 
        EXPECT_EQ(mapMat(1,21),-coeffs[1]*mat(1,2)) << "Unexpected {0,1,2} value in " << map;
122
 
        EXPECT_EQ(mapMat(1,22),-coeffs[1]*mat(1,1)) << "Unexpected {0,1,3} value in " << map;
123
 
        EXPECT_EQ(mapMat(1,24),coeffs[1]*mat(1,4)) << "Unexpected {0,2,3} value in " << map;
124
 
        EXPECT_EQ(mapMat(1,30),coeffs[1]*mat(1,3)) << "Unexpected {1,2,3} value in " << map;
125
261
}
126
262
 
127
263
} //namespace MAUS