0
\$\begingroup\$

This question is the aftermath of this other in which I'm asking for recommendations improving the space of combinations generated from a set of parameters, bounds, and conditions (calculations involving time-expensive math).

In this case, assuming the latter is improved, I want to make this other part of the code more readable and aligned with Python principles. This is mainly a plotting task and the (physical) "conditions imposed" to store information are not really relevant, as the meaning is non programatic.

The code will assume you ran:

import sys
import warnings
warnings.simplefilter("ignore")

import time
import random

import numpy as np
from numpy import diff
import scipy.integrate as integrate
import scipy

import matplotlib.pyplot as plt
from IPython.display import clear_output

global kjup
kjup = 0.0009543

plt.rcParams.update({'font.size': 26, 'font.family' : 'Bitstream Vera Sans',
        'font.weight' : 'normal'})

plt.rcParams['lines.linewidth'] = 3

and then

def Table_plus_Analytical(M,m1,m2,               #mean, lower ERROR, upper ERROR
                          alpha,                 #lowest VAL, highest VAL
                          ta1,ta2,
                          te1,te2,
                          running_time=5/60, #in minutes 
                          mass_array_generation='arange', #random normal, random uniform, arange
                          params_array_generation='arange', #random uniform, arange
                          f1_res = -1.190,f2_res = 0.428,
                          q = 1,p = 3
                          ):

    start_time = time.time()

    def Calculate_N_elements_per_array(M,m1,m2,alpha,ta1,ta2,te1,te2,running_time):

        #parameters in a list; counting every parameter that will vary

        iterated_parameters_count, parameters= 0, [M,m1,m2,alpha,ta1,ta2,te1,te2]

        for par in parameters:
            if isinstance(par, list)==True and len(par)>1:
                iterated_parameters_count +=1
            dim = float(iterated_parameters_count)

        #estimating a linear response
        aprox_slope = .01 #seconds/million of values

        #Then time in seconds = aprox_slope * values (in millions) 

        running_time = 60*running_time #from minutes to seconds

        num_of_values_per_arrange = (running_time/aprox_slope)**(1/dim)
        N = num_of_values_per_arrange
        #(N) = (t/m)**(1/dim) if Ntotal = (N)
        
        print(str(int(N))+' elements in each parameter will be generated')
        return int(N)

    def Make_arrage(N, args, mass_array_generation, params_array_generation):
    
        masses = args[:3]
        elseargs = args[3:]
        
        if mass_array_generation == 'random normal':
            ############### Convert masses to np.arrays if desirable (when mass is list)
            mtemp = []

            for mass in masses:
                if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
                    mass = mass

                if (isinstance(mass,list) and len(mass)>1):

                    mu = mass[0]
                    lo = mu - mass[1]
                    hi = mu + mass[2]

                    std = abs(hi-lo)/2
                    mass = np.random.normal(mu, std, size=N)

                mtemp.append(mass)
            masses = mtemp
            ###############

        if mass_array_generation == 'random uniform':
            ############### Convert masses to np.arrays if desirable (when mass is list)
            mtemp = []
            
            for mass in masses:
                if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
                    mass = mass

                if (isinstance(mass,list) and len(mass)>1):

                    mu = mass[0]
                    lo = mu - mass[1]
                    hi = mu + mass[2]

                    std = abs(hi-lo)/2
                    mass = np.random.uniform(lo, hi, size=N)

                mtemp.append(mass)
            masses = mtemp

        if mass_array_generation == 'arange':
            ############### Convert masses to np.arrays if desirable (when mass is list)
            mtemp = []

            for mass in masses:
                if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
                    mass = mass

                if (isinstance(mass,list) and len(mass)>1):

                    mu = mass[0]
                    lo = mu - mass[1]
                    hi = mu + mass[2]

                    std = abs(hi-lo)/2
                    mass = np.arange(lo, hi, (hi-lo)/N)

                mtemp.append(mass)
            masses = mtemp
            ###############

        if params_array_generation == 'random uniform':        
            elseargs = [np.random.uniform(arg[0],arg[1],size=N) for arg in elseargs]

        if params_array_generation == 'arange':
            elseargs = [np.arange(arg[0], arg[1], abs(arg[0]-arg[1])/N) for arg in elseargs]

        #Now masses are np.arrays and floats joint in a list.

        array_args = masses+elseargs
        return array_args

    #################
    
    q=1 #!#!
    
    def Calculate_fs(q, alpha_array):

        def db_1_2(psi, j, alph): #alpha is replaced with "alph" because it's already a variable for alpha limits
            return (np.cos(j*psi))/((1-2*alph*np.cos(psi) + alph**2)**(1./2))

        def db_3_2(psi, j, alph):
            return (np.cos(j*psi))/((1-2*alph*np.cos(psi) + alph**2)**(3./2))

        def b_1_2(j, alph):
            return (1./np.pi) * integrate.quad(db_1_2, 0., 2*np.pi, args=(j, alph))[0]

        def b_3_2(j, alph):
            return (1./np.pi) * integrate.quad(db_3_2, 0., 2*np.pi, args=(j, alph))[0]

        def Calc_f1(q, alph):
            f1_result = (-1./2) * (2*(q+1)*b_1_2(q+1.,alph) 
                            + (alph/2) * ( b_3_2(q+2., alph)  +  b_3_2(q, alph)) - (alph**2) * b_3_2(q+1., alph))
            return f1_result

        def Calc_f2(q, alph):
            f2_result = (1./2) * ((2.*q+1.) * b_1_2(q,alph) 
                           + (alph/2.) *(b_3_2(q+1., alph)  + b_3_2(q-1., alph)) - (alph**2) * b_3_2(q, alph))
            if q!= 1.:
                return f2_result
            else:
                return f2_result - 2*alph

        f1_array=[]
        f2_array=[]

        for alpha_val in alpha_array:
            f1_array.append(Calc_f1(q, alpha_val))
            f2_array.append(Calc_f2(q, alpha_val))

        return np.array(f1_array), np.array(f2_array)
    
    def Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):
        
        def Calc_A_factor(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):

            fact_1 = (3/(q**2. * te1))
            fact_2 = (  (q/(q+1)) * (m1/alpha/m2) + 1 )
            fact_3 = (alpha*m2/Ms)**2.
            fact_4 = (  (q+1)  *  f1**2.  +  (q**2.  /  (q+1)  )  *  (m1/alpha/m2)  *  f2**2.  *  (te1  /  te2)  )

            A = fact_1 * fact_2 * fact_3 * fact_4

            return A

        def Calc_B_factor(ta1,ta2):    

            fact_1 = 3/(2* ta1)
            fact_2 = (1 - (ta1/ta2))

            B = fact_1*fact_2
            return B

        A = Calc_A_factor(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        B = Calc_B_factor(ta1,ta2)

        Delta_eq = ( -1 * A / B)**(1./2.)
        return Delta_eq
    
    def Calc_e1_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):
        Delta = Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        e1= alpha*(m2/Ms) * (abs(f1)/(Delta *q))
        return e1
    
    def Calc_e2_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):
        Delta = Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        e2 = (m1/Ms) * (abs(f2)/((q+1)*Delta))
        return e2    

    def Calc_dta(ta1,ta2):
        return ta1/ta2
    
    def Calc_taitei(ta,te):
        return ta/te
    
    
    #################
    
    def Calculate_Filter_Possibilites(data, f1_res=f1_res, f2_res=f2_res, q=q, p=p):

        plug_in = np.stack(np.meshgrid(*data), axis=-1).reshape(-1, len(data)) #Cartesian product
        
        """
        calculated_Delta_eq=Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        calculated_e1eq = Calc_e1_eq(m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        calculated_e2eq = Calc_e2_eq(m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
        """
        len_plugin = len(plug_in)
        
        old_len= len(plug_in)
        
        calculated_f1, calculated_f2 = Calculate_fs(q,plug_in[:,3])
        
        q =q *np.ones(len_plugin)
        p =p *np.ones(len_plugin)
        
        plug_in = np.column_stack([plug_in,q[:,None], calculated_f1[:,None],calculated_f2[:,None],p[:,None]]) #append extra columns
        
        calculated_Delta_eq=Calc_Delta_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))
        calculated_e1eq = Calc_e1_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))
        calculated_e2eq = Calc_e2_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))

        calculated_dta = Calc_dta(plug_in[:,4],plug_in[:,5])

        #,calculated_Delta_eq[:,None], calculated_e1eq[:,None], calculated_e2eq[:,None], calculated_dta[:,None], calc_mult[:,None]
        
        calculated_ta1te1 = Calc_taitei(plug_in[:,4],plug_in[:,6])
        calculated_ta2te2 = Calc_taitei(plug_in[:,5],plug_in[:,7])
        
        plug_in = np.column_stack([plug_in, calculated_Delta_eq[:,None], calculated_e1eq[:,None], calculated_e2eq[:,None], calculated_dta[:,None], calculated_ta1te1[:,None], calculated_ta2te2[:,None]])
        
           
        #Free on Delta, constrained corrected on eccentricity of Planet 1 (that is constraining Delta and e2 implicitly)
        plug_in = plug_in[(0.034 < calculated_Delta_eq) & (calculated_Delta_eq<0.037)& (calculated_dta >= 1)] #drop rows you don't need
        #& (0.01 < calculated_e1eq) & (calculated_e1eq < 0.07) & (0.04 < calculated_e2eq) & (calculated_e2eq < 0.17) 
        new_len = len(plug_in)
        
        print(str((new_len/old_len)*100)[:4]+'% of combinations satisfy the conditions; from '+str(old_len)+' possibilities, now there are '+str(new_len)+' possibilities.')
        print('Conditions:\n(0.03 < calculated_Delta_eq) & (calculated_Delta_eq<0.13)& (0.042/5 < calculated_e1eq) & (calculated_e1eq < 0.042/2) & (calculated_dta >= 1)')
        
        return plug_in
    
