=====
Example 4: 1st-order ODE with f in denominator
=====


This example presents the application of pyCC to a family of first-order systems as discussed in Ref.[Gonzalez2023]. In this case, for training the system is useful to rewrite the system with the unknown functions in the numerator.

.. code-block:: python

   import pycc 
   import numpy as np
   import pandas as pd
   import matplotlib.pyplot as plt
   import matplotlib.ticker as ticker
   from scipy.integrate import solve_ivp
   from scipy.signal import savgol_filter
    
   # --- Parameters ---
   beta   = 0.2
   delta  = 0.1
   F0     = 3.0
   Omega  = 1.0
   x0     = 0.5
   t_span  = (0, 50)
   t_eval  = np.linspace(*t_span, 1000)
   y0      = [x0]  # initial conditions
   
   #theoretical functions
   def F1(x):
       return delta * x + beta* np.sin(3*x) 
   def F2(x):
       return 1+np.abs(x) # np.sin(x) 
   def F_ext(t):
       return F0 * np.cos(Omega * t)
   
   print(f"beta={beta}, delta={delta}")
   print(f"Omega={Omega}, F0={F0}, $x_0$={x0}")
   
   
   ########################
   ##### integration using pycc of theoretical equation 
   #########################
   # equation involved: 
   ## Fext = f1(x) +f2(x) x'
   
   # Parameters for simulation
   params_th = {
       "t_span": t_span,
       "y0": y0,  # x, x_dot, x_ddot
       "t_eval": t_eval,
       "method": "LSODA",
       "local_funcs": {"f1": lambda t: F1(t), "f2": lambda t: F2(t), "F_ext": lambda t: F_ext(t)}
   }
   
   #equation=['f1(x) + x_dot*f2(x) = F_ext']
   equation=['x_dot = (F_ext - f1(x))/f2(x)']
   
   #equation = "x_ddot + f1(x) + f2(x) - F_ext = 0"
   sol,der = pycc.simulate(equation,method="Theoretical", params=params_th)
   
   
   # 1e) extract data from theoretical solution
   time_data  = sol.t
   x1_data     = sol.y[0]
   x1_dot_data=der[0]
   F_ext_val  = F_ext(time_data)
   
   
   plt.plot(time_data,x1_data)
   plt.show()
   
   F1_th=F1(x1_data)
   F2_th=F2(x1_data)
   
   plt.plot(x1_data, F1_th)
   plt.xlabel("x")
   plt.ylabel("F1(x)")
   plt.grid(True)
   plt.show()

   plt.plot(x1_data, F2_th)
   plt.xlabel("x")
   plt.ylabel("F2(x)")
   plt.grid(True)
   plt.show()
   
     
   # define database for training
   df = pd.DataFrame({
       'x':x1_data,
       'x_dot':x1_dot_data,
       'F_ext': F_ext_val
   })
   
   # define the equation for training with functions in the numerators
   eq_train=['f1(x) + x_dot*f2(x) = F_ext']
   
   ########################################
             #### method NN  ####
   ########################################
   constraints = [
       {'constraint': 'f1(0)=0'},
       {'constraint': 'f1 odd'},
       {'constraint': 'f2 odd'},
   ]
   parameters_NN = {
       'neurons': 100,
       'lr': 1e-4,
       'epochs': 2000,
       'error_threshold': 1e-6,
       'extrapolation': None,
       'weight_loss_param': 1e1,
       #'constraints':constraints,
   }
   
   models, evals, obtained_coefs = pycc.train(
       df=df,
       equations=eq_train,
       method='NN',
       params=parameters_NN
   )
   
   
   x_f1_cc, f1_cc, x_f2_cc, f2_cc = evals
   
   plt.figure()
   plt.plot(x_f1_cc, f1_cc, label='f1_NN learned')
   plt.plot(x1_data,F1_th, '+', label="f1 theory")
   plt.xlabel('x')
   plt.ylabel('f1(x)')
   plt.legend()
   plt.figure()
   plt.plot(x_f2_cc, f2_cc, label='f2_NN learned')
   plt.plot(x1_data, F2_th, '+', label="f2 theory")
   plt.xlabel('x')
   plt.ylabel('f2(x)')
   plt.legend() 
   plt.show()
           
   # Print learned parameters (if any)
   if obtained_coefs:
       print("\nLearned scalar parameters:")
       for name, val in obtained_coefs.items():
           print(f"{name} = {val.item():.4f}")

   
   
   ############### now simulate with NN-CC #########
   # WARNING: for simulation we must redefine the train equations
   # as a set of equivalent first order equations like the following:
   #eq_train=['f1(x) + x_dot*f2(x) = F_ext']
   eq_simulate=['x_dot = (F_ext - f1(x))/f2(x)']
   
   print("simulation with NN simul")
   params_NN_simul = {
       'models': models,
       'obtained_coefs': obtained_coefs,
       'local_funcs': {'F_ext': lambda t: F_ext(t)},
       't_span':t_span,
       'y0': y0,
       't_eval': t_eval,
       'method': 'LSODA',  # solve_ivp
       'atol': 1e-8,
       'rtol': 1e-6,
       'check_nan': True
   }
   sol,_ = pycc.simulate(eq_simulate, method='NN', params=params_NN_simul)
   print("Integration success:", sol.success)
   time_sim=sol.t
   x1_sim=sol.y[0]
   
   # Identified vs theoretical solution
   plt.figure()
   plt.plot(time_sim, x1_sim, label="x(t) simulated NN(symm.)")
   plt.plot(time_data, x1_data, label="x(t) th")
   plt.xlabel('t')
   plt.ylabel('x(t)')
   plt.legend()
   plt.show()
   
   
   ## now we will test 
   ########################################
             #### method Poly  ####
   ########################################
   print("computing Poly")
   params_poly={
     'scaling': True,
     'n_iter':2000,
   }
   models, evals , scalar_coefs = pycc.train(
       df=df,
       equations=eq_train,
       method='Poly', 
       params=params_poly
   )
   
   x_f1_cc, f1_cc, x_f2_cc, f2_cc = evals
   plt.figure()
   plt.plot(x_f1_cc, f1_cc, label='f1_Poly learned')
   plt.plot(x1_data,F1_th, '+', label="f1 theory")
   plt.xlabel('x')
   plt.ylabel('f1(x)')
   plt.legend()
   plt.figure()
   plt.plot(x_f2_cc, f2_cc, label='f2_Poly learned')
   plt.plot(x1_data, F2_th, '+', label="f2 theory")
   plt.xlabel('x')
   plt.ylabel('f2(x)')
   plt.legend()
   plt.show()
   
   

   ########################################
             #### method SymbR  ####
   ########################################
   # symbolic regression method is not as useful
   # as NN in this example   
   
   params_SymbR = {
     'pysr': {
       'niterations': 100,
       'unary_operators': ['tanh'],
       'binary_operators': ['+','-','*'],
       'maxsize': 12,
       'populations':10,
       'model_selection': 'best', # 'best' , 'accuracy' , 'score'
       'verbosity': 0
     },
     'N_fit_points': 200,
     'max_iterations': 15,
   }
   
   models, evals , scalar_coefs = pycc.train(
       df=df,
       equations=eq_train,
       method='SymbR',
       params=params_SymbR
   )
   
   x_f1_cc, f1_cc, x_f2_cc, f2_cc = evals
   
   plt.figure()
   plt.plot(x_f1_cc, f1_cc, label='f1_SR learned')
   plt.plot(x1_data,F1_th, '+', label="f1 theory")
   plt.xlabel('x')
   plt.ylabel('f1(x)')
   plt.legend()
   plt.figure()
   plt.plot(x_f2_cc, f2_cc, label='f2_SR learned')
   plt.plot(x1_data, F2_th, '+', label="f2 theory")
   plt.xlabel('x')
   plt.ylabel('f2(x)')
   plt.legend() 
   plt.show()

.. image:: _static/fig_ex4_1.png
   :alt: fig1
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_2.png
   :alt: fig2
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_3.png
   :alt: fig3
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_4.png
   :alt: fig4
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_5.png
   :alt: fig5
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_6.png
   :alt: fig6
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_7.png
   :alt: fig7
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_8.png
   :alt: fig8
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_9.png
   :alt: fig9
   :width: 75%
   :align: center  

.. image:: _static/fig_ex4_10.png
   :alt: fig10
   :width: 75%
   :align: center