<jupyter_start><jupyter_text>Fixed parameters (initial model)<jupyter_code>bin_filepath = "/home/nienke/Documents/Research/SS_MTI/External_packages/reflectivity_Mars/SRC/test/crfl_sac"
save_path_OG = "/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/Gradient_descent_classic"
if not exists(join(save_path_OG, "start_v")):
    makedirs(join(save_path_OG, "start_v"))
f_start = join(save_path_OG, "start_v")

## Fixed parameters:
src_depth = 20.0
epi_in_km = 1774.7380
epi = kilometer2degrees(epi_in_km, radius=3389.5)
baz = 0.0

dt = 0.025

phases = ["P", "S", "P", "S", "S"]
comps = ["Z", "T", "R", "Z", "R"]
t_pres = [1, 1, 1, 1, 1]
t_posts = [30, 30, 30, 30, 30]
ylims = [2e-10, 2e-10, 1e-10, 3e-10, 2e-10]

fmin = 0.2
fmax = 0.6
zerophase = False


## Start parameters:
bm_start_model = "/home/nienke/Documents/Research/Data/MTI/MT_vs_STR/bm_models/TAYAK.bm"
m_rr = 0.3000
m_tt = 0.1000 
m_pp = 0.2000
m_rt = 0.2000
m_rp = 0.5000
m_tp = 0.2000
# mtt,mtp,mrt,mpp,mrp,mrr = 0.1000    0.2000    0.2000    0.1000    0.1000    0.2000 # TRUE MODEL


focal_mech = [m_rr, m_tt, m_pp, m_rt, m_rp, m_tp]
Moho_d = 75.0
m0 = np.hstack((focal_mech, Moho_d))
sigmas = np.ones(len(phases)) * 1e-10

Create_Vmod.create_dat_file(
    src_depth, epi_in_km, baz, focal_mech, dt, f_start, bm_start_model,
)
Create_Vmod.update_dat_file(
    dat_folder=f_start,
    m=m0,
    vpvs=False,
    depth=True,
    produce_tvel= True,
    tvel_name = "Init",
)
m0<jupyter_output><empty_output><jupyter_text>Observed data<jupyter_code>path_observed = (
    "/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/obs_2/"
)
npz_file = "/home/nienke/Documents/Research/Data/npz_files/TAYAK.npz"
st_obs = Gradient.read_refl_mseeds(path=path_observed)
Taup = TauPyModel(npz_file)
obs_tts = [PhaseTracer.get_traveltime(Taup, phase, src_depth, epi) for phase in phases]
st_obs_w, st_obs_full, s_obs = Gradient.window(
    st_obs, phases, comps, obs_tts, t_pres, t_posts, fmin, fmax, zerophase,
)<jupyter_output><empty_output><jupyter_text>Classic gradient descent<jupyter_code>update_nr = 25
current_update = 15

epsilon = 0.001
prior_crfl_filepath = join(f_start, "crfl.dat")
Inversion.gradient_descent(
    bin_path=bin_filepath,
    save_path=save_path_OG,
    epsilon=epsilon,
    update_nr=update_nr,
    dt=dt,
    sigmas=sigmas,
    st_obs_w=st_obs_w,
    current_update=current_update,
    prior_crfl_filepath=None,
    alphas=[1e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 1e-2, 1e-1],
    fmin=fmin,
    fmax=fmax,
    phases=phases,
    comps=comps,
    t_pres=t_pres,
    t_posts=t_posts,
)
Gradient.plot_updates(
    save_path=save_path_OG,
    st_obs_full=st_obs_full,
    st_obs_w=st_obs_w,
    obs_tts=obs_tts,
    phases=phases,
    comps=comps,
    t_pres=t_pres,
    t_posts=t_posts,
    fmin=fmin,
    fmax=fmax,
    zerophase=zerophase,
    ylims=ylims,
)
Gradient.plot_misfits(save_path_OG,epsilon=0.001)<jupyter_output>[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]<jupyter_text>Gaus-Newton method   Jacobian calculation<jupyter_code>def Exact_gradient(
    thread_folder: str,
    thread_nr: int,
    dat_folder: str,
    unit_v: np.array,
    bin_filepath: str,
):
    print(thread_nr)
    print(unit_v)
    """ Thread worker function """

    """ Copy .dat file into new folder """
    dat_filepath = join(dat_folder, "crfl.dat")
    subprocess.call(f"scp {dat_filepath} .", shell=True, cwd=thread_folder)
    """ Copy binary file into new folder """
    subprocess.call(f"scp {bin_filepath} .", shell=True, cwd=thread_folder)

    """ Plug unit vector into .dat file """
    Create_Vmod.update_dat_file(
        dat_folder=thread_folder,
        m=unit_v,
        vpvs=False,
        depth=False,
        produce_tvel=True,
        tvel_name="exact",
    )

    """ Run the reflectivity code """
    print(f"Running exact gradient on thread: {thread_nr}\n")
    subprocess.call("./crfl_sac", shell=True, cwd=thread_folder)
    print(f"thread: {thread_nr} is done\n")
    np.save(join(thread_folder, "m.npy"), unit_v)


