.gitignore

+# Ignore Jupyter Notebook checkpoints
+.ipynb_checkpoints/
+
+#Ignore pycache
+*.pyc
+__pycache__/
+
+#And all figures
+fig*/
+
+func_bk/
+
+#Also no pkl files
+*.pkl
+
+plot-data/
+theresa-data/
+obsolete/
+

LinAdvFunc/LinAdvFunc.egg-info/PKG-INFO

+Metadata-Version: 2.1
+Name: LinAdvFunc
+Version: 0.1

LinAdvFunc/LinAdvFunc.egg-info/SOURCES.txt

+setup.py
+LinAdvFunc.egg-info/PKG-INFO
+LinAdvFunc.egg-info/SOURCES.txt
+LinAdvFunc.egg-info/dependency_links.txt
+LinAdvFunc.egg-info/top_level.txt
\ No newline at end of file

LinAdvFunc/LinAdvFunc.egg-info/dependency_links.txt

+

LinAdvFunc/LinAdvFunc.egg-info/top_level.txt

+

LinAdvFunc/__init__.py

+from LinAdvFunc.da_functions import *
+from LinAdvFunc.misc_functions import *    
+from LinAdvFunc.model_functions import *
+from LinAdvFunc.plot_functions import *
+from LinAdvFunc.efsoi_functions import *
\ No newline at end of file

da_functions.py → LinAdvFunc/da_functions.py

@@ -3,8 +3,9 @@
 #Contains all functions related to setting up and running the ensemble filter. There is quite a bit of overlap with the model functions, but I felt it was a bit cleaner to try to separate them a bit. 
 
 import numpy as np
-from model_functions import *
-from misc_functions import *
+from LinAdvFunc.model_functions import *
+from LinAdvFunc.misc_functions import *
+from numba import jit
 
 
 
@@ -556,7 +557,8 @@ def run_linear_advection_KF_22(m_const,da_const,sat_operator):
                 states[j]["an"] = states[j]["an"]+[analysis]
             
             if da_const['method'] == 'LETKF':
-                analysis,bla = LETKF_analysis_22(background,obs,obs_sat,m_const,da_const,sat_operator)
+                # analysis,bla = LETKF_analysis_22(background,obs,obs_sat,m_const,da_const,sat_operator)
+                analysis,bla,W_a = LETKF_analysis_23(background,obs,obs_sat,m_const,da_const,sat_operator)
                 states[j]["an"] = states[j]["an"]+[analysis]
             
             if da_const['method'] == 'sqEnKF':
@@ -759,6 +761,14 @@ def sum_mid_tri(x):
     idx_end = int(2*nx/3.)
     return np.sum(x[idx_str:idx_end])
 
+def mid_hh(x):
+    """
+    For Herbert Hartl, the response function is simply the value at the middle of the state vector
+    """
+    nx = len(x)
+    idx_str = int(nx/2.)
+    return x[idx_str]
+
 def sum_mid_tri_abs(x):
     """
     Is a sum over the absolute values over the middle of the state 
@@ -1084,8 +1094,12 @@ def loc_matrix(da_const,m_const):
     return C
 
 def gaspari_cohn(x,loc_length):
-    """Gaspari-Cohn function intended to be applied to the x grid, which mirrors the matrix to account for periodic boundary domains."""
+    """
+    Gaspari-Cohn function intended to be applied to the x grid, which mirrors the matrix to account for periodic boundary domains.
+    2024-02-07: fixed the index in the mirroring that mean one side of the mirroring was one index off
+    """
     
+
     ra = x/loc_length
     gp = np.zeros_like(ra)
     i=np.where(ra<=1.)[0]
@@ -1095,7 +1109,8 @@ def gaspari_cohn(x,loc_length):
         
     #Now we mirror things to zero for the periodic boundary domain
     half_idx = int(ra.shape[0]/2)
-    gp[-half_idx:] = gp[half_idx-1::-1]
+    gp[-half_idx+1:] = gp[half_idx-1:0:-1]
+    # gp[-half_idx:] = gp[half_idx-1::-1]
     return gp
 
 def gaspari_cohn_non_mirrored(x,loc_length):
@@ -2081,6 +2096,134 @@ def LETKF_analysis_22(bg,obs,obs_sat,m_const,da_const,sat_operator):
     return x_a, x_ol_a
 
 
+def LETKF_analysis_23(bg,obs,obs_sat,m_const,da_const,sat_operator):
+    """
+    Same as LETKF_analysis_22 but speed things up by putting the middle loop into a separate numba function.
+
+    It seems roughly 10 times faster than the 22 version.
+
+    Further speed up could probably be achieved by reducing the size of the state arrays to not include areas outside the localized observations influence. 
+    
+    Follows the recipe and notation of Hunt 2007, e.g. what is normally called B is now P_a, and _ol refers to over line, so the mean.
+    P_tilde_a refers to B but in ensemble space.
+
+    Has an additional outputs needed for the EFSOI, namely the full weight matrix W_a
+    """
+
+
+
+    n_obs_h  =len(da_const["obs_loc"])
+    n_obs_sat =len(da_const["obs_loc_sat"])
+    n_obs = n_obs_h + n_obs_sat
+    obs_loc = np.hstack([da_const["obs_loc"],da_const["obs_loc_sat"]]).astype(int)
+    
+    
+    L = da_const['nens']
+    x_b = bg
+    x_ol_b = np.mean(x_b,axis=1)
+    X_b = x_b.T-x_ol_b
+    X_b = X_b.T
+
+    R = da_const['R'] 
+   
+    # this is the function that ports the background ensemble state to pertubations in observation space
+    Y_b, y_ol_b = state_to_observation_space(bg,m_const,da_const,sat_operator)
+    
+    # ineligant way to merge the obs vector depending on which observations occur
+    if n_obs_h>0:  y_obs = obs[da_const["obs_loc"]]
+    if n_obs_sat>0: sat_obs = obs_sat[da_const["obs_loc_sat"]]
+    if n_obs_h>0 and n_obs_sat>0 :  y_obs = np.hstack([y_obs,sat_obs])
+    if n_obs_sat>0 and n_obs_h==0 : y_obs   = sat_obs
+    
+    delta_y = y_obs-y_ol_b
+        
+    if da_const['loc']==False:
+        """ Now that all the variables are set, we start by computing the covariance matrix in ensemble state """
+        YRY = np.dot(Y_b.T,np.dot(np.linalg.inv(R),Y_b))
+        P_tilde_a = np.linalg.inv((L-1)*np.identity(L)+YRY)
+
+        """Next step, computing the enesemble mean analysis via the weighting vector w_ol_a"""
+        w_ol_a = np.dot(P_tilde_a,np.dot(Y_b.T,np.dot(np.linalg.inv(R),delta_y)))
+        x_ol_a = x_ol_b+np.dot(X_b,w_ol_a)
+
+        """We now get the ensemble by calculating the weighting matrix through a square root of the error covariance matrix, and adding the mean values to the ensemble deviations"""
+        W_a =  np.real(sqrtm((L-1)*P_tilde_a))
+        w_a = W_a+w_ol_a
+
+        x_a =  np.dot(X_b,w_a).T+ x_ol_a
+        x_a = x_a.T
+
+    else:
+        # print('making sure that the fix is being used')
+        # print('using the old version')
+        # print('trying the third version')
+        # print('sigh, what is not working now?')
+        # bla, blub,bleurgh=LETKF_numba_loop(m_const['x_grid'],da_const['loc_length'],obs_loc,x_ol_b,x_b,R,Y_b,L,delta_y,X_b,m_const['nx'],m_const['dx'])
+        # print('this is annoying, I hate numba')
+        x_a, x_ol_a,W_a=LETKF_numba_loop(m_const['x_grid'],da_const['loc_length'],obs_loc,x_ol_b,x_b,R,Y_b,L,delta_y,X_b,m_const['nx'],m_const['dx'])
+    return x_a, x_ol_a,W_a
+
+# @jit
+def LETKF_numba_loop(x_grid,loc_length,obs_loc,x_ol_b,x_b,R,Y_b,L,delta_y,X_b,N,dx):
+    x_grid_neg = -x_grid[::-1]-dx
+    x_grid_ext = x_grid+dx*N
+    x_a_loc = np.zeros((N,L))
+    W_a_loc = np.zeros((N,L,L))
+    x_ol_a_loc = np.zeros((N))
+    
+    for g in range(N):
+        dist_reg = np.abs(x_grid[obs_loc]-x_grid[g])
+        dist_neg = x_grid[g]-x_grid_neg[obs_loc]
+        dist_ext = x_grid_ext[obs_loc]-x_grid[g]
+        dist = np.minimum(dist_reg,dist_ext)#apparently minimum doesn't like 3 variables   ,dist_neg)
+        dist = np.minimum(dist,dist_neg)
+        #And now we calculate the gaspari cohn weighting and apply it to inverse R 
+        """Gaspari-Cohn function, with no mirroring."""
+    
+        ra = np.abs(dist)/loc_length
+    
+        gp = np.zeros_like(ra)
+        i=np.where(ra<=1.)[0]
+        gp[i]=-0.25*ra[i]**5+0.5*ra[i]**4+0.625*ra[i]**3-5./3.*ra[i]**2+1.
+        i=np.where((ra>1.)*(ra<=2.))[0]
+        gp[i]=1./12.*ra[i]**5-0.5*ra[i]**4+0.625*ra[i]**3+5./3.*ra[i]**2-5.*ra[i]+4.-2./3./ra[i]
+
+        gc_loc=gp 
+
+        if np.max(gc_loc) == 0.:
+            #If no observations are within the double gc radius no need to bother computing things
+            x_ol_a_loc[g] = x_ol_b[g]
+            x_a_loc[g,:]  = x_b[g,:]
+            W_a_loc[g,:,:] = np.identity(L)
+        else:
+            R_inv =  np.linalg.inv(R) * np.diag(gc_loc)
+
+            YRY = np.dot(Y_b.T,np.dot(R_inv,Y_b))
+            """ Now that all the variables are set, we start by computing the covariance matrix in ensemble state """
+            P_tilde_a = np.linalg.inv((L-1)*np.identity(L)+YRY)
+
+            """Next step, computing the enesemble mean analysis via the weighting vector w_ol_a"""
+            w_ol_a = np.dot(P_tilde_a,np.dot(Y_b.T,np.dot(R_inv,delta_y)))
+            x_ol_a = x_ol_b[g]+np.dot(X_b[g],w_ol_a)
+            # Trying to replace sqrtm with numpy eigenvectors so jit can work.
+            evalues, evectors = np.linalg.eigh((L-1.)*P_tilde_a)
+            W_a = evectors * np.sqrt(evalues) @ np.linalg.inv(evectors)
+            # w_a = W_a+w_ol_a
+            # x_a =  np.dot(X_b[g,:],w_a).T+ x_ol_a #this seems to be wrong
+            # x_a =  np.dot(X_b[g,:],w_a).T+ x_ol_b[g] # this should be right but isn't returning the right stuff
+            # This fix somehow makes things worse. I'll try something slightly different
+            # Well try Xa = Xb Wa, and then add x_ol_a to that
+            X_a = np.dot(X_b[g,:],W_a)
+            x_a = X_a+x_ol_a
+            x_ol_a_loc[g]   = x_ol_a
+            x_a_loc[g,:]    = x_a
+            W_a_loc[g,:,:]  = W_a
+            # print(W_a_loc.shape)
+            # print(x_a_loc.shape)
+            # print(x_ol_a_loc.shape)
+
+    return x_a_loc,x_ol_a_loc, W_a_loc
+
 
 def state_to_observation_space(X,m_const,da_const,sat_operator):
     """