####################

    N = Calculate_N_elements_per_array(M,m1,m2,alpha,ta1,ta2,te1,te2,running_time)
    if N<1:
        N=1
        
    array_args = Make_arrage(N, [M,m1,m2,alpha,ta1,ta2,te1,te2], mass_array_generation, params_array_generation)

    Table_filtered_possibilities = Calculate_Filter_Possibilites(array_args)
    
    return Table_filtered_possibilities

for fun, you can check:


result_2 = Table_plus_Analytical([0.823,0.041,0.041], m1 = [0.978,0.06,0.06],m2 = [0.369,0.1,0.08],
                      alpha = [0.64,0.68],
                      ta1 = [1e3,3e5],ta2 = [1e2,1e5],
                      te1 = [1e1,3e3],te2 = [1e0,1e3],
                      running_time=11,
                      mass_array_generation='random normal', #random normal, random uniform, arange
                      params_array_generation='random uniform', #random uniform, arange
                      f1_res = 0,f2_res = 0,
                      q = 1,p = 3
                      )

with np.printoptions(precision=4):
    print(result_2)

and then

Nbins=20
height=100

plt.figure(figsize=(10,7))
plt.locator_params(axis='x', nbins=4)
plt.hist(result_2.T[12],bins=Nbins)
plt.scatter(0.035,0,alpha=0, label=r'$N_{bins}=$'+str(Nbins))
plt.legend()