def Get_J_moment(
    save_path: str,
    prior_dat_folder: str,
    bin_file_path: str,
    phases: [str],
    comps: [str],
    t_pres: [float],
    t_posts: [float],
    fmin: float,
    fmax: float,
    zerophase: bool,
):

    unit_vs = np.array(
        (
            [1, 0, 0, 0, 0, 0],
            [0, 1, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0],
            [0, 0, 0, 1, 0, 0],
            [0, 0, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 1],
        )
    )
    J_approx = []

    threads = []
    dat_folder = prior_dat_folder  # f_start
    """ Creating the exact seismograms: """
    for i in range(6):
        """ Create folder for each thread """
        thread_folder = join(save_path, f"Exact_{i}")
        if not exists(thread_folder):
            makedirs(thread_folder)
        """ Run forward model with each unit vector """
        exact_st = [
            f
            for f in listdir(thread_folder)
            if f.startswith("st00")
            if isfile(join(thread_folder, f))
        ]
        if not exact_st:
            t = threading.Thread(
                target=Exact_gradient,
                args=[thread_folder, i, dat_folder, unit_vs[:, i], bin_filepath],
            )
            threads.append(t)
            t.start()
    if not exact_st:
        for thread in threads:
            thread.join()

    """ Reading in seismograms and substituting into approximate Green's functions: """
    for i in range(6):
        thread_folder = join(save_path, f"Exact_{i}")
        st = Gradient.read_refl_mseeds(path=thread_folder, stack=False)

        """ Window the data """
        npz_name = [
            f
            for f in listdir(thread_folder)
            if f.endswith(".npz")
            if isfile(join(thread_folder, f))
        ]
        if npz_name:
            npz_file = join(thread_folder, npz_name[0],)
            dat_file = join(thread_folder)

            Taup = TauPyModel(npz_file)
            depth = Create_Vmod.read_depth_from_dat(dat_file)
            epi = Create_Vmod.read_epi_from_dat(dat_file)
            syn_tts = []
            for j, phase in enumerate(phases):
                syn_tts.append(PhaseTracer.get_traveltime(Taup, phase, depth, epi))
        else:
            syn_tts = Gradient.get_tt_from_dat_file(phases, thread_folder, "exact")
        st_syn_w, st_syn_full, s_syn = Gradient.window(
            st, phases, comps, syn_tts, t_pres, t_posts, fmin, fmax, zerophase,
        )

        if i == 0:
            J_approx = np.zeros((len(s_syn), 6))
            J_approx[:, i] = s_syn
        else:
            J_approx[:, i] = s_syn
    return J_approx


