problem using linear trend in stwave example

[mfarthin]

@jonghyunharrylee when I try to specify a linear trend using the 'drift':'linear' flag in example_inv_stwave.py I get the following traceback

***** Iteration 1 ******
computed Jacobian-Matrix products in 189.261032 secs
solve saddle point (co-kriging) systems with Levenberg-Marquardt
Traceback (most recent call last):
  File "example_inv_stwave.py", line 90, in <module>
    s_hat, simul_obs, post_diagv, iter_best = prob.Run()
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 1254, in Run
    s_hat, simul_obs, post_diagv, iter_best = self.GaussNewton()
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 1157, in GaussNewton
    s_cur, beta_cur, simul_obs_cur, obj = self.LinearIteration(s_past, simul_obs)
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 1071, in LinearIteration
    s_hat, beta, simul_obs_new = self.IterativeSolve(s_cur, simul_obs, precond = precond)
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 902, in IterativeSolve
    x, info = gmres(Afun, b, tol=itertol, restart=restart, maxiter=solver_maxiter, callback=callback, M=P)
  File "<decorator-gen-5>", line 2, in gmres
  File "/Users/rdchlmwf/anaconda2/envs/enkf/lib/python2.7/site-packages/scipy/_lib/_threadsafety.py", line 46, in caller
    return func(*a, **kw)
  File "/Users/rdchlmwf/anaconda2/envs/enkf/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/iterative.py", line 460, in gmres
    work[slice1] = psolve(work[slice2])
  File "/Users/rdchlmwf/anaconda2/envs/enkf/lib/python2.7/site-packages/scipy/sparse/linalg/interface.py", line 219, in matvec
    y = self._matvec(x)
  File "/Users/rdchlmwf/anaconda2/envs/enkf/lib/python2.7/site-packages/scipy/sparse/linalg/interface.py", line 471, in _matvec
    return self.__matvec_impl(x)
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 892, in Pmv
    S = np.dot(HX.T, invPsi(HX)) # p by p matrix
  File "/Users/rdchlmwf/Public/code/pyPCGA/pyPCGA/pcga.py", line 878, in invPsi
    return np.multiply(np.multiply((1./alpha[i]),self.invR),v) - np.dot(Psi_U_i, np.multiply(Dvec[:Psi_U_i.shape[1]].reshape(Psi_UTv.shape), Psi_UTv))
ValueError: cannot reshape array of size 69 into shape (69,3)

I get a similar complaint execept the size is (69,2) when I try to manually set X. Do I need to set another flag? Thanks

[mfarthin]

It looks like this is the command that is failing Dvec[:Psi_U_i.shape[1]].reshape(Psi_UTv.shape) . during the call S = np.dot(HX.T, invPsi(HX)) Is the goal to replicate Dvec[:Psi_U_i.shape[1]] into an array of p columns?

[mfarthin]

@jonghyunharrylee I pushed up a branch mwf/frf_uas_example where I tried a quick fix. The code runs at least and the linear trend seems to help

[jonghyunharrylee]

@mfarthin Thanks for reporting the error! The latest update to support heterogeneous data sets (by allowing measurement error matrix R as a vector) caused this error. I fixed this based on your implementation and now it is working for stwave duck site example.

[jonghyunharrylee]

@mfarthin I also checked frf_uas_20160722.py in mwf/frf_uas_example worked well. Note that when testing linear/quadratic trend (drift) prior, usually I would set prior variance small compared to a case with no prior.