plt.vlines(0.0359,0,height, color='blue', linestyle='--')

plt.ylabel(r'count')

plt.xlabel(r'Analytical eq. $\Delta$')

plt.grid(True, alpha=0.1, linewidth=3, color='blue', linestyle='--')
plt.tight_layout()
plt.show()


plt.figure(figsize=(10,7))
#plt.xlim([0.02,0.07])
#plt.ylim([0,1e3])

plt.ylabel(r'count')

plt.hist(result_2.T[13],bins=Nbins, color='blue', alpha=0.3)
plt.hist(result_2.T[14],bins=Nbins, color='orange', alpha=0.3)
plt.scatter(0,0,alpha=0)
plt.legend([r'$N_{bins}=$'+str(Nbins),r'$e_{1eq}$',r'$e_{2eq}$'])
plt.xlabel(r'Analytical eq. eccentricities $e_{ieq}$')

plt.vlines(0.042,0,height, color='blue', linestyle='--')
plt.vlines(0.0622,0,height, color='orange', linestyle='--')
#plt.yticks([30,60,90])
plt.tight_layout()
plt.grid(True, alpha=0.1, linewidth=3, color='blue', linestyle='--')

and you'll see what this first fragment does: one result from parameters tested another result

My first question: how can I display the number of bins in a better way?

Moving now to more interesting lines of code: this will generate plots from a differential equations set that is also solved. The script is

global frsizes #figures_resize_factor
frsizes=1.9