def proposal(
    m:[float],
    save_path_OG: str,
    new_folder_name: str,
    dat_folder: str,
    phases: [str],
    comps: [str],
    t_pres: [float],
    t_posts: [float],
    dt: float,
    sigmas: [float],
    fmin: float,
    fmax:float,
    zerophase:bool,
):
    save_folder = join(save_path_OG, new_folder_name)
    if not exists(save_folder):
        makedirs(save_folder)

    if exists(join(save_folder, "It_0")):
        st_file = [
            f
            for f in listdir(join(save_folder, "It_0"))
            if f.startswith("st00")
            if isfile(join(save_folder, "It_0", f))
        ]
    else:
        st_file = []
    if st_file:
        st_syn = Gradient.read_refl_mseeds(path=join(save_folder, "It_0"), stack=False)
    else:
        src_str = Gradient.SRC_STR(
            binary_file_path=bin_filepath,
            prior_dat_filepath=join(dat_folder, "crfl.dat"),
            save_folder=save_folder,
            phases=phases,
            components=comps,
            t_pres=t_pres,
            t_posts=t_posts,
            vpvs=False,
            depth=True,
            dt=dt,
            sigmas=sigmas,
            tstars=None,
            fmin=fmin,
            fmax=fmax,
            zerophase=zerophase,
            start_it=0,
        )

        st_syn = src_str.forward(m)
        
        
    """ Window the data """
    if exists(join(save_folder, "It_0")):
        npz_name = [
            f
            for f in listdir(join(save_folder, "It_0"))
            if f.endswith(".npz")
            if isfile(join(save_folder, "It_0", f))
        ]
    else:
        npz_name = []

    if npz_name:
        npz_file = join(save_folder, "It_0", npz_name[0],)
        dat_file = join(save_folder, "It_0")

        Taup = TauPyModel(npz_file)
        depth = Create_Vmod.read_depth_from_dat(dat_file)
        epi = Create_Vmod.read_epi_from_dat(dat_file)
        syn_tts = []
        for i, phase in enumerate(phases):
            syn_tts.append(PhaseTracer.get_traveltime(Taup, phase, depth, epi))
    else:
        syn_tts = Gradient.get_tt_from_dat_file(phases, join(save_folder, "It_0"), m[-1])
    st_syn_w, st_syn_full, s_syn = Gradient.window(
        st_syn, phases, comps, syn_tts, t_pres, t_posts, fmin, fmax, zerophase,
    )
    
    
    return st_syn_w, st_syn_full, s_syn<jupyter_output><empty_output><jupyter_text>Start the inversion:<jupyter_code>save_path_OG = "/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/Gauss_newton/"
update_nr = 4
current_update = 2
update = 0

    
lambd_0 = 100. # Initial lambda value
nu = 1.5

misfits = np.zeros(update_nr + 1)
while update < update_nr:
    if update == 0:
        accepted_folder = f_start
        
