Solar Wind Models

[Hazboun6]

This is a new version of the solar wind model, that I am planning on writing up a methods paper on. This pull request includes a number of changes. The first is actually a change to model_utils.py.

1. I have included a new jump proposal that will propose jumps from the prior of any "non-standard" signals. I had wanted to included a flag where one could supply the name of a signal, trying to avoid the proliferation of jump proposals, but PTMCMCSampler only takes a standard set of object inputs. If someone can think of a cleaner, or more specific way to do this please make suggestions. @paulthebaker or @svigeland Any thoughts on this?

    def draw_from_signal_prior(self, x, iter, beta):

        q = x.copy()
        lqxy = 0
        std = ['linear timing model',
               'red noise',
               'phys_ephem',
               'gw',
               'cw',
               'bwm',
               'gp_sw',
               'ecorr_sherman-morrison',
               'ecorr',
               'efac',
               'equad',
               ]
        non_std = [nm for nm in self.snames.keys() if nm not in std]
        # draw parameter from signal model
        signal_name = np.random.choice(non_std)
        param = np.random.choice(self.snames[signal_name])
        if param.size:
            idx2 = np.random.randint(0, param.size)
            q[self.pmap[str(param)]][idx2] = param.sample()[idx2]

        # scalar parameter
        else:
            q[self.pmap[str(param)]] = param.sample()

        # forward-backward jump probability
        lqxy = (param.get_logpdf(x[self.pmap[str(param)]]) -
                param.get_logpdf(q[self.pmap[str(param)]]))

        return q, float(lqxy)

[Hazboun6]

I have removed all of the references to the solar wind code, except in the /electromagnetic/solar_wind.py file. This file includes:

1. A deterministic signal for the solar wind that has a common mean solar electron density for all of the pulsars. One can choose one value for the entire data set, or partition into bins.

#Here there are 11 bins
n_earth = parameter.Uniform(0,20,size=11)('n_earth')
sw = SW.solar_wind(n_earth=n_earth,n_earth_bins=11,t_init=tmin,t_final=tmax)
mean_sw = deterministic_signals.Deterministic(sw, name='mean_sw')

2. I have also included an astrophysical prior for the electron density:
   Here is a single value example using the ACE prior:

n_earth = SW.ACE_SWEPAM_Parameter()('n_earth')
sw = SW.solar_wind(n_earth=n_earth)
mean_sw = deterministic_signals.Deterministic(sw, name='mean_sw')

3. There is also a function to create a Fourier design matrix for constructing a Gaussian process. This can be used as desired, but works best as a perturbation to the mean signal defined above.

dm_sw_basis = SW.createfourierdesignmatrix_solar_dm(nmodes=15,Tspan=Tspan)
dm_sw_prior = utils.powerlaw(log10_A=log10_A_dm_sw, gamma=gamma_dm_sw)
gp_sw = gp_signals.BasisGP(priorFunction=dm_sw_prior, basisFunction=dm_sw_basis, name='gp_sw')

4. There is also a convenience function for making the "best" version of the SW model. The user can choose to include a simple power-law GP for the DM as well. This is again, tuned to the best version from a lot of testing that I have done. The issue is that frequencies higher than 1/yr in the DM GP are strongly covariant with the cusps of the solar wind. If one includes these frequencies the two models give and take giving funky results, like negative electron densities. There are options to provide various priors and bases, but here is a simple example:

dm_block = SW.solar_wind_block(Tspan=Tspan, include_dmgp=True)

[paulthebaker]

Re: point 1, see issue #19