Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Mpp
MLUQ
Commits
5d3c196e
Commit
5d3c196e
authored
Aug 26, 2019
by
niklas.baumgarten
Browse files
naming
parent
c116980f
Changes
1
Hide whitespace changes
Inline
Side-by-side
mlmc/src/CirculantEmbedding.C
View file @
5d3c196e
#include
"CirculantEmbedding.h"
vector
<
double
>
CirculantEmbedding1D
::
get
_f
ieldf
()
{
vector
<
double
>
CirculantEmbedding1D
::
get
F
ieldf
()
{
vector
<
double
>
fieldf
;
fieldf
.
resize
(
N
);
if
(
internal
_c
ounter
==
0
)
Xf
=
generate
_f
ield
();
if
(
internal
C
ounter
==
0
)
Xf
=
generate
F
ield
();
for
(
int
i
=
0
;
i
<
N
;
i
++
)
{
if
(
field
_t
ype
==
"Gaussian"
)
{
if
(
internal
_c
ounter
==
0
)
fieldf
[
i
]
=
Xf
[
i
].
real
()
+
mean
;
if
(
field
T
ype
==
"Gaussian"
)
{
if
(
internal
C
ounter
==
0
)
fieldf
[
i
]
=
Xf
[
i
].
real
()
+
mean
;
else
fieldf
[
i
]
=
Xf
[
i
].
imag
()
+
mean
;
}
else
if
(
field
_t
ype
==
"LogNormal"
)
{
if
(
internal
_c
ounter
==
0
)
fieldf
[
i
]
=
exp
(
Xf
[
i
].
real
()
+
mean
);
}
else
if
(
field
T
ype
==
"LogNormal"
)
{
if
(
internal
C
ounter
==
0
)
fieldf
[
i
]
=
exp
(
Xf
[
i
].
real
()
+
mean
);
else
fieldf
[
i
]
=
exp
(
Xf
[
i
].
imag
()
+
mean
);
}
else
Exit
(
"Has to be Gaussian or LogNormal"
)
}
internal
_c
ounter
+=
1
;
internal
_c
ounter
%=
2
;
internal
C
ounter
+=
1
;
internal
C
ounter
%=
2
;
return
fieldf
;
}
vector
<
complex
<
double
>>
CirculantEmbedding1D
::
generate
_f
ield
()
{
vector
<
complex
<
double
>>
CirculantEmbedding1D
::
generate
F
ield
()
{
vector
<
complex
<
double
>>
Xf_
;
vector
<
complex
<
double
>>
Z
;
Z
.
resize
(
N_embedded
);
...
...
@@ -34,7 +34,7 @@ vector<complex<double>> CirculantEmbedding1D::generate_field() {
vector
<
complex
<
double
>>
product
;
product
.
resize
(
N_embedded
);
for
(
int
i
=
0
;
i
<
N_embedded
;
i
++
)
{
product
[
i
]
=
sqrt
_ev
[
i
]
*
Z
[
i
];
product
[
i
]
=
sqrt
EV
[
i
]
*
Z
[
i
];
}
fft
(
product
,
Xf_
);
double
divisor
=
sqrt
(
N_embedded
);
...
...
@@ -43,7 +43,7 @@ vector<complex<double>> CirculantEmbedding1D::generate_field() {
return
Xf_
;
}
vector
<
double
>
CirculantEmbedding1D
::
generate
_i
nterpol
_f
ield
(
const
vector
<
double
>
&
fieldf
)
{
vector
<
double
>
CirculantEmbedding1D
::
generate
I
nterpol
F
ield
(
const
vector
<
double
>
&
fieldf
)
{
vector
<
double
>
fieldc
;
if
(
N
%
2
!=
0
)
Exit
(
"Interpolation of field does only work for evenly spaced grids."
)
int
N_interpol
=
N
/
2
;
...
...
@@ -54,7 +54,7 @@ vector<double> CirculantEmbedding1D::generate_interpol_field(const vector<double
return
fieldc
;
}
void
CirculantEmbedding1D
::
create
_c
ovariance
_m
atrix
(
vector
<
double
>
&
ar
,
vector
<
double
>
&
ac
)
{
void
CirculantEmbedding1D
::
create
C
ovariance
M
atrix
(
vector
<
double
>
&
ar
,
vector
<
double
>
&
ac
)
{
const
vector
<
double
>
&
coordinate1
=
linspace
(
domain_start
,
domain_end
,
N
);
ar
.
resize
(
N
);
ac
.
resize
(
N
);
...
...
@@ -74,13 +74,13 @@ void CirculantEmbedding1D::create_covariance_matrix(vector<double> &ar, vector<d
}
}
void
CirculantEmbedding1D
::
embed
_c
ovariance
_m
atrix
(
vector
<
double
>
&
ar
,
vector
<
double
>
&
ac
,
vector
<
double
>
&
br
)
{
void
CirculantEmbedding1D
::
embed
C
ovariance
M
atrix
(
vector
<
double
>
&
ar
,
vector
<
double
>
&
ac
,
vector
<
double
>
&
br
)
{
br
.
reserve
(
2
*
ar
.
size
()
-
1
);
br
.
insert
(
br
.
end
(),
ar
.
begin
(),
ar
.
end
());
br
.
insert
(
br
.
end
(),
ac
.
rbegin
(),
ac
.
rend
()
-
1
);
}
void
CirculantEmbedding1D
::
compute
_ev
(
vector
<
double
>
&
br
,
vector
<
double
>
&
ev
)
{
void
CirculantEmbedding1D
::
compute
EV
(
vector
<
double
>
&
br
,
vector
<
double
>
&
ev
)
{
fft
(
br
,
ev
);
}
...
...
@@ -112,16 +112,16 @@ void CirculantEmbedding1D::padding(vector<double> &ar, vector<double> &ac) {
}
}
vector
<
double
>
CirculantEmbedding1D
::
get
_s
qrt
_ev
()
{
vector
<
double
>
CirculantEmbedding1D
::
get
S
qrt
EV
()
{
vector
<
double
>
sqrt_ev_return
,
ar
,
ac
,
br
,
ev
;
create
_c
ovariance
_m
atrix
(
ar
,
ac
);
embed
_c
ovariance
_m
atrix
(
ar
,
ac
,
br
);
compute
_ev
(
br
,
ev
);
create
C
ovariance
M
atrix
(
ar
,
ac
);
embed
C
ovariance
M
atrix
(
ar
,
ac
,
br
);
compute
EV
(
br
,
ev
);
for
(
int
i
=
0
;
i
<
ev
.
size
();
i
++
)
{
if
(
ev
[
i
]
>
0
)
continue
;
else
if
(
ev
[
i
]
<
0
&&
ev
[
i
]
>
e
ps
)
ev
[
i
]
=
0
.
0
;
else
if
(
ev
[
i
]
<
0
&&
ev
[
i
]
>
e
vtol
)
ev
[
i
]
=
0
.
0
;
else
padding
(
ar
,
ac
);
}
...
...
@@ -131,29 +131,29 @@ vector<double> CirculantEmbedding1D::get_sqrt_ev() {
return
sqrt_ev_return
;
}
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
get
_f
ieldf
()
{
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
get
F
ieldf
()
{
vector
<
vector
<
double
>>
fieldf
;
fieldf
.
resize
(
N
[
1
]);
if
(
internal
_c
ounter
f
==
0
)
Xf
=
generate
_f
ield
();
if
(
internal
C
ounter
==
0
)
Xf
=
generate
F
ield
();
for
(
int
i
=
0
;
i
<
N
[
1
];
i
++
)
{
fieldf
[
i
].
resize
(
N
[
0
]);
for
(
int
j
=
0
;
j
<
N
[
0
];
j
++
)
{
if
(
field
_t
ype
==
"Gaussian"
)
{
if
(
internal
_c
ounter
f
==
0
)
fieldf
[
i
][
j
]
=
Xf
[
i
][
j
].
real
()
+
mean
;
if
(
field
T
ype
==
"Gaussian"
)
{
if
(
internal
C
ounter
==
0
)
fieldf
[
i
][
j
]
=
Xf
[
i
][
j
].
real
()
+
mean
;
else
fieldf
[
i
][
j
]
=
Xf
[
i
][
j
].
imag
()
+
mean
;
}
else
if
(
field
_t
ype
==
"LogNormal"
)
{
if
(
internal
_c
ounter
f
==
0
)
fieldf
[
i
][
j
]
=
exp
(
Xf
[
i
][
j
].
real
()
+
mean
);
}
else
if
(
field
T
ype
==
"LogNormal"
)
{
if
(
internal
C
ounter
==
0
)
fieldf
[
i
][
j
]
=
exp
(
Xf
[
i
][
j
].
real
()
+
mean
);
else
fieldf
[
i
][
j
]
=
exp
(
Xf
[
i
][
j
].
imag
()
+
mean
);
}
else
Exit
(
"Has to be Gaussian or LogNormal"
)
}
}
internal
_c
ounter
f
+=
1
;
internal
_c
ounter
f
%=
2
;
internal
C
ounter
+=
1
;
internal
C
ounter
%=
2
;
return
fieldf
;
}
vector
<
vector
<
complex
<
double
>>>
CirculantEmbedding2D
::
generate
_f
ield
()
{
vector
<
vector
<
complex
<
double
>>>
CirculantEmbedding2D
::
generate
F
ield
()
{
vector
<
vector
<
complex
<
double
>>>
Xf_
;
vector
<
vector
<
complex
<
double
>>>
Z
;
Z
.
resize
(
N_embedded
[
1
]);
...
...
@@ -174,7 +174,7 @@ vector<vector<complex<double>>> CirculantEmbedding2D::generate_field() {
for
(
int
i
=
0
;
i
<
N_embedded
[
1
];
i
++
)
{
product
[
i
].
resize
(
N_embedded
[
0
]);
for
(
int
j
=
0
;
j
<
N_embedded
[
0
];
j
++
)
{
product
[
i
][
j
]
=
sqrt
_ev
[
i
][
j
]
*
Z
[
i
][
j
];
product
[
i
][
j
]
=
sqrt
EV
[
i
][
j
]
*
Z
[
i
][
j
];
}
}
fft2
(
product
,
Xf_
);
...
...
@@ -185,7 +185,7 @@ vector<vector<complex<double>>> CirculantEmbedding2D::generate_field() {
return
Xf_
;
}
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
generate
_i
nterpol
_f
ield
(
vector
<
vector
<
double
>>
&
fieldf
)
{
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
generate
I
nterpol
F
ield
(
vector
<
vector
<
double
>>
&
fieldf
)
{
vector
<
vector
<
double
>>
fieldc
;
if
(
N
[
0
]
%
2
!=
0
||
N
[
1
]
%
2
!=
0
)
Exit
(
"Interpolation of field does only work for evenly spaced grids."
)
int
N_interpol
[
2
]
=
{
N
[
0
]
/
2
,
N
[
1
]
/
2
};
...
...
@@ -201,7 +201,7 @@ vector<vector<double>> CirculantEmbedding2D::generate_interpol_field(vector<vect
return
fieldc
;
}
void
CirculantEmbedding2D
::
create
_c
ovariance
_m
atrix
(
vector
<
vector
<
double
>>
&
ar
,
vector
<
vector
<
double
>>
&
ac
)
{
void
CirculantEmbedding2D
::
create
C
ovariance
M
atrix
(
vector
<
vector
<
double
>>
&
ar
,
vector
<
vector
<
double
>>
&
ac
)
{
const
vector
<
double
>
&
coordinate1
=
linspace
(
domain1_start
,
domain1_end
,
N
[
0
]);
const
vector
<
double
>
&
coordinate2
=
linspace
(
domain2_start
,
domain2_end
,
N
[
1
]);
...
...
@@ -230,8 +230,8 @@ void CirculantEmbedding2D::create_covariance_matrix(vector<vector<double>> &ar,
}
}
void
CirculantEmbedding2D
::
embed
_c
ovariance
_m
atrix
(
vector
<
vector
<
double
>>
&
ar
,
vector
<
vector
<
double
>>
&
ac
,
vector
<
vector
<
double
>>
&
br
)
{
void
CirculantEmbedding2D
::
embed
C
ovariance
M
atrix
(
vector
<
vector
<
double
>>
&
ar
,
vector
<
vector
<
double
>>
&
ac
,
vector
<
vector
<
double
>>
&
br
)
{
int
br_size
=
2
*
ar
.
size
()
-
1
;
br
.
resize
(
br_size
);
for
(
int
i
=
0
;
i
<
br
.
size
();
i
++
)
{
...
...
@@ -245,7 +245,7 @@ void CirculantEmbedding2D::embed_covariance_matrix(vector<vector<double>> &ar, v
}
}
void
CirculantEmbedding2D
::
compute
_ev
(
vector
<
vector
<
double
>>
&
br
,
vector
<
vector
<
double
>>
&
ev
)
{
void
CirculantEmbedding2D
::
compute
EV
(
vector
<
vector
<
double
>>
&
br
,
vector
<
vector
<
double
>>
&
ev
)
{
fft2
(
br
,
ev
);
}
...
...
@@ -308,17 +308,17 @@ void CirculantEmbedding2D::padding(vector<vector<double>> &ar, vector<vector<dou
}
}
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
get
_s
qrt
_ev
()
{
vector
<
vector
<
double
>>
CirculantEmbedding2D
::
get
S
qrt
EV
()
{
vector
<
vector
<
double
>>
sqrt_ev_return
,
ar
,
ac
,
br
,
ev
;
create
_c
ovariance
_m
atrix
(
ar
,
ac
);
embed
_c
ovariance
_m
atrix
(
ar
,
ac
,
br
);
compute
_ev
(
br
,
ev
);
create
C
ovariance
M
atrix
(
ar
,
ac
);
embed
C
ovariance
M
atrix
(
ar
,
ac
,
br
);
compute
EV
(
br
,
ev
);
for
(
int
i
=
0
;
i
<
ev
[
0
].
size
();
i
++
)
{
for
(
int
j
=
0
;
j
<
ev
[
1
].
size
();
j
++
)
{
if
(
ev
[
i
][
j
]
>
0
)
continue
;
else
if
(
ev
[
i
][
j
]
<
0
&&
ev
[
i
][
j
]
>
e
ps
)
ev
[
i
][
j
]
=
0
.
0
;
else
if
(
ev
[
i
][
j
]
<
0
&&
ev
[
i
][
j
]
>
e
vtol
)
ev
[
i
][
j
]
=
0
.
0
;
else
padding
(
ar
,
ac
);
}
}
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment