<jupyter_start><jupyter_text>Raw data<jupyter_code>pynumdiff.utils.simulate.lorenz_x?
timeseries_length = 4
r = pynumdiff.utils.simulate.lorenz_x(timeseries_length, noise_parameters=[0, 0.2])
x, x_truth, dxdt_truth, _ = r
t = np.linspace(0, timeseries_length, len(x))
layout = fifi.svg_to_axes.FigureLayout('wind_speed.svg', autogenlayers=True, 
                                           make_mplfigures=True, hide_layers=[])

ax = layout.axes[('raw', 'wind_speed')]

ax.plot(t, x, '.', color='blue', markersize=1)
ax.plot(t, x_truth, '--', color='black', markersize=1)
ax.set_xlim(0, 4)
ax.set_ylim(-2,2)
fifi.mpl_functions.adjust_spines(ax, ['left', 'bottom'], 
                                 tick_length=2.5,
                                 yticks = [-2, 0, 2],
                                 spine_locations={'left': 3, 'bottom': 3},
                                 linewidth=0.5)
fifi.mpl_functions.set_fontsize(ax, 6)
    
layout.append_figure_to_layer(layout.figures['raw'], 'raw', cleartarget=True)
layout.write_svg('wind_speed.svg')<jupyter_output><empty_output><jupyter_text>Power spectrakeep 20% of the power in this case<jupyter_code>wind_speed = x
yf = scipy.fftpack.fft(wind_speed)
N = len(t)
dt = np.mean(np.diff(t))
xf = np.linspace(0.0, 1.0/(2.0*dt), N/2)
layout = fifi.svg_to_axes.FigureLayout('wind_speed.svg', autogenlayers=True, 
                                           make_mplfigures=True, hide_layers=[])
ax = layout.axes[('spectra', 'spectra')]