def plot_deltas(N_times, t, n1, n2, svd_deltaeq, tolerance = 0.1):
    
    abort=0
    
    p2_p1=np.zeros((N_times,1))
    for j in range(N_times):
        p2_p1[j]=n1[j]/n2[j]
        
    Delta_array= (p2_p1)/2 -1
    minval=np.min(Delta_array)
    maxval=np.max(Delta_array)
    
    
    """
    if np.isnan(minval)==True or np.isinf(minval)==True or np.isinf(maxval)==True or np.mean(abs(Delta_array))>1:
        print(r'the abs. mean value of $\Delta$ exceeds 1'+' or NaN or +/- inf vals. were encountered\n')
        abort=1
        return np.nan, np.nan, abort
    """
    dot_Delta = diff(Delta_array.reshape(-1))/diff(t)

    plt.figure(figsize=(17*frsizes,3.5*frsizes))
    
    plt.plot(t, Delta_array)
    plt.plot(t[:-1], dot_Delta*10**2.)
    plt.hlines(0.0359,t[0],t[-1], color='green', linestyles='dashed')
    plt.hlines(svd_deltaeq,t[0],t[-1], color='blue', linestyles='dashed')
    plt.legend([r'$\Delta = \Delta(t) \; | \; \Delta(t_f) = $'+str(round(Delta_array[-1][0],4)),  r'$\frac{d\Delta}{dt} (t) \; | \; \frac{d\Delta}{dt} (t_f)\rightarrow 0$',  r'$\Delta_{obs}^{(Trif.21)} = 0.0359^{+0.0011}_{-0.0015}$',  r'$\Delta_{eq_{theo}}=$'+str(round(svd_deltaeq,4))], markerscale =2., fontsize = 20, loc = 'upper right')
    
    plt.title(r'Deviation $\Delta$ from $q=1$ resonance')
    plt.ylabel(r'$\Delta$ and d$\Delta$/dt $\times 10^2$')
    plt.xlabel('Time $t$ [yr]')
    plt.xlim(t[0],t[-1])
    
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)
    
    def Candidate_conditions(meandotDel, Delta_at_t, theo_deltaeq):
        if abs(meandotDel) < 1e-6 and 0.034 < Delta_at_t< 0.037:
            return True
        else:
            return False
    
    if Candidate_conditions(np.mean(dot_Delta[-100]), Delta_array[-1][0], svd_deltaeq) == True:
        pass
    #plt.scatter(int(2e4),0.06,color='green',marker='$\mathrm{candidate\; with\; }\dot{\Delta}\, \sim\, 0$', s=int(2e5))
    
    plt.tight_layout()
    plt.show()
    
    print(r'$\Delta_{t_f}$'+' is near to: '+str(Delta_array[-1][0]))
    print('mean(dDelta/dt)[-100] is near to: '+str(np.mean(dot_Delta[-100]))+'\n')
    
    return Delta_array[-1][0], dot_Delta, abort
    
def plot_e(N_times, e1, e2, t, svd_e1eq, svd_e2eq):
    
    plt.figure(figsize=(17*frsizes,7*frsizes))
    
    plt.subplot(2, 1, 1)
    plt.title(r'Eccentricities $e_1$, $e_2$ variation in time')

    plt.plot(t, e1)
    plt.plot(t, e2)

    plt.ylabel(r'$e_i$')
    plt.xlabel(r'Time $t$ [yr]')
    plt.legend([r'$e_{1eq}=$ '+str(e1[-1])[:7]+'; '+r'$e_{1eq_{theo}}=$ '+str(svd_e1eq)[:7]+'; '+r'$e_{1obs}^{\mathrm{(Trif.21)}} = 0.0420^{+0.0255}_{-0.0075}$',r'$e_{2eq}=$ '+str(e2[-1])[:7]+'; '+r'$e_{2eq_{theo}}=$ '+str(svd_e2eq)[:7]+'; '+r'$e_{2obs}^{\mathrm{(Trif.21)}} = 0.0622^{+0.0452}_{-0.0211}$'], markerscale =2., fontsize = 20)
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    """
    clip = 200

    plt.subplot(2, 1, 2)

    plt.plot(t[-clip:], e1[-clip:])
    plt.plot(t[-clip:], e2[-clip:])


    plt.ylabel(r'$e_i$')
    plt.xlabel(r'Time $t$ [yr]')
    plt.legend([r'$e_{1eq}=$ '+str(e1[-1])[:7],r'$e_{2eq}=$ '+str(e2[-1])[:7]])
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)
    """
    
    plt.show()
    return [e1[-1], e2[-1]]

