dneuman · February 25, 2022 21:25
diff --git a/Impulse_Response.py b/Impulse_Response.py
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
 Derive  Earth system response from atmospheric composition and
 temperature record using Fourier Ttansform to perform a deconvolution.

 Created on Thu Feb 24 18:42:31 2022

 @author: dan
 """
 import matplotlib.pyplot as plt
 import pandas as pd
 import numpy as np
 import datastore as ds  # private data support
 import tools as tls  # private misc tools
 from scipy.stats import pearsonr


 #%precision 2
 pd.options.display.float_format = '{:.2f}'.format  # change print format
 yn, mn, dn, en, nn, sn = ['Year', 'Month', 'Data',
                          'Error', 'Normalized', 'Smooth']
 emn, epn = ['Err-', 'Err+']  # difference from mean/median

 dst = ds.dst  # Data store object for getting data.

 def qplot(name='quick'):
    fig, axs = plt.subplots(3, 1, num=name, clear=True, sharex=True,
                            figsize=(6.75, 9))
    return axs

 def work(remove_mei=True):
    # relative strengths for best fit, from previous work.
    rc = np.array([0.79008,    # forcing
                   0.50020,    # so2
                   0.78961,    # vol
                   -0.15628])  # offset
    suffix = '' # for file names
    tsuffix = '' # for titles
    frc = dst.forcing().loc[1860:]
    frc.fillna(value=0.0, inplace=True)
    frc['so2'] *= rc[1]/rc[0]
    frc['vol'] *= rc[2]/rc[0]
    f = frc.sum(axis=1)
    t = dst.best().loc[1860:, dn]
    pie = list(range(1860, 1900))  # pre-industrial era
    t -= t.loc[pie].mean()  # remove PIE average
    n = len(t)//2

    if remove_mei:
        m = dst.mei()
        mei = pd.Series(index=t.index, data=0)
        mei[m.index] = m[dn]
        st = tls.loess(t, 12)[sn]  # smoothed temperatue
        # the optimal window value above (12) was found by iteration
        dt = t - st  # residual
        A = np.vstack([mei.values.T, np.ones(len(mei))]).T
        c = np.linalg.lstsq(A, dt.values, rcond=None)[0]  # get best fit
        mei = pd.Series(index=mei.index, data=A @ c, name='MEI')
        vt = dt.std()
        t -= mei
        suffix = '-mei'
        tsuffix = ' (MEI removed)'
        # plot comparison of mei vs temperature residual
        fig, mx = plt.subplots(num='MEI', clear=True, figsize=(12, 6.75))
        mx.plot(dt, label='Temperature Index', alpha=.25)
        mx.plot(mei, label='MEI', alpha=.25)
        mx.plot(dt-mei, color='r', label='Difference')
        scale = 1. - (dt-mei).std()/vt
        text = f'Standard Deviation improved by {scale*100:.1f}%'
        mx.text(.02, .98, text, transform=mx.transAxes, va='top')
        mx.legend()


    # Now take FFTs of t and f to put in frequency domain.
    ff = np.fft.fft(f)
    tf = np.fft.fft(t)
    kf = tf/ff
    labels = 'Temperature,Forcing,Impulse Response (T/F)'.split(',')
    #plot ff, tf, and kf (frequency domains)
    fx, tx, kx = qplot('FFT'+suffix)
    for ax, y, lbl in zip([tx, fx, kx], [tf, ff, kf], labels):
        ax.plot(np.abs(y[:n+1]), label=lbl)
        ax.legend(loc='upper right')
    fx.set_title('Fourier Transforms'+tsuffix, loc='left',
             size='xx-large', weight='bold')

    k = np.fft.ifft(kf)
    x = t.index
    # plot f, t, and k
    fx, tx, kx = qplot('IFFT'+suffix)
    for ax, y, lbl in zip([tx, fx, kx], [t, f, k.real], labels):
        ax.plot(x, y, label=lbl)
        ax.legend(loc='upper right')
    fx.set_title('Global Values'+tsuffix, loc='left',
             size='xx-large', weight='bold')
    # plot (F*K) against T
    fig, ax = plt.subplots(num='Modeled'+suffix, clear=True,
                        figsize=(12., 6.75))
    n = len(k)
    # make kernel with past and future influence, but future zeroed
    kernel = np.zeros(2*n-1, dtype=np.complex128)
    kernel[-n:] = k
    fk = np.convolve(f, kernel, mode='valid')
    fk = fk.real
    # fk = np.abs(fk)
    ax.plot(t, label='Temperature')
    ax.plot(x, fk, label='Impulse Model')
    # Get correlations
    R2 = tls.R2(fk, t)
    P, p = pearsonr(fk, t)
    text = f'R²: {R2:0.3f}\nPearson: {P:0.3f}, p={p:0.5f}'
    ax.text(0.1, 0.8, text, transform=ax.transAxes, size='large')
    ax.set_title('Impulse Model'+tsuffix, loc='left',
             size='xx-large', weight='bold')
    ax.legend(loc='upper right')

    plt.show()
	#!/usr/bin/env python3
	# -- coding: utf-8 --
	"""
	Derive Earth system response from atmospheric composition and
	temperature record using Fourier Ttansform to perform a deconvolution.

	Created on Thu Feb 24 18:42:31 2022

	@author: dan
	"""
	import matplotlib.pyplot as plt
	import pandas as pd
	import numpy as np
	import datastore as ds # private data support
	import tools as tls # private misc tools
	from scipy.stats import pearsonr


	#%precision 2
	pd.options.display.float_format = '{:.2f}'.format # change print format
	yn, mn, dn, en, nn, sn = ['Year', 'Month', 'Data',
	'Error', 'Normalized', 'Smooth']
	emn, epn = ['Err-', 'Err+'] # difference from mean/median

	dst = ds.dst # Data store object for getting data.

	def qplot(name='quick'):
	fig, axs = plt.subplots(3, 1, num=name, clear=True, sharex=True,
	figsize=(6.75, 9))
	return axs

	def work(remove_mei=True):
	# relative strengths for best fit, from previous work.
	rc = np.array([0.79008, # forcing
	0.50020, # so2
	0.78961, # vol
	-0.15628]) # offset
	suffix = '' # for file names
	tsuffix = '' # for titles
	frc = dst.forcing().loc[1860:]
	frc.fillna(value=0.0, inplace=True)
	frc['so2'] *= rc[1]/rc[0]
	frc['vol'] *= rc[2]/rc[0]
	f = frc.sum(axis=1)
	t = dst.best().loc[1860:, dn]
	pie = list(range(1860, 1900)) # pre-industrial era
	t -= t.loc[pie].mean() # remove PIE average
	n = len(t)//2

	if remove_mei:
	m = dst.mei()
	mei = pd.Series(index=t.index, data=0)
	mei[m.index] = m[dn]
	st = tls.loess(t, 12)[sn] # smoothed temperatue
	# the optimal window value above (12) was found by iteration
	dt = t - st # residual
	A = np.vstack([mei.values.T, np.ones(len(mei))]).T
	c = np.linalg.lstsq(A, dt.values, rcond=None)[0] # get best fit
	mei = pd.Series(index=mei.index, data=A @ c, name='MEI')
	vt = dt.std()
	t -= mei
	suffix = '-mei'
	tsuffix = ' (MEI removed)'
	# plot comparison of mei vs temperature residual
	fig, mx = plt.subplots(num='MEI', clear=True, figsize=(12, 6.75))
	mx.plot(dt, label='Temperature Index', alpha=.25)
	mx.plot(mei, label='MEI', alpha=.25)
	mx.plot(dt-mei, color='r', label='Difference')
	scale = 1. - (dt-mei).std()/vt
	text = f'Standard Deviation improved by {scale*100:.1f}%'
	mx.text(.02, .98, text, transform=mx.transAxes, va='top')
	mx.legend()


	# Now take FFTs of t and f to put in frequency domain.
	ff = np.fft.fft(f)
	tf = np.fft.fft(t)
	kf = tf/ff
	labels = 'Temperature,Forcing,Impulse Response (T/F)'.split(',')
	#plot ff, tf, and kf (frequency domains)
	fx, tx, kx = qplot('FFT'+suffix)
	for ax, y, lbl in zip([tx, fx, kx], [tf, ff, kf], labels):
	ax.plot(np.abs(y[:n+1]), label=lbl)
	ax.legend(loc='upper right')
	fx.set_title('Fourier Transforms'+tsuffix, loc='left',
	size='xx-large', weight='bold')

	k = np.fft.ifft(kf)
	x = t.index
	# plot f, t, and k
	fx, tx, kx = qplot('IFFT'+suffix)
	for ax, y, lbl in zip([tx, fx, kx], [t, f, k.real], labels):
	ax.plot(x, y, label=lbl)
	ax.legend(loc='upper right')
	fx.set_title('Global Values'+tsuffix, loc='left',
	size='xx-large', weight='bold')
	# plot (F*K) against T
	fig, ax = plt.subplots(num='Modeled'+suffix, clear=True,
	figsize=(12., 6.75))
	n = len(k)
	# make kernel with past and future influence, but future zeroed
	kernel = np.zeros(2*n-1, dtype=np.complex128)
	kernel[-n:] = k
	fk = np.convolve(f, kernel, mode='valid')
	fk = fk.real
	# fk = np.abs(fk)
	ax.plot(t, label='Temperature')
	ax.plot(x, fk, label='Impulse Model')
	# Get correlations
	R2 = tls.R2(fk, t)
	P, p = pearsonr(fk, t)
	text = f'R²: {R2:0.3f}\nPearson: {P:0.3f}, p={p:0.5f}'
	ax.text(0.1, 0.8, text, transform=ax.transAxes, size='large')
	ax.set_title('Impulse Model'+tsuffix, loc='left',
	size='xx-large', weight='bold')
	ax.legend(loc='upper right')

	plt.show()