#     if current_update != 0:
#         """ Check where the previous update ended and take this crfl.dat file as prior file"""
#         prev_update = current_update - 1

    """ 
    Step 1: Get the Gradient of seismogram w.r.t moment tensor paramters based on unit vectors: 
    """
    J_m = Get_J_moment(
        save_path=join(save_path_OG,f"Update_{update}"),
        prior_dat_folder=accepted_folder,
        bin_file_path=bin_filepath,
        phases=phases,
        comps=comps,
        t_pres=t_pres,
        t_posts=t_posts,
        fmin=fmin,
        fmax=fmax,
        zerophase=zerophase,
    )
    """ 
    Step 2: Forward run with the actual model parameters
    - This is necessary for the Jacobian of the structural parameters
    - This is necessary for the gradient of the misfit w.r.t. the seismometer
    - This is necessary for the misfit calculation
    """


    st_syn_w0, st_syn_full0, s_syn0 = proposal(
        m=m0,
        save_path_OG=save_path_OG,
        new_folder_name=f"Update_{update}",
        dat_folder=accepted_folder,
        phases=phases,
        comps=comps,
        t_pres=t_pres,
        t_posts=t_posts,
        dt=dt,
        sigmas=sigmas,
        fmin=fmin,
        fmax=fmax,
        zerophase = False,
    )

    """ Calculate misfit """
    xi0 = np.sum(Misfit.L2().run_misfit(phases, st_obs_w, st_syn_w0, sigmas ** 2))

    
    """ Save initial guess parameters: """
    if update == 0:
        if not exists(join(save_path_OG, "start_v")):
            makedirs(join(save_path_OG, "start_v"))
        np.save(join(save_path_OG, "start_v","m1_initial.npy"),m0)
        st = Gradient.read_refl_mseeds(path = join(save_path_OG,f"Update_{update}","It_0"),stack = False)
        st.write(join(save_path_OG, "start_v", "st_m1.mseed"), format="MSEED")
        np.save(join(save_path_OG, "start_v", "misfit.npy"),xi0)
        

    """ 
    Step 3: Get the Jacobian of the structural parameters
    """

    """ Then we need S1 and epsilon """
    epsilon = 0.001
    m1 = m0
    m1[-1] += epsilon * m1[-1]

    st_syn_w1, st_syn_full1, s_syn1 = proposal(
        m=m1,
        save_path_OG=join(save_path_OG,f"Update_{update}"),
        new_folder_name="Depth_update",
        dat_folder=accepted_folder,
        phases=phases,
        comps=comps,
        t_pres=t_pres,
        t_posts=t_posts,
        dt=dt,
        sigmas=sigmas,
        fmin=fmin,
        fmax=fmax,
        zerophase = False,
    )

    J_str = np.expand_dims((s_syn1 - s_syn0) / epsilon, axis=1)
    """ 
    Step 4: Combine the two Jacobians (moment tensor + structure)
    """
    J_total = np.hstack((J_m, J_str)) / np.mean(sigmas)**2

    """ 
    Step 5: Get the gradient of the misfit w.r.t. the seismogram 
    (i.e., derivative of the L2-norm) 
    """
    dxi_ds = np.expand_dims(-(s_obs - s_syn0) / (np.mean(sigmas) ** 2), axis=1)


    """ 
    Step 6: Do update
    NOTE: we have a singular matrix unfortunately, so therefore we will make use of the Levenberg-Marquardt algorithm,
    which is basically adding damping.
    Instead of the identity matrix we use diag(J.T@J): 
    Fletcher(1971)-A modified Marquardt subroutine for non-linear least squares, 
    zero values will be filled with average logaritmic values.
    """

#     I = np.eye(J_total.shape[1]) # Levenberg-Marquardt algorithm
    
    diag = np.diag(J_total.T @ J_total)
    diag.setflags(write=1)
    # Fill zero-values with logaritmic average values:
    zero_diag = np.where(diag == 0)
    non_zero_inds = np.where(diag != 0)
    if zero_diag:
        for z in zero_diag[0]:
            print(np.mean(diag[diag != 0]))
            diag[z] =np.exp(np.mean(np.log(diag[diag != 0])))
    diag[non_zero_inds]=0.
    I = np.diag(diag) 
    
    