@@ -2125,7 +2268,9 @@ def state_to_observation_space(X,m_const,da_const,sat_operator):
 def Kalman_gain_observation_deviations(bg,m_const,da_const,sat_operator):
     """
     Calculates Kalman gain, but using YYT instead of HBH, and XYT intead of BHT,
-    This is important for satellite data because we don't need a linearized version of the satellite operator 
+    This is important for satellite data because we don't need a linearized version of the satellite operator
+
+    Important fix! Localization was applied incorrectly. In stead of localizing BH the whole Kalman gain was localized :( 
     """
     n_obs_h  =len(da_const["obs_loc"])
     n_obs_sat =len(da_const["obs_loc_sat"])
@@ -2143,7 +2288,8 @@ def Kalman_gain_observation_deviations(bg,m_const,da_const,sat_operator):
     if da_const['loc']:
         L_obs, L_obs_state = localization_matrices_observation_space(m_const,da_const)
         YYlocR_inv = np.linalg.inv(L_obs*np.dot(dY_b,dY_b.T)/(L-1)+R)
-        K = L_obs_state*np.dot(X_b,np.dot(dY_b.T,YYlocR_inv))/(L-1)
+        # K = L_obs_state*np.dot(X_b,np.dot(dY_b.T,YYlocR_inv))/(L-1)
+        K = np.dot(L_obs_state*np.dot(X_b,dY_b.T),YYlocR_inv)/(L-1)
     else:
         YYR_inv = np.linalg.inv(np.dot(dY_b,dY_b.T)/(L-1)+R)
         K = np.dot(X_b,np.dot(dY_b.T,YYR_inv))/(L-1)