P = np.log(2.0/N * np.abs(yf[:N//2]))
ax.plot(np.log(xf), P, color='black', zorder=-10)
#ax.vlines(np.log(0.005), np.log(1e-4), np.log(1), linestyle='--', color='red')
ax.vlines(np.log(2), np.log(1e-4), np.log(1), linestyle='--', color='red')

frequencies = np.array([1e-2, 1e-1, 1e0, 1e1, 1e2])
power = np.array([1e-4, 1e-3, 1e-2, 1e-1, 1e0])

ax.set_ylim(1e-4, 1)

fifi.mpl_functions.adjust_spines(ax, ['left', 'bottom'], 
                                 xticks=np.log(frequencies),
                                 yticks=np.log(power),
                                 tick_length=2.5,
                                 spine_locations={'left': 3, 'bottom': 3},
                                 linewidth=0.5)

#ax.set_xticklabels([])
#ax.set_yticklabels([])

fifi.mpl_functions.set_fontsize(ax, 6)
    
layout.append_figure_to_layer(layout.figures['spectra'], 'spectra', cleartarget=True)
layout.write_svg('wind_speed.svg')<jupyter_output>/home/caveman/BMD2/lib/python3.5/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in log<jupyter_text>choose gamma based on frequency<jupyter_code>frequency = 2# 0.005
dt
def plot_method(method_name, method_parent, method, wind_speed, frequency, dt):
    
    if method_name == 'savgol':
        gamma_errcorr = np.exp(0.303*np.log(frequency) + 0.61*np.log(dt) - 2.41)
        gamma_rmse = np.exp(-1.1*np.log(frequency) -0.9*np.log(dt) - 5.3)
        spines = ['left']
    elif method_name == 'tvrj':
        gamma_errcorr = np.exp(0.33*np.log(frequency) + 1.13*np.log(dt) - 0.44)
        gamma_rmse = np.exp(-1.23*np.log(frequency) -0.43*np.log(dt) -1.81)
        spines = []
    elif method_name == 'butter':
        gamma_errcorr = np.exp(0.14*np.log(frequency) + 1.11*np.log(dt) - 0.82)
        gamma_rmse = np.exp(-0.88*np.log(frequency) + 0.26*np.log(dt) - 0.59)
        spines = []
    elif method_name == 'kalman':
        gamma_errcorr = np.exp(0.07*np.log(frequency) + 0.46*np.log(dt) - 8.01)
        gamma_rmse = np.exp(-1.12*np.log(frequency) -0.97*np.log(dt) - 6.47)
        spines = []
    
    
    gamma_general = gamma_rmse/5.
    #gamma_general = np.mean([(gamma_rmse), (gamma_errcorr)])
        
    
    
    print('______________________________________________________')
    print(method_name)
    print('gamma_rmse: ', gamma_rmse)
    print('gamma_errcorr: ', gamma_errcorr)
    print('gamma_general: ', gamma_general)


    layout = fifi.svg_to_axes.FigureLayout('wind_speed.svg', autogenlayers=True, 
                                               make_mplfigures=True, hide_layers=[])
    
    ylim = [-15,15]
    xlim = [0,4]
    
    
    ## model discovery
    ax = layout.axes[(method_name, 'model_discovery')]
    params, v = pynumdiff.optimize.__dict__[method_parent].__dict__[method](wind_speed, dt, tvgamma=gamma_errcorr)
    x_smooth, xdot_smooth = pynumdiff.__dict__[method_parent].__dict__[method](wind_speed, dt, params)

    ax.plot(t, xdot_smooth, color='red')
    ax.plot(t, dxdt_truth, '--', color='black')
    ax.set_ylim(ylim[0], ylim[-1])
    ax.set_xlim(xlim[0], xlim[-1])
    fifi.mpl_functions.adjust_spines(ax, spines, 
                                     tick_length=2.5,
                                     #yticks=[-0.2, 0, 0.2, 0.4],
                                     #xticks=[0,20,40,60,80,100],
                                     spine_locations={'left': 3, 'bottom': 3},
                                     linewidth=0.5)
    print('___________model discovery_____________________')
    rmse = np.sqrt(np.mean((xdot_smooth - dxdt_truth)**2))
    errcorr = pynumdiff.utils.evaluate.error_correlation(xdot_smooth, dxdt_truth)
    print('error correlation: ', errcorr)
    print('rmse: ', rmse)
    
    
    ## general
    ax = layout.axes[(method_name, 'general')]
    params, v = pynumdiff.optimize.__dict__[method_parent].__dict__[method](wind_speed, dt, tvgamma=gamma_general)
    x_smooth, xdot_smooth = pynumdiff.__dict__[method_parent].__dict__[method](wind_speed, dt, params)

    ax.plot(t, xdot_smooth, color='red')
    ax.plot(t, dxdt_truth, '--', color='black')
    ax.set_ylim(ylim[0], ylim[-1])
    ax.set_xlim(xlim[0], xlim[-1])
    fifi.mpl_functions.adjust_spines(ax, spines, 
                                     tick_length=2.5,
                                     #yticks=[-0.2, 0, 0.2, 0.4],
                                     #xticks=[0,20,40,60,80,100],
                                     spine_locations={'left': 3, 'bottom': 3},
                                     linewidth=0.5)
    
    print('___________________general_____________________')
    rmse = np.sqrt(np.mean((xdot_smooth - dxdt_truth)**2))
    errcorr = pynumdiff.utils.evaluate.error_correlation(xdot_smooth, dxdt_truth)
    print('error correlation: ', errcorr)
    print('rmse: ', rmse)
    
    
    
    
    ## control
    ax = layout.axes[(method_name, 'control')]
    params, v = pynumdiff.optimize.__dict__[method_parent].__dict__[method](wind_speed, dt, tvgamma=gamma_rmse)
    x_smooth, xdot_smooth = pynumdiff.__dict__[method_parent].__dict__[method](wind_speed, dt, params)

    ax.plot(t, xdot_smooth, color='red')
    ax.plot(t, dxdt_truth, '--', color='black')
    ax.set_ylim(ylim[0], ylim[-1])
    ax.set_xlim(xlim[0], xlim[-1])
    spines.append('bottom')
    fifi.mpl_functions.adjust_spines(ax, spines, 
                                     tick_length=2.5,
                                     #yticks=[-0.2, 0, 0.2, 0.4],
                                     #xticks=[0,20,40,60,80,100],
                                     spine_locations={'left': 3, 'bottom': 3},
                                     linewidth=0.5)
    
    print('___________________control_____________________')
    rmse = np.sqrt(np.mean((xdot_smooth - dxdt_truth)**2))
    errcorr = pynumdiff.utils.evaluate.error_correlation(xdot_smooth, dxdt_truth)
    print('error correlation: ', errcorr)
    print('rmse: ', rmse)
    
    
    

    fifi.mpl_functions.set_fontsize(ax, 6)
    layout.append_figure_to_layer(layout.figures[method_name], method_name, cleartarget=True)
    layout.write_svg('wind_speed.svg')
    
    
    
method_name = 'savgol'
method_parent = 'linear_model'
method = 'savgoldiff'
plot_method(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'tvrj'
method_parent = 'total_variation_regularization'
method = 'jerk_sliding'
plot_method(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'butter'
method_parent = 'smooth_finite_difference'
method = 'butterdiff'
plot_method(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'kalman'
method_parent = 'kalman_smooth'
method = 'constant_acceleration'
plot_method(method_name, method_parent, method, wind_speed, frequency, dt)
def pareto_plot(method_name, method_parent, method, wind_speed, frequency, dt):
    
    if method_name == 'savgol':
        gamma_errcorr = np.exp(0.303*np.log(frequency) + 0.61*np.log(dt) - 2.41)
        gamma_rmse = np.exp(-1.1*np.log(frequency) -0.9*np.log(dt) - 5.3)
        spines = ['left']
    elif method_name == 'tvrj':
        gamma_errcorr = np.exp(0.33*np.log(frequency) + 1.13*np.log(dt) - 0.44)
        gamma_rmse = np.exp(-1.23*np.log(frequency) -0.43*np.log(dt) -1.81)
        spines = []
    elif method_name == 'butter':
        gamma_errcorr = np.exp(0.14*np.log(frequency) + 1.11*np.log(dt) - 0.82)
        gamma_rmse = np.exp(-0.88*np.log(frequency) + 0.26*np.log(dt) - 0.59)
        spines = []
    elif method_name == 'kalman':
        gamma_errcorr = np.exp(0.07*np.log(frequency) + 0.46*np.log(dt) - 8.01)
        gamma_rmse = np.exp(-1.12*np.log(frequency) -0.97*np.log(dt) - 6.47)
        spines = []
    
    #if method_name == 'tvrj':
    #    print(gamma_errcorr, gamma_rmse)
    #    gammas = np.linspace(gamma_errcorr, gamma_rmse/2, 10)
    #else:
    gammas = np.linspace(gamma_errcorr, gamma_rmse, 20)
    errcorrs = []
    rmses = []
    
    
    ## model discovery
    for gamma in gammas:
        params, v = pynumdiff.optimize.__dict__[method_parent].__dict__[method](wind_speed, dt, tvgamma=gamma)
        x_smooth, xdot_smooth = pynumdiff.__dict__[method_parent].__dict__[method](wind_speed, dt, params)

        rmse = np.sqrt(np.mean((xdot_smooth - dxdt_truth)**2))
        errcorr = pynumdiff.utils.evaluate.error_correlation(xdot_smooth, dxdt_truth)
        rmses.append(rmse)
        errcorrs.append(errcorr)
        
    layout = fifi.svg_to_axes.FigureLayout('wind_speed.svg', autogenlayers=True, 
                                               make_mplfigures=True, hide_layers=[])
    
    ax = layout.axes[(method_name+'_pareto', 'pareto')]
    ax.plot(np.log(rmses), errcorrs)
    ax.plot(np.log(rmses), errcorrs, '.')
    ax.plot(np.log(rmses[int(len(gammas)/4)]), errcorrs[int(len(gammas)/4)], '*')
    ax.set_xlim(0.1, 2.5)
    ax.set_ylim(0,1)
    
    fifi.mpl_functions.set_fontsize(ax, 6)
    layout.append_figure_to_layer(layout.figures[method_name+'_pareto'], method_name+'_pareto', cleartarget=True)
    layout.write_svg('wind_speed.svg')
    
    
    
    
    return gammas, errcorrs, rmses
method_name = 'savgol'
method_parent = 'linear_model'
method = 'savgoldiff'
gammas, errcorrs, rmses = pareto_plot(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'butter'
method_parent = 'smooth_finite_difference'
method = 'butterdiff'
gammas, errcorrs, rmses = pareto_plot(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'tvrj'
method_parent = 'total_variation_regularization'
method = 'jerk'
gammas, errcorrs, rmses = pareto_plot(method_name, method_parent, method, wind_speed, frequency, dt)
method_name = 'kalman'
method_parent = 'kalman_smooth'
method = 'constant_acceleration'
gammas, errcorrs, rmses = pareto_plot(method_name, method_parent, method, wind_speed, frequency, dt)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14823900473103208, 'fraction_numerical': 0.029612756264236904, 'mean_word_length': 2.8408479138627185, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '531013', 'license': 'MIT', 'n_tokens_mistral': 4442, 'n_tokens_neox': 4132, 'n_words': 610, 'path': 'notebooks/paper_figures/misc/wind_speed_derivative.ipynb', 'repo_name': 'fossabot/PyNumDiff', 'size': '11414'}