#     diag = np.diag(J_total.T @ J_total)
#     diag.setflags(write=1)
#     # Fill zero-values with logaritmic average values:
#     zero_diag = np.where(diag == 0)
#     if zero_diag:
#         for z in zero_diag[0]:
#             print(z)
#             diag[z] = np.exp(np.mean(np.log(diag[diag != 0])))
#     I = np.diag(diag)  # <NAME>
   
    """ 
    Step 7: Test two scenarios:
    1. Lambda
    2. Lambda/nu
    """
    update_current = update
    while update_current == update:
        print(f"lamda 0 = {lambd_0}")
        

        xi1s=np.ones(2)*9e9
        m1_props = np.zeros((len(m0),2))
        for i in range(2):    
            """ 1. lambda_0"""
            if i == 0:
                lambd = lambd_0
                prop_folder = join(save_path_OG,f"Update_{update}","lambd_0")

            elif i == 1:
                lambd = lambd_0/nu
                prop_folder = join(save_path_OG,f"Update_{update}","lambd_0_v")

            J_inv = np.linalg.inv(J_total.T @ J_total + lambd * I)
            J_d = J_total.T @ dxi_ds

            m1_props[:,i] = m0 + (J_inv @ J_d)[:,0]
            xi1s[i] = Gradient.SRC_STR(
                    binary_file_path=bin_filepath,
                    prior_dat_filepath=join(accepted_folder,"crfl.dat"),
                    save_folder=prop_folder,
                    phases=phases,
                    components=comps,
                    t_pres=t_pres,
                    t_posts=t_posts,
                    vpvs=False,
                    depth=True,
                    dt=dt,
                    sigmas=sigmas,
                    tstars=None,
                    fmin=fmin,
                    fmax=fmax,
                    zerophase=zerophase,
                    start_it=0,
                ).misfit(m1_props[:,i],st_obs_w)

        ## Make a choice for updating 
        if xi1s[0] < xi0 and xi1s[1] < xi0: 
            """
            Both proposals are smaller in misfit,
            we then choose lambda/nu because we aim for least damping
            """
            print("update 1")
            lambd_0 /= nu 
            # Save Everything in current update:

            np.save(join(save_path_OG,f"Update_{update}",f"m1_.npy"),m1_props[:,1])
            st = Gradient.read_refl_mseeds(path = join(save_path_OG,f"Update_{update}","lambd_0_v","It_0"),stack = False)
            st.write(join(save_path_OG,f"Update_{update}", "st_m1.mseed"), format="MSEED")
            np.save(join(save_path_OG,f"Update_{update}", "misfit.npy"),xi1s[1])
            np.save(join(save_path_OG,f"Update_{update}", "lambda_0_v.npy"),lambd_0 / nu)


            #DO UPTDATE
            accepted_folder = join(save_path_OG,f"Update_{update}","lambd_0_v","It_0")
            update +=1
        elif xi1s[0] < xi0 and xi1s[1] > xi0: 
            print("update 2")
            lambd_0 = lambd_0  
            #DO UPTDATE

            np.save(join(save_path_OG,f"Update_{update}",f"m1_.npy"),m1_props[:,0])
            st = Gradient.read_refl_mseeds(path = join(save_path_OG,f"Update_{update}","lambd_0","It_0"),stack = False)
            st.write(join(save_path_OG,f"Update_{update}", "st_m1.mseed"), format="MSEED")
            np.save(join(save_path_OG,f"Update_{update}", "misfit.npy"),xi1s[0])
            np.save(join(save_path_OG,f"Update_{update}", "lambda_0.npy"),lambd_0)

            accepted_folder = join(save_path_OG,f"Update_{update}","lambd_0","It_0")

            update+= 1
        else: # No update, because both new misfits are smaller
            print("update 3 (no update)")
            lambd_0 *= nu


diag = np.diag(J_total.T @ J_total)
diag.setflags(write=1)
# Fill zero-values with logaritmic average values:
zero_diag = np.where(diag == 0)
non_zero_inds = np.where(diag != 0)
if zero_diag:
    for z in zero_diag[0]:
        print(np.mean(diag[diag != 0]))
        diag[z] =100#np.exp(np.mean(np.log(diag[diag != 0])))
diag[non_zero_inds]=0.
I = np.diag(diag)  
print(I)
J_inv = np.linalg.inv(J_total.T @ J_total + 1e2 * I)
J_d = J_total.T @ dxi_ds

prop = m0 + (J_inv @ J_d)[:,0]
J_inv @ J_total.T
J_total[:,2]
Gradient.plot_updates(
    save_path=save_path_OG,
    st_obs_full=st_obs_full,
    st_obs_w=st_obs_w,
    obs_tts=obs_tts,
    phases=phases,
    comps=comps,
    t_pres=t_pres,
    t_posts=t_posts,
    fmin=fmin,
    fmax=fmax,
    zerophase=zerophase,
    ylims=ylims,
)
fig, ax = plt.subplots(nrows=1, ncols=1, sharex="all", figsize=(8, 8))

