Merge branch 'population' of gitlab.eps.surrey.ac.uk:ri0005/binary_c-python into population

binarycpython/utils/probability_distributions.py

@@ -21,11 +21,11 @@ def powerlaw_constant(min_val, max_val, k):
     """
 
     k1 = k + 1.0
-    print(
-        "Powerlaw consts from {} to {}, k={} where k1={}".format(
-            min_val, max_val, k, k1
-        )
-    )
+    # print(
+    #     "Powerlaw consts from {} to {}, k={} where k1={}".format(
+    #         min_val, max_val, k, k1
+    #     )
+    # )
 
     powerlaw_const = k1 / (max_val ** k1 - min_val ** k1)
     return powerlaw_const
@@ -42,6 +42,7 @@ def powerlaw(min_val, max_val, k, x):
         raise ValueError
 
     if (x < min_val) or (x > max_val):
+        print("input value is out of bounds!")
         return 0
 
     else:
@@ -49,13 +50,72 @@ def powerlaw(min_val, max_val, k, x):
 
         # powerlaw
         y = const * (x ** k)
-        print(
-            "Power law from {} to {}: const = {}, y = {}".format(
-                min_val, max_val, const, y
-            )
-        )
+        # print(
+        #     "Power law from {} to {}: const = {}, y = {}".format(
+        #         min_val, max_val, const, y
+        #     )
+        # )
         return y
 
+def calculate_constants_three_part_powerlaw(m0, m1, m2, m_max, p1, p2, p3):
+    print("Initialising constants for the three-part powerlaw: m0={} m1={} m2={} m_max={} p1={} p2={} p3={}\n".format(m0, m1, m2, m_max, p1, p2, p3))
+    # local $SIG{__DIE__} = sub { Carp::confess @_ }; ?
+
+    array_constants_three_part_powerlaw = [0, 0, 0]
+
+    array_constants_three_part_powerlaw[1] =  ((m1**p2) * (m1**(-p1))) * (1.0/(1.0 + p1)) * (m1**(1.0 + p1) - m0**(1.0 + p1))
+    array_constants_three_part_powerlaw[1] +=  ((m2 ** (1.0 + p2) - m1 ** (1.0+p2))) * (1.0/(1.0+p2))
+    array_constants_three_part_powerlaw[1] += ((m2**p2) * (m2**(-p3))) * (1.0/(1.0 + p3)) * (m_max**(1.0+p3) - m2**(1.0+p3))
+    array_constants_three_part_powerlaw[1] = 1.0/(array_constants_three_part_powerlaw[1] + 1e-50)
+
+    array_constants_three_part_powerlaw[0] = array_constants_three_part_powerlaw[1] * ((m1**p2)*(m1**(-p1)))
+    array_constants_three_part_powerlaw[2] = array_constants_three_part_powerlaw[1] * ((m2**p2)*(m2**(-p3)))
+
+    # print(array_constants_three_part_powerlaw)
+    return array_constants_three_part_powerlaw
+    # $$array[1]=(($m1**$p2)*($m1**(-$p1)))*
+    # (1.0/(1.0+$p1))*
+    # ($m1**(1.0+$p1)-$m0**(1.0+$p1))+
+    # (($m2**(1.0+$p2)-$m1**(1.0+$p2)))*
+    # (1.0/(1.0+$p2))+
+    # (($m2**$p2)*($m2**(-$p3)))*
+    # (1.0/(1.0+$p3))*
+    # ($mmax**(1.0+$p3)-$m2**(1.0+$p3));
+    # $$array[1]=1.0/($$array[1]+1e-50);
+    # $$array[0]=$$array[1]*$m1**$p2*$m1**(-$p1);
+    # $$array[2]=$$array[1]*$m2**$p2*$m2**(-$p3);
+    # #print "ARRAY SET @_ => @$array\n";
+    # $threepart_powerlaw_consts{"@_"}=[@$array];
+
+def three_part_powerlaw(M, M0, M1, M2, M_MAX, P1, P2, P3):
+    """
+    Generalized three-part power law, usually used for mass distributions
+    """
+
+    three_part_powerlaw_constants = calculate_constants_three_part_powerlaw(M0, M1, M2, M_MAX, P1, P2, P3)
+
+    # Check trick with using the global value
+
+    if (M < M1): 
+        if M < M0: # Below lower bound
+            p = 0
+        else: 
+            p = three_part_powerlaw_constants[0] * (M ** P1) # Between M0 and M1
+    elif (M < M_MAX): # Below M_MAX
+        if (M < M2): # Between M1 and M2
+            p = three_part_powerlaw_constants[1] * (M ** P2)
+        else: # Between M2 and M_MAX
+            p = three_part_powerlaw_constants[2] * (M ** P3)            
+    else: # Above M_MAX
+        p = 0
+    return p
+
+
+
+
+# calculate_constants_three_part_powerlaw(0.08, 0.5, 1, 100, -1.3, -2.3, -2.3)
+# print(three_part_powerlaw(10, 0.08, 0.5, 1, 100, -1.3, -2.3, -2.3))
+
 
 # sub const
 # {
@@ -67,29 +127,6 @@ def powerlaw(min_val, max_val, k, x):
 #     return $C;
 # }
 
-# sub initialize_three_part_power_law_consts
-# {
-#     local $SIG{__DIE__} = sub { Carp::confess @_ };
-#     # calculate constants for this power law
-#     my ($m0,$m1,$m2,$mmax,$p1,$p2,$p3)=@_;
-#     #print "INIT three-part power law m0=$m0 m1=$m1 m2=$m2 mmax=$mmax p1=$p1 p2=$p2 p3=$p3\n";
-#     my $array=[];
-
-#     $$array[1]=(($m1**$p2)*($m1**(-$p1)))*
-#     (1.0/(1.0+$p1))*
-#     ($m1**(1.0+$p1)-$m0**(1.0+$p1))+
-#     (($m2**(1.0+$p2)-$m1**(1.0+$p2)))*
-#     (1.0/(1.0+$p2))+
-#     (($m2**$p2)*($m2**(-$p3)))*
-#     (1.0/(1.0+$p3))*
-#     ($mmax**(1.0+$p3)-$m2**(1.0+$p3));
-#     $$array[1]=1.0/($$array[1]+1e-50);
-#     $$array[0]=$$array[1]*$m1**$p2*$m1**(-$p1);
-#     $$array[2]=$$array[1]*$m2**$p2*$m2**(-$p3);
-
-#     #print "ARRAY SET @_ => @$array\n";
-#     $threepart_powerlaw_consts{"@_"}=[@$array];
-# }
 
 # sub ktg93_lnspace
 # {
@@ -191,36 +228,6 @@ def powerlaw(min_val, max_val, k, x):
 #                 $$opts{p3});
 # }
 
-# sub three_part_power_law
-# {
-#     # generalized three part power law, usually used for mass distributions
-#     # args are M, M0, M1, M2, MMAX, P1, P2, P3 (powers)
-#     my $m=shift;
-
-#     # use pre-calculated consts if possible, otherwise calculate them
-#     my $consts=$threepart_powerlaw_consts{"@_"} //
-#     initialize_three_part_power_law_consts(@_);
-
-#     my $p;
-#     if($m<$_[1])
-#     {
-#     $p = $m<$_[0]
-#        ? 0.0
-#        : $$consts[0]*($m**$_[4]);
-#     }
-#     elsif($m<$_[3])
-#     {
-#     $p = $m<$_[2]
-#         ? $$consts[1]*($m**$_[5])
-#         : $$consts[2]*($m**$_[6]);
-#     }
-#     else
-#     {
-#     $p=0.0;
-#     }
-#     return $p;
-# }
-
 # sub gaussian
 # {
 #     # Gaussian distribution function e.g. for Duquennoy + Mayor 1991

tests/population/global_variable_for_distributions.py

+import time
+
+powerlaw_const = None
+
+###
+# Function/test to see the speed up of making values global so we dont have to calculate them each time
+# Theres only a factor ~2 increase here.
+
+def with_glob():
+    global powerlaw_const
+    if not powerlaw_const:
+        # print('not defined')
+        powerlaw_const = powerlaw_constant(10, 100, -2)
+    else:
+        # print('defined')
+        return powerlaw_const
+
+def without_glob():
+    powerlaw_const = powerlaw_constant(10, 100, -2) 
+    return powerlaw_const
+
+
+def powerlaw_constant(min_val, max_val, k):
+    """
+    Function that returns the constant to normalise a powerlaw
+    """
+
+    k1 = k + 1.0
+    # print(
+    #     "Powerlaw consts from {} to {}, k={} where k1={}".format(
+    #         min_val, max_val, k, k1
+    #     )
+    # )
+
+    global powerlaw_const
+    powerlaw_const = k1 / (max_val ** k1 - min_val ** k1)
+    return powerlaw_const
+
+def powerlaw(min_val, max_val, k, x):
+    """
+    Single powerlaw with index k at x from min to max
+    """
+
+    # Handle faulty value
+    if k == -1:
+        print("wrong value for k")
+        raise ValueError
+
+    if (x < min_val) or (x > max_val):
+        print('value is out of bounds')
+        return 0
+
+    else:
+        powerlaw_const = powerlaw_constant(min_val, max_val, k)
+
+        # powerlaw
+        y = powerlaw_const * (x ** k)
+        # print(y)
+        # print(
+        #     "Power law from {} to {}: const = {}, y = {}".format(
+        #         min_val, max_val, const, y
+        #     )
+        # )
+        return y
+
+def powerlaw_with(min_val, max_val, k, x):
+    """
+    Single powerlaw with index k at x from min to max
+    """
+
+    # Handle faulty value
+    if k == -1:
+        print("wrong value for k")
+        raise ValueError
+
+    if (x < min_val) or (x > max_val):
+        return 0
+
+    else:
+        global powerlaw_const
+        if not powerlaw_const:
+            powerlaw_const = powerlaw_constant(min_val, max_val, k)
+
+        # powerlaw
+        y = powerlaw_const * (x ** k)
+        # print(
+        #     "Power law from {} to {}: const = {}, y = {}".format(
+        #         min_val, max_val, const, y
+        #     )
+        # )
+        return y
+
+steps = 1000000
+
+start_time_without_glob = time.time()
+for i in range(steps):
+    powerlaw(10, 100, -2, 20)
+stop_time_without_glob = time.time()
+
+start_time_with_glob = time.time()
+for i in range(steps):
+    powerlaw_with(10, 100, -2, 20)
+stop_time_with_glob = time.time()
+
+total_time_without = stop_time_without_glob-start_time_without_glob
+total_time_with = stop_time_with_glob - start_time_with_glob
+
+print("without: {}\nwith: {}\nRatio: {}".format(total_time_without, total_time_with, total_time_without/total_time_with))

tests/population/plot_scaling.py

+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+
+
+def calc_mean_and_std(arr):
+    return np.mean(arr), np.std(arr)
+
+# Readout
+result_file = 'comparison_result.dat'
+results = []
+name_testcase = 'laptop David'
+
+with open(result_file, 'r') as f:
+    for line in f:
+        res = list(eval(line.strip()))
+        res.append(res[-2]/res[-1])
+        results.append(res)
+
+# make dataframe
+headers = ['cores', 'total_systems', 'total_time_sequentially', 'total_time_multiprocessing', 'ratio']
+df = pd.DataFrame(results)
+df.columns = headers
+
+# Select unique amounts
+unique_amt_cores = df['cores'].unique()
+unique_amt_systems = df['total_systems'].unique()
+
+# Create dictionary with calculated means and stdevs etc
+calculated_results = []
+for i in unique_amt_cores:
+    for j in unique_amt_systems:
+        total_time_sequential_list, total_time_multiprocessing_list, total_ratio_list = [], [], []
+
+        for res in results:
+            if (res[0]==i) & (res[1]==j):
+                total_time_sequential_list.append(res[2])
+                total_time_multiprocessing_list.append(res[3])
+                total_ratio_list.append(res[4])
+
+        if (total_time_sequential_list) and (total_time_multiprocessing_list) and (total_ratio_list):
+            # calculate stuff
+            mean_time_sequential, std_sequential = calc_mean_and_std(np.array(total_time_sequential_list))
+            mean_time_multiprocessing, std_multiprocessing = calc_mean_and_std(np.array(total_time_multiprocessing_list))
+            mean_ratio, std_ratio = calc_mean_and_std(np.array(total_ratio_list))
+
+            # make dict
+            res_dict = {
+                    'cores': i, 
+                    'systems': j, 
+                    'mean_time_sequential': mean_time_sequential, 
+                    'mean_time_multiprocessing':mean_time_multiprocessing,
+                    'mean_ratio': mean_ratio,
+                    'std_sequential': std_sequential,
+                    'std_multiprocessing': std_multiprocessing,
+                    'std_ratio': std_ratio,
+                    'total_runs': len(total_time_sequential_list)
+                }
+
+            calculated_results.append(res_dict)
+
+
+# Plot
+x_position_shift = 0
+y_position_shift = -0.05
+for amt_systems in unique_amt_systems:
+    cores = []
+    speedup = []
+    std = [] 
+
+    for el in calculated_results:
+        if el['systems']==amt_systems:
+
+            cores.append(el['cores'] + x_position_shift)
+            speedup.append(el['mean_ratio'])
+            std.append(el['std_ratio'])
+
+            # add number of runs its based on
+            plt.text(el['cores'] + x_position_shift+0.01, el['mean_ratio']+y_position_shift, el['total_runs'])
+
+    plt.errorbar(cores, speedup, std, linestyle='None', marker='^', label='{} systems'.format(amt_systems))
+    x_position_shift += 0.04
+
+plt.title("Speed up ratio vs amount of cores for different amounts of systems on {}".format(name_testcase))
+plt.xlabel("Amount of cores used")
+plt.ylabel("Speed up ratio (time_linear/time_parallel)")
+plt.legend()
+plt.show()
+
+
+
+
+
+