yuditsabet commented 4 years ago

Hello, jajcayn I see that you have done a job similar to that of Martín Vejmelka. I am trying to adapt some of your scripts for my work with causal links but I encounter this error which I have not been able to solve for two days. I do not know if you had any similar error during the development of your work. sgf.prepare_surrogates(pool) File "/2020/Era5New/ndw-climate-master/scripts/surr_geo_field_ar.py", line 97, in prepare_surrogates job_results = map_func(_prepare_surrogates, job_data) File "/usr/lib/python2.7/multiprocessing/pool.py", line 253, in map return self.map_async(func, iterable, chunksize).get() File "/usr/lib/python2.7/multiprocessing/pool.py", line 572, in get raise self._value TypeError: tuple indices must be integers, not tuple

Creo que el error está asociado a estas líneas de código:

run the job in parallel

    job_data = [(i, j, self.order_range, self.crit, self.d[:, i, j]) for i in range(num_lats) for j in range(num_lons)]             
    job_results = map_func(_prepare_surrogates, job_data)

Regards

yuditsabet commented 4 years ago

The problem is here:

def estimate(self, time_series, prng, fit_intercept = False, crit_type = 'sbc', ndx = None): """Stepwise selection of the model using QR decomposition according to Neumaier and Schneider across model orders [prng[0], prng[1]]. The criterion 'crit' is the model order selection criterion. Optionally the (integer) index array indicates which parts of the time series are contiguous and which are not. Code is mostly a port from the MATLAB ARFIT toolbox by the same authors."""

    p_min, p_max = prng[0], prng[1]
    fi = 1 if fit_intercept else 0

    ts = time_series[:, np.newaxis] if time_series.ndim == 1 else time_series

    N, m = ts.shape
    ****print N,m** in my case N=156 and m=1**
    n_p = np.zeros(shape=(p_max+1,), dtype = np.int)
    n_p[p_max] = fi + p_max * m

    # remove "presample" data (p_max values from start) 
    N = N - p_max

    # construct matrix K (add row space for regularization matrix deltaD)
    K = np.zeros((N + n_p[p_max] + m, n_p[p_max] + m))

    # add intercept if required
    if(fit_intercept):
        K[:N, 0] = 1.0

    # set predictors u
    for j in range(1, p_max + 1):
        K[:N, fi+(j-1)*m:fi+j*m] = ts[p_max-j:N+p_max-j, :]
        **print (ts[p_max-j:N+p_max-j, :])   all 0 that is a problem??**

jajcayn commented 4 years ago

Dear @yuditsabet ,

unfortunately, I do not work in this field anymore and furthermore, this code is unmaintained and 5 years old. I can give you only couple of pointers:

VAR model was written indeed by Martin Vejmelka, not by me I just used it in this project
the whole codebase is python2.7 which is discontinued and not advisable to use anyway, I strongly advise you to port the parts of the code you need to python3
when I was using AR surrogates, Martin's port of Matlab's ARFIT toolbox was the only way how to do AR fits in python. However, a lot changed since then and e.g. statsmodels offers perfect functionality for you: they have ARMA models, and also AR models: for model selection (i.e. estimating model order p) you can use ar_select_order and then for fit their AutoReg class
all 0s btw seem to be a problem - it is like the model cannot be fitted - if you haven't done it before, I'd strongly advise doing some basic exploratory analysis on your timeseries, maybe they cannot be fitted using AR(p) model hence you cannot use AR surrogates

You can easily rewrite parts of SurrogateField class using statsmodels - I am pretty sure it will work.

Best, N.

yuditsabet commented 4 years ago

Thank you so much. I wrote Martin but does not reply to me. Sorry I am reviewing your propose Regards

jajcayn / multi-scale

job_results = map_func(_prepare_surrogates, job_data) #1

run the job in parallel