Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Removes use of statsmodels.tsa.arima_model.ARMA #1896

Merged
merged 9 commits into from
Aug 2, 2022
Merged

Removes use of statsmodels.tsa.arima_model.ARMA #1896

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a547424
Moves to ARIMA class in statsmodels
Jul 18, 2022
4d72d07
regold for new ARIMA class used in ARMA ROM
Jul 18, 2022
968008d
formatting and comments
Jul 18, 2022
0034272
Merge branch 'devel' into arima
Jul 22, 2022
5879bb4
handles storage of ARIMAResults
j-bryan Jul 25, 2022
1b155b5
corrects variable names
j-bryan Jul 27, 2022
579cca0
Merge branch 'arima' of https://github.com/j-bryan/raven into arima
Jul 27, 2022
9b2450d
fixes setting ARMA parameters for interpolated ROMs
Jul 27, 2022
3d9d72a
loosens Basic test tolerance
Aug 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 31 additions & 24 deletions ravenframework/SupervisedLearning/ARMA.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -915,15 +915,15 @@ def _generateARMASignal(self, model, numSamples=None,randEngine=None):
numSamples = len(self.pivotParameterValues)
if randEngine is None:
randEngine=self.randomEng
import statsmodels.tsa
hist = statsmodels.tsa.arima_process.arma_generate_sample(ar = np.append(1., -model.arparams),
ma = np.append(1., model.maparams),
nsample = numSamples,
distrvs = functools.partial(randomUtils.randomNormal,engine=randEngine),
# functool.partial provide the random number generator as a function
# with normal distribution and take engine as the positional arguments keywords.
scale = np.sqrt(model.sigma2),
burnin = 2*max(self.P,self.Q)) # @alfoa, 2020
import statsmodels.api
hist = statsmodels.tsa.arima_process.arma_generate_sample(ar=model.polyar,
ma=model.polyma,
nsample=numSamples,
distrvs=functools.partial(randomUtils.randomNormal,engine=randEngine),
# functool.partial provide the random number generator as a function
# with normal distribution and take engine as the positional arguments keywords.
scale=model.sigma,
burnin=2*max(self.P,self.Q)) # @alfoa, 2020
return hist

def _generateFourierSignal(self, pivots, periods):
Expand Down Expand Up @@ -1094,7 +1094,13 @@ def _trainARMA(self, data, masks=None):
if masks is not None:
data = data[masks]
import statsmodels.api
results = statsmodels.tsa.arima_model.ARMA(data, order = (self.P, self.Q)).fit(disp = False)
results = statsmodels.tsa.arima.model.ARIMA(data, order=(self.P, 0, self.Q), trend='c').fit()
# The ARIMAResults object here can cause problems with ray when running in parallel. Dropping it
# in exchange for the armaResultsProxy class avoids the issue while preserving what we really
# care out from the ARIMAResults object.
results = armaResultsProxy(results.polynomial_ar,
results.polynomial_ma,
np.sqrt(results.params[results.param_names.index('sigma2')]))
return results

def _trainCDF(self, data, binOps=None):
Expand Down Expand Up @@ -1382,7 +1388,7 @@ def writeXML(self, writeTo, targets=None, skip=None):
root.append(targetNode)
armaNode = xmlUtils.newNode('ARMA_params')
targetNode.append(armaNode)
armaNode.append(xmlUtils.newNode('std', text=np.sqrt(arma.sigma2)))
armaNode.append(xmlUtils.newNode('std', text=arma.sigma))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line is failing in HERON tests: https://civet.inl.gov/job/1125079/

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is because we don't have a way to update the trained ARMAs that we use for testing as prerequisites to running the HERON cases, so if ever the ARMA ROM gets updated, we have to manually update those serialized testing ARMAs. This may not be avoidable, but we should make sure that @j-bryan can update those trained ARMAs in HERON immediately after merging this. Will that work?

Copy link
Contributor

