In [486]:

sc = 1e6
for m in range(M):
    plot(P[:,m]/sc,altvec,'k-',linewidth=0.2,alpha=0.4)
xlabel('$n_e [cm^{-3} \cdot 10^{-6}]$')
ylabel('Altitude [km]')
title('500 Random Nighttime IRI profiles')

Out[486]:

<matplotlib.text.Text at 0x7f127facadd0>

In [144]:

Out[144]:

array([  3.50003717e+12,   3.04689693e+11,   4.31592528e+10,
         7.89405079e+09,   1.80227116e+09,   8.76517173e+08,
         2.14227350e+08,   9.99631976e+07,   3.43183064e+07,
         1.42262266e+07,   9.69662774e+06,   4.91740582e+06,
         2.94734872e+06,   1.83720305e+06,   1.52815088e+06,
         1.00301485e+06,   5.96784664e+05,   4.72746655e+05,
         4.07735601e+05,   2.80916201e+05,   2.03274607e+05,
         1.36518344e+05,   9.29780345e+04,   7.56312707e+04,
         5.56431559e+04,   2.47717542e+04,   1.73758168e+04,
         1.33858788e+04,   9.18796520e+03,   6.42127385e+03,
         4.80505347e+03,   4.13195143e+03,   2.06284388e+03,
         1.69968067e+03,   1.29322269e+03,   9.47274982e+02,
         5.88812346e+02,   5.64390077e+02,   4.89124809e+02,
         3.44086780e+02,   2.79918760e+02,   2.32660246e+02,
         1.44986673e+02,   2.05367864e+01,   3.34320428e+00,
         4.47809235e-01,   6.46653843e-02,   3.66110307e-03,
         9.26529921e-04,   7.21857762e-04,   6.50804925e-04,
         5.68078117e-04,   3.76073493e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   3.58489040e-04,
         3.58489040e-04,   3.58489040e-04,   2.40322805e-04,
         8.30298751e-05])

In [472]:

sigma = sqrt(diag(Sigx))
errorbar(mu/sc,altvec,xerr=sigma/sc,fmt='k-',alpha=0.5)
xlabel('$n_e [cm^{-3} \cdot 10^{-6}]$')
ylabel('Altitude [km]')
title('Mean and Cov of ne from IRI')

Out[472]:

<matplotlib.text.Text at 0x7f128168abd0>

Pick a single profile, do forward model, then invert.

This is the measurement, with noise:

$$ z = Ax $$

In [490]:

plot(ztrue,altvec,'r')
plot(z,altvec,'k')
xlabel('Measurement [arb]')
ylabel('Tangent Altitude [km]')

Out[490]:

<matplotlib.text.Text at 0x7f127ef8ba50>

Solution to Linear MMSE problem: $$ x_{LMMSE} = \Sigma_x A^T (A \Sigma_x A^T + \Sigma_z)^{-1} (z - A \mu_x) + \mu_x$$

Note similarity with Bayesian MAP estimator using Gaussian prior:

$$ x \sim \mathcal{N}(\mu_x, \Sigma_x)$$$$ z \sim \mathcal{N}(A x, \Sigma_z)$$

Also note that Tikhonov regularization is of an similar form, except $\Sigma_x$ is replaced with $ (\alpha^2 L^T L)^{-1}$. So the difference here is that:

We're using IRI to create $\Sigma_x$ instead of pre-supposing a form for it.
We don't have to choose $\alpha$.

In [492]:

sc = 1e6
plot(xtrue/sc,altvec,'k-',linewidth=2,label='Truth')
plot(x_lmmse/sc,altvec,'r-',label='Rec. IRI prior')
#plot(x_t/sc,altvec,'g-',label='Rec. Tikh 2nd order, Lcurve')
xlabel('$n_e [cm^{-3} \cdot 10^{-6}]$')
ylabel('Altitude [km]')
legend(loc='best',fancybox=True)

Out[492]:

<matplotlib.legend.Legend at 0x7f127ea94950>

In [504]:

alphastar = alphavec[istar]
x_t,_,_ = get_x_t(alphastar)

Here is the comparison with 2nd-order Tikhonov regularization using L-curve:

In [507]:

sc = 1e6
plot(xtrue/sc,altvec,'k-',linewidth=2,label='Truth')
plot(x_lmmse/sc,altvec,'r-',label='Rec. IRI prior')
plot(x_t/sc,altvec,'g-',label='Rec. Tikh 2nd order, Lcurve')
xlabel('$n_e [cm^{-3} \cdot 10^{-6}]$')
ylabel('Altitude [km]')
legend(loc='best',fancybox=True,prop={'size':8})