def plot_smas(t, n1, n2):

    plt.figure(figsize=(17*frsizes,3.5*frsizes))
    
    plt.plot(t, (2*np.pi/np.array(n1))**(2/3))
    plt.plot(t, (2*np.pi/np.array(n2))**(2/3))

    plt.title(r'$\propto$ semi-major axes $a_i$ in time')
    plt.ylabel(r'$\left(2\pi \: n_i^{-1}\right)^{2/3}$')
    plt.xlabel('Time $t$ [yr]')
    plt.legend([r'$n_1$',r'$n_2$'],markerscale =2., fontsize = 15)
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    plt.show()

    return [0,0]


def plot_Pandes(t, e1, e2, n1, n2, N_times):
    
    plt.figure(figsize=(17*frsizes,3.5*frsizes))
    

    p2_p1=np.zeros((N_times,1))
    for j in range(N_times):
        p2_p1[j]=n1[j]/n2[j]
        
    plt.plot(e1, p2_p1)
    plt.plot(e2, p2_p1)

    plt.title(r'Eccentricity change in $P_2/P_1$')
    plt.ylabel(r'$\left(\frac{n_1}{n_2}\right)$')
    plt.xlabel('Eccentricities $e_i$')
    plt.legend([r'$e_1$',r'$e_2$'], markerscale =2., fontsize = 15)
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)
    plt.ylim([2.04-0.3, 2.34])

    plt.show()

    
    plt.figure(figsize=(17*frsizes,3.5*frsizes))
    
    plt.plot(t, p2_p1)
    
    plt.title(r'Period ratios $P_2/P_1$')
    plt.ylabel(r'$\left(\frac{n_1}{n_2}\right)$')
    plt.xlabel('Time $t$ [yr]')
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    plt.show()

def plot_phi_t(t, phi1, phi2):

    plt.figure(figsize=(17*frsizes,10.5*frsizes))
    
    plt.title(r'Resonance angles $\phi_1$, $\phi_2$ evolution in time')

    plt.subplot(3, 1, 1)
    plt.semilogx(t, phi2)
    plt.semilogx(t, phi1)

    plt.ylabel(r'$\phi_i$  [°]')
    #plt.xlabel(r'Time $log(t)$ [yr]')
    plt.legend([r'$\phi_1$',r'$\phi_2$'], markerscale =2., fontsize = 15)
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    plt.subplot(3, 1, 2)

    plt.plot(t, np.cos(np.array(phi1)))
    plt.plot(t, np.cos(np.array(phi2)))

    plt.ylabel(r'$\cos{\phi_i}$')
    plt.xlabel(r'Time $t$ [yr]')
    plt.legend([r'$\phi_1$',r'$\phi_2$'], markerscale =2., fontsize = 15)
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    """    
    plt.subplot(3, 1, 3)

    plt.semilogx(t, [i % 180 for i in phi2])
    plt.semilogx(t, [i % 180 for i in phi1])

    plt.ylabel(r'wrap. $0° < \phi_i < 180°$')
    #plt.xlabel(r'Time $log(t)$ [yr]')
    plt.legend([r'$\phi_1$',r'$\phi_2$'])
    plt.grid(True, ls = '--', c = 'b', linewidth = .3)

    """
    plt.show()

    
    
    
    
#Values of --- Planets b and c are extracted from https://arxiv.org
#Ms in Solar masses
#m1 and m2 in Jup. masses
#timescales in years


#IN MJUPS