update_folders = np.sort(
    np.asarray(
        [int(f.strip("Update_")) for f in listdir(save_path_OG) if f.startswith("misfit")]
    )
)
print(update_folders)

Xs = np.array(
    [
        np.load(join(save_path, f"Update_{update_folders[up_nr]}", f"misfit.npy"))
        for up_nr in update_folders
    ]
)
ms = np.arange(0, len(Xs))
ax.semilogy(ms, Xs)
ax.tick_params(axis="both", which="major", labelsize=15)
ax.tick_params(axis="both", which="minor", labelsize=15)
ax.set_xlabel("Update nr", fontsize=20)
ax.set_ylabel("Misfit", fontsize=20)
ax.set_xticks(ms)
ax.set_xticklabels(ms)
(xi1s < xi0).any()
dxi_ds
J_d
diag
dxi_ds.T.shape
J_total.T @ dxi_ds.T<jupyter_output><empty_output><jupyter_text>Two-step inversion:  Get misfit w.r.t. seismogramSeismogram is computed using m0 (so initial guess)<jupyter_code>""" Forward run with the actual model parameters """
save_folder = join(save_path_OG, "Update_1")
if not exists(save_folder):
    makedirs(save_folder)

st_file = [
    f
    for f in listdir(join(save_folder, "It_0"))
    if f.startswith("st00")
    if isfile(join(save_folder, "It_0", f))
]
if st_file:
    st_syn = Gradient.read_refl_mseeds(path=st_file, stack=False)
else:
    src_str = Gradient.SRC_STR(
        binary_file_path=bin_filepath,
        prior_dat_filepath=join(f_start, "crfl.dat"),
        save_folder=save_folder,
        phases=phases,
        components=comps,
        t_pres=t_pres,
        t_posts=t_posts,
        vpvs=False,
        depth=True,
        dt=dt,
        sigmas=sigmas,
        tstars=None,
        fmin=fmin,
        fmax=fmax,
        zerophase=False,
        start_it=0,
    )

    st_syn = src_str.forward(m0)

""" Window the data """
npz_name = [
    f
    for f in listdir(join(save_folder, "It_0"))
    if f.endswith(".npz")
    if isfile(join(save_folder, "It_0", f))
]
if npz_name:
    npz_file = join(save_folder, "It_0", npz_name[0],)
    dat_file = join(save_folder, "It_0")

    Taup = TauPyModel(npz_file)
    depth = Create_Vmod.read_depth_from_dat(dat_file)
    epi = Create_Vmod.read_epi_from_dat(dat_file)
    syn_tts = []
    for i, phase in enumerate(phases):
        syn_tts.append(PhaseTracer.get_traveltime(Taup, phase, depth, epi))
else:
    syn_tts = Gradient.get_tt_from_dat_file(phases, join(save_folder, "It_0"), m0[-1])
st_syn_w, st_syn_full, s_syn = Gradient.window(
    st_syn, phases, comps, syn_tts, t_pres, t_posts, fmin, fmax, zerophase,
)<jupyter_output><empty_output><jupyter_text>Gradient of seismogram w.r.t. structural update<jupyter_code>""" We have the seismogram of m0 (from previous cell)"""
s0 = s_syn

""" Then we need S1 and epsilon """
epsilon = 0.001
m1 = m0
m1[-1] += epsilon * m1[-1]

save_folder = join(save_path_OG, "Depth_update")
if not exists(save_folder):
    makedirs(save_folder)

st_file = [
    f
    for f in listdir(join(save_folder, "It_0"))
    if f.startswith("st00")
    if isfile(join(save_folder, "It_0", f))
]
if st_file:
    st_syn = Gradient.read_refl_mseeds(path=st_file, stack=False)
else:
    src_str = Gradient.SRC_STR(
        binary_file_path=bin_filepath,
        prior_dat_filepath=join(f_start, "crfl.dat"),
        save_folder=save_folder,
        phases=phases,
        components=comps,
        t_pres=t_pres,
        t_posts=t_posts,
        vpvs=False,
        depth=True,
        dt=dt,
        sigmas=sigmas,
        tstars=None,
        fmin=fmin,
        fmax=fmax,
        zerophase=False,
        start_it=0,
    )

    st_syn = src_str.forward(m1)

