diff --git a/binarycpython/utils/probability_distributions.py b/binarycpython/utils/probability_distributions.py
index 5506cf4eba7dd6e11a0aaa3b6bfcbe98ff2661c9..c546dd66b5f62bb0a6f489a91854781639a6493f 100644
--- a/binarycpython/utils/probability_distributions.py
+++ b/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
diff --git a/tests/population/global_variable_for_distributions.py b/tests/population/global_variable_for_distributions.py
new file mode 100644
index 0000000000000000000000000000000000000000..47b5d3f8aa8cdb531b6ee04469d15750eec842e5
--- /dev/null
+++ b/tests/population/global_variable_for_distributions.py
@@ -0,0 +1,108 @@
+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))
diff --git a/tests/population/plot_scaling.py b/tests/population/plot_scaling.py
new file mode 100644
index 0000000000000000000000000000000000000000..c819a207347dcb2263748cc14f42df932293caa7
--- /dev/null
+++ b/tests/population/plot_scaling.py
@@ -0,0 +1,94 @@
+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()
+
+
+
+
+
+