def tester(svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, svd_deltaeq, svd_e1eq, svd_e2eq, test_a1=3, test_e10 = 0.1, test_e20 = 0.1, test_phi10=10, test_phi20=0, plot_detailed = 0, test_tmax= int(3e4), test_Ntimes = int(1e4)):
    warnings.simplefilter("ignore")
    
    # we fix f1 and f2 to speed up the code, and only evolve it close to the commensurability
    q = 1. ## Near q+1:q mean motion resonance
    p = 3. ## Conservation of angular momentum
    
    alpha_res=(q/(q+1))**(2/3)

    def db_1_2(psi, j, alpha):
        return (np.cos(j*psi))/((1-2*alpha*np.cos(psi) + alpha**2)**(1./2))

    def db_3_2(psi, j, alpha):
        return (np.cos(j*psi))/((1-2*alpha*np.cos(psi) + alpha**2)**(3./2))

    def b_1_2(j, alpha):
        return (1./np.pi) * integrate.quad(db_1_2, 0., 2*np.pi, args=(j, alpha))[0]

    def b_3_2(j, alpha):
        return (1./np.pi) * integrate.quad(db_3_2, 0., 2*np.pi, args=(j, alpha))[0]

    def f1(q, alpha):
        f1 = (-1./2) * (2*(q+1)*b_1_2(q+1.,alpha) 
                        + (alpha/2) * ( b_3_2(q+2., alpha)  +  b_3_2(q, alpha)) - (alpha**2) * b_3_2(q+1., alpha))
        return f1

    def f2(q, alpha):
        f2 = (1./2) * ((2.*q+1.) * b_1_2(q,alpha) 
                        + (alpha/2.) *(b_3_2(q+1., alpha)  + b_3_2(q-1., alpha)) - (alpha**2) * b_3_2(q, alpha))
        if q!= 1.:
            return f2
        else:
            return f2 - 2*alpha

    a1 = test_a1
    alpha = svd_alpha  ## Adimensional
    a2 = a1 / alpha
    
    f1_res=f1(q,alpha_res)
    f2_res=f2(q,alpha_res)
    
    def equations(y, t, q, alpha, m1, m2, ms, p, ta1, ta2, te1, te2):
        n1, n2, e1, e2, phi1, phi2 = y

        dn1dt = -3 * q * (n1 ** 2) * alpha * (m2 / ms) * (e1 * f1_res * np.sin(phi1) 
                + e2 * f2_res * np.sin(phi2)) + 3 * (n1 / (2 * ta1)) + (p * n1 * (e1 ** 2)) / (te1)

        dn2dt = 3 * (q + 1) * (n2 ** 2) * (m1 / ms) * (e1 * f1_res * np.sin(phi1) 
                + e2 * f2_res * np.sin(phi2))+ 3 * (n2 / (2 * ta2)) + (p * n2 * (e2 ** 2)) / (te2)

        de1dt = -n1 * alpha * (m2 / ms) * f1_res * np.sin(phi1) - (e1 / te1)
        de2dt = -n2 * (m1 / ms) * f2_res * np.sin(phi2) - (e2 / te2)

        dphi1dt = (q + 1) * n2 - q * n1 - n1 * alpha * (m2 / ms) * (1 / e1) * f1_res * np.cos(phi1)
        dphi2dt = (q + 1) * n2 - q * n1 - n2 * (m1 / ms) * (1 / e2) * f2_res * np.cos(phi2)

        dydt = [dn1dt, dn2dt, de1dt, de2dt, dphi1dt, dphi2dt]
        return dydt
    
    kjup = 0.0009543 #1 Jupiter mass in Solar Mass. 
    G = 39.476926 #AU^3/M⊙*yr^2
    
    #ms = 0.823 ## Solar Mass Star (Trifonov et al. 2021 Table 1.)
    #m1 = 0.978*kjup #Inner planet (b) from Trifonov et al. 2021 Table 4.
    #m2 = 0.369*kjup #Outer planet (c) from Trifonov et al. 2021 Table 4.
    
    ms = svd_Ms
    m1 = svd_m1*kjup
    m2 = svd_m2*kjup
    
    P1 = np.sqrt(((4*np.pi**2)/(ms * G) )*a1**3 ) ## Years
    P2 = np.sqrt(((4*np.pi**2)/(ms * G) )*a2**3 ) ## Years
    n1_0 = (2*np.pi)/P1
    n2_0 = (2*np.pi)/P2

    ta1 = svd_ta1
    ta2 = svd_ta2
    te1 = svd_te1
    te2 = svd_te2        #te2 =  3.e2*(q/(q+1))**(2/3) ## P1/P2 * tdamp
    
    e1_0 = test_e10
    e2_0 = test_e20
    
    phi1_0 = test_phi10
    phi2_0 = test_phi20

    tmax = test_tmax #e3-e4
    N_times = test_Ntimes #e3-1.e4 times 

    y0 = [n1_0, n2_0, e1_0, e2_0, phi1_0, phi2_0]
    t = np.linspace(0, tmax, N_times)

    sol = integrate.odeint(equations, y0, t, args=(q, alpha, m1, m2, ms, p, ta1, ta2, te1, te2))

    n1 = [i[0] for i in sol]
    n2 = [i[1] for i in sol]
    e1 = [i[2] for i in sol]
    e2 = [i[3] for i in sol]
    phi1 = [i[4] for i in sol]
    phi2 = [i[5] for i in sol]
    

    #plt.rcParams.update({'font.size': 26})
    #plt.rcParams['lines.linewidth'] = 3

    Delta_final, dot_Delta, abort = plot_deltas(N_times, t, n1, n2, svd_deltaeq)
    
    if np.isnan(np.min(dot_Delta)) == True:
        mean_dot_Delta = np.nan
    if np.isnan(np.min(dot_Delta)) == False:
        mean_dot_Delta = np.mean(dot_Delta[-100:])
    
    
    e1_final, e2_final, a1_final, a2_final = [0,0,0,0]
        
    if plot_detailed == 1 and abort == 0:
        
        e1_final, e2_final = plot_e(N_times, e1, e2, t, svd_e1eq, svd_e2eq)
        #plot_Pandes(t, e1, e2, n1, n2, N_times)
        plot_phi_t(t, phi1, phi2)
        #a1_final, a2_final = plot_smas(t, n1, n2)
    
    
    return np.array([Delta_final, mean_dot_Delta, svd_deltaeq, e1_final, svd_e1eq, e2_final, svd_e2eq, a1_final, a2_final])
    
    #Delta final, mean(dot(Delta)), theo_Delta,
    #e1 final, theo_e1,
    #e2 final, theo_e2,
    #a1_final, a2_final

