From c6fdcf844a194896065777818786e38200c84417 Mon Sep 17 00:00:00 2001 From: berberto Date: Fri, 15 Oct 2021 13:00:03 +0100 Subject: [PATCH 1/4] packaged in sindyae --- examples/lorenz/analyze_lorenz_model1.py | 317 +++++++++++++++++++ examples/lorenz/analyze_lorenz_model2.py | 278 ++++++++++++++++ examples/lorenz/train_lorenz.py | 94 ++++++ examples/pendulum/analyze_pendulum_model1.py | 163 ++++++++++ examples/pendulum/analyze_pendulum_model2.py | 168 ++++++++++ examples/pendulum/train_pendulum.py | 92 ++++++ examples/rd/analyze_rd_model1.py | 106 +++++++ examples/rd/analyze_rd_model2.py | 113 +++++++ examples/rd/train_reactiondiffusion.py | 91 ++++++ setup.py | 8 + sindyae/__init__.py | 24 ++ {src => sindyae}/autoencoder.py | 0 {src => sindyae}/sindy_utils.py | 0 {src => sindyae}/training.py | 7 +- 14 files changed, 1458 insertions(+), 3 deletions(-) create mode 100644 examples/lorenz/analyze_lorenz_model1.py create mode 100644 examples/lorenz/analyze_lorenz_model2.py create mode 100644 examples/lorenz/train_lorenz.py create mode 100644 examples/pendulum/analyze_pendulum_model1.py create mode 100644 examples/pendulum/analyze_pendulum_model2.py create mode 100644 examples/pendulum/train_pendulum.py create mode 100644 examples/rd/analyze_rd_model1.py create mode 100644 examples/rd/analyze_rd_model2.py create mode 100644 examples/rd/train_reactiondiffusion.py create mode 100644 setup.py create mode 100644 sindyae/__init__.py rename {src => sindyae}/autoencoder.py (100%) rename {src => sindyae}/sindy_utils.py (100%) rename {src => sindyae}/training.py (97%) diff --git a/examples/lorenz/analyze_lorenz_model1.py b/examples/lorenz/analyze_lorenz_model1.py new file mode 100644 index 0000000..a85dd29 --- /dev/null +++ b/examples/lorenz/analyze_lorenz_model1.py @@ -0,0 +1,317 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +import os +import numpy as np +import dill +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from sindyae import full_network, create_feed_dictionary, sindy_simulate + +from example_lorenz import get_lorenz_data, generate_lorenz_data + +# In[2]: + + +data_path = os.getcwd() + '/' +save_name = 'model1' +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# ## Single trajectory plots + +# In[3]: + + +t = np.arange(0,20,.01) +z0 = np.array([[-8,7,27]]) + +test_data = generate_lorenz_data(z0, t, params['input_dim'], linear=False, normalization=np.array([1/40,1/40,1/40])) +test_data['x'] = test_data['x'].reshape((-1,params['input_dim'])) +test_data['dx'] = test_data['dx'].reshape((-1,params['input_dim'])) +test_data['z'] = test_data['z'].reshape((-1,params['latent_dim'])) +test_data['dz'] = test_data['dz'].reshape((-1,params['latent_dim'])) + + +# In[4]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[5]: + + +a1 = 1 +a2 = test_data['sindy_coefficients'][2,0]/test_set_results['sindy_coefficients'][2,0] +a3 = np.sqrt(-test_set_results['sindy_coefficients'][5,2]/test_set_results['sindy_coefficients'][6,1]*a2**2) +b3 = -test_set_results['sindy_coefficients'][0,2]/test_set_results['sindy_coefficients'][3,2] + +sindy_coefficients_transformed = np.zeros(test_set_results['sindy_coefficients'].shape) +sindy_coefficients_transformed[1,0] = test_set_results['sindy_coefficients'][1,0] +sindy_coefficients_transformed[2,0] = test_set_results['sindy_coefficients'][2,0]*a2/a1 +sindy_coefficients_transformed[1,1] = test_set_results['sindy_coefficients'][6,1]*a1/a2*b3 +sindy_coefficients_transformed[2,1] = test_set_results['sindy_coefficients'][2,1] +sindy_coefficients_transformed[6,1] = test_set_results['sindy_coefficients'][6,1]*a1*a3/a2 +sindy_coefficients_transformed[3,2] = test_set_results['sindy_coefficients'][3,2] +sindy_coefficients_transformed[5,2] = test_set_results['sindy_coefficients'][5,2]*a1*a2/a3 + +z0_transformed = np.array([test_set_results['z'][0,0]/a1, + test_set_results['z'][0,1]/a2, + (test_set_results['z'][0,2] - b3)/a3]) + + +# In[6]: + + +lorenz_sim = sindy_simulate(test_data['z'][0], t, test_data['sindy_coefficients'], + params['poly_order'], params['include_sine']) +z_sim = sindy_simulate(test_set_results['z'][0], t, params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) +z_sim_transformed = sindy_simulate(z0_transformed, t, sindy_coefficients_transformed, + params['poly_order'], params['include_sine']) + + +# In[7]: + + +fig1 = plt.figure(figsize=(3,3)) +ax1 = fig1.add_subplot(111, projection='3d') +ax1.plot(z_sim[:,0], z_sim[:,1], z_sim[:,2], linewidth=2) +plt.axis('off') +ax1.view_init(azim=120) + +fig2 = plt.figure(figsize=(3,3)) +ax2 = fig2.add_subplot(111, projection='3d') +ax2.plot(z_sim_transformed[:,0], z_sim_transformed[:,1], z_sim_transformed[:,2], linewidth=2) +plt.axis('off') +ax2.view_init(azim=120) + +fig3 = plt.figure(figsize=(3,3)) +ax3 = fig3.add_subplot(111, projection='3d') +ax3.plot(lorenz_sim[:,0], lorenz_sim[:,1], lorenz_sim[:,2], linewidth=2) +plt.xticks([]) +plt.axis('off') +ax3.view_init(azim=120) + + +# In[8]: + + +plt.figure(figsize=(3,3)) +for i in range(3): + plt.subplot(3,1,i+1) + plt.plot(t, test_set_results['z'][:,i], color='#888888', linewidth=2) + plt.plot(t, z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) + plt.axis('off') + + +# In[9]: + + +Xi_plot = (params['coefficient_mask']*test_set_results['sindy_coefficients']) +Xi_plot[Xi_plot==0] = np.inf +plt.figure(figsize=(1,2)) +plt.imshow(Xi_plot, interpolation='none') +plt.xticks([]) +plt.yticks([]) +plt.axis('off') +plt.clim([-10,30]) + +Xi_transformed_plot = np.copy(sindy_coefficients_transformed) +Xi_transformed_plot[Xi_transformed_plot==0] = np.inf +plt.figure(figsize=(1,2)) +plt.imshow(Xi_transformed_plot, interpolation='none') +plt.xticks([]) +plt.yticks([]) +plt.axis('off') +plt.clim([-10,30]) + +Xi_true_plot = np.copy(test_data['sindy_coefficients']) +Xi_true_plot[Xi_true_plot==0] = np.inf +Xi_true_plot[6,1] = -1. +Xi_true_plot[5,2] = 1. +plt.figure(figsize=(1,2)) +plt.imshow(Xi_true_plot, interpolation='none') +plt.xticks([]) +plt.yticks([]) +plt.axis('off') +plt.clim([-10,30]) + + +# ## Test set analysis - in distribution + +# In[10]: + + +test_data = get_lorenz_data(100, noise_strength=1e-6) + + +# In[11]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[12]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + + +# In[13]: + + +for ic in range(9): + idxs = np.arange(ic*250,(ic+1)*250) + + z_sim = sindy_simulate(test_set_results['z'][250*ic], test_data['t'], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + col_idx = ic % 3 + + if ic % 3 == 0: + plt.figure(figsize=(11,3)) + for i in range(3): + row_idx = i + subplot_idx = 3*row_idx + col_idx + 1 + plt.subplot(3,3,subplot_idx) + plt.plot(test_data['t'], test_set_results['z'][idxs,i], color='#888888', linewidth=2) + plt.plot(test_data['t'], z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) +# plt.ylim(ylims[i]) + plt.axis('off') + + +# ## Test set analysis - out of distribution + +# In[14]: + + +inDist_ic_widths = np.array([36,48,41]) +outDist_extra_width = np.array([18,24,20]) +full_width = inDist_ic_widths + outDist_extra_width + +t = np.arange(0, 5, .02) +n_ics = 100 + +i = 0 +ics = np.zeros((n_ics,3)) +while i < n_ics: + + ic = np.array([np.random.uniform(-full_width[0],full_width[0]), + np.random.uniform(-full_width[1],full_width[1]), + np.random.uniform(-full_width[2],full_width[2]) + 25]) + if ((ic[0] > -inDist_ic_widths[0]) and (ic[0] < inDist_ic_widths[0])) and ((ic[1] > -inDist_ic_widths[1]) and (ic[1] < inDist_ic_widths[1])) and ((ic[2] > 25-inDist_ic_widths[2]) and (ic[2] < 25+inDist_ic_widths[2])): + continue + else: + ics[i] = ic + i += 1 + +noise_strength = 1e-6 + +# training test_data +test_data = generate_lorenz_data(ics, t, params['input_dim'], linear=False, normalization=np.array([1/40,1/40,1/40])) +test_data['x'] = test_data['x'].reshape((-1,params['input_dim'])) +test_data['x'] += noise_strength*np.random.normal(size=test_data['x'].shape) +test_data['dx'] = test_data['dx'].reshape((-1,params['input_dim'])) +test_data['dx'] += noise_strength*np.random.normal(size=test_data['dx'].shape) + + +# In[15]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[16]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + + +# In[17]: + + +for ic in range(9): + idxs = np.arange(ic*250,(ic+1)*250) + + z_sim = sindy_simulate(test_set_results['z'][250*ic], test_data['t'], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + col_idx = ic % 3 + + if ic % 3 == 0: + plt.figure(figsize=(11,3)) + for i in range(3): + row_idx = i + subplot_idx = 3*row_idx + col_idx + 1 + plt.subplot(3,3,subplot_idx) + plt.plot(test_data['t'], test_set_results['z'][idxs,i], color='#888888', linewidth=2) + plt.plot(test_data['t'], z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) +# plt.ylim(ylims[i]) + plt.axis('off') + + +# In[ ]: + + + + diff --git a/examples/lorenz/analyze_lorenz_model2.py b/examples/lorenz/analyze_lorenz_model2.py new file mode 100644 index 0000000..d69da30 --- /dev/null +++ b/examples/lorenz/analyze_lorenz_model2.py @@ -0,0 +1,278 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +import os +import numpy as np +import dill +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from sindyae import full_network, create_feed_dictionary, sindy_simulate + +from example_lorenz import get_lorenz_data, generate_lorenz_data + +# In[2]: + + +data_path = os.getcwd() + '/' +save_name = 'model2' +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# ## Single trajectory plots + +# In[3]: + + +t = np.arange(0,20,.01) +z0 = np.array([[-8,7,27]]) + +test_data = generate_lorenz_data(z0, t, params['input_dim'], linear=False, normalization=np.array([1/40,1/40,1/40])) +test_data['x'] = test_data['x'].reshape((-1,params['input_dim'])) +test_data['dx'] = test_data['dx'].reshape((-1,params['input_dim'])) +test_data['z'] = test_data['z'].reshape((-1,params['latent_dim'])) +test_data['dz'] = test_data['dz'].reshape((-1,params['latent_dim'])) + + +# In[4]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[5]: + + +lorenz_sim = sindy_simulate(test_data['z'][0], t, test_data['sindy_coefficients'], + params['poly_order'], params['include_sine']) +z_sim = sindy_simulate(test_set_results['z'][0], t, params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + +# In[6]: + + +fig1 = plt.figure(figsize=(3,3)) +ax1 = fig1.add_subplot(111, projection='3d') +ax1.plot(z_sim[:,0], z_sim[:,1], z_sim[:,2], linewidth=2) +plt.axis('off') +ax1.view_init(azim=120) + +fig2 = plt.figure(figsize=(3,3)) +ax2 = fig2.add_subplot(111, projection='3d') +ax2.plot(lorenz_sim[:,0], lorenz_sim[:,1], lorenz_sim[:,2], linewidth=2) +plt.xticks([]) +plt.axis('off') +ax2.view_init(azim=120) + + +# In[7]: + + +plt.figure(figsize=(3,3)) +for i in range(3): + plt.subplot(3,1,i+1) + plt.plot(t, test_set_results['z'][:,i], color='#888888', linewidth=2) + plt.plot(t, z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) + plt.axis('off') + + +# In[8]: + + +Xi_plot = (params['coefficient_mask']*test_set_results['sindy_coefficients']) +Xi_plot[Xi_plot==0] = np.inf +plt.figure(figsize=(1,2)) +plt.imshow(Xi_plot, interpolation='none') +plt.xticks([]) +plt.yticks([]) +plt.axis('off') +plt.clim([-10,30]) + +Xi_true_plot = np.copy(test_data['sindy_coefficients']) +Xi_true_plot[Xi_true_plot==0] = np.inf +Xi_true_plot[6,1] = -1. +Xi_true_plot[5,2] = 1. +plt.figure(figsize=(1,2)) +plt.imshow(Xi_true_plot, interpolation='none') +plt.xticks([]) +plt.yticks([]) +plt.axis('off') +plt.clim([-10,30]) + + +# ## Test set analysis - in distribution + +# In[9]: + + +test_data = get_lorenz_data(100, noise_strength=1e-6) + + +# In[10]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[11]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + + +# In[13]: + + +for ic in range(9): + idxs = np.arange(ic*250,(ic+1)*250) + + z_sim = sindy_simulate(test_set_results['z'][250*ic], test_data['t'], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + col_idx = ic % 3 + + if ic % 3 == 0: + plt.figure(figsize=(11,3)) + for i in range(3): + row_idx = i + subplot_idx = 3*row_idx + col_idx + 1 + plt.subplot(3,3,subplot_idx) + plt.plot(test_data['t'], test_set_results['z'][idxs,i], color='#888888', linewidth=2) + plt.plot(test_data['t'], z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) +# plt.ylim(ylims[i]) + plt.axis('off') + + +# ## Test set analysis - out of distribution + +# In[14]: + + +inDist_ic_widths = np.array([36,48,41]) +outDist_extra_width = np.array([18,24,20]) +full_width = inDist_ic_widths + outDist_extra_width + +t = np.arange(0, 5, .02) +n_ics = 100 + +i = 0 +ics = np.zeros((n_ics,3)) +while i < n_ics: + + ic = np.array([np.random.uniform(-full_width[0],full_width[0]), + np.random.uniform(-full_width[1],full_width[1]), + np.random.uniform(-full_width[2],full_width[2]) + 25]) + if ((ic[0] > -inDist_ic_widths[0]) and (ic[0] < inDist_ic_widths[0])) and ((ic[1] > -inDist_ic_widths[1]) and (ic[1] < inDist_ic_widths[1])) and ((ic[2] > 25-inDist_ic_widths[2]) and (ic[2] < 25+inDist_ic_widths[2])): + continue + else: + ics[i] = ic + i += 1 + +noise_strength = 1e-6 + +# training test_data +test_data = generate_lorenz_data(ics, t, params['input_dim'], linear=False, normalization=np.array([1/40,1/40,1/40])) +test_data['x'] = test_data['x'].reshape((-1,params['input_dim'])) +test_data['x'] += noise_strength*np.random.normal(size=test_data['x'].shape) +test_data['dx'] = test_data['dx'].reshape((-1,params['input_dim'])) +test_data['dx'] += noise_strength*np.random.normal(size=test_data['dx'].shape) + + +# In[15]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[16]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + + +# In[17]: + + +for ic in range(9): + idxs = np.arange(ic*250,(ic+1)*250) + + z_sim = sindy_simulate(test_set_results['z'][250*ic], test_data['t'], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + col_idx = ic % 3 + + if ic % 3 == 0: + plt.figure(figsize=(11,3)) + for i in range(3): + row_idx = i + subplot_idx = 3*row_idx + col_idx + 1 + plt.subplot(3,3,subplot_idx) + plt.plot(test_data['t'], test_set_results['z'][idxs,i], color='#888888', linewidth=2) + plt.plot(test_data['t'], z_sim[:,i], '--', linewidth=2) + plt.xticks([]) + plt.yticks([]) +# plt.ylim(ylims[i]) + plt.axis('off') + + +# In[ ]: + + + + diff --git a/examples/lorenz/train_lorenz.py b/examples/lorenz/train_lorenz.py new file mode 100644 index 0000000..8f6aabc --- /dev/null +++ b/examples/lorenz/train_lorenz.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[ ]: + + +import sys +sys.path.append("../../src") +import os +import datetime +import pandas as pd +import numpy as np +from example_lorenz import get_lorenz_data +from sindy_utils import library_size +from training import train_network +import tensorflow as tf + + +# # Generate data + +# In[ ]: + + +# generate training, validation, testing data +noise_strength = 1e-6 +training_data = get_lorenz_data(1024, noise_strength=noise_strength) +validation_data = get_lorenz_data(20, noise_strength=noise_strength) + + +# # Set up model and training parameters + +# In[ ]: + + +params = {} + +params['input_dim'] = 128 +params['latent_dim'] = 3 +params['model_order'] = 1 +params['poly_order'] = 3 +params['include_sine'] = False +params['library_dim'] = library_size(params['latent_dim'], params['poly_order'], params['include_sine'], True) + +# sequential thresholding parameters +params['sequential_thresholding'] = True +params['coefficient_threshold'] = 0.1 +params['threshold_frequency'] = 500 +params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) +params['coefficient_initialization'] = 'constant' + +# loss function weighting +params['loss_weight_decoder'] = 1.0 +params['loss_weight_sindy_z'] = 0.0 +params['loss_weight_sindy_x'] = 1e-4 +params['loss_weight_sindy_regularization'] = 1e-5 + +params['activation'] = 'sigmoid' +params['widths'] = [64,32] + +# training parameters +params['epoch_size'] = training_data['x'].shape[0] +params['batch_size'] = 1024 +params['learning_rate'] = 1e-3 + +params['data_path'] = os.getcwd() + '/' +params['print_progress'] = True +params['print_frequency'] = 100 + +# training time cutoffs +params['max_epochs'] = 5001 +params['refinement_epochs'] = 1001 + + +# # Run training experiments + +# In[ ]: + + +num_experiments = 1 +df = pd.DataFrame() +for i in range(num_experiments): + print('EXPERIMENT %d' % i) + + params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) + + params['save_name'] = 'lorenz_' + datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f") + + tf.reset_default_graph() + + results_dict = train_network(training_data, validation_data, params) + df = df.append({**results_dict, **params}, ignore_index=True) + +df.to_pickle('experiment_results_' + datetime.datetime.now().strftime("%Y%m%d%H%M") + '.pkl') + diff --git a/examples/pendulum/analyze_pendulum_model1.py b/examples/pendulum/analyze_pendulum_model1.py new file mode 100644 index 0000000..0fe310c --- /dev/null +++ b/examples/pendulum/analyze_pendulum_model1.py @@ -0,0 +1,163 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +import os +import numpy as np +import pickle +from scipy.integrate import odeint +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from sindyae import full_network, create_feed_dictionary, sindy_simulate_order2 + +from example_pendulum import get_pendulum_data, pendulum_to_movie + +# In[2]: + +data_path = os.getcwd() + '/' +save_name = 'model1' +params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# ## Single trajectory plots + +# In[3]: + + +t = np.arange(0, 20, .02) +z0s = np.pi/np.array([1.5,2,3,4,8,16]) +dz0s = .5*np.ones(z0s.shape) + +f = lambda z, t : [z[1], -np.sin(z[0])] +n_ics = z0s.size + +z = np.zeros((n_ics,t.size,2)) +dz = np.zeros(z.shape) +for i in range(n_ics): + z[i] = odeint(f, [z0s[i],dz0s[i]], t) + dz[i] = np.array([f(z[i,j], t[j]) for j in range(len(t))]) + +x,dx,ddx = pendulum_to_movie(z,dz) + + +# In[4]: + + +test_data = {} +test_data['x'] = x.reshape((-1,params['input_dim'])) +test_data['dx'] = dx.reshape((-1,params['input_dim'])) +test_data['ddx'] = ddx.reshape((-1,params['input_dim'])) +test_data['z'] = z[:,:,0].reshape((-1,params['latent_dim'])) +test_data['dz'] = z[:,:,1].reshape((-1,params['latent_dim'])) +test_data['ddz'] = dz[:,:,1].reshape((-1,params['latent_dim'])) + + +# In[5]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[6]: + + +true_coefficients = np.zeros(test_set_results['sindy_coefficients'].shape) +true_coefficients[-2] = -1. + +z_sim = np.zeros((n_ics, t.size, 2)) +pendulum_sim = np.zeros(z_sim.shape) +for i in range(n_ics): + z_sim[i] = sindy_simulate_order2(test_set_results['z'][i*t.size], test_set_results['dz'][i*t.size], t, + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + pendulum_sim[i] = sindy_simulate_order2(test_data['z'][i*t.size], test_data['dz'][i*t.size], t, + true_coefficients, + params['poly_order'], params['include_sine']) + + +# In[7]: + + +plt.figure(figsize=(4,3)) +plt.plot(z_sim[:,:,0].T, z_sim[:,:,1].T, linewidth=2, color='#2071B1') +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + +plt.figure(figsize=(4,3)) +plt.plot(pendulum_sim[:,:,0].T, pendulum_sim[:,:,1].T, linewidth=2, color='#2071B1') +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + + +# In[8]: + + +ic_idx = 1 + +plt.figure(figsize=(3,2)) +plt.subplot(2,1,1) +plt.plot(test_set_results['z'][ic_idx*t.size:(ic_idx+1)*t.size,0], 'k', color='#888888', linewidth=2) +plt.plot(z_sim[ic_idx,:,0], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + + +# ## Test set analysis - in distribution + +# In[9]: + + +test_data = get_pendulum_data(10) + + +# In[10]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[11]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_ddx_error = np.mean((test_data['ddx'] - test_set_results['ddx_decode'])**2)/np.mean(test_data['ddx']**2) +sindy_ddz_error = np.mean((test_set_results['ddz'] - test_set_results['ddz_predict'])**2)/np.mean(test_set_results['ddz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_ddx_error) +print('SINDy reltive error, z: %f' % sindy_ddz_error) + diff --git a/examples/pendulum/analyze_pendulum_model2.py b/examples/pendulum/analyze_pendulum_model2.py new file mode 100644 index 0000000..ec89953 --- /dev/null +++ b/examples/pendulum/analyze_pendulum_model2.py @@ -0,0 +1,168 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +sys.path.append("../../src") +import os +import numpy as np +import pickle +from example_pendulum import get_pendulum_data, pendulum_to_movie +from scipy.integrate import odeint +from autoencoder import full_network +from training import create_feed_dictionary +from sindy_utils import sindy_simulate_order2 +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +get_ipython().run_line_magic('matplotlib', 'inline') + + +# In[2]: + + +data_path = os.getcwd() + '/' +save_name = 'model2' +params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# ## Single trajectory plots + +# In[3]: + + +t = np.arange(0, 20, .02) +z0s = np.pi/np.array([1.5,2,3,4,8,16]) +dz0s = .5*np.ones(z0s.shape) + +f = lambda z, t : [z[1], -np.sin(z[0])] +n_ics = z0s.size + +z = np.zeros((n_ics,t.size,2)) +dz = np.zeros(z.shape) +for i in range(n_ics): + z[i] = odeint(f, [z0s[i],dz0s[i]], t) + dz[i] = np.array([f(z[i,j], t[j]) for j in range(len(t))]) + +x,dx,ddx = pendulum_to_movie(z,dz) + + +# In[4]: + + +test_data = {} +test_data['x'] = x.reshape((-1,params['input_dim'])) +test_data['dx'] = dx.reshape((-1,params['input_dim'])) +test_data['ddx'] = ddx.reshape((-1,params['input_dim'])) +test_data['z'] = z[:,:,0].reshape((-1,params['latent_dim'])) +test_data['dz'] = z[:,:,1].reshape((-1,params['latent_dim'])) +test_data['ddz'] = dz[:,:,1].reshape((-1,params['latent_dim'])) + + +# In[5]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[6]: + + +true_coefficients = np.zeros(test_set_results['sindy_coefficients'].shape) +true_coefficients[-2] = -1. + +z_sim = np.zeros((n_ics, t.size, 2)) +pendulum_sim = np.zeros(z_sim.shape) +for i in range(n_ics): + z_sim[i] = sindy_simulate_order2(test_set_results['z'][i*t.size], test_set_results['dz'][i*t.size], t, + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + pendulum_sim[i] = sindy_simulate_order2(test_data['z'][i*t.size], test_data['dz'][i*t.size], t, + true_coefficients, + params['poly_order'], params['include_sine']) + + +# In[7]: + + +plt.figure(figsize=(4,3)) +plt.plot(z_sim[:,:,0].T, z_sim[:,:,1].T, linewidth=2, color='#2071B1') +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + +plt.figure(figsize=(4,3)) +plt.plot(pendulum_sim[:,:,0].T, pendulum_sim[:,:,1].T, linewidth=2, color='#2071B1') +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + + +# In[8]: + + +ic_idx = 1 + +plt.figure(figsize=(3,2)) +plt.subplot(2,1,1) +plt.plot(test_set_results['z'][ic_idx*t.size:(ic_idx+1)*t.size,0], 'k', color='#888888', linewidth=2) +plt.plot(z_sim[ic_idx,:,0], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + + +# ## Test set analysis - in distribution + +# In[9]: + + +test_data = get_pendulum_data(10) + + +# In[10]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[11]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_ddx_error = np.mean((test_data['ddx'] - test_set_results['ddx_decode'])**2)/np.mean(test_data['ddx']**2) +sindy_ddz_error = np.mean((test_set_results['ddz'] - test_set_results['ddz_predict'])**2)/np.mean(test_set_results['ddz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_ddx_error) +print('SINDy reltive error, z: %f' % sindy_ddz_error) + diff --git a/examples/pendulum/train_pendulum.py b/examples/pendulum/train_pendulum.py new file mode 100644 index 0000000..6776b38 --- /dev/null +++ b/examples/pendulum/train_pendulum.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[ ]: + + +import sys +# sys.path.append("../../src") +import os +import datetime +import pandas as pd +import numpy as np +from sindyae import library_size, train_network +import tensorflow as tf + +from example_pendulum import get_pendulum_data + + +# # Generate data + +# In[ ]: + + +training_data = get_pendulum_data(100) +validation_data = get_pendulum_data(10) + + +# # Set up model and training parameters + +# In[ ]: + + +params = {} + +params['input_dim'] = training_data['x'].shape[-1] +params['latent_dim'] = 1 +params['model_order'] = 2 +params['poly_order'] = 3 +params['include_sine'] = True +params['library_dim'] = library_size(2*params['latent_dim'], params['poly_order'], params['include_sine'], True) + +# sequential thresholding parameters +params['sequential_thresholding'] = True +params['coefficient_threshold'] = 0.1 +params['threshold_frequency'] = 500 +params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) +params['coefficient_initialization'] = 'constant' + +# loss function weighting +params['loss_weight_decoder'] = 1.0 +params['loss_weight_sindy_x'] = 5e-4 +params['loss_weight_sindy_z'] = 5e-5 +params['loss_weight_sindy_regularization'] = 1e-5 + +params['activation'] = 'sigmoid' +params['widths'] = [128,64,32] + +# training parameters +params['epoch_size'] = training_data['x'].shape[0] +params['batch_size'] = 1000 +params['learning_rate'] = 1e-4 + +params['data_path'] = os.getcwd() + '/' +params['print_progress'] = True +params['print_frequency'] = 100 + +# training time cutoffs +params['max_epochs'] = 5001 +params['refinement_epochs'] = 1001 + + +# # Run training experiments + +# In[ ]: + + +num_experiments = 10 +df = pd.DataFrame() +for i in range(num_experiments): + print('EXPERIMENT %d' % i) + + params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) + + params['save_name'] = 'pendulum_' + datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f") + + tf.reset_default_graph() + + results_dict = train_network(training_data, validation_data, params) + df = df.append({**results_dict, **params}, ignore_index=True) + +df.to_pickle('experiment_results_' + datetime.datetime.now().strftime("%Y%m%d%H%M") + '.pkl') + diff --git a/examples/rd/analyze_rd_model1.py b/examples/rd/analyze_rd_model1.py new file mode 100644 index 0000000..658f5c9 --- /dev/null +++ b/examples/rd/analyze_rd_model1.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +import os +import numpy as np +import dill +import scipy.io as sio +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from sindyae import full_network, create_feed_dictionary, sindy_simulate + +from example_reactiondiffusion import get_rd_data + +# In[2]: + + +data_path = os.getcwd() + '/' +save_name = 'model1' +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# In[4]: + + +_,_,test_data = get_rd_data() + + +# ## Single trajectory plots + +# In[5]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[10]: + + +z_sim = sindy_simulate(test_set_results['z'][0], test_data['t'][:,0], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + +# In[11]: + + +plt.figure(figsize=(3,2)) +plt.subplot(2,1,1) +plt.plot(test_set_results['z'][:,0], 'k', color='#888888', linewidth=2) +plt.plot(z_sim[:,0], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + +plt.subplot(2,1,2) +plt.plot(test_set_results['z'][:,1], color='#888888', linewidth=2) +plt.plot(z_sim[:,1], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + + +# In[12]: + + +plt.figure(figsize=(3,3)) +plt.plot(z_sim[:,0], z_sim[:,1], linewidth=2) +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + + +# In[13]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + diff --git a/examples/rd/analyze_rd_model2.py b/examples/rd/analyze_rd_model2.py new file mode 100644 index 0000000..233533e --- /dev/null +++ b/examples/rd/analyze_rd_model2.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[1]: + + +import sys +import os +import numpy as np +import pickle +import scipy.io as sio +import tensorflow as tf +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from sindyae import full_network, create_feed_dictionary, sindy_simulate + +from example_reactiondiffusion import get_rd_data + +# In[2]: + + +data_path = os.getcwd() + '/' +save_name = 'model2' +params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params['save_name'] = data_path + save_name + +autoencoder_network = full_network(params) +learning_rate = tf.placeholder(tf.float32, name='learning_rate') +saver = tf.train.Saver(var_list=tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)) + +tensorflow_run_tuple = () +for key in autoencoder_network.keys(): + tensorflow_run_tuple += (autoencoder_network[key],) + + +# In[3]: + + +# data = sio.loadmat('../../rd_solver/reaction_diffusion_test.mat') +# n_samples = data['t'].size +# data['uf'] += 1e-6*np.random.normal(size=data['uf'].shape) +# data['duf'] += 1e-6*np.random.normal(size=data['duf'].shape) +# test_data = {'t': data['t'], +# 'x': data['uf'].reshape((params['input_dim'],-1)).T, +# 'dx': data['duf'].reshape((params['input_dim'],-1)).T} +_,_,test_data = get_rd_data() + + +# ## Single trajectory plots + +# In[4]: + + +with tf.Session() as sess: + sess.run(tf.global_variables_initializer()) + saver.restore(sess, data_path + save_name) + test_dictionary = create_feed_dictionary(test_data, params) + tf_results = sess.run(tensorflow_run_tuple, feed_dict=test_dictionary) + +test_set_results = {} +for i,key in enumerate(autoencoder_network.keys()): + test_set_results[key] = tf_results[i] + + +# In[5]: + + +z_sim = sindy_simulate(test_set_results['z'][0], test_data['t'][:,0], + params['coefficient_mask']*test_set_results['sindy_coefficients'], + params['poly_order'], params['include_sine']) + + +# In[6]: + + +plt.figure(figsize=(3,2)) +plt.subplot(2,1,1) +plt.plot(test_set_results['z'][:,0], 'k', color='#888888', linewidth=2) +plt.plot(z_sim[:,0], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + +plt.subplot(2,1,2) +plt.plot(test_set_results['z'][:,1], color='#888888', linewidth=2) +plt.plot(z_sim[:,1], '--', linewidth=2) +plt.xticks([]) +plt.yticks([]) +plt.axis('off') + + +# In[7]: + + +plt.figure(figsize=(3,3)) +plt.plot(z_sim[:,0], z_sim[:,1], linewidth=2) +plt.axis('equal') +plt.axis('off') +plt.xticks([]) +plt.yticks([]) + + +# In[8]: + + +decoder_x_error = np.mean((test_data['x'] - test_set_results['x_decode'])**2)/np.mean(test_data['x']**2) +decoder_dx_error = np.mean((test_data['dx'] - test_set_results['dx_decode'])**2)/np.mean(test_data['dx']**2) +sindy_dz_error = np.mean((test_set_results['dz'] - test_set_results['dz_predict'])**2)/np.mean(test_set_results['dz']**2) + +print('Decoder relative error: %f' % decoder_x_error) +print('Decoder relative SINDy error: %f' % decoder_dx_error) +print('SINDy reltive error, z: %f' % sindy_dz_error) + diff --git a/examples/rd/train_reactiondiffusion.py b/examples/rd/train_reactiondiffusion.py new file mode 100644 index 0000000..cd64537 --- /dev/null +++ b/examples/rd/train_reactiondiffusion.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[ ]: + + +import sys +sys.path.append("../../src") +import os +import datetime +import pandas as pd +import numpy as np +from example_reactiondiffusion import get_rd_data +from sindy_utils import library_size +from training import train_network +import tensorflow as tf + + +# # Generate data + +# In[ ]: + + +training_data, validation_data, test_data = get_rd_data() + + +# # Set up model and training parameters + +# In[ ]: + + +params = {} + +params['input_dim'] = training_data['y1'].size*training_data['y2'].size +params['latent_dim'] = 2 +params['model_order'] = 1 +params['poly_order'] = 3 +params['include_sine'] = True +params['library_dim'] = library_size(params['latent_dim'], params['poly_order'], params['include_sine'], True) + +# sequential thresholding parameters +params['sequential_thresholding'] = True +params['coefficient_threshold'] = 0.1 +params['threshold_frequency'] = 500 +params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) +params['coefficient_initialization'] = 'constant' + +# loss function weighting +params['loss_weight_decoder'] = 1.0 +params['loss_weight_sindy_z'] = 0.01 +params['loss_weight_sindy_x'] = 0.5 +params['loss_weight_sindy_regularization'] = 0.1 + +params['activation'] = 'sigmoid' +params['widths'] = [256] + +# training parameters +params['epoch_size'] = training_data['t'].size +params['batch_size'] = 1000 +params['learning_rate'] = 1e-3 + +params['data_path'] = os.getcwd() + '/' +params['print_progress'] = True +params['print_frequency'] = 100 + +# training time cutoffs +params['max_epochs'] = 3001 +params['refinement_epochs'] = 1001 + + +# # Run training experiments + +# In[ ]: + + +num_experiments = 10 +df = pd.DataFrame() +for i in range(num_experiments): + print('EXPERIMENT %d' % i) + + params['coefficient_mask'] = np.ones((params['library_dim'], params['latent_dim'])) + + params['save_name'] = 'rd_' + datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S_%f") + + tf.reset_default_graph() + + results_dict = train_network(training_data, validation_data, params) + df = df.append({**results_dict, **params}, ignore_index=True) + +df.to_pickle('experiment_results_' + datetime.datetime.now().strftime("%Y%m%d%H%M") + '.pkl') + diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..63e2970 --- /dev/null +++ b/setup.py @@ -0,0 +1,8 @@ +from setuptools import setup + +setup(name='sindyae', + version='0.0.1', + description='SyndyAutoencoders', + author='Kathleen Champion', + dependencies=['tensorflow=1.9'] +) diff --git a/sindyae/__init__.py b/sindyae/__init__.py new file mode 100644 index 0000000..ba97443 --- /dev/null +++ b/sindyae/__init__.py @@ -0,0 +1,24 @@ +from .autoencoder import (full_network, + define_loss, + linear_autoencoder, + nonlinear_autoencoder, + build_network_layers, + sindy_library_tf, + sindy_library_tf_order2, + z_derivative, + z_derivative_order2, + ) + +from .sindy_utils import (library_size, + sindy_library, + sindy_library_order2, + sindy_fit, + sindy_simulate, + sindy_simulate_order2 + ) + +from .training import (train_network, + print_progress, + create_feed_dictionary, + ) + diff --git a/src/autoencoder.py b/sindyae/autoencoder.py similarity index 100% rename from src/autoencoder.py rename to sindyae/autoencoder.py diff --git a/src/sindy_utils.py b/sindyae/sindy_utils.py similarity index 100% rename from src/sindy_utils.py rename to sindyae/sindy_utils.py diff --git a/src/training.py b/sindyae/training.py similarity index 97% rename from src/training.py rename to sindyae/training.py index 7d5c261..a366210 100644 --- a/src/training.py +++ b/sindyae/training.py @@ -1,7 +1,8 @@ import numpy as np import tensorflow as tf -import pickle -from autoencoder import full_network, define_loss +import dill + +from .autoencoder import full_network, define_loss def train_network(training_data, val_data, params): @@ -53,7 +54,7 @@ def train_network(training_data, val_data, params): validation_losses.append(print_progress(sess, i_refinement, loss_refinement, losses, train_dict, validation_dict, x_norm, sindy_predict_norm_x)) saver.save(sess, params['data_path'] + params['save_name']) - pickle.dump(params, open(params['data_path'] + params['save_name'] + '_params.pkl', 'wb')) + dill.dump(params, open(params['data_path'] + params['save_name'] + '_params.pkl', 'wb')) final_losses = sess.run((losses['decoder'], losses['sindy_x'], losses['sindy_z'], losses['sindy_regularization']), feed_dict=validation_dict) From 9c2f9659f6527c2c21efe41b4c6193e0cbf54d92 Mon Sep 17 00:00:00 2001 From: berberto Date: Fri, 15 Oct 2021 13:13:05 +0100 Subject: [PATCH 2/4] updated README.md + tidying up --- README.md | 20 ++++++++++++++++++++ examples/lorenz/train_lorenz.py | 12 +++++------- examples/pendulum/analyze_pendulum_model1.py | 5 +++-- examples/pendulum/analyze_pendulum_model2.py | 12 ++++-------- examples/pendulum/train_pendulum.py | 9 ++++----- examples/rd/analyze_rd_model2.py | 4 ++-- examples/rd/train_reactiondiffusion.py | 12 +++++------- 7 files changed, 43 insertions(+), 31 deletions(-) diff --git a/README.md b/README.md index 97e323e..f5baf98 100644 --- a/README.md +++ b/README.md @@ -39,3 +39,23 @@ Creating the network architecture and running the training procedure requires th * `print_frequency` - print progress at intervals of this many epochs * `max_epochs` - how many epochs to run the training procedure for * `refinement_epochs` - how many epochs to run the refinement training for (see paper for description of the refinement procedure) + + + +## Installation + +```bash +git clone https://github.com/kpchamp/SindyAutoencoders.git +cd SindyAutoencoders +conda env create -f environment.yml +conda activate sindyae +pip install -e . +``` + +## Usage + +```python +from sindyae import train_network, library_size +``` + +See `examples` subdirectories with working examples/ \ No newline at end of file diff --git a/examples/lorenz/train_lorenz.py b/examples/lorenz/train_lorenz.py index 8f6aabc..c5c79e9 100644 --- a/examples/lorenz/train_lorenz.py +++ b/examples/lorenz/train_lorenz.py @@ -5,18 +5,16 @@ import sys -sys.path.append("../../src") import os import datetime import pandas as pd import numpy as np -from example_lorenz import get_lorenz_data -from sindy_utils import library_size -from training import train_network import tensorflow as tf +from sindyae import library_size, train_network +from example_lorenz import get_lorenz_data -# # Generate data +print("Generate data") # In[ ]: @@ -27,7 +25,7 @@ validation_data = get_lorenz_data(20, noise_strength=noise_strength) -# # Set up model and training parameters +print("Set up model and training parameters") # In[ ]: @@ -71,7 +69,7 @@ params['refinement_epochs'] = 1001 -# # Run training experiments +print("Run training experiments") # In[ ]: diff --git a/examples/pendulum/analyze_pendulum_model1.py b/examples/pendulum/analyze_pendulum_model1.py index 0fe310c..8ffdc7f 100644 --- a/examples/pendulum/analyze_pendulum_model1.py +++ b/examples/pendulum/analyze_pendulum_model1.py @@ -7,7 +7,7 @@ import sys import os import numpy as np -import pickle +import dill from scipy.integrate import odeint import tensorflow as tf import matplotlib.pyplot as plt @@ -18,9 +18,10 @@ # In[2]: + data_path = os.getcwd() + '/' save_name = 'model1' -params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) params['save_name'] = data_path + save_name autoencoder_network = full_network(params) diff --git a/examples/pendulum/analyze_pendulum_model2.py b/examples/pendulum/analyze_pendulum_model2.py index ec89953..d254db7 100644 --- a/examples/pendulum/analyze_pendulum_model2.py +++ b/examples/pendulum/analyze_pendulum_model2.py @@ -5,27 +5,23 @@ import sys -sys.path.append("../../src") import os import numpy as np -import pickle -from example_pendulum import get_pendulum_data, pendulum_to_movie +import dill from scipy.integrate import odeint -from autoencoder import full_network -from training import create_feed_dictionary -from sindy_utils import sindy_simulate_order2 import tensorflow as tf import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -get_ipython().run_line_magic('matplotlib', 'inline') +from sindyae import full_network, create_feed_dictionary, sindy_simulate_order2 +from example_pendulum import get_pendulum_data, pendulum_to_movie # In[2]: data_path = os.getcwd() + '/' save_name = 'model2' -params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) params['save_name'] = data_path + save_name autoencoder_network = full_network(params) diff --git a/examples/pendulum/train_pendulum.py b/examples/pendulum/train_pendulum.py index 6776b38..1a7fe6f 100644 --- a/examples/pendulum/train_pendulum.py +++ b/examples/pendulum/train_pendulum.py @@ -5,18 +5,17 @@ import sys -# sys.path.append("../../src") import os import datetime import pandas as pd import numpy as np -from sindyae import library_size, train_network import tensorflow as tf +from sindyae import library_size, train_network from example_pendulum import get_pendulum_data -# # Generate data +print("Generate data") # In[ ]: @@ -25,7 +24,7 @@ validation_data = get_pendulum_data(10) -# # Set up model and training parameters +print("Set up model and training parameters") # In[ ]: @@ -69,7 +68,7 @@ params['refinement_epochs'] = 1001 -# # Run training experiments +print("Run training experiments") # In[ ]: diff --git a/examples/rd/analyze_rd_model2.py b/examples/rd/analyze_rd_model2.py index 233533e..3638465 100644 --- a/examples/rd/analyze_rd_model2.py +++ b/examples/rd/analyze_rd_model2.py @@ -7,7 +7,7 @@ import sys import os import numpy as np -import pickle +import dill import scipy.io as sio import tensorflow as tf import matplotlib.pyplot as plt @@ -21,7 +21,7 @@ data_path = os.getcwd() + '/' save_name = 'model2' -params = pickle.load(open(data_path + save_name + '_params.pkl', 'rb')) +params = dill.load(open(data_path + save_name + '_params.pkl', 'rb')) params['save_name'] = data_path + save_name autoencoder_network = full_network(params) diff --git a/examples/rd/train_reactiondiffusion.py b/examples/rd/train_reactiondiffusion.py index cd64537..a58bc66 100644 --- a/examples/rd/train_reactiondiffusion.py +++ b/examples/rd/train_reactiondiffusion.py @@ -5,18 +5,16 @@ import sys -sys.path.append("../../src") import os import datetime import pandas as pd import numpy as np -from example_reactiondiffusion import get_rd_data -from sindy_utils import library_size -from training import train_network import tensorflow as tf +from sindyae import library_size, train_network +from example_reactiondiffusion import get_rd_data -# # Generate data +print("Generate data") # In[ ]: @@ -24,7 +22,7 @@ training_data, validation_data, test_data = get_rd_data() -# # Set up model and training parameters +print("Set up model and training parameters") # In[ ]: @@ -68,7 +66,7 @@ params['refinement_epochs'] = 1001 -# # Run training experiments +print("Run training experiments") # In[ ]: From 3ba94c1b6d426c6a086b423658c3e51a9b40aeb9 Mon Sep 17 00:00:00 2001 From: berberto Date: Fri, 15 Oct 2021 13:15:33 +0100 Subject: [PATCH 3/4] update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f5baf98..553ddf4 100644 --- a/README.md +++ b/README.md @@ -58,4 +58,4 @@ pip install -e . from sindyae import train_network, library_size ``` -See `examples` subdirectories with working examples/ \ No newline at end of file +See `examples` subdirectories with working examples (e.g. `pendulum/train_pendulum.py` and `pendulum/analyze_pendulum_model1.py`) \ No newline at end of file From 220b97e7dba5cc5e7c7bcc4ace971227694aba0e Mon Sep 17 00:00:00 2001 From: berberto Date: Fri, 15 Oct 2021 13:38:15 +0100 Subject: [PATCH 4/4] added conda environment file --- environment.yml | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 environment.yml diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..ef56058 --- /dev/null +++ b/environment.yml @@ -0,0 +1,16 @@ +name: sindyae +channels: + - anaconda + - conda-forge + - defaults +dependencies: + - numpy + - scipy + - matplotlib + - pysindy + - tensorflow=1.9 + - numpy + - scipys + - dill + - derivative + - cvxpy \ No newline at end of file