Skip to content

Commit 3dcf587

Browse files
lambdacdmlambdacdm
authored
Add files via upload
1 parent 74dd329 commit 3dcf587

File tree

2 files changed

+31
-25
lines changed

2 files changed

+31
-25
lines changed

‎numerical.h‎

Lines changed: 13 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,9 @@ class Matrix
5252

5353
Matrix<DM> &operator=(const Matrix<DM> &);
5454
DM &operator()(int,int);
55+
DM operator()(int, int) const;
5556
DM &operator()(int);
57+
DM operator()(int) const;
5658
const vector<DM> &operator[](int);
5759
template <class DB> friend bool operator==(const Matrix<DB> &, const Matrix<DB> &);
5860
template<class DB> friend ostream &operator<< (ostream &,const Matrix<DB> &);
@@ -330,18 +332,22 @@ template<class DB> DB& Matrix<DB>::operator()(int r,int c)
330332
}
331333
return value[r][c];
332334
}
335+
template<class DB> DB Matrix<DB>::operator()(int r,int c) const
336+
{
337+
if(r>=row_num || c>=column_num)
338+
{
339+
cerr << "错误：矩阵下标越界" << '\n';
340+
return value[0][0];
341+
}
342+
return value[r][c];
343+
}
333344
template<class DB> DB& Matrix<DB>::operator()(int r)
334345
{
335346
return (*this)(r,0);
336347
}
337-
template<class DB> const vector<DB>& Matrix<DB>::operator[](int r)
348+
template<class DB> DB Matrix<DB>::operator()(int r) const
338349
{
339-
if(r>=row_num)
340-
{
341-
cerr << "错误：矩阵下标越界" << '\n';
342-
return value[0];
343-
}
344-
return value[r];
350+
return (*this)(r,0);
345351
}
346352
template <class DB> bool operator==(const Matrix<DB> &A, const Matrix<DB> &B)
347353
{

‎optimization.h‎

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -70,8 +70,8 @@ template<class OP> Matrix<OP> Gradient(const function<OP(Matrix<OP>)> &f, const
7070
{
7171
Matrix<OP> x1 = x;
7272
Matrix<OP> x2 = x;
73-
x1(i, 0) = Get(x, i, 0) - h;
74-
x2(i, 0) = Get(x, i, 0) + h;
73+
x1(i, 0) = x(i, 0) - h;
74+
x2(i, 0) = x(i, 0) + h;
7575
df(i, 0) = (f(x2) - f(x1)) / (2*h);
7676
}
7777
return df;
@@ -86,8 +86,8 @@ template<class OP> Matrix<OP> GradientT(const function<OP(Matrix<OP>)> &f, const
8686
{
8787
Matrix<OP> x1 = x;
8888
Matrix<OP> x2 = x;
89-
x1(i, 0) = Get(x, i, 0) - h;
90-
x2(i, 0) = Get(x, i, 0) + h;
89+
x1(i, 0) = x(i, 0) - h;
90+
x2(i, 0) = x(i, 0) + h;
9191
df(0, i) = (f(x2) - f(x1)) / twoh;
9292
}
9393
return df;
@@ -102,25 +102,25 @@ template<class OP> Matrix<OP> Hessian(const function<OP(Matrix<OP>)> &f, const M
102102
{
103103
Matrix<OP> x1 = x;
104104
Matrix<OP> x2 = x;
105-
x1(i, 0) = Get(x, i, 0) + twoh;
106-
x2(i, 0) = Get(x, i, 0) - twoh;
105+
x1(i, 0) = x(i, 0) + twoh;
106+
x2(i, 0) = x(i, 0) - twoh;
107107
H(i, i) = ((f(x1) - f(x)) / twoh + (f(x2) - f(x)) / twoh) / twoh;
108108
for (int j = i + 1; j < n;++j)
109109
{
110110
Matrix<OP> x1 = x;
111111
Matrix<OP> x2 = x;
112112
Matrix<OP> x3 = x;
113113
Matrix<OP> x4 = x;
114-
x1(i, 0) = Get(x, i, 0) + h; x1(j, 0) = Get(x, j, 0) + h;
115-
x2(i, 0) = x1(i,0); x2(j, 0) = Get(x, j, 0) - h;
116-
x3(i, 0) = Get(x, i, 0) - h; x3(j, 0) = x1(j, 0);
114+
x1(i, 0) = x(i, 0) + h; x1(j, 0) = x(j, 0) + h;
115+
x2(i, 0) = x1(i,0); x2(j, 0) = x(j, 0) - h;
116+
x3(i, 0) = x(i, 0) - h; x3(j, 0) = x1(j, 0);
117117
x4(i, 0) = x3(i,0); x4(j, 0) = x2(j, 0);
118118
H(i, j) = ((f(x1) - f(x2)) / twoh + (f(x4) - f(x3)) / twoh) / twoh;
119119
}
120120
}
121121
for (int i = 1; i < n;++i)
122122
for (int j = 0; j < i;++j)
123-
H(i, j) = Get(H, j, i);
123+
H(i, j) = H(j, i);
124124
return H;
125125
}
126126

@@ -135,7 +135,7 @@ const Matrix<OP> &x, const Matrix<OP> &direction, OP stepsize, OP gamma, OP c)
135135
{
136136
const uint8_t times = 10;
137137
uint8_t cnt = 0;
138-
while(f(x + (stepsize * direction)) > f(x)+c*stepsize*Get(GradientT<OP>(f,x)*direction,0,0)
138+
while(f(x + (stepsize * direction)) > f(x)+c*stepsize*(GradientT<OP>(f,x)*direction)(0,0)
139139
&& cnt<times)
140140
{
141141
stepsize = gamma * stepsize;
@@ -154,7 +154,7 @@ const Matrix<OP> &x, const Matrix<OP> &direction, OP stepsize, OP c)
154154
const uint8_t times = 10;
155155
uint8_t cnt = 0;
156156
OP fx = f(x);
157-
OP fd = Get(GradientT<OP>(f, x) * direction, 0, 0);
157+
OP fd = (GradientT<OP>(f, x) * direction)(0, 0);
158158
OP lhs = f(x + stepsize * direction);
159159
OP afd = stepsize * fd;
160160
OP left_end = OP(0);
@@ -189,7 +189,7 @@ const Matrix<OP> &x, const Matrix<OP> &direction, OP stepsize, OP c1, OP c2, boo
189189
{
190190
const uint8_t times = 10;
191191
uint8_t cnt = 0;
192-
OP fd = Get(GradientT<OP>(f, x) * direction, 0, 0);
192+
OP fd = (GradientT<OP>(f, x) * direction)(0, 0);
193193
Matrix<OP> xad = x + stepsize * direction;
194194
OP left_end = OP(0);
195195
OP right_end = OP(4) * stepsize;
@@ -200,12 +200,12 @@ const Matrix<OP> &x, const Matrix<OP> &direction, OP stepsize, OP c1, OP c2, boo
200200
right_end = stepsize;
201201
stepsize = (left_end + right_end) / OP(2);
202202
}
203-
else if (!strong && Get(GradientT<OP>(f, xad)*direction,0,0) < c2 * fd)
203+
else if (!strong && Abs((GradientT<OP>(f, xad)*direction)(0,0)) < c2 * fd)
204204
{
205205
left_end = stepsize;
206206
stepsize = std::min(OP(2) * stepsize, (left_end + right_end) / OP(2));
207207
}
208-
else if (strong && Abs(Get(GradientT<OP>(f, xad)*direction,0,0)) > -c2 * fd)
208+
else if (strong && Abs((GradientT<OP>(f, xad)*direction)(0,0)) > -c2 * fd)
209209
{
210210
left_end = stepsize;
211211
stepsize = std::min(OP(2) * stepsize, (left_end + right_end) / OP(2));
@@ -526,7 +526,7 @@ const Matrix<OP> &x_pre, const Matrix<OP> &gradient, const Matrix<OP> &gradient_
526526
Matrix<OP> y = gradient - gradient_pre;
527527
Matrix<OP> M = s - H * y;
528528
Matrix<OP> MT = Transpose(M);
529-
return H + (OP(1) / Get(MT * y, 0, 0)) * (M * MT);
529+
return H + (OP(1) / (MT * y)(0, 0)) * (M * MT);
530530
}
531531
template<class OP> Matrix<OP> SR1(const function<OP(Matrix<OP>)> &f,
532532
const Matrix<OP> &initial_x, OP epsilon)
@@ -540,8 +540,8 @@ const Matrix<OP> &x_pre, const Matrix<OP> &gradient, const Matrix<OP> &gradient_
540540
Matrix<OP> y = gradient - gradient_pre;
541541
Matrix<OP> Hy = H * y;
542542
Matrix<OP> yT = Transpose(y);
543-
return H - (OP(1) / Get(yT * Hy, 0, 0)) * Hy * Transpose(Hy) +
544-
(OP(1) / Get(yT * s, 0, 0)) * s * Transpose(s);
543+
return H - (OP(1) / (yT * Hy)(0, 0)) * Hy * Transpose(Hy) +
544+
(OP(1) / (yT * s)(0, 0)) * s * Transpose(s);
545545
}
546546
template<class OP> Matrix<OP> DFP(const function<OP(Matrix<OP>)> &f,
547547
const Matrix<OP> &initial_x, OP epsilon)

0 commit comments

Comments
 (0)