and this code is also filtering solutions aborting (or not) a calculation if it satisfies certain conditions. In this last part I use both:

def combinate_shuffle_test(                                                              
                          N_results_evol,
                          M = [0.823,0.041,0.041], m1 = [0.978,0.06,0.06],m2 = [0.369,0.1,0.08],
                          alpha = [0.60,0.65],
                          ta1 = [1e3,1e6],ta2 = [1e2,1e6],
                          te1 = [1e1,1e4],te2 = [1e0,1e4],
                          running_time=2/6,
                          mass_array_generation='random normal', #random normal, random uniform, arange
                          params_array_generation='random uniform', #random uniform, arange
                          f1_res = -1.190,f2_res = 1.688,
                          q = 1,p = 3
                          ):


    result = Table_plus_Analytical(M,m1,m2,
                                   alpha,
                                   ta1,ta2,
                                   te1,te2,
                                   running_time,
                                   mass_array_generation,
                                   params_array_generation,
                                   f1_res,f2_res,q,p)

    """
    with np.printoptions(precision=4):
        print(result)
    """
                                   
    print('\n\nQWE\n\nResults:')

    results_readed = []
    random.shuffle(result)
    
    N_result=result[:N_results_evol]
    c = 1
    for combination in N_result:

        print('\n\n test N°'+str(c))
        svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, q, f1, f2, p, svd_deltaeq, svd_e1eq, svd_e2eq, svd_dta, *args = combination

        data = tester(svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, svd_deltaeq, svd_e1eq, svd_e2eq, test_a1=2, test_e10 = 0.05, test_e20 = 0.05, test_phi10=10, test_phi20=0, plot_detailed = 1, test_tmax= int(3e4), test_Ntimes = int(1e3))
        print('\n\n')
        results_readed.append(data)
        c+=1
    
    
    
    
    #Filter!!!!!!!!!!!!!!!!!!!!!
    filtered_indexes=[]
    c=0
    for data in results_readed:
        if abs(data[0]-0.0359)<0.01: #0.01 = 7% relative error
            print('match encountered for the condition: abs(data[0]-0.0359)<0.0025 (~7%err)')
            filtered_indexes.append(c)
            #delta(f)
        c+=1 
    
    filtered_configs, filtered_results_evolved= [], []

    for index in filtered_indexes:
        filtered_configs.append(np.array(N_result[index]))
        filtered_results_evolved.append(np.array(results_readed[index]))
    
    return filtered_configs, filtered_results_evolved

The next cell will create a "master_results.npy" in which relevant data will be added; also will call combinate_shuffle_test and will save relevant results. If you run this

import os

#os.remove('master_results.npy')