""" Window the data """
npz_name = [
    f
    for f in listdir(join(save_folder, "It_0"))
    if f.endswith(".npz")
    if isfile(join(save_folder, "It_0", f))
]
if npz_name:
    npz_file = join(save_folder, "It_0", npz_name[0],)
    dat_file = join(save_folder, "It_0")

    Taup = TauPyModel(npz_file)
    depth = Create_Vmod.read_depth_from_dat(dat_file)
    epi = Create_Vmod.read_epi_from_dat(dat_file)
    syn_tts = []
    for i, phase in enumerate(phases):
        syn_tts.append(PhaseTracer.get_traveltime(Taup, phase, depth, epi))
else:
    syn_tts = Gradient.get_tt_from_dat_file(phases, join(save_folder, "It_0"), m0[-1])
st_syn_w1, st_syn_full1, s1 = Gradient.window(
    st_syn, phases, comps, syn_tts, t_pres, t_posts, fmin, fmax, zerophase,
)
""" Get the approximate gradient of the depth """
J_depth = np.expand_dims((s1 - s0) / epsilon, axis=1)
J_depth.shape<jupyter_output><empty_output><jupyter_text>Add depth gradient to exact gradient<jupyter_code>J_total = np.hstack((J_approx, J_depth))<jupyter_output><empty_output><jupyter_text>Calculate derivative of misfit w.r.t. seismogram (using L2)<jupyter_code># Considering L2:
dxi_ds = np.expand_dims(-(s_obs - s_syn) / (np.mean(sigmas) ** 2), axis=0)
dxi_ds.shape

G_total.shape<jupyter_output><empty_output><jupyter_text>Determine the gradient of the misfit w.r.t. model params<jupyter_code>dxi_dm = dxi_ds @ G_total
dxi_dm
plt.semilogy(dxi_dm[0, :])<jupyter_output><empty_output><jupyter_text>Two step inversion: 1. Invert the moment tensor using the exact gradients<jupyter_code>""" 1.1 Get the Green's function for the moment tensor:"""
save_path_OG = "/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/Exact_Gradient"
G_approx = Get_G(
    save_path=save_path_OG,
    prior_dat_folder=f_start,
    bin_file_path=bin_filepath,
    phases=phases,
    comps=comps,
    t_pres=t_pres,
    t_posts=t_posts,
    fmin=fmin,
    fmax=fmax,
    zerophase=False
)

np.diag(1/(sigmas**2)).shape
Wd = np.diag(np.ones(G_approx.shape[0])) * (1/sigmas[0]**2)
s_obs_dims = np.expand_dims(s_obs,axis = 1)
A = G_approx.T @ G_approx
B = G_approx.T @ s_obs_dims


M = np.linalg.solve(A, B)
# M = np.linalg.lstsq(A, B)[0]
M
plt.figure(figsize=(16,8))

# for i in range(6):
plt.plot(G_approx[:,2])
st_file = "/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/Exact_Gradient/Exact_5/"
st_syn = Gradient.read_refl_mseeds(path=st_file, stack=False)
plt.figure(figsize=(16,8))

for tr in st_syn:
    plt.plot(tr.times(),tr.data)<jupyter_output><empty_output><jupyter_text>2. Invert the structure using approximate gradient<jupyter_code>np.load("/home/nienke/Documents/Research/SS_MTI/External_packages/Test_reflectivity/Gradient_descent_classic/Update_6/X1s_0.001.npy")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1264929424538545, 'fraction_numerical': 0.02121439906456193, 'mean_word_length': 2.7341337907375642, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '7816', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 9613, 'n_tokens_neox': 9163, 'n_words': 1689, 'path': 'Scripts_MT_Structure/.ipynb_checkpoints/Exact_gradients-checkpoint.ipynb', 'repo_name': 'nienkebrinkman/SS_MTI', 'size': '23968'}