@joshua-cogliati-inl joshua-cogliati-inl Aug 2, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that is fine with me. (Is there written documentation on how to do this? I checked Troubleshooting and some other pages in the HERON wiki and didn't see how to update those ARMAs.)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It probably should be documented better. There's RAVEN XML files in the ARMA test folder along with the associated training data, so it should be a simple matter of running those training XML files with the new changes in RAVEN present.

# TODO covariances, P and Q, etc
for target,peakInfo in self.peaks.items():
targetNode = root.find(target)
Expand Down Expand Up @@ -1564,13 +1570,13 @@ def getFundamentalFeatures(self, requestedFeatures, featureTemplate=None):
for target, arma in self.armaResult.items():
# sigma
feature = featureTemplate.format(target=target, metric='arma', id='std')
features[feature] = np.sqrt(arma.sigma2)
features[feature] = arma.sigma
# autoregression
for p, val in enumerate(arma.arparams):
for p, val in enumerate(-arma.polyar[1:]): # The AR coefficients are stored in polynomial form here (flipped sign and with a term in the zero position of the array for lag=0)
feature = featureTemplate.format(target=target, metric='arma', id='AR_{}'.format(p))
features[feature] = val
# moving average
for q, val in enumerate(arma.maparams):
for q, val in enumerate(arma.polyma[1:]): # keep only the terms for lag>0
feature = featureTemplate.format(target=target, metric='arma', id='MA_{}'.format(q))
features[feature] = val
for target, cdfParam in self.cdfParams.items():
Expand Down Expand Up @@ -1864,15 +1870,15 @@ def _setArmaResults(self, paramDict):
if 'AR' in info:
AR_keys, AR_vals = zip(*list(info['AR'].items()))
AR_keys, AR_vals = zip(*sorted(zip(AR_keys, AR_vals), key=lambda x:x[0]))
AR_vals = np.asarray(AR_vals)
AR_vals = np.concatenate(([1], -np.asarray(AR_vals))) # convert the AR params into polynomial form
else:
AR_vals = np.array([])
AR_vals = np.array([1]) # must include a 1 for the zero lag term
if 'MA' in info:
MA_keys, MA_vals = zip(*list(info['MA'].items()))
MA_keys, MA_vals = zip(*sorted(zip(MA_keys, MA_vals), key=lambda x:x[0]))
MA_vals = np.asarray(MA_vals)
MA_vals = np.concatenate(([1], np.asarray(MA_vals))) # converts the MA params into polynomial form
else:
MA_vals = np.array([])
MA_vals = np.array([1])
if 'bin' in info:
bin_keys, bin_vals = zip(*list(info['bin'].items()))
bin_keys, bin_vals = zip(*sorted(zip(bin_keys, bin_vals), key=lambda x:x[0]))
Expand Down Expand Up @@ -2163,6 +2169,7 @@ def adjustLocalRomSegment(self, settings, picker):
Call this before training the subspace segment ROMs
Note this is called on the LOCAL subsegment ROMs, NOT on the GLOBAL templateROM from the ROMcollection!
@ In, settings, object, arbitrary information about ROM clustering settings from getGlobalRomSegmentSettings
@ In, picker, slice, slice object for selecting the desired segment
@ Out, None
"""
if self.zeroFilterTarget:
Expand Down Expand Up @@ -2537,14 +2544,14 @@ class armaResultsProxy:
Class that can be used to artifically construct ARMA information
from pre-determined values
"""
def __init__(self, arparams, maparams, sigma):
def __init__(self, polyar, polyma, sigma):
"""
Constructor.
@ In, arparams, np.array(float), autoregressive coefficients
@ In, maparams, np.array(float), moving average coefficients
@ In, polyar, np.array(float), autoregressive coefficients
@ In, polyma, np.array(float), moving average coefficients
@ In, sigma, float, standard deviation of ARMA residual noise
@ Out, None
"""
self.arparams = np.atleast_1d(arparams)
self.maparams = np.atleast_1d(maparams)
self.sigma2 = sigma**2
self.polyar = polyar
self.polyma = polyma
self.sigma = sigma
2 changes: 1 addition & 1 deletion tests/framework/ROM/TimeSeries/ARMA/gold/ARMAreseed/rommeta.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
</period>
</Fourier>
<ARMA_params>
<std>0.258086103993</std>
<std>0.258080826596</std>
</ARMA_params>
</Speed>
</ARMA>
Expand Down
2 changes: 1 addition & 1 deletion tests/framework/ROM/TimeSeries/ARMA/gold/Basic/romMeta.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@
</period>
</Fourier>
<ARMA_params>
<std>0.258086103993</std>
<std>0.258080826596</std>
</ARMA_params>
</Speed>
</arma>
Expand Down
35 changes: 18 additions & 17 deletions tests/framework/ROM/TimeSeries/ARMA/gold/Peaks/romMeta.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
<seg_index_start>segment_number</seg_index_start>
</dims>
<general>
<datasetName>romMeta</datasetName>
<pointwise_meta>seg_Time_end,seg_Time_start,seg_index_end,seg_index_start</pointwise_meta>
<sampleTag>RAVEN_sample_ID</sampleTag>
</general>
Expand All @@ -24,7 +25,7 @@
<SegmentROM segment="0" type="Static">
<Signal>
<ARMA_params>
<std>1.39905186671</std>
<std>1.39906272172</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -41,7 +42,7 @@
<SegmentROM segment="1" type="Static">
<Signal>
<ARMA_params>
<std>1.40116930574</std>
<std>1.40117407499</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -58,7 +59,7 @@
<SegmentROM segment="2" type="Static">
<Signal>
<ARMA_params>
<std>1.44350608151</std>
<std>1.44351127777</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -75,7 +76,7 @@
<SegmentROM segment="3" type="Static">
<Signal>
<ARMA_params>
<std>1.40741081917</std>
<std>1.4074139545</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -92,7 +93,7 @@
<SegmentROM segment="4" type="Static">
<Signal>
<ARMA_params>
<std>1.40295178957</std>
<std>1.402957792</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -109,7 +110,7 @@
<SegmentROM segment="5" type="Static">
<Signal>
<ARMA_params>
<std>1.37645138378</std>
<std>1.37645633643</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -126,7 +127,7 @@
<SegmentROM segment="6" type="Static">
<Signal>
<ARMA_params>
<std>1.4551579758</std>
<std>1.45516109199</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -143,7 +144,7 @@
<SegmentROM segment="7" type="Static">
<Signal>
<ARMA_params>
<std>1.40077421912</std>
<std>1.40077918937</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -160,7 +161,7 @@
<SegmentROM segment="8" type="Static">
<Signal>
<ARMA_params>
<std>1.39729542433</std>
<std>1.3972982007</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -177,7 +178,7 @@
<SegmentROM segment="9" type="Static">
<Signal>
<ARMA_params>
<std>1.39409898092</std>
<std>1.3941043431</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -194,7 +195,7 @@
<SegmentROM segment="10" type="Static">
<Signal>
<ARMA_params>
<std>1.40451155501</std>
<std>1.40451697137</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -211,7 +212,7 @@
<SegmentROM segment="11" type="Static">
<Signal>
<ARMA_params>
<std>1.41542838444</std>
<std>1.41542097263</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -228,7 +229,7 @@
<SegmentROM segment="12" type="Static">
<Signal>
<ARMA_params>
<std>1.39571008147</std>
<std>1.39571506343</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -245,7 +246,7 @@
<SegmentROM segment="13" type="Static">
<Signal>
<ARMA_params>
<std>1.42764320201</std>
<std>1.42763966138</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -262,7 +263,7 @@
<SegmentROM segment="14" type="Static">
<Signal>
<ARMA_params>
<std>1.39176698371</std>
<std>1.3917729782</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -279,7 +280,7 @@
<SegmentROM segment="15" type="Static">
<Signal>
<ARMA_params>
<std>1.38512558105</std>
<std>1.38513462051</std>
</ARMA_params>
<Peak_params>
<peak>
Expand All @@ -296,7 +297,7 @@
<SegmentROM segment="16" type="Static">
<Signal>
<ARMA_params>
<std>1.38110648914</std>
<std>1.38110661531</std>
</ARMA_params>
<Peak_params>
<peak>
Expand Down
2 changes: 1 addition & 1 deletion tests/framework/ROM/TimeSeries/ARMA/tests
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
input = 'basic.xml'
output = 'Basic/romMeta.csv'
UnorderedXml = 'Basic/romMeta.xml'
rel_err = 1e-6
rel_err = 7e-6
[../]

[./ARMAparallel]
Expand Down
4 changes: 2 additions & 2 deletions tests/framework/unit_tests/SupervisedLearning/testARMA.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -467,8 +467,8 @@ def plotPDF(edges, bins, ax, label, color, s='.', alpha=1.0):
arma.setEngine(eng, seed=901017)
signal7 = arma._generateARMASignal(testVal)

sig7 = [0.39975177, -0.14531468, 0.13138866, -0.56565224, 0.06020252,
0.60752306, -0.29076173, -1.1758456, 0.41108591, -0.05735384]
sig7=[0.39974990, -0.14531400, 0.1313880, -0.565649598, 0.06020223,
0.60752023, -0.29076037, -1.1758401, 0.41108399, -0.05735358]
for n in range(10):
checkFloat('signal 7, evaluation ind{}'.format(n), signal7[n], sig7[n], tol=1e-7)

Expand Down