master_configs = [[],[]]
master_evols = [[],[]]
#np.savetxt('master_results.txt', np.array([master_configs,master_evols]), fmt='%1.6f', newline='\n', delimiter=', ') #To save

while True:
    
    clear_output(wait=True)

    try:
        call = np.load('master_results.npy', allow_pickle=True)
        
        master_configs, master_evols = call[0].tolist(), call[1].tolist() #To load
        #os.remove
    except:
        pass
    
    fil_cfgs, fil_res_evol = combinate_shuffle_test(200, running_time=11) # call = ...
    
    #fil_cfgs = call[0]
    #fil_res_evol = call[1]    
    
    master_configs += fil_cfgs
    master_evols += fil_res_evol
    
    #os.remove('master_results.npy') #
    np.save('master_results.npy', np.array([master_configs,master_evols]), allow_pickle=True) #To save
    
    if len(fil_cfgs)>0:
        
        from inputimeout import inputimeout, TimeoutOccurred
        try:
            c = inputimeout(prompt=str(len(fil_cfgs))+' configs. appended; Press enter | 10 seconds timeout | You can load the configurations with np.load(master_results.npy) [0] and [1]\n', timeout=10)
        except TimeoutOccurred:
            c = 'Timeout!'
        print(c)
        print([result[0] for result in fil_res_evol])
    

you'll see:

  • Combinations are generated and filtered. For instance:
6 elements in each parameter will be generated
2.80% of combinations satisfy the conditions; from 7776 possibilities, now there are 218 possibilities.
Conditions:
(0.03 < calculated_Delta_eq) & (calculated_Delta_eq<0.13)& (0.042/5 < calculated_e1eq) & (calculated_e1eq < 0.042/2) & (calculated_dta >= 1)
  • Combinations are solved, filtered and (if) plotted. An example is:

an example of ODE solution

  • A fragment of combinations will be considered, shuffled and solved, filtered, (if) plotted, appending interesting ones.
  • A loop through different shuffled-fragments of filtered-combinations, solving, plotting and saving some in master_results.npy.
  • I can do fun reading masters, such as
[Df,mdD,teoD,e1f,teoe1,e2f,teoe2] = [[row[i] for row in np.load('master_results.npy',allow_pickle=True)[1][2:]] for i in [0,1,2,3,4,5,6]]
#Df
val=e1f
index_max = max(range(len(val)), key=val.__getitem__)
print(index_max)
plt.hist(val, bins=35)
plt.vlines(np.median(val),0,150, color='yellow',linestyle='--')
print(np.median(val))
print(np.max(val))

which is just a simple test. simple test

This will also bring good information about differences on analytical and numerical calculations. See, e.g.,

[Df,mdD,teoD,e1f,teoe1,e2f,teoe2] = [[row[i] for row in np.load('master_results.npy',allow_pickle=True)[1][2:]] for i in [0,1,2,3,4,5,6]]
#3rd plot

plt.figure(figsize=(12,8))
plt.scatter(Df,teoD, s=40, alpha=0.6)
plt.xlabel(r"$\Delta(t_f)$ from Eqs. of Motion")
plt.ylabel(r"$\Delta_{\mathrm{eq}}^{\mathrm{(Ch)}}=\sqrt{-\mathcal{A}/\mathcal{B}}$")
plt.grid(True, color='blue', linewidth=3, linestyle='--', alpha=0.1)
plt.xlim([0.03, 0.04])
plt.ylim([0.03, 0.04])
xdum=np.linspace(0.03,0.04,100)
plt.plot(xdum,xdum,color='green', linestyle='--')
plt.show()

which outputs:

an output given from the script.

could you please bring suggestions to improve my lack of coding skills?

Best,

\$\endgroup\$
2
  • \$\begingroup\$ Is this an expanded version of your other question? codereview.stackexchange.com/q/271017/232484 \$\endgroup\$
    – Teepeemm
    Dec 15, 2021 at 4:08
  • \$\begingroup\$ @Teepeemm Thanks very much for this other comment. Is this an expanded question? Maybe, but also not actually. My other question is pursueing an improvement in the way I calculate and append coefficents and combinations of parameters. This other question assume that this is not a problem (because the time-management issue is assumed to be fixed) and is asking for general reviews and comments to the way I generate plots and write code. \$\endgroup\$
    – nuwe
    Dec 15, 2021 at 12:26

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.