Out[507]:

<matplotlib.legend.Legend at 0x7f127ec01a90>

In [ ]:

figure(figsize=(3,3))
plot(Xvec,Yvec,'k-')
plot(Xvec[istar],Yvec[istar],'ro')
title('L-curve')

Out[ ]:

<matplotlib.text.Text at 0x7f127e4326d0>

In [500]:

idx = randint(M)
print idx
xtrue = P[:,idx] # pick a "random" column

ztrue = A.dot(xtrue)

Sign = diag(ztrue)
n = random.multivariate_normal(zeros(N), Sign)

z = ztrue + n


R = A.dot(Sigx.dot(A.T)) + Sign
x_lmmse = Sigx.dot(A.T).dot(np.linalg.solve(R, z-A.dot(mux))) + mux




L = -2*diag(ones(N)) + diag(ones(N-1),1) + diag(ones(N-1),-1)
L = L[1:-1,:]
Signinv = np.linalg.inv(Sign)

    
def get_x_t(alpha):
    R = A.T.dot(Signinv.dot(A)) + alpha* L.T.dot(L)
    x_t = np.linalg.solve(R, A.T.dot(Signinv.dot(z)))
    res = A.dot(x_t) - z
    #Y = log10(res.dot(Siginv.dot(res)))
    Y = log10(res.dot(res))
    X = log10(norm(L.dot(x_t))**2)
    return x_t,X,Y

Nalpha = 200
alphavec = logspace(-10,0,Nalpha)
Xvec = zeros(Nalpha)
Yvec = zeros(Nalpha)
for i in range(Nalpha):
    x_t,X,Y = get_x_t(alphavec[i])
    Xvec[i] = X
    Yvec[i] = Y
    
# Evaluate second derivative at every point and find max
Cvec = zeros(Nalpha)
for i in range(1,Nalpha-1):
    da = alphavec[i+1]-alphavec[i-1]
    dxda = (Xvec[i+1]-Xvec[i-1])/da
    dyda = (Yvec[i+1]-Yvec[i-1])/da
    da2 = ((alphavec[i+1]-alphavec[i-1])/2)**2
    d2yda2 = (Yvec[i-1] - 2*Yvec[i] + Yvec[i+1])/da2
    d2xda2 = (Xvec[i-1] - 2*Xvec[i] + Xvec[i+1])/da2
    Cvec[i] = -(dxda*d2yda2 - dyda*d2xda2)/(dxda**2 + dyda**2)**(1.5)
    
    #da2 = (alphavec[i+1] - alphavec[i-1])/2
    #d2y = (Yvec[i+1] - 2*Yvec[i] + Yvec[i-1])/da2
    #d2x = (Xvec[i+1] - 2*Xvec[i] + Xvec[i-1])/da2
    #Cvec[i] = d2y/d2x
    
    
istar = argmax(Cvec)
alphastar = alphavec[istar]
x_t,_,_ = get_x_t(alphastar)


plot(Xvec,Yvec,'k.-')
plot(Xvec[istar],Yvec[istar],'ro')
#xlim((Xvec[istar]-0.1,Xvec[istar]+0.1))
#ylim((Yvec[istar]-0.1,Yvec[istar]+0.1))

figure()
sc = 1e6
plot(xtrue/sc,altvec,'k-',linewidth=2,label='Truth')
plot(x_lmmse/sc,altvec,'r-',label='Rec. IRI prior')
plot(x_t/sc,altvec,'g-',label='Rec. Tikh 2nd order, Lcurve')
xlabel('$n_e [cm^{-3} \cdot 10^{-6}]$')
ylabel('Altitude [km]')
legend(loc='best',fancybox=True)

Out[500]:

<matplotlib.legend.Legend at 0x7f127ebe65d0>

What if ne isn't close to IRI?¶

In [ ]:

idx = 67
xtrue = P[:,idx] # pick a "random" column
xtrue[]

ztrue = A.dot(xtrue)

Sign = diag(ztrue)
n = random.multivariate_normal(zeros(N), Sign)

z = ztrue + n

In [ ]:

Use LMMSE estimator by constructing mean and cov matrix using randomized IRI runs¶

Simulate inversion¶

Solve LMMSE problem with IRI-generated statistics¶

Solve using 2nd order Tikhonov with L-curve¶

Put it all together¶

What if ne isn't close to IRI?¶