@@ -2176,7 +2322,7 @@ def square_root_Kalman_gain_observation_deviations(bg,m_const,da_const,sat_opera
         L_obs, L_obs_state = localization_matrices_observation_space(m_const,da_const)
     else:
         L_obs = np.ones([n_obs,n_obs])
-        L_obs_state = np.ones([n_x,n_obs])
+        L_obs_state = np.ones([m_const['nx'],n_obs])
         
     YYR = L_obs*np.dot(dY_b,dY_b.T)/(L-1)+R
     XY  = L_obs_state*np.dot(X_b,dY_b.T)/(L-1)
@@ -2246,13 +2392,20 @@ def single_step_analysis_forecast_22(background,truth,da_const,m_const,sat_opera
     """
     obs, obs_sat = generate_obs_22_single(truth,m_const,da_const,sat_operator,obs_seed) 
     
+    """
+    create dictionary to store this single step 
+    """
+    quad_state = {}
+    
     #Getting the analysis
     if da_const['method'] == 'EnKF':
         np.random.seed(obs_seed+100000)
         an = ENKF_analysis_22(background,obs,obs_sat,da_const,m_const,sat_operator)
     
     if da_const['method'] == 'LETKF':
-        an,bla = LETKF_analysis_22(background,obs,obs_sat,m_const,da_const,sat_operator)
+        # an,bla = LETKF_analysis_22(background,obs,obs_sat,m_const,da_const,sat_operator)
+        an,bla,W_a = LETKF_analysis_23(background,obs,obs_sat,m_const,da_const,sat_operator)
+        quad_state['W_a'] = W_a
     
     if da_const['method'] == 'sqEnKF':
         an = sqEnKF_analysis_22(background,obs,obs_sat,da_const,m_const,sat_operator)
@@ -2275,10 +2428,6 @@ def single_step_analysis_forecast_22(background,truth,da_const,m_const,sat_opera
         fc[:,i]    = linear_advection_model(an[:,i],u_ens,dhdt_ens,m_const["dx"],da_const["dt"],da_const["nt"])
 
     
-    """
-    create dictionary to store this single step 
-    """
-    quad_state = {}
     quad_state['bg'] = background
     quad_state['an'] = an
     quad_state['bf'] = bf

LinAdvFunc/efsoi_functions.py

+#!/usr/bin/env python
+
+#Contains all functions generated for the EFSOI project.  
+
+from LinAdvFunc.model_functions import *
+from LinAdvFunc.misc_functions import *
+from LinAdvFunc.da_functions import *
+import numpy as np
+import seaborn as sns
+import matplotlib.pyplot as plt
+
+
+def set_obs_v_dict(
+    obs_loc         = np.arange(10,299,10),
+    obs_loc_sat     = np.array([]),
+    True_std_obs_h  = 0.01,
+    True_std_obs_sat= 0.01,
+    loc =False,
+    loc_length=2000):
+    
+    """
+    This function sets up the dictionary for the observations.
+    obs_loc is the location of the observations
+    obs_loc_sat is the location of the satellite observations
+    True_std_obs_h is the standard deviation of the observations
+    True_std_obs_sat is the standard deviation of the satellite observations
+    loc is a boolean that determines whether or not to use localization
+    loc_length is the length of the localization
+    loc length stuff is for localizing the explicit sensitivity
+
+    """
+
+    
+    obs_v_dict = {}
+    # obs_loc_sat     = np.arange(30,299,40)
+
+    obs_v_dict["obs_loc"]     = obs_loc   
+    obs_v_dict["obs_loc_sat"] = obs_loc_sat   
+    obs_v_dict["True_std_obs"] = True_std_obs_h
+    obs_v_dict["True_std_obs_sat"] = True_std_obs_sat   
+    obs_v_dict['n_obs_h']   = len(obs_loc)
+    obs_v_dict['n_obs_sat'] = len(obs_loc_sat)
+    obs_v_dict['loc_length'] = loc_length
+    obs_v_dict['loc'] = loc
+    obs_v_dict['loc_type'] = 'gaspari_cohn'
+    return obs_v_dict
+
+
+def localization_matrices_obs_obs(m_const, da_const, obs_v_dict, sigprop_error=1.0):
+    """
+    Creates the localization matrix needed to localize the YaYf covariances used in the EFSOI.
+
+    _a is for analysis, _v for verification (forecast)
+    
+    The distance calculation is still far from elegant.
+
+    Annoyingly, this version doesn't work with the advection routine, because that only works for vectors defined on the regular grid.
+    So in a break with tradition we achieve things by moving the x_grid first
+
+    2024-02-07: fixed something in the x_grid_neg definition, messed up an indes
+    """
+    n_obs_a       = da_const['n_obs_h'] + da_const['n_obs_sat']
+    obs_loc_a     = np.hstack([da_const["obs_loc"], da_const["obs_loc_sat"]]).astype(int)
+    n_obs_v       = obs_v_dict['n_obs_h'] + obs_v_dict['n_obs_sat']
+    obs_loc_v     = np.hstack([obs_v_dict["obs_loc"], obs_v_dict["obs_loc_sat"]]).astype(int)
+    L             = np.zeros([n_obs_a, n_obs_v])
+
+    # Now we need the positions of the verification observations
+    x_grid_v_obs  = m_const['x_grid'][obs_loc_v]
+    x_grid_neg    = -m_const['x_grid'][::-1][obs_loc_v] - m_const['dx']
+    # x_grid_neg    = -x_grid_v_obs[::-1] - m_const['dx'] # This is wrong
+    x_grid_ext    = x_grid_v_obs + (m_const['dx'] * m_const['nx'])
+
+    # Add on advection here with the option to add a multiplicative factor to the advection as a signal propagation error
+    x_grid_v_obs  = x_grid_v_obs - m_const['u_ref'] *sigprop_error* da_const['dt']
+    x_grid_neg    = x_grid_neg   - m_const['u_ref'] *sigprop_error* da_const['dt']
+    x_grid_ext    = x_grid_ext   - m_const['u_ref'] *sigprop_error* da_const['dt']
+
+    for o in range(n_obs_a):
+        g          = obs_loc_a[o] 
+        # First the one using the distance between observations and the state 
+        dist_reg   = np.abs(x_grid_v_obs - m_const['x_grid'][g])
+        dist_neg   = m_const['x_grid'][g] - x_grid_neg
+        dist_ext   = x_grid_ext - m_const['x_grid'][g]
+        dist       = np.minimum(dist_reg, dist_ext)  # apparently minimum doesn't like 3 variables, dist_neg)
+        dist       = np.minimum(dist, dist_neg)
+        gc_loc     = gaspari_cohn_non_mirrored(dist, da_const['loc_length'])
+        L[o, :]    = gc_loc
+
+       
+    return L
+
+
+def quad_plotter_paper_efsoi(quad_state,m_const,da_const,obs_v_dict):
+    """
+    Plots the initial background and blind forecast, as well as the analysis and forecast.
+    Both the assimilated and verification observations are plotted with a standard deviation error bar.
+
+
+    Returns: 
+    figure and axes
+    """
+
+
+    sns.set()
+    sns.set_style('whitegrid')
+
+
+    alpha = np.sqrt(1/da_const['nens'])+0.1
+    
+    fig, ax = plt.subplots(2,2,figsize=(10,6),sharex='all',sharey='all')
+  
+    
+    
+    for i in range(da_const["nens"]):
+        ax[0,0].plot(m_const['x_grid']/1000.,quad_state['bg'][:,i],'r',alpha =alpha,zorder=1)
+        ax[1,0].plot(m_const['x_grid']/1000.,quad_state['bf'][:,i],'b',alpha =alpha,zorder=1)
+        ax[0,1].plot(m_const['x_grid']/1000.,quad_state['an'][:,i],'magenta',alpha =alpha,zorder=1)
+        ax[1,1].plot(m_const['x_grid']/1000.,quad_state['fc'][:,i],'c',alpha =alpha,zorder=1)
+    
+
+    
+    # ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'w',alpha =1,zorder=1,lw=8,label='ens mean')
+    # ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'w',alpha =1,zorder=1,lw=8,label='ens mean')
+    # ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'w',alpha =1,zorder=1,lw=8,label='ens mean')
+    # ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'w',alpha =1,zorder=1,lw=8,label='ens mean')
+    # ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'lightgrey',alpha =1,zorder=2,lw=6,label='ens mean')
+    # ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'lightgrey',alpha =1,zorder=2,lw=6,label='ens mean')
+    # ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'lightgrey',alpha =1,zorder=2,lw=6,label='ens mean')
+    # ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'lightgrey',alpha =1,zorder=2,lw=6,label='ens mean')
+    # ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'r'      ,alpha =1,zorder=2,lw=3,ls='--',label='ens mean')
+    # ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'magenta',alpha =1,zorder=2,lw=3,ls='--',label='ens mean')
+    # ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'b'      ,alpha =1,zorder=2,lw=3,ls='--',label='ens mean')
+    # ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'c'      ,alpha =1,zorder=2,lw=3,ls='--',label='ens mean')
+    ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[0,1].legend(loc='upper right')
+
+    if da_const['n_obs_h']:
+        ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[0,0].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+
+     
+    
+    ax[0,0].set_title('background')
+    ax[1,0].set_title('free-forecast')
+    ax[0,1].set_title('analysis')
+    ax[1,1].set_title('forecast')
+    ax[1,0].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]/1000])
+    plt.subplots_adjust(wspace=0.03,hspace=0.20)
+    ax[1,1].set_xlabel('x [km]')
+    ax[1,0].set_xlabel('x [km]')
+    ax[0,0].set_ylabel(r'$\phi$')
+    ax[1,0].set_ylabel(r'$\phi$')
+    
+    
+    if obs_v_dict['n_obs_h']:
+        ax[1,0].scatter(m_const["x_grid"][obs_v_dict["obs_loc"]] /1000.,quad_state['obs_h_v']      ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].scatter(m_const["x_grid"][obs_v_dict["obs_loc"]] /1000.,quad_state['obs_h_v']      ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].errorbar(m_const["x_grid"][obs_v_dict["obs_loc"]]/1000.,quad_state['obs_h_v']   ,yerr=obs_v_dict['True_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[1,0].errorbar(m_const["x_grid"][obs_v_dict["obs_loc"]]/1000.,quad_state['obs_h_v']   ,yerr=obs_v_dict['True_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+    
+    
+    return fig,ax
+
+def single_step_analysis_forecast_23(background,truth,da_const,m_const,sat_operator,model_seed=0,obs_seed=0,flag_mean=True):
+    """
+    Revisited version of single step analysis forecast to deal with sat data as well.
+    
+
+
+    If the LETKF is used, the ensemble weighting matrix W_a is also returned.
+    """
+    obs, obs_sat = generate_obs_22_single(truth,m_const,da_const,sat_operator,obs_seed) 
+    
+    """
+    create dictionary to store this single step 
+    """
+    quad_state = {}
+    
+    #Getting the analysis
+    if da_const['method'] == 'EnKF':
+        np.random.seed(obs_seed+100000)
+        an = ENKF_analysis_22(background,obs,obs_sat,da_const,m_const,sat_operator)
+    
+    if da_const['method'] == 'LETKF':
+        an,bla,W_a = LETKF_analysis_23(background,obs,obs_sat,m_const,da_const,sat_operator)
+        quad_state['W_a'] = W_a
+    
+    if da_const['method'] == 'sqEnKF':
+        an = sqEnKF_analysis_22(background,obs,obs_sat,da_const,m_const,sat_operator)
+
+    
+    
+
+
+    """
+    Predict blind forecast and forecast
+    """
+    bf = np.zeros_like(background)
+    fc = np.zeros_like(background)
+    
+    # In this version the mean of the ensemble deviations is not forced to be zero
+    if flag_mean==False:
+        for i in range(da_const["nens"]):
+            np.random.seed(i+model_seed)
+            u_ens      = np.random.normal(m_const["u_ref"],da_const["u_std_ens"])
+            dhdt_ens   = np.random.normal(m_const["dhdt_ref"],da_const["dhdt_std_ens"])
+            bf[:,i]    = linear_advection_model(background[:,i],u_ens,dhdt_ens,m_const["dx"],da_const["dt"],da_const["nt"])
+            fc[:,i]    = linear_advection_model(an[:,i],u_ens,dhdt_ens,m_const["dx"],da_const["dt"],da_const["nt"])
+    
+    # This version normalizes the mean ensemble deviations to zero
+    if flag_mean==True:
+        u_ens    = np.zeros(da_const["nens"])
+        dhdt_ens = np.zeros(da_const["nens"])
+        for i in range(da_const["nens"]):
+            np.random.seed(i+model_seed)
+            u_ens[i]      = np.random.normal(m_const["u_ref"],da_const["u_std_ens"])
+            dhdt_ens[i]   = np.random.normal(m_const["dhdt_ref"],da_const["dhdt_std_ens"])
+        # print('dhdt_ens mean before',dhdt_ens.mean())
+        # print('u_ens mean',u_ens.mean())
+        u_ens    = u_ens - u_ens.mean()+m_const["u_ref"] 
+        dhdt_ens = dhdt_ens - dhdt_ens.mean() + m_const["dhdt_ref"]
+        # print('u_ens mean after',u_ens.mean())
+        # print('dhdt_ens mean',dhdt_ens.mean())
+        for i in range(da_const["nens"]):
+            bf[:,i]    = linear_advection_model(background[:,i],u_ens[i],dhdt_ens[i],m_const["dx"],da_const["dt"],da_const["nt"])
+            fc[:,i]    = linear_advection_model(an[:,i],u_ens[i],dhdt_ens[i],m_const["dx"],da_const["dt"],da_const["nt"])
+
+    
+    quad_state['bg'] = background
+    quad_state['an'] = an
+    quad_state['bf'] = bf
+    quad_state['fc'] = fc
+    quad_state['tr_bg'] = truth
+    quad_state['obs'] = obs
+    quad_state['obs_sat'] = obs_sat
+    return quad_state
+
+
+def quad_plotter_sat_a_efsoi(quad_state,m_const,da_const,obs_v_dict,sat_operator):
+    """
+    Not sure how to do this, but will start off trying to plot the assimilated non-linear observations in reflectance space, and the linear verification observations.  
+    Both the assimilated and verification observations are plotted with a standard deviation error bar.
+
+
+    Returns: 
+    figure and axes
+    """
+
+
+    sns.set()
+    sns.set_style('whitegrid')
+
+
+    alpha = np.sqrt(1/da_const['nens'])+0.1
+    
+    # for plotting reflectance
+    sat_an    = sat_operator(quad_state['an'])  
+    sat_bg    = sat_operator(quad_state['bg'])
+    n_obs_h  =len(da_const["obs_loc"])
+    n_obs_sat =len(da_const["obs_loc_sat"])
+
+    if n_obs_sat==0: print('No satellite observations, wrong function buddy')
+     
+    dY_b, y_ol_b = state_to_observation_space(quad_state['bg'],m_const,da_const,sat_operator)
+    y_b = dY_b.T + y_ol_b
+    y_b = y_b.T
+    
+    
+    # for plotting the pixels
+    window = m_const['x_grid'][da_const["obs_loc_sat"][1]]-m_const['x_grid'][da_const["obs_loc_sat"][0]]
+    xpixmin = m_const['x_grid'][da_const["obs_loc_sat"]]-window/2 
+    xpixmax = m_const['x_grid'][da_const["obs_loc_sat"]]+window/2 
+    
+    
+    
+    fig, ax = plt.subplots(2,2,figsize=(10,6),sharex='all',sharey='row')
+  
+    
+    
+    for i in range(da_const["nens"]):
+        # ax[0,0].plot(m_const['x_grid']/1000.,quad_state['bg'][:,i],'r',alpha =alpha,zorder=1)
+        ax[1,0].plot(m_const['x_grid']/1000.,quad_state['bf'][:,i],'b',alpha =alpha,zorder=1)
+        # ax[0,1].plot(m_const['x_grid']/1000.,quad_state['an'][:,i],'magenta',alpha =alpha,zorder=1)
+        ax[1,1].plot(m_const['x_grid']/1000.,quad_state['fc'][:,i],'c',alpha =alpha,zorder=1)
+    
+        ax[0,1].hlines(
+                     sat_an [da_const["obs_loc_sat"],i],
+                     xpixmin/1000,xpixmax/1000,
+                     color='magenta',alpha=0.2,lw=3)
+        ax[0,0].hlines(
+                     sat_bg [da_const["obs_loc_sat"],i],
+                     xpixmin/1000,xpixmax/1000,
+                     color='r',alpha=0.2,lw=3)
+    ax[0,0].hlines(sat_bg [da_const["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='white',alpha=1,lw=4,zorder=1)
+    ax[0,0].hlines(sat_bg [da_const["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='black',alpha=1,lw=2,zorder=2)
+    ax[0,1].hlines(sat_an [da_const["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='white',alpha=1,lw=4,zorder=1)
+    ax[0,1].hlines(sat_an [da_const["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='black',alpha=1,lw=2,zorder=2)
+
+    
+    ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'k'      ,alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'k'      ,alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+
+    if da_const['n_obs_sat']:
+        # ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,zorder=3,s=30,label='verification obs',color='k',marker='_')
+        # ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,zorder=3,s=30,label='verification obs',color='k',marker='_')
+        ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].errorbar(m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,yerr=da_const['used_std_obs_sat'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[0,0].errorbar(m_const["x_grid"][da_const["obs_loc_sat"]]/1000.,quad_state['obs_sat'][da_const["obs_loc_sat"]]   ,yerr=da_const['used_std_obs_sat'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+    
+    
+    ax[0,0].set_title('background')
+    ax[1,0].set_title('free-forecast')
+    ax[0,1].set_title('analysis')
+    ax[1,1].set_title('forecast')
+    ax[1,0].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]/1000])
+    plt.subplots_adjust(wspace=0.03,hspace=0.20)
+    ax[1,1].set_xlabel('x [km]')
+    ax[1,0].set_xlabel('x [km]')
+    ax[0,0].set_ylabel('reflectance')
+    ax[1,0].set_ylabel(r'$\phi$')
+    
+    ax[0,0].set_ylim(0,1)
+    ax[0,0].set_yticks([0,0.3,0.7,1.0])
+    # ax[0,1].yaxis.tick_right()
+    # ax[0,1].tick_params(labelright=True)
+    # ax[0,0].set_yticklabels(['','no cloud','cloudy',''],rotation=70)
+    
+    if obs_v_dict['n_obs_h']:
+        ax[1,0].scatter(m_const["x_grid"][obs_v_dict["obs_loc"]] /1000.,quad_state['obs_h_v']      ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].scatter(m_const["x_grid"][obs_v_dict["obs_loc"]] /1000.,quad_state['obs_h_v']      ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].errorbar(m_const["x_grid"][obs_v_dict["obs_loc"]]/1000.,quad_state['obs_h_v']   ,yerr=obs_v_dict['True_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[1,0].errorbar(m_const["x_grid"][obs_v_dict["obs_loc"]]/1000.,quad_state['obs_h_v']   ,yerr=obs_v_dict['True_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+    
+    
+    return fig,ax
+
+
+def quad_plotter_sat_v_efsoi(quad_state,m_const,da_const,obs_v_dict,sat_operator):
+    """
+    For experiment with non-linear verification observations
+    Both the assimilated and verification observations are plotted with a standard deviation error bar.
+
+
+    Returns: 
+    figure and axes
+    """
+
+
+    sns.set()
+    sns.set_style('whitegrid')
+
+
+    alpha = np.sqrt(1/da_const['nens'])+0.1
+    
+    # for plotting reflectance
+    sat_fc    = sat_operator(quad_state['fc'])  
+    sat_ff    = sat_operator(quad_state['bf'])
+    n_obs_h  = len(obs_v_dict["obs_loc"])
+    n_obs_sat =len(obs_v_dict["obs_loc_sat"])
+
+    if n_obs_sat==0: print('No satellite observations, wrong function buddy')
+     
+    dY_ff, y_ol_ff = state_to_observation_space(quad_state['bf'],m_const,da_const,sat_operator)
+    y_ff = dY_ff.T + y_ol_ff
+    y_ff = y_ff.T
+    
+    
+    # for plotting the pixels
+    window  = m_const['x_grid'][obs_v_dict["obs_loc_sat"][1]]-m_const['x_grid'][obs_v_dict["obs_loc_sat"][0]]
+    xpixmin = m_const['x_grid'][obs_v_dict["obs_loc_sat"]]-window/2 
+    xpixmax = m_const['x_grid'][obs_v_dict["obs_loc_sat"]]+window/2 
+    
+    
+    
+    fig, ax = plt.subplots(2,2,figsize=(10,6),sharex='all',sharey='row')
+  
+    
+    
+    for i in range(da_const["nens"]):
+        ax[0,0].plot(m_const['x_grid']/1000.,quad_state['bg'][:,i],'r',alpha =alpha,zorder=1)
+        # ax[1,0].plot(m_const['x_grid']/1000.,quad_state['bf'][:,i],'b',alpha =alpha,zorder=1)
+        ax[0,1].plot(m_const['x_grid']/1000.,quad_state['an'][:,i],'magenta',alpha =alpha,zorder=1)
+        # ax[1,1].plot(m_const['x_grid']/1000.,quad_state['fc'][:,i],'c',alpha =alpha,zorder=1)
+    
+        ax[1,1].hlines(
+                     sat_fc [obs_v_dict["obs_loc_sat"],i],
+                     xpixmin/1000,xpixmax/1000,
+                     color='c',alpha=0.2,lw=3)
+        ax[1,0].hlines(
+                     sat_ff [obs_v_dict["obs_loc_sat"],i],
+                     xpixmin/1000,xpixmax/1000,
+                     color='b',alpha=0.2,lw=3)
+
+    ax[1,0].hlines(sat_ff [obs_v_dict["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='white',alpha=1,lw=4,zorder=1)
+    ax[1,0].hlines(sat_ff [obs_v_dict["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='black',alpha=1,lw=2,zorder=2)
+    ax[1,1].hlines(sat_fc [obs_v_dict["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='white',alpha=1,lw=4,zorder=1)
+    ax[1,1].hlines(sat_fc [obs_v_dict["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='black',alpha=1,lw=2,zorder=2)
+    # ax[1,0].hlines(sat_ff [obs_v_dict["obs_loc_sat"],:].mean(axis=1),xpixmin/1000,xpixmax/1000,color='blue',alpha=1,lw=2,zorder=3)
+    
+    ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+    ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'k',alpha =1,zorder=2,lw=2,ls='-',label='ens mean')
+
+    if obs_v_dict['n_obs_sat']:
+        ax[1,0].scatter( m_const["x_grid"][obs_v_dict["obs_loc_sat"]]/1000.,quad_state['obs_sat_v']   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].scatter( m_const["x_grid"][obs_v_dict["obs_loc_sat"]]/1000.,quad_state['obs_sat_v']   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[1,1].errorbar(m_const["x_grid"][obs_v_dict["obs_loc_sat"]]/1000.,quad_state['obs_sat_v']   ,yerr=obs_v_dict['True_std_obs_sat'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[1,0].errorbar(m_const["x_grid"][obs_v_dict["obs_loc_sat"]]/1000.,quad_state['obs_sat_v']   ,yerr=obs_v_dict['True_std_obs_sat'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+    
+    if da_const['n_obs_h']:
+        ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='verification obs',color='k')
+        ax[0,1].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[0,0].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+    
+    
+    ax[0,0].set_title('background')
+    ax[1,0].set_title('free-forecast')
+    ax[0,1].set_title('analysis')
+    ax[1,1].set_title('forecast')
+    ax[1,0].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]/1000])
+    plt.subplots_adjust(wspace=0.03,hspace=0.20)
+    ax[1,1].set_xlabel('x [km]')
+    ax[1,0].set_xlabel('x [km]')
+    ax[1,0].set_ylabel('reflectance')
+    ax[0,0].set_ylabel(r'$\phi$')
+    
+    # ax[1,0].set_ylim(0,1)
+    ax[1,0].set_yticks([0.3,0.7])
+    
+    
+    
+    return fig,ax
+
+
+def EFSOI_LETKF(background,truth,m_const,da_const,obs_v_dict,sat_operator,
+                quad_state = None,exp_num=0,obs_seed=0,model_seed=0,obs_v_seed=1000,alpha=0.1,sigprop_error=1.0,flag_mean=False):
+
+    """
+    EFSOI, specifically for the LETKF. It has 3 different versions of computing Ya
+    
+    1. Ya = H Xa   (HXa, not great for non-linear H)
+    2. Ya = Wa_l Yb_l (WlYbl, notation is a bit strange, but the _l indicates the l gridpoint of where the observation is. ) 
+    3. Ya_l = Wa_l Yb (WYb, this Ya needs to calculated separately for each observation point.) 
+ 
+    Shows how saving the ensemble weights leads to the correct answer for dt=0, because the Kalman gain is correctly captured. Is like calculating the Kalman 
+    from the background instead of from the analysis.
+
+    If the observation operator is non-linear, option 2 is better than option 1, option 3 is always better than 2, but not feasible for large systems.
+
+    For now I've thrown out all the localization stuff for the explicit sensitivity. Not sure if or how it should be included in the paper
+    
+    Finally we also have the Sommer formulation, which is the same as the WYb formulation, but with a different error in front. Should be closer to a data denial experiment.
+
+    """
+
+    if obs_v_dict['n_obs_h']*obs_v_dict['n_obs_sat']>0: 
+        print('This function was not tested for both h and sat obs, only for one of them. It could work, but it was not tested')
+        return
+    ###########################################################################
+    #First, need we calculate the basic quad (analysis, forecast, blind forecast)
+    #A precomputed quad_state can be supplied to quickly use different response functions
+    ###########################################################################
+    if quad_state==None:
+        quad_state = single_step_analysis_forecast_23(background,truth,da_const,m_const,sat_operator,obs_seed=obs_seed,model_seed=model_seed,flag_mean=flag_mean)
+
+    # if random_truth_u_flag:
+    np.random.seed(model_seed*da_const['nens'])
+    u_tr      = np.random.normal(m_const["u_ref"],da_const["u_std_ens"])
+    dhdt_tr   = np.random.normal(m_const["dhdt_ref"],da_const["dhdt_std_ens"])
+    tr_ff     = linear_advection_model(truth,u_tr,dhdt_tr,m_const["dx"],da_const["dt"],da_const["nt"])
+    # print('warning,truth uses random speed')
+    # else:
+    # tr_ff     = linear_advection_model(truth,m_const['u_ref'],m_const['dhdt_ref'],m_const["dx"],da_const["dt"],da_const["nt"])
+    quad_state['tr_ff'] = tr_ff
+
+    ###########################################################################
+    # generating verification obs, and model equivalents
+    ###########################################################################
+    obs_h_v, obs_sat_v = generate_obs_22_single(quad_state['tr_ff'],m_const,obs_v_dict,sat_operator,obs_v_seed)
+    if obs_v_dict['n_obs_h']*obs_v_dict['n_obs_sat']>0:
+        obs_v =  np.hstack([obs_h_v[obs_v_dict['obs_loc']],obs_sat_v[obs_v_dict['obs_loc_sat']]])
+        quad_state['obs_h_v']       =  obs_h_v[obs_v_dict['obs_loc']]
+        quad_state['obs_sat_v']     =  obs_sat_v[obs_v_dict['obs_loc_sat']]
+    if obs_v_dict['n_obs_h']>1 and obs_v_dict['n_obs_sat']==0:
+        quad_state['obs_h_v']       =  obs_h_v[obs_v_dict['obs_loc']]
+        quad_state['obs_sat_v']     =  [] 
+        obs_v =  obs_h_v[obs_v_dict['obs_loc']]
+    if obs_v_dict['n_obs_sat']>1 and obs_v_dict['n_obs_h']==0:
+        quad_state['obs_sat_v']     =  obs_sat_v[obs_v_dict['obs_loc_sat']]
+        quad_state['obs_h_v']     =  [] 
+        obs_v =  obs_sat_v[obs_v_dict['obs_loc_sat']]
+
+    
+    #model equivalents of verification observation
+    bf_obs_v_deviations,bf_obs_v_ol = state_to_observation_space(quad_state['bf'],m_const,obs_v_dict,sat_operator)
+    bf_obs_v                        = (bf_obs_v_deviations.T+bf_obs_v_ol).T
+    fc_obs_v_deviations,fc_obs_v_ol = state_to_observation_space(quad_state['fc'],m_const,obs_v_dict,sat_operator)
+    fc_obs_v                        = (fc_obs_v_deviations.T+fc_obs_v_ol).T
+    quad_state['y_bf']        = bf_obs_v
+    quad_state['y_bf_ol']     = bf_obs_v_ol
+    quad_state['y_fc']        = fc_obs_v
+    quad_state['y_fc_ol']     = fc_obs_v_ol
+    quad_state['y_obs_v']     = obs_v 
+
+
+    #And observation equivalent of the observation that is being denied 
+    an_obs_deviations,an_obs_ol = state_to_observation_space(quad_state['an'],m_const,da_const,sat_operator)
+    an_obs                      = (an_obs_deviations.T+an_obs_ol).T
+    bg_obs_deviations,bg_obs_ol = state_to_observation_space(quad_state['bg'],m_const,da_const,sat_operator)
+    
+    #calculating error of mean, and for now we use as error simply the difference to the obs
+    quad_state['e_bf']    = bf_obs_v_ol-obs_v
+    quad_state['e_fc']    = fc_obs_v_ol-obs_v
+    quad_state['e_sq_bf']    = (bf_obs_v_ol-obs_v)**2
+    quad_state['e_sq_fc']    = (fc_obs_v_ol-obs_v)**2
+    
+    quad_state['Jall'] = quad_state['e_sq_fc']-quad_state['e_sq_bf'] 
+
+     
+    
+    # Ok, this was all just foreplay. Now the real fun starts, estimating Jd' via efsoi, and then looping over the observations
+
+    n_obs_v =obs_v_dict['n_obs_h']+obs_v_dict['n_obs_sat']
+    n_obs_a =da_const['n_obs_h']+da_const['n_obs_sat']
+    
+    # Initializing the matrices that will be used to calculate the explicit sensitivity of the observation
+    efsoi_Jd_WYb   = np.zeros([n_obs_a,n_obs_v])
+    efsoi_Jd_som   = np.zeros([n_obs_a,n_obs_v])
+    efsoi_Jd_WlYbl = np.zeros([n_obs_a,n_obs_v])
+    efsoi_Jd_HXa   = np.zeros([n_obs_a,n_obs_v])
+    
+    
+    if da_const['loc']: L= localization_matrices_obs_obs(m_const,da_const,obs_v_dict,sigprop_error=sigprop_error)#Already includes advection
+    else: L = np.ones([n_obs_a,n_obs_v]) 
+
+
+    # Just saving the locactions into a single array
+    obs_loc_a = np.hstack([da_const["obs_loc"],da_const["obs_loc_sat"]]).astype(int)
+    obs_loc_v = np.hstack([obs_v_dict["obs_loc"],obs_v_dict["obs_loc_sat"]]).astype(int)
+    
+    
+    Yf = fc_obs_v_deviations
+    
+    #Now the 3 different ways of computing computing Ya
+    W_a = quad_state['W_a'] 
+    
+    #1. Directly from the analysis state stuff for Ya = H Xa
+    HXa = an_obs_deviations
+    #2. From W_a_ll. This could be done more efficiently by including it in the LETKF analysis function, but for now we do it here so we don't have to change the functions to also pass Ya
+    WlYbl = bg_obs_deviations*0.
+    for a in range(n_obs_a):
+        g = obs_loc_a[a]
+        WlYbl[a,:] = np.dot(bg_obs_deviations[a,:],W_a[g,:,:])
+    
+
+    
+    for o in range(n_obs_v):
+        
+        
+        #LETKF version, which useses the ya vector, and a localized R^-1
+        # first order calculate localization matrix for R
+        g = obs_loc_v[o]
+
+
+        # We need the Wa where the observation was assimilated
+        # Finding the upstream index of the verification observation
+        g_upstream = (g-m_const['u_ref']*da_const['dt']/m_const['dx']).astype(int) 
+        WYb = np.dot(bg_obs_deviations,W_a[g_upstream,:])
+        # print(g,g_upstream)
+
+        # We also need the advected localization matrix
+        gc_loc=L[:,o] 
+        
+        
+        WYbRY_loc   = np.dot(np.dot(np.diag(gc_loc)*np.linalg.inv(da_const['R']),WYb),Yf[o,:])
+        WlYblRY_loc = np.dot(np.dot(np.diag(gc_loc)*np.linalg.inv(da_const['R']),WlYbl),Yf[o,:])
+        HXaRY_loc   = np.dot(np.dot(np.diag(gc_loc)*np.linalg.inv(da_const['R']),HXa),Yf[o,:])
+
+        init_error =quad_state['e_bf'][o]+quad_state['e_fc'][o]
+        sommer_error = 2*quad_state['e_fc'][o]
+        for a in range(n_obs_a):
+            if a < da_const['n_obs_h']:
+                dbg = quad_state['obs'][da_const['obs_loc'][a]]-bg_obs_ol[a]
+            else:
+                dbg = quad_state['obs_sat'][da_const['obs_loc_sat'][a-da_const['n_obs_h']]]-bg_obs_ol[a]
+
+            efsoi_Jd_WYb[a,o]     = init_error  *WYbRY_loc[a]  *dbg/(da_const['nens']-1)
+            efsoi_Jd_som[a,o]     = sommer_error*WYbRY_loc[a]  *dbg/(da_const['nens']-1)
+            efsoi_Jd_WlYbl[a,o]   = init_error  *WlYblRY_loc[a]*dbg/(da_const['nens']-1)
+            efsoi_Jd_HXa[a,o]     = init_error  *HXaRY_loc[a]  *dbg/(da_const['nens']-1)
+
+
+    quad_state['efsoi_WYb']            = efsoi_Jd_WYb
+    quad_state['efsoi_som']            = efsoi_Jd_som
+    quad_state['efsoi_WlYbl']          = efsoi_Jd_WlYbl
+    quad_state['efsoi_HXa']            = efsoi_Jd_HXa
+
+    
+    return quad_state
\ No newline at end of file

plot_functions.py → LinAdvFunc/plot_functions.py

@@ -6,7 +6,7 @@ import seaborn as sns
 import matplotlib.pyplot as plt
 import numpy as np
 
-from da_functions import *
+from LinAdvFunc.da_functions import *
 
 def ensemble_plotter(states,m_const,da_const,j=0,t_start=0,t_end=0):
     """
@@ -253,8 +253,8 @@ def B_plotter(states,j=0,ncyc=0,matrix='bg'):
     
     fig, ax = plt.subplots(1,1,figsize=(5,5))
 
-    ax.set_ylabel('[h1,h2,h3,...]')
-    ax.set_xlabel('[h1,h2,h3,...]')
+    ax.set_ylabel(r'$[\phi_1,\phi_2,\phi_3,...]$')
+    ax.set_xlabel(r'$[\phi_1,\phi_2,\phi_3,...]$')
 
     #B = np.cov(ensemble_diff.transpose())
     B = np.cov(states[j][matrix][ncyc],ddof=1)
@@ -510,35 +510,34 @@ def ensemble_plotter_22(states,m_const,da_const,j=0,t_start=0,t_end=0,h_c=None):
     #for c in range(nc):
     for c in range(na):
         for i in range(da_const["nens"]):
-            ax[c,0].plot(m_const['x_grid'],states[j]['bg'][c+t_start][:,i],'r',alpha =alpha,zorder=1)
-            ax[c,1].plot(m_const['x_grid'],states[j]['an'][c+t_start][:,i],'magenta',alpha =alpha,zorder=1)
-        ax[c,0].plot(m_const['x_grid'],np.mean(states[j]['bg'][c+t_start][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
-        ax[c,0].plot(m_const["x_grid"],states[j]['truth'][c+t_start],'k',zorder=10,label='truth')
+            ax[c,0].plot(m_const['x_grid']/1000.,states[j]['bg'][c+t_start][:,i],'r',alpha =alpha,zorder=1)
+            ax[c,1].plot(m_const['x_grid']/1000.,states[j]['an'][c+t_start][:,i],'magenta',alpha =alpha,zorder=1)
+        ax[c,0].plot(m_const['x_grid']/1000.,np.mean(states[j]['bg'][c+t_start][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+        ax[c,0].plot(m_const["x_grid"]/1000.,states[j]['truth'][c+t_start],'k',zorder=10,label='truth')
         
-        ax[c,1].plot(m_const['x_grid'],np.mean(states[j]['an'][c+t_start][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
-        ax[c,1].plot(m_const["x_grid"],states[j]['truth'][c+t_start],'k',zorder=10,label='truth')
+        ax[c,1].plot(m_const['x_grid']/1000.,np.mean(states[j]['an'][c+t_start][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+        ax[c,1].plot(m_const["x_grid"]/1000.,states[j]['truth'][c+t_start],'k',zorder=10,label='truth')
         
-        if da_const['n_obs_h']>0:
-            ax[c,0].scatter(m_const["x_grid"][da_const["obs_loc"]],states[j]['obs'][c+t_start][da_const["obs_loc"]],zorder=3,label='observations',c='k')
-            ax[c,1].scatter(m_const["x_grid"][da_const["obs_loc"]],states[j]['obs'][c+t_start][da_const["obs_loc"]],zorder=3,label='observations',c='k')
+        # if da_const['n_obs_h']>0:
+        if da_const['n_obs_h']:
+            ax[c,0].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,states[j]['obs'][c+t_start][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+            ax[c,1].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,states[j]['obs'][c+t_start][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+            ax[c,1].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,states[j]['obs'][c+t_start][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+            ax[c,0].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,states[j]['obs'][c+t_start][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
         if h_c is not None:
-            ax[c,0].hlines(h_c,m_const['x_grid'][0],m_const['x_grid'][-1],color='k',ls=':',label='sat threshold')
-            ax[c,1].hlines(h_c,m_const['x_grid'][0],m_const['x_grid'][-1],color='k',ls=':',label='sat threshold')
+            ax[c,0].hlines(h_c,m_const['x_grid'][0]/1000.,m_const['x_grid'][-1]/1000.,color='k',ls=':',label='sat threshold')
+            ax[c,1].hlines(h_c,m_const['x_grid'][0]/1000.,m_const['x_grid'][-1]/1000.,color='k',ls=':',label='sat threshold')
 
-        #if da_const['n_obs_sat']>0:
-        #    #sat obs location will be added as well, I'll add them in with their value for now Maybe I'll get fancy and use a full circle where cloudy and empty circle where empty, because why not
-        #    ax[c,0].scatter(m_const["x_grid"][da_const["obs_loc_sat"]],states[j]['obs_sat'][c+t_start][da_const["obs_loc_sat"]],zorder=3,label='sat obs loc',marker='x',c='k')
-        #    ax[c,1].scatter(m_const["x_grid"][da_const["obs_loc_sat"]],states[j]['obs_sat'][c+t_start][da_const["obs_loc_sat"]],zorder=3,label='sat obs loc',marker='x',c='k')
-        ax[c,0].set_ylabel('assim step '+str(c+t_start) + '\n h [m]')
-        
+        # ax[c,0].set_ylabel('assim step: ' + str(c+t_start)+ '\n' +r'$\phi$') 
+        ax[c,0].set_ylabel('step: ' + str(c+t_start)+ '\n' +r'$\phi$') 
     
-    ax[0,0].set_title('background ensemble and observations')
-    ax[0,1].set_title('analysis ensemble after update')
-    ax[0,1].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]])
+    ax[0,0].set_title('background ensemble \n and observations')
+    ax[0,1].set_title('analysis ensemble \n after assimilation')
+    ax[0,1].set_xlim([m_const["x_grid"][0]/1000.,m_const["x_grid"][-1]/1000.])
     plt.subplots_adjust(wspace=0.01,hspace=0.01)
 #     ax[0].set_xticklabels([])
-    ax[na-1,1].set_xlabel('x [m]')
-    ax[na-1,0].set_xlabel('x [m]')
+    ax[na-1,1].set_xlabel('x [km]')
+    ax[na-1,0].set_xlabel('x [km]')
     plt.legend()
     return fig,ax
 
@@ -608,6 +607,134 @@ def quad_plotter_22(quad_state,m_const,da_const):
     return fig,ax
 
 
+def quad_plotter_25(quad_state,m_const,da_const):
+    """
+    Plots the initial background and blind forecast, as well as the analysis and forecast with observations.
+    
+    Difference to earlier _22 is the arrangement of the subplots, and the addition of uncertainity bars for the observations. 
+ 
+
+    Returns: 
+    figure and axes
+    """
+    sns.set()
+    sns.set_style('whitegrid')
+
+
+    alpha = np.sqrt(1/da_const['nens'])+0.1
+    
+    fig, ax = plt.subplots(2,2,figsize=(7.5,4),sharex='all',sharey='all')
+  
+    
+    
+    for i in range(da_const["nens"]):
+        ax[0,0].plot(m_const['x_grid']/1000.,quad_state['bg'][:,i],'r',alpha =alpha,zorder=1)
+        ax[1,0].plot(m_const['x_grid']/1000.,quad_state['bf'][:,i],'b',alpha =alpha,zorder=1)
+        ax[0,1].plot(m_const['x_grid']/1000.,quad_state['an'][:,i],'magenta',alpha =alpha,zorder=1)
+        ax[1,1].plot(m_const['x_grid']/1000.,quad_state['fc'][:,i],'c',alpha =alpha,zorder=1)
+    
+
+    ax[0,0].plot(m_const["x_grid"]/1000.,quad_state['tr_bg'],'k',zorder=10,label='truth')
+    ax[0,1].plot(m_const["x_grid"]/1000.,quad_state['tr_bg'],'k',zorder=10,label='truth')
+    #ax[0,1].plot(m_const["x_grid"],quad_stat/1000.e['tr_fc'],'k',zorder=10,label='truth')
+    #ax[1,1].plot(m_const["x_grid"],quad_stat/1000.e['tr_fc'],'k',zorder=10,label='truth')
+    
+    ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+
+    if da_const['n_obs_h']:
+        ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+        ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+        ax[0,1].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[0,0].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        # ax[0,0].scatter(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]],zorder=3,label='observations',color='k')
+        # ax[1,0].scatter(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]],zorder=3,label='observations',color='k')
+    
+    
+    ax[0,0].set_title('background')
+    ax[1,0].set_title('free-forecast')
+    ax[0,1].set_title('analysis')
+    ax[1,1].set_title('forecast')
+    ax[0,1].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]/1000])
+    plt.subplots_adjust(wspace=0.03,hspace=0.20)
+    ax[1,1].set_xlabel('x [km]')
+    ax[1,0].set_xlabel('x [km]')
+    ax[0,0].set_ylabel(r'$\phi$')
+    ax[1,0].set_ylabel(r'$\phi$')
+    ax[1,0].legend()
+    
+    # now to add in shading for response domain
+    ylimits = ax[0,0].get_ylim()
+    ax[1,0].fill_between([10,20], ylimits[0], y2=ylimits[1],color='grey',alpha=0.2,zorder=-1)
+    ax[1,1].fill_between([10,20], ylimits[0], y2=ylimits[1],color='grey',alpha=0.2,zorder=-1)
+    ax[1,1].set_ylim(ylimits)
+    
+    return fig,ax
+
+def quad_plotter_hh(quad_state,m_const,da_const):
+    """
+    Plots the initial background and blind forecast, as well as the analysis and forecast with observations.
+    
+    Intended for Herbert Hartl, but only difference to _25 is that the grey shading used to show the middle of the domain is not used. 
+ 
+
+    Returns: 
+    figure and axes
+    """
+    sns.set()
+    sns.set_style('whitegrid')
+
+
+    alpha = np.sqrt(1/da_const['nens'])+0.1
+    
+    fig, ax = plt.subplots(2,2,figsize=(7.5,4),sharex='all',sharey='all')
+  
+    
+    
+    for i in range(da_const["nens"]):
+        ax[0,0].plot(m_const['x_grid']/1000.,quad_state['bg'][:,i],'r',alpha =alpha,zorder=1)
+        ax[1,0].plot(m_const['x_grid']/1000.,quad_state['bf'][:,i],'b',alpha =alpha,zorder=1)
+        ax[0,1].plot(m_const['x_grid']/1000.,quad_state['an'][:,i],'magenta',alpha =alpha,zorder=1)
+        ax[1,1].plot(m_const['x_grid']/1000.,quad_state['fc'][:,i],'c',alpha =alpha,zorder=1)
+    
+
+    ax[0,0].plot(m_const["x_grid"]/1000.,quad_state['tr_bg'],'k',zorder=10,label='truth')
+    ax[0,1].plot(m_const["x_grid"]/1000.,quad_state['tr_bg'],'k',zorder=10,label='truth')
+    #ax[0,1].plot(m_const["x_grid"],quad_stat/1000.e['tr_fc'],'k',zorder=10,label='truth')
+    #ax[1,1].plot(m_const["x_grid"],quad_stat/1000.e['tr_fc'],'k',zorder=10,label='truth')
+    
+    ax[0,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bg'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[0,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['an'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[1,0].plot(m_const['x_grid']/1000.,np.mean(quad_state['bf'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+    ax[1,1].plot(m_const['x_grid']/1000.,np.mean(quad_state['fc'][:,:],axis=1),'k--',alpha =1,zorder=2,label='ens mean')
+
+    if da_const['n_obs_h']:
+        ax[0,0].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+        ax[0,1].scatter( m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,zorder=3,s=10,label='obs',color='k')
+        ax[0,1].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        ax[0,0].errorbar(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]]   ,yerr=da_const['used_std_obs'] ,zorder=3,color='k',capsize=3,capthick=2,ls='none')
+        # ax[0,0].scatter(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]],zorder=3,label='observations',color='k')
+        # ax[1,0].scatter(m_const["x_grid"][da_const["obs_loc"]]/1000.,quad_state['obs'][da_const["obs_loc"]],zorder=3,label='observations',color='k')
+    
+    
+    ax[0,0].set_title('background')
+    ax[1,0].set_title('free-forecast')
+    ax[0,1].set_title('analysis')
+    ax[1,1].set_title('forecast')
+    ax[0,1].set_xlim([m_const["x_grid"][0],m_const["x_grid"][-1]/1000])
+    plt.subplots_adjust(wspace=0.03,hspace=0.20)
+    ax[1,1].set_xlabel('x [km]')
+    ax[1,0].set_xlabel('x [km]')
+    ax[0,0].set_ylabel(r'$\phi$')
+    ax[1,0].set_ylabel(r'$\phi$')
+    ax[1,0].legend()
+    
+
+    
+    return fig,ax
+
 def plot_ensemble_sat_analysis(bg,an,obs,obs_sat,truth,sat_operator,m_const,da_const,h_c=None):
     fig ,ax = plt.subplots(2,2,figsize=(10,6),sharey='col',sharex='all')
     """
@@ -909,8 +1036,8 @@ def plot_J_quad_paper(J_dict,quad,sens,dx,bw=0.3,dJ=True):
     dX_an=dx
     dJ_ff=np.dot(sens,dX_bg)
     dJ_fc=np.dot(sens,dX_an)
-    print('vr_reductions:',np.var(dJ_fc,ddof=1)-np.var(dJ_ff,ddof=1 ),np.var(J_dict['fc'],ddof=1)-np.var(J_dict['bf'],ddof=1))
-    print('variance:',np.var(J_dict['bf'],ddof=1),np.var(dJ_ff,ddof=1),np.var(J_dict['fc'],ddof=1),np.var(dJ_fc,ddof=1 ))
+    print('estimated and real var reductions:',np.var(dJ_fc,ddof=1)-np.var(dJ_ff,ddof=1 ),np.var(J_dict['fc'],ddof=1)-np.var(J_dict['bf'],ddof=1))
+    print('variance of free forecast, estimated variance of ff, and the same for the forecast:',np.var(J_dict['bf'],ddof=1),np.var(dJ_ff,ddof=1),np.var(J_dict['fc'],ddof=1),np.var(dJ_fc,ddof=1 ))
             #'response' : np.hstack([J_dict['bf']-np.mean(J_dict['bf']),dJ_ff,J_dict['fc']-np.mean(J_dict['fc']),J_dict['es']-np.mean(J_dict['es'])]),
     if dJ:
         plot_data = {
@@ -929,7 +1056,8 @@ def plot_J_quad_paper(J_dict,quad,sens,dx,bw=0.3,dJ=True):
     'boxprops':{'facecolor':'none', 'edgecolor':'black'},
     }
     #ax = sns.violinplot(data=plot_data, inner='quartile', orient="v",cut=0,bw=bw,y='response',x='x_pos',palette=my_pal)#sns.color_palette('cool',n_colors=3))#,x='type')#,y='response',x='cyc',hue='type',,split=True,palette={dict_label[left_var]:dict_color[left_var],dict_label[right_var]:dict_color[right_var]}
-    ax = sns.stripplot(data=plot_data, y='response',x='x_pos',alpha=0.7,jitter=0.15,size=5,palette=my_pal)#color='0.0')#
+    ax = sns.stripplot(data=plot_data, y='response', x='x_pos', hue='x_pos', alpha=0.7, jitter=0.15, size=5, palette=my_pal, legend=False)#color='0.0')#
+    # ax = sns.stripplot(data=plot_data, y='response',x='x_pos',alpha=0.7,jitter=0.15,size=5,palette=my_pal)#color='0.0')#
     #ax = sns.boxplot(data=plot_data, y='response',x='x_pos',showfliers=False,**PROPS)#,patch_artist=False)#color='0.0')#,palette=my_pal
     #plot errorbars
     plt.errorbar(np.arange(4),np.zeros(4),[np.std(J_dict['bf'],ddof=1),np.std(dJ_ff,ddof=1),np.std(J_dict['fc'],ddof=1),np.std(dJ_fc,ddof=1 )],fmt='.',capsize=15,lw=3,color='k') 

LinAdvFunc/setup.py

+from setuptools import setup, find_packages
+
+setup(
+    name="LinAdvFunc",
+    version="0.1",
+    packages=find_packages(),
+)