-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
veenstrajelmer
committed
Jul 22, 2022
1 parent
dc87d49
commit a744c98
Showing
3 changed files
with
148 additions
and
96 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -73,7 +73,7 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
|
||
from hatyan.hatyan_core import get_const_list_hatyan, sort_const_list, robust_timedelta_sec, robust_daterange_fromtimesextfreq | ||
from hatyan.hatyan_core import get_freqv0_generic, get_uf_generic | ||
from hatyan.timeseries import check_ts, check_rayleigh, Timeseries_Statistics | ||
from hatyan.timeseries import check_ts, nyquist_folding, check_rayleigh, Timeseries_Statistics | ||
|
||
|
||
class HatyanSettings: | ||
|
@@ -303,7 +303,7 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
elif len(kwargs)>0: | ||
raise Exception('both arguments hatyan_settings and other settings (e.g. nodalfactors) are provided, this is not valid') | ||
|
||
print('-'*50) | ||
#print('-'*50) | ||
print('ANALYSIS initializing') | ||
print(hatyan_settings) | ||
|
||
|
@@ -331,6 +331,10 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
const_list_counts = pd.DataFrame({'constituent':const_list_uniq,'occurences':const_list_uniq_counts}) | ||
raise Exception('remove duplicate constituents from const_list:\n%s'%(const_list_counts.loc[const_list_counts['occurences']>1])) | ||
|
||
#check for length | ||
if len(ts_pd)<2: | ||
raise Exception('provided timeseries is less than 2 timesteps long, analysis not possible') | ||
|
||
#remove nans | ||
ts_pd_nonan = ts_pd[~ts_pd['values'].isna()] | ||
if len(ts_pd_nonan)==0: | ||
|
@@ -356,24 +360,9 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
u_i_rad, f_i = get_uf_generic(hatyan_settings, const_list, dood_date_fu) | ||
v_u = v_0i_rad.values + u_i_rad.values | ||
|
||
check_rayleigh(ts_pd,t_const_freq_pd) | ||
|
||
#TODO: nyquist stuk code generieker maken met Henrique | ||
stats = Timeseries_Statistics(ts=ts_pd) | ||
unique_timesteps = stats.stats['timeseries unique timesteps (minutes)'] | ||
if len(unique_timesteps)==1: #constant freq, so folding is possible | ||
constant_freq_phr = list(unique_timesteps)[0]/60 | ||
fs = 1/constant_freq_phr | ||
nyquist_freq = 0.5*fs | ||
bool_isnyquist = t_const_freq_pd['freq']==nyquist_freq | ||
if bool_isnyquist.any(): | ||
raise Exception(f'there is a component on the Nyquist frequency ({nyquist_freq} [1/hr]), this not possible:\n{t_const_freq_pd.loc[bool_isnyquist]}') | ||
bool_freqtoohigh = t_const_freq_pd['freq']>nyquist_freq | ||
t_const_freq_pd_folded = t_const_freq_pd.copy() | ||
t_const_freq_pd_folded.loc[bool_freqtoohigh,'freq'] = np.abs(fs - t_const_freq_pd_folded.loc[bool_freqtoohigh,'freq'] ) #remainder is better: 0.408%(1/3) | ||
t_const_freq_pd_folded = t_const_freq_pd_folded.sort_values('freq') #TODO: sorting also done in rayleigh check | ||
print('check folded rayleigh') | ||
check_rayleigh(ts_pd,t_const_freq_pd_folded) | ||
#check rayleigh frequency after nyquist frequency folding process | ||
freq_rem = nyquist_folding(ts_pd,t_const_freq_pd) | ||
check_rayleigh(ts_pd,freq_rem) #TODO: maybe sometimes valuable to not fold with nyquist (eg with strongly varying time interval), in that case: check_rayleigh(ts_pd,t_const_freq_pd) | ||
|
||
#### TIMESERIES ANALYSIS | ||
N = len(const_list) | ||
|
@@ -389,18 +378,6 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
tic = dt.datetime.now() | ||
xTxmat = np.dot(xTmat,xmat) | ||
|
||
if 0: | ||
#TODO | ||
xTxmat2 = xTxmat/xTxmat[0,0] | ||
xTxmat2[np.abs(xTxmat2)<0.05] = 0 | ||
xTxmat2_condition = np.linalg.cond(xTxmat2) | ||
svd_u,svd_s,svd_vh = np.linalg.svd(xTxmat2) # s zijn singuliere waarden van groot naar klein, u en v zijn gespiegeld over diagonaal | ||
xTxmat2_check = [email protected](svd_s)@svd_vh #matrix multiplicatie geeft weer originele matrix (ongeveer) | ||
xTxmat2_check = svd_u[:,:2]@np.diag(svd_s[:2])@svd_vh[:2,:] #helft van info weggooien, geeft approximatie van originele matrix | ||
|
||
#samp freq is 3hr, freq [1/hr],dus fs is 1/3, nyquist freq is 1/6 (0.5*fs) | ||
#probleem is freq van 3M2S10 (0.408201 [1/hr]) is groter dan 1/6 | ||
|
||
print('xTx matrix calculated') | ||
if 'A0' in const_list: #correct center value [N,N] for better matrix condition | ||
xTxmat_condition = np.linalg.cond(xTxmat) | ||
|
@@ -529,7 +506,7 @@ <h1 class="title">Module <code>hatyan.analysis_prediction</code></h1> | |
elif len(kwargs)>0: | ||
raise Exception('both arguments hatyan_settings and other settings (e.g. nodalfactors) are provided, this is not valid') | ||
|
||
print('-'*50) | ||
#print('-'*50) | ||
print('PREDICTION initializing') | ||
print(hatyan_settings) | ||
|
||
|
@@ -857,7 +834,7 @@ <h2 id="returns">Returns</h2> | |
elif len(kwargs)>0: | ||
raise Exception('both arguments hatyan_settings and other settings (e.g. nodalfactors) are provided, this is not valid') | ||
|
||
print('-'*50) | ||
#print('-'*50) | ||
print('ANALYSIS initializing') | ||
print(hatyan_settings) | ||
|
||
|
@@ -885,6 +862,10 @@ <h2 id="returns">Returns</h2> | |
const_list_counts = pd.DataFrame({'constituent':const_list_uniq,'occurences':const_list_uniq_counts}) | ||
raise Exception('remove duplicate constituents from const_list:\n%s'%(const_list_counts.loc[const_list_counts['occurences']>1])) | ||
|
||
#check for length | ||
if len(ts_pd)<2: | ||
raise Exception('provided timeseries is less than 2 timesteps long, analysis not possible') | ||
|
||
#remove nans | ||
ts_pd_nonan = ts_pd[~ts_pd['values'].isna()] | ||
if len(ts_pd_nonan)==0: | ||
|
@@ -910,24 +891,9 @@ <h2 id="returns">Returns</h2> | |
u_i_rad, f_i = get_uf_generic(hatyan_settings, const_list, dood_date_fu) | ||
v_u = v_0i_rad.values + u_i_rad.values | ||
|
||
check_rayleigh(ts_pd,t_const_freq_pd) | ||
|
||
#TODO: nyquist stuk code generieker maken met Henrique | ||
stats = Timeseries_Statistics(ts=ts_pd) | ||
unique_timesteps = stats.stats['timeseries unique timesteps (minutes)'] | ||
if len(unique_timesteps)==1: #constant freq, so folding is possible | ||
constant_freq_phr = list(unique_timesteps)[0]/60 | ||
fs = 1/constant_freq_phr | ||
nyquist_freq = 0.5*fs | ||
bool_isnyquist = t_const_freq_pd['freq']==nyquist_freq | ||
if bool_isnyquist.any(): | ||
raise Exception(f'there is a component on the Nyquist frequency ({nyquist_freq} [1/hr]), this not possible:\n{t_const_freq_pd.loc[bool_isnyquist]}') | ||
bool_freqtoohigh = t_const_freq_pd['freq']>nyquist_freq | ||
t_const_freq_pd_folded = t_const_freq_pd.copy() | ||
t_const_freq_pd_folded.loc[bool_freqtoohigh,'freq'] = np.abs(fs - t_const_freq_pd_folded.loc[bool_freqtoohigh,'freq'] ) #remainder is better: 0.408%(1/3) | ||
t_const_freq_pd_folded = t_const_freq_pd_folded.sort_values('freq') #TODO: sorting also done in rayleigh check | ||
print('check folded rayleigh') | ||
check_rayleigh(ts_pd,t_const_freq_pd_folded) | ||
#check rayleigh frequency after nyquist frequency folding process | ||
freq_rem = nyquist_folding(ts_pd,t_const_freq_pd) | ||
check_rayleigh(ts_pd,freq_rem) #TODO: maybe sometimes valuable to not fold with nyquist (eg with strongly varying time interval), in that case: check_rayleigh(ts_pd,t_const_freq_pd) | ||
|
||
#### TIMESERIES ANALYSIS | ||
N = len(const_list) | ||
|
@@ -943,18 +909,6 @@ <h2 id="returns">Returns</h2> | |
tic = dt.datetime.now() | ||
xTxmat = np.dot(xTmat,xmat) | ||
|
||
if 0: | ||
#TODO | ||
xTxmat2 = xTxmat/xTxmat[0,0] | ||
xTxmat2[np.abs(xTxmat2)<0.05] = 0 | ||
xTxmat2_condition = np.linalg.cond(xTxmat2) | ||
svd_u,svd_s,svd_vh = np.linalg.svd(xTxmat2) # s zijn singuliere waarden van groot naar klein, u en v zijn gespiegeld over diagonaal | ||
xTxmat2_check = [email protected](svd_s)@svd_vh #matrix multiplicatie geeft weer originele matrix (ongeveer) | ||
xTxmat2_check = svd_u[:,:2]@np.diag(svd_s[:2])@svd_vh[:2,:] #helft van info weggooien, geeft approximatie van originele matrix | ||
|
||
#samp freq is 3hr, freq [1/hr],dus fs is 1/3, nyquist freq is 1/6 (0.5*fs) | ||
#probleem is freq van 3M2S10 (0.408201 [1/hr]) is groter dan 1/6 | ||
|
||
print('xTx matrix calculated') | ||
if 'A0' in const_list: #correct center value [N,N] for better matrix condition | ||
xTxmat_condition = np.linalg.cond(xTxmat) | ||
|
@@ -1126,7 +1080,7 @@ <h2 id="returns">Returns</h2> | |
elif len(kwargs)>0: | ||
raise Exception('both arguments hatyan_settings and other settings (e.g. nodalfactors) are provided, this is not valid') | ||
|
||
print('-'*50) | ||
#print('-'*50) | ||
print('PREDICTION initializing') | ||
print(hatyan_settings) | ||
|
||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -128,7 +128,7 @@ <h2 id="getting-started">Getting Started</h2> | |
|
||
__author__ = """Jelmer Veenstra""" | ||
__email__ = '[email protected]' | ||
__version__ = '2.5.62' | ||
__version__ = '2.5.64' | ||
|
||
from hatyan.analysis_prediction import * | ||
from hatyan.astrog import * | ||
|
Oops, something went wrong.