Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80226569
pod_engine.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 30, 13:50
Size
5 KB
Mime Type
text/x-python
Expires
Sun, Sep 1, 13:50 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20326806
Attached To
R6746 RationalROMPy
pod_engine.py
View Options
# Copyright (C) 2018-2020 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import
numpy
as
np
from
copy
import
deepcopy
as
copy
from
warnings
import
catch_warnings
from
rrompy.utilities.base.types
import
Np1D
,
Np2D
,
Tuple
,
HFEng
,
sampList
from
rrompy.sampling
import
sampleList
__all__
=
[
'PODEngine'
]
class
PODEngine
:
"""
POD engine for general matrix orthogonalization.
"""
def
__init__
(
self
,
HFEngine
:
HFEng
):
self
.
HFEngine
=
HFEngine
def
name
(
self
)
->
str
:
return
self
.
__class__
.
__name__
def
__str__
(
self
)
->
str
:
return
self
.
name
()
def
__repr__
(
self
)
->
str
:
return
self
.
__str__
()
+
" at "
+
hex
(
id
(
self
))
def
normalize
(
self
,
A
:
Np2D
,
is_state
:
bool
=
True
)
->
Tuple
[
Np1D
,
Np1D
]:
"""
Normalize column-wise by norm.
Args:
A: matrix to be normalized;
is_state: whether state-norm should be used.
Returns:
Resulting normalized matrix, column-wise norm.
"""
r
=
self
.
HFEngine
.
norm
(
A
,
is_state
=
is_state
)
return
A
/
r
,
r
def
GS
(
self
,
a
:
Np1D
,
Q
:
sampList
,
n
:
int
=
-
1
,
is_state
:
bool
=
True
)
->
Tuple
[
Np1D
,
Np1D
,
bool
]:
"""
Compute 1 Gram-Schmidt step with given projector.
Args:
a: vector to be projected;
Q: orthogonal projection matrix;
n: number of columns of Q to be considered;
is_state: whether state-norm should be used.
Returns:
Resulting normalized vector, coefficients of a wrt the updated
basis, whether computation is ill-conditioned.
"""
if
n
==
-
1
:
n
=
Q
.
shape
[
1
]
r
=
np
.
zeros
((
n
+
1
,),
dtype
=
Q
.
dtype
)
if
n
>
0
:
from
rrompy.utilities.numerical
import
dot
Q
=
Q
[:
n
]
for
j
in
range
(
2
):
# twice is enough!
nu
=
self
.
HFEngine
.
innerProduct
(
a
,
Q
,
is_state
=
is_state
)
a
=
a
-
dot
(
Q
,
nu
)
r
[:
-
1
]
=
r
[:
-
1
]
+
nu
.
flatten
()
r
[
-
1
]
=
self
.
HFEngine
.
norm
(
a
,
is_state
=
is_state
)
ill_cond
=
False
with
catch_warnings
(
record
=
True
)
as
w
:
snr
=
np
.
abs
(
r
[
-
1
])
/
np
.
linalg
.
norm
(
r
)
if
len
(
w
)
>
0
or
snr
<
np
.
finfo
(
np
.
complex
)
.
eps
*
len
(
r
):
ill_cond
=
True
r
[
-
1
]
=
1.
a
=
a
/
r
[
-
1
]
return
a
,
r
,
ill_cond
def
generalizedQR
(
self
,
A
:
sampList
,
Q0
:
sampList
=
None
,
only_R
:
bool
=
False
,
genTrials
:
int
=
10
,
is_state
:
bool
=
True
)
->
Tuple
[
sampList
,
Np2D
]:
"""
Compute generalized QR decomposition of a matrix through Householder
method; see Trefethen, IMA J.N.A., 2010.
Args:
A: matrix to be decomposed;
Q0(optional): initial orthogonal guess for Q; defaults to random;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
genTrials(optional): number of trials of generation of linearly
independent vector; defaults to 10.
is_state(optional): whether state-norm should be used; defaults to
True.
Returns:
Resulting (orthogonal and )upper-triangular factor(s).
"""
Nh
,
N
=
A
.
shape
B
=
copy
(
A
)
V
=
sampleList
(
np
.
zeros
(
A
.
shape
,
dtype
=
A
.
dtype
))
R
=
np
.
zeros
((
N
,
N
),
dtype
=
A
.
dtype
)
Q
=
copy
(
V
)
if
Q0
is
None
else
sampleList
(
Q0
)
for
k
in
range
(
N
):
a
=
B
[
k
]
R
[
k
,
k
]
=
self
.
HFEngine
.
norm
(
a
,
is_state
=
is_state
)
if
Q0
is
None
and
k
<
Nh
:
for
_
in
range
(
genTrials
):
Q
[
k
],
_
,
illC
=
self
.
GS
(
np
.
random
.
randn
(
Nh
),
Q
,
k
,
is_state
)
if
not
illC
:
break
else
:
illC
=
k
>=
Nh
if
illC
:
if
Q0
is
not
None
or
k
<
Nh
:
Q
[
k
]
=
0.
else
:
alpha
=
self
.
HFEngine
.
innerProduct
(
a
,
Q
[
k
],
is_state
=
is_state
)
if
np
.
isclose
(
np
.
abs
(
alpha
),
0.
):
s
=
1.
else
:
s
=
-
alpha
/
np
.
abs
(
alpha
)
Q
[
k
]
=
s
*
Q
[
k
]
V
[
k
],
_
,
_
=
self
.
GS
(
R
[
k
,
k
]
*
Q
[
k
]
-
a
,
Q
,
k
,
is_state
)
J
=
np
.
arange
(
k
+
1
,
N
)
vtB
=
self
.
HFEngine
.
innerProduct
(
B
[
J
],
V
[
k
],
is_state
=
is_state
)
B
[
J
]
=
(
B
[
J
]
-
2
*
np
.
outer
(
V
[
k
],
vtB
))
.
T
if
not
illC
:
R
[
k
,
J
]
=
self
.
HFEngine
.
innerProduct
(
B
[
J
],
Q
[
k
],
is_state
=
is_state
)
B
[
J
]
=
(
B
[
J
]
-
np
.
outer
(
Q
[
k
],
R
[
k
,
J
]))
.
T
if
only_R
:
return
R
for
k
in
range
(
min
(
Nh
,
N
)
-
1
,
-
1
,
-
1
):
J
=
np
.
arange
(
k
,
min
(
Nh
,
N
))
vtQ
=
self
.
HFEngine
.
innerProduct
(
Q
[
J
],
V
[
k
],
is_state
=
is_state
)
Q
[
J
]
=
(
Q
[
J
]
-
2
*
np
.
outer
(
V
[
k
],
vtQ
))
.
T
return
Q
,
R
Event Timeline
Log In to Comment