def sde(self):
<EOL>if self.approx_order is not None:<EOL><INDENT>N = int(self.approx_order)<EOL><DEDENT>else:<EOL><INDENT>N = <NUM_LIT:7> <EOL><DEDENT>p_period = float(self.period)        <EOL>p_lengthscale = <NUM_LIT:2>*float(self.lengthscale)<EOL>p_variance = float(self.variance)        <EOL>w0 = <NUM_LIT:2>*np.pi/p_period <EOL>[q2,dq2l] = seriescoeff(N, p_lengthscale, p_variance)        <EOL>dq2l = <NUM_LIT:2>*dq2l  <EOL>eps = <NUM_LIT><EOL>if np.any( np.isfinite(q2) == False) or np.any( np.abs(q2) > <NUM_LIT:1.0>/eps) or np.any( np.abs(q2) < eps):<EOL><INDENT>warnings.warn("<STR_LIT>".format(eps) + q2.__format__("<STR_LIT>") )<EOL><DEDENT>if np.any( np.isfinite(dq2l) == False) or np.any( np.abs(dq2l) > <NUM_LIT:1.0>/eps) or np.any( np.abs(dq2l) < eps):<EOL><INDENT>warnings.warn("<STR_LIT>".format(eps) + q2.__format__("<STR_LIT>") )<EOL><DEDENT>F    = np.kron(np.diag(range(<NUM_LIT:0>,N+<NUM_LIT:1>)),np.array( ((<NUM_LIT:0>, -w0), (w0, <NUM_LIT:0>)) ) )<EOL>L    = np.eye(<NUM_LIT:2>*(N+<NUM_LIT:1>))<EOL>Qc   = np.zeros((<NUM_LIT:2>*(N+<NUM_LIT:1>), <NUM_LIT:2>*(N+<NUM_LIT:1>)))<EOL>P_inf = np.kron(np.diag(q2),np.eye(<NUM_LIT:2>))<EOL>H    = np.kron(np.ones((<NUM_LIT:1>,N+<NUM_LIT:1>)),np.array((<NUM_LIT:1>,<NUM_LIT:0>)) )<EOL>P0 = P_inf.copy()<EOL>dF = np.empty((F.shape[<NUM_LIT:0>], F.shape[<NUM_LIT:1>], <NUM_LIT:3>))<EOL>dQc = np.empty((Qc.shape[<NUM_LIT:0>], Qc.shape[<NUM_LIT:1>], <NUM_LIT:3>))<EOL>dP_inf = np.empty((P_inf.shape[<NUM_LIT:0>], P_inf.shape[<NUM_LIT:1>], <NUM_LIT:3>))         <EOL>dF[:,:,<NUM_LIT:0>] = np.zeros(F.shape)<EOL>dQc[:,:,<NUM_LIT:0>] = np.zeros(Qc.shape)<EOL>dP_inf[:,:,<NUM_LIT:0>] = P_inf / p_variance<EOL>dF[:,:,<NUM_LIT:1>] = np.kron(np.diag(range(<NUM_LIT:0>,N+<NUM_LIT:1>)),np.array( ((<NUM_LIT:0>,  w0), (-w0, <NUM_LIT:0>)) ) / p_period );<EOL>dQc[:,:,<NUM_LIT:1>] = np.zeros(Qc.shape)<EOL>dP_inf[:,:,<NUM_LIT:1>] = np.zeros(P_inf.shape)      <EOL>dF[:,:,<NUM_LIT:2>] = np.zeros(F.shape)<EOL>dQc[:,:,<NUM_LIT:2>] = np.zeros(Qc.shape)<EOL>dP_inf[:,:,<NUM_LIT:2>] = np.kron(np.diag(dq2l),np.eye(<NUM_LIT:2>))<EOL>dP0 = dP_inf.copy()<EOL>if self.balance:<EOL><INDENT>import GPy.models.state_space_main as ssm<EOL>(F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf,dP0) = ssm.balance_ss_model(F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0 )<EOL><DEDENT>return (F, L, Qc, H, P_inf, P0, dF, dQc, dP_inf, dP0)<EOL>
Return the state space representation of the standard periodic covariance.


! Note: one must constrain lengthscale not to drop below 0.2. (independently of approximation order)
After this Bessel functions of the first becomes NaN. Rescaling
time variable might help.

! Note: one must keep period also not very low. Because then
the gradients wrt wavelength become ustable. 
However this might depend on the data. For test example with
300 data points the low limit is 0.15.