diff --git a/.gitignore b/.gitignore
index 600d2d33..700d0ba8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,3 @@
-.vscode
\ No newline at end of file
+.vscode
+.idea
+*.pyc
\ No newline at end of file
diff --git a/ramp.py b/ramp.py
index eb8852a0..e7cb311f 100644
--- a/ramp.py
+++ b/ramp.py
@@ -39,7 +39,8 @@
                 num_profiles = num_profiles * len(input_files_to_run)
             else:
                 if len(num_profiles) != len(input_files_to_run):
-                    raise ValueError("The number of profiles parameters  should match the number of input files provided")
+                    raise ValueError(
+                        "The number of profiles parameters  should match the number of input files provided")
         else:
             num_profiles = [None] * len(input_files_to_run)
 
diff --git a/ramp/core/constants.py b/ramp/core/constants.py
index c83cf9cf..016fc6fd 100644
--- a/ramp/core/constants.py
+++ b/ramp/core/constants.py
@@ -1,3 +1,25 @@
+def switch_on_parameters():
+    """
+    Calibration parameters. These can be changed in case the user has some real data against which the model can be calibrated
+    They regulate the probability of coincident switch-on within the peak window
+
+    mu_peak corresponds to \mu_{%} in [1], p.8
+    s_peak corresponds to \sigma_{%} in [1], p.8
+
+    Notes
+    -----
+    [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+        Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+        Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+    """
+
+    mu_peak = 0.5  # median value of gaussian distribution [0,1] by which the number of coincident switch_ons is randomly selected
+    s_peak = 0.5  # standard deviation (as percentage of the median value) of the gaussian distribution [0,1] above mentioned
+    op_factor = 0.5  # off-peak coincidence calculation parameter
+
+    return mu_peak, s_peak, op_factor
+
+
 OLD_TO_NEW_MAPPING = {
     "name": "user_name",
     "n_users": "num_users",
diff --git a/ramp/core/core.py b/ramp/core/core.py
index 323205c2..0ab77ced 100644
--- a/ramp/core/core.py
+++ b/ramp/core/core.py
@@ -5,8 +5,10 @@
 import numpy.ma as ma
 import pandas as pd
 import warnings
-from ramp.core.constants import NEW_TO_OLD_MAPPING, APPLIANCE_ATTRIBUTES, APPLIANCE_ARGS, WINDOWS_PARAMETERS, MAX_WINDOWS, DUTY_CYCLE_PARAMETERS
-from ramp.core.utils import read_input_file
+import random
+import math
+from ramp.core.constants import NEW_TO_OLD_MAPPING, APPLIANCE_ATTRIBUTES, APPLIANCE_ARGS, WINDOWS_PARAMETERS, DUTY_CYCLE_PARAMETERS, switch_on_parameters
+from ramp.core.utils import random_variation, duty_cycle, random_choice, read_input_file
 
 #%% Definition of Python classes that constitute the model architecture
 """
@@ -101,15 +103,13 @@ def load(self, filename):
                 user.add_appliance(**appliance_parameters)
 
 
-
 class User:
     def __init__(self, user_name="", num_users=1, user_preference=0):
         self.user_name = user_name
         self.num_users = num_users  # specifies the number of users within the class
         self.user_preference = user_preference  # allows to check if random number coincides with user preference, to distinguish between various appliance_use options (e.g. different cooking options)
-        self.App_list = (
-            []
-        )  # each instance of User (i.e. each user class) has its own list of Appliances
+        self.load = None
+        self.App_list = []  # each instance of User (i.e. each user class) has its own list of Appliances
 
     def add_appliance(self, *args, **kwargs):
 
@@ -143,6 +143,16 @@ def add_appliance(self, *args, **kwargs):
 
         return app
 
+    @property
+    def maximum_profile(self):
+        """Aggregate the theoretical maximal profiles of each appliance of the user by switching the appliance always on"""
+        user_max_profile = np.zeros(1440)
+        for appliance in self.App_list:
+            # Calculate windows curve, i.e. the theoretical maximum curve that can be obtained, for each app, by switching-on always all the 'n' apps altogether in any time-step of the functioning windows
+            app_max_profile = appliance.maximum_profile  # this computes the curve for the specific App
+            user_max_profile = np.vstack([user_max_profile, app_max_profile])  # this stacks the specific App curve in an overall curve comprising all the Apps within a User class
+        return np.transpose(np.sum(user_max_profile, axis=0)) * self.num_users
+
     def save(self, filename=None):
         answer = pd.concat([app.save() for app in self.App_list], ignore_index=True)
         if filename is not None:
@@ -195,6 +205,7 @@ def __eq__(self, other_user):
     def export_to_dataframe(self):
         return self.save()
 
+
     def Appliance(
         self,
         user,
@@ -233,6 +244,100 @@ def Appliance(
             name=name,
         )
 
+    def generate_single_load_profile(self, prof_i, peak_time_range, Year_behaviour):
+        """Generates a load profile for a single user taking all its appliances into consideration
+
+        Parameters
+        ----------
+        prof_i: int
+            ith profile requested by the user
+        peak_time_range: numpy array
+            randomised peak time range calculated using calc_peak_time_range function
+        Year_behaviour: numpy array
+            array consisting of a yearly pattern of weekends and weekdays peak_time_range
+        """
+
+        single_load = np.zeros(1440)
+        rand_daily_pref = 0 if self.user_preference == 0 else random.randint(1, self.user_preference)
+
+        for App in self.App_list:  # iterates for all the App types in the given User class
+            # initialises variables for the cycle
+            App.daily_use = np.zeros(1440)
+
+            # skip this appliance in any of the following applies
+            if (
+                    # evaluates if occasional use happens or not
+                    (random.uniform(0, 1) > App.occasional_use
+                     # evaluates if daily preference coincides with the randomised daily preference number
+                     or (App.pref_index != 0 and rand_daily_pref != App.pref_index)
+                     # checks if the app is allowed in the given yearly behaviour pattern
+                     or App.wd_we_type not in [Year_behaviour[prof_i], 2])
+            ):
+                continue
+
+            # recalculate windows start and ending times randomly, based on the inputs
+            rand_window_1 = App.calc_rand_window(window_idx=1)
+            rand_window_2 = App.calc_rand_window(window_idx=2)
+            rand_window_3 = App.calc_rand_window(window_idx=3)
+            rand_windows = [rand_window_1, rand_window_2, rand_window_3]
+
+            # random variability is applied to the total functioning time and to the duration
+            # of the duty cycles provided they have been specified
+            # step 2a of [1]
+            rand_time = App.rand_total_time_of_use(rand_window_1, rand_window_2, rand_window_3)
+
+            # redefines functioning windows based on the previous randomisation of the boundaries
+            # step 2b of [1]
+            if App.flat == 'yes':
+                # for "flat" appliances the algorithm stops right after filling the newly
+                # created windows without applying any further stochasticity
+                total_power_value = App.power[prof_i] * App.number
+                for rand_window in rand_windows:
+                    App.daily_use[rand_window[0]:rand_window[1]] = np.full(np.diff(rand_window),
+                                                                           total_power_value)
+                single_load = single_load + App.daily_use
+                continue
+            else:
+                # "non-flat" appliances a mask is applied on the newly defined windows and
+                # the algorithm goes further on
+                for rand_window in rand_windows:
+                    App.daily_use[rand_window[0]:rand_window[1]] = np.full(np.diff(rand_window), 0.001)
+            App.daily_use_masked = np.zeros_like(ma.masked_not_equal(App.daily_use, 0.001))
+
+            # calculates randomised cycles taking the random variability in the duty cycle duration
+            App.assign_random_cycles()
+
+            # steps 2c-2e repeated until the sum of the durations of all the switch-on events equals rand_time
+            App.generate_load_profile(rand_time, peak_time_range, rand_window_1, rand_window_2, rand_window_3, power=App.power[prof_i])
+
+            single_load = single_load + App.daily_use  # adds the Appliance load profile to the single User load profile
+        return single_load
+
+    def generate_aggregated_load_profile(self, prof_i, peak_time_range, Year_behaviour):
+        """Generates an aggregated load profile from single load profile of each user
+
+            Each single load profile has its own separate randomisation
+
+        Parameters
+        ----------
+        prof_i: int
+            ith profile requested by the user
+        peak_time_range: numpy array
+            randomised peak time range calculated using calc_peak_time_range function
+        Year_behaviour: numpy array
+            array consisting of a yearly pattern of weekends and weekdays peak_time_range
+
+        Returns
+        -------
+        User.load : numpy array
+            the aggregate load profile of all the users within a user class
+        """
+
+        self.load = np.zeros(1440)  # initialise empty load for User instance
+        for _ in range(self.num_users):
+            # iterates for every single user within a User class.
+            self.load = self.load + self.generate_single_load_profile(prof_i, peak_time_range, Year_behaviour)
+
 
 class Appliance:
     def __init__(
@@ -415,7 +520,7 @@ def windows(self, window_1=None, window_2=None, random_var_w=0, window_3=None):
                 raise ValueError("Windows 3 is not provided although 3 windows were declared")
         else:
             self.window_3 = window_3
-        print(self.window_1, self.window_2, self.window_3)
+
         self.random_var_w = random_var_w #percentage of variability in the start and ending times of the windows
         self.daily_use = np.zeros(1440) #create an empty daily use profile
         self.daily_use[self.window_1[0]:(self.window_1[1])] = np.full(np.diff(self.window_1),0.001) #fills the daily use profile with infinitesimal values that are just used to identify the functioning windows
@@ -431,6 +536,82 @@ def windows(self, window_1=None, window_2=None, random_var_w=0, window_3=None):
             self.cw11 = self.window_1
             self.cw12 = self.window_2
 
+    def assign_random_cycles(self):
+        """Calculates randomised cycles taking the random variability in the duty cycle duration"""
+        if self.fixed_cycle >= 1:
+            p_11 = random_variation(var=self.thermal_p_var, norm=self.p_11) #randomly variates the power of thermal apps, otherwise variability is 0
+            p_12 = random_variation(var=self.thermal_p_var, norm=self.p_12) #randomly variates the power of thermal apps, otherwise variability is 0
+            self.random_cycle1 = duty_cycle(var=self.r_c1, t1=self.t_11, p1=p_11, t2=self.t_12, p2=p_12) #randomise also the fixed cycle
+            self.random_cycle2 = self.random_cycle1
+            self.random_cycle3 = self.random_cycle1
+            if self.fixed_cycle >= 2:
+                p_21 = random_variation(var=self.thermal_p_var, norm=self.p_21) #randomly variates the power of thermal apps, otherwise variability is 0
+                p_22 = random_variation(var=self.thermal_p_var, norm=self.p_22) #randomly variates the power of thermal apps, otherwise variability is 0
+                self.random_cycle2 = duty_cycle(var=self.r_c2, t1=self.t_21, p1=p_21, t2=self.t_22, p2=p_22) #randomise also the fixed cycle
+
+                if self.fixed_cycle >= 3:
+                    p_31 = random_variation(var=self.thermal_p_var, norm=self.p_31) #randomly variates the power of thermal apps, otherwise variability is 0
+                    p_32 = random_variation(var=self.thermal_p_var, norm=self.p_32) #randomly variates the power of thermal apps, otherwise variability is 0
+                    self.random_cycle1 = random_choice(self.r_c1, t1=self.t_11, p1=p_11, t2=self.t_12, p2=p_12)
+
+                    self.random_cycle2 = random_choice(self.r_c2, t1=self.t_21, p1=p_21, t2=self.t_22, p2=p_22)
+
+                    self.random_cycle3 = random_choice(self.r_c3, t1=self.t_31, p1=p_31, t2=self.t_32, p2=p_32)
+
+
+    def update_daily_use(self, coincidence, power, indexes):
+        """Update the daily use depending on existence of duty cycles of the Appliance instance
+
+            This corresponds to step 2d. and 2e. of [1]
+
+        Notes
+        -----
+        [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+            Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+            Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+
+        """
+
+        if self.fixed_cycle > 0:  # evaluates if the app has some duty cycles to be considered
+            evaluate = np.round(np.mean(indexes)) if indexes.size > 0 else 0
+            # selects the proper duty cycle and puts the corresponding power values in the indexes range
+            if evaluate in range(self.cw11[0], self.cw11[1]) or evaluate in range(self.cw12[0], self.cw12[1]):
+                np.put(self.daily_use, indexes, (self.random_cycle1 * coincidence))
+                np.put(self.daily_use_masked, indexes, (self.random_cycle1 * coincidence), mode='clip')
+            elif evaluate in range(self.cw21[0], self.cw21[1]) or evaluate in range(self.cw22[0], self.cw22[1]):
+                np.put(self.daily_use, indexes, (self.random_cycle2 * coincidence))
+                np.put(self.daily_use_masked, indexes, (self.random_cycle2 * coincidence), mode='clip')
+            else:
+                np.put(self.daily_use, indexes, (self.random_cycle3 * coincidence))
+                np.put(self.daily_use_masked, indexes, (self.random_cycle3 * coincidence), mode='clip')
+        else:  # if no duty cycles are specified, a regular switch_on event is modelled
+            # randomises also the App Power if thermal_p_var is on
+            np.put(self.daily_use, indexes, (random_variation(var=self.thermal_p_var, norm=coincidence * power)))
+            np.put(self.daily_use_masked, indexes, (random_variation(var=self.thermal_p_var, norm=coincidence * power)), mode='clip')
+        # updates the mask excluding the current switch_on event to identify the free_spots for the next iteration
+        self.daily_use_masked = np.zeros_like(ma.masked_greater_equal(self.daily_use_masked, 0.001))
+
+    def calc_rand_window(self, window_idx=1, window_range_limits=np.array([0, 1440])):
+        _window = self.__getattribute__(f'window_{window_idx}')
+        _random_var = self.__getattribute__(f'random_var_{window_idx}')
+        rand_window = np.array([int(random.uniform(_window[0] - _random_var, _window[0] + _random_var)),
+                                int(random.uniform(_window[1] - _random_var, _window[1] + _random_var))])
+        if rand_window[0] < window_range_limits[0]:
+            rand_window[0] = window_range_limits[0]
+        if rand_window[1] > window_range_limits[1]:
+            rand_window[1] = window_range_limits[1]
+
+        return rand_window
+
+    @property
+    def maximum_profile(self):
+        """Virtual maximum load profile of the appliance
+
+            It assumes the appliance is always switched-on with maximum power and
+            numerosity during all of its potential windows of use
+        """
+        return self.daily_use * np.mean(self.power) * self.number
+
     def specific_cycle(self, cycle_num, **kwargs):
         if cycle_num == 1:
             self.specific_cycle_1(**kwargs)
@@ -439,7 +620,8 @@ def specific_cycle(self, cycle_num, **kwargs):
         elif cycle_num == 3:
             self.specific_cycle_3(**kwargs)
 
-        #if needed, specific duty cycles can be defined for each Appliance, for a maximum of three different ones
+        # if needed, specific duty cycles can be defined for each Appliance, for a maximum of three different ones
+
     def specific_cycle_1(self, p_11 = 0, t_11 = 0, p_12 = 0, t_12 = 0, r_c1 = 0, cw11=None, cw12=None):
         self.p_11 = p_11 #power absorbed during first part of the duty cycle
         self.t_11 = int(t_11) #duration of first part of the duty cycle
@@ -490,3 +672,185 @@ def cycle_behaviour(self, cw11 = np.array([0,0]), cw12 = np.array([0,0]), cw21 =
         self.cw31 = cw31 #same for cycle 3
         self.cw32 = cw32
 
+    def rand_total_time_of_use(self, rand_window_1, rand_window_2, rand_window_3):
+        """Randomised total time of use of the Appliance instance
+
+            This corresponds to step 2a. of [1]
+
+        Notes
+        -----
+        [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+            Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+            Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+        """
+        random_var_t = random_variation(var=self.time_fraction_random_variability)
+
+        rand_time = round(random.uniform(self.func_time, int(self.func_time * random_var_t)))
+
+        # check that the total randomised time of use does not exceed the total space available in the windows
+        if rand_time > 0.99 * (np.diff(rand_window_1) + np.diff(rand_window_2) + np.diff(rand_window_3)):
+            rand_time = int(0.99 * (np.diff(rand_window_1) + np.diff(rand_window_2) + np.diff(rand_window_3)))
+        return rand_time
+
+    def switch_on(self, rand_window_1, rand_window_2, rand_window_3):
+        """Return a random switch-on time of the Appliance instance
+
+            This corresponds to step 2c. of [1]
+
+        Notes
+        -----
+        [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+            Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+            Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+        """
+        # check how many windows to consider
+        if self.num_windows == 1:
+            return int(random.choice([random.uniform(rand_window_1[0], (rand_window_1[1]))]))
+        elif self.num_windows == 2:
+            return int(random.choice(np.concatenate((np.arange(rand_window_1[0], rand_window_1[1]), np.arange(rand_window_2[0], rand_window_2[1])), axis=0)))
+        else:
+            return int(random.choice(np.concatenate((np.arange(rand_window_1[0], rand_window_1[1]), np.arange(rand_window_2[0], rand_window_2[1]), np.arange(rand_window_3[0], rand_window_3[1]),), axis=0)))
+
+    def calc_indexes_for_rand_switch_on(self, switch_on, rand_time, max_free_spot, rand_window):
+        """Identifies a random switch on time within the available functioning windows"""
+        if np.any(self.daily_use[switch_on:rand_window[1]] != 0.001):  # control to check if there are any other switch on times after the current one
+            next_switch = [switch_on + k[0] for k in np.where(self.daily_use[switch_on:] != 0.001)]  # identifies the position of next switch on time and sets it as a limit for the duration of the current switch on
+            if (next_switch[0] - switch_on) >= self.func_cycle and max_free_spot >= self.func_cycle:
+                upper_limit = min((next_switch[0] - switch_on), min(rand_time, rand_window[1] - switch_on))
+            elif (next_switch[0] - switch_on) < self.func_cycle and max_free_spot >= self.func_cycle:  # if next switch_on event does not allow for a minimum functioning cycle without overlapping, but there are other larger free spots, the cycle tries again from the beginning
+                return None # this will trigger a continue in the loop
+            else:
+                upper_limit = next_switch[0] - switch_on  # if there are no other options to reach the total time of use, empty spaces are filled without minimum cycle restrictions until reaching the limit
+        else:
+            upper_limit = min(rand_time, rand_window[1] - switch_on)  # if there are no other switch-on events after the current one, the upper duration limit is set this way
+        if upper_limit >= self.func_cycle:
+            indexes = np.arange(switch_on, switch_on + (int(random.uniform(self.func_cycle, upper_limit))))
+        else:
+            indexes = np.arange(switch_on, switch_on + upper_limit)
+        return indexes
+
+    def calc_coincident_switch_on(self, inside_peak_window=True):
+        """Computes how many of the 'n' Appliance instance are switched on simultaneously
+
+            Implement eqs. 3 and 4 of [1]
+
+        Notes
+        -----
+        [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+            Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+            Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+        """
+        s_peak, mu_peak, op_factor = switch_on_parameters()
+
+        # check if indexes are within peak window
+        if inside_peak_window is True and self.fixed == 'no':
+            # calculates coincident behaviour within the peak time range
+            # eq. 4 of [1]
+            coincidence = min(self.number, max(1, math.ceil(random.gauss(mu=(self.number * mu_peak + 0.5), sigma=(s_peak * self.number * mu_peak)))))
+        # check if indexes are off-peak
+        elif inside_peak_window is False and self.fixed == 'no':
+            # calculates probability of coincident switch_ons off-peak
+            # eq. 3 of [1]
+            prob = random.uniform(0, (self.number - op_factor) / self.number)
+
+            # randomly selects how many appliances are on at the same time
+            array = np.arange(0, self.number) / self.number
+            try:
+                on_number = np.max(np.where(prob >= array)) + 1
+            except ValueError:
+                on_number = 1
+
+            coincidence = on_number
+        else:
+            # All 'n' copies of an Appliance instance are switched on altogether
+            coincidence = self.number
+        return coincidence
+
+    def generate_load_profile(self, rand_time, peak_time_range, rand_window_1, rand_window_2, rand_window_3, power):
+        """Generate load profile of the Appliance instance by updating its daily_use attribute
+
+            Repeat steps 2c – 2e of [1] until the sum of the durations of all the switch-on events equals
+            the randomised total time of use of the Appliance
+
+        Notes
+        -----
+        [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+            Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+            Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+        """
+        max_free_spot = rand_time  # free spots are used to detect if there's still space for switch_ons. Before calculating actual free spots, the max free spot is set equal to the entire randomised func_time
+        tot_time = 0
+        while tot_time <= rand_time:
+            # Identifies a random switch on time within the available functioning windows
+            # step 2c of [1]
+            switch_on = self.switch_on(rand_window_1, rand_window_2, rand_window_3)
+
+            if self.daily_use[switch_on] == 0.001:
+                # control to check if the Appliance instance is not already on
+                # at the randomly selected switch-on time
+                if switch_on in range(rand_window_1[0], rand_window_1[1]):
+                    # random switch_on happens in rand_windows_1
+                    rand_window = rand_window_1
+                elif switch_on in range(rand_window_2[0], rand_window_2[1]):
+                    # random switch_on happens in rand_windows_2
+                    rand_window = rand_window_2
+                else:
+                    # if switch_on is not in rand_window_1 nor in rand_window_2, it shall be in rand_window_3.
+                    rand_window = rand_window_3
+
+                indexes = self.calc_indexes_for_rand_switch_on(
+                    switch_on=switch_on,
+                    rand_time=rand_time,
+                    max_free_spot=max_free_spot,
+                    rand_window=rand_window
+                )
+                if indexes is None:
+                    continue
+
+                # the count of total time is updated with the size of the indexes array
+                tot_time = tot_time + indexes.size
+
+                if tot_time > rand_time:
+                    # the total functioning time is reached, a correction is applied to avoid overflow of indexes
+                    indexes_adj = indexes[:-(tot_time - rand_time)]
+                    inside_peak_window = np.in1d(peak_time_range, indexes_adj).any()
+                    # Computes how many of the 'n' of the Appliance instance are switched on simultaneously
+                    coincidence = self.calc_coincident_switch_on(
+                        inside_peak_window
+                    )
+                    # Update the daily use depending on existence of duty cycles of the Appliance instance
+                    self.update_daily_use(
+                        coincidence,
+                        power=power,
+                        indexes=indexes_adj
+                    )
+                    break  # exit cycle and go to next Appliance
+
+                else:
+                    # the total functioning time has not yet exceeded the rand_time
+                    # Computes how many of the 'n' of the Appliance instance are switched on simultaneously
+                    inside_peak_window = np.in1d(peak_time_range, indexes).any()
+                    coincidence = self.calc_coincident_switch_on(
+                        inside_peak_window
+                    )
+                    # Update the daily use depending on existence of duty cycles of the Appliance instance
+                    self.update_daily_use(
+                        coincidence,
+                        power=power,
+                        indexes=indexes
+                    )
+
+                free_spots = []  # calculate how many free spots remain for further switch_ons
+                try:
+                    for j in ma.notmasked_contiguous(self.daily_use_masked):
+                        free_spots.append(j.stop - j.start)
+                except TypeError:
+                    free_spots = [0]
+                max_free_spot = max(free_spots)
+
+            else:
+                # if the random switch_on falls somewhere where the Appliance instance
+                # has been already turned on, tries again from beginning of the while cycle
+                # TODO not sure this code is useful as the while loop would continue anyway in that case
+                continue
+
diff --git a/ramp/core/initialise.py b/ramp/core/initialise.py
index 09456083..a47e01bd 100644
--- a/ramp/core/initialise.py
+++ b/ramp/core/initialise.py
@@ -2,65 +2,78 @@
 
 #%% Initialisation of a model instance
 
-import numpy as np 
+import numpy as np
 import importlib
 from ramp.core.core import UseCase
 
 
 def yearly_pattern():
-    '''
-    Definition of a yearly pattern of weekends and weekdays, in case some appliances have specific wd/we behaviour
-    '''
-    #Yearly behaviour pattern
-    Year_behaviour = np.zeros(365)
-    Year_behaviour[5:365:7] = 1
-    Year_behaviour[6:365:7] = 1
-    
-    return(Year_behaviour)
+    """Definition of a yearly pattern of weekends and weekdays, in case some appliances have specific wd/we behaviour"""
+
+    year_behaviour = np.zeros(365)
+    year_behaviour[5:365:7] = 1
+    year_behaviour[6:365:7] = 1
+
+    return year_behaviour
 
 
 def user_defined_inputs(j=None, fname=None):
-    '''
-    Imports an input file and returns a processed User_list
-    '''
+    """Imports an input file and returns a processed user_list
+
+    Parameters
+    ----------
+    j: int
+        the index of a .py input file that is in the format input_files/input_file_j.py
+    fname: path
+        path to a .xlsx input file
+        if provided, overrides the loading of input_files/input_file_j.py
+
+    Returns
+    -------
+    A list of User instances
+    """
+
     # Back compatibility with old code
     if j is not None:
-        file_module = importlib.import_module(f'ramp.input_files.input_file_{j}')
-        User_list = file_module.User_list
+        file_module = importlib.import_module(f"ramp.input_files.input_file_{j}")
+        user_list = file_module.User_list
 
     if fname is not None:
         usecase = UseCase()
         usecase.load(fname)
-        User_list = usecase.users
+        user_list = usecase.users
 
-    return(User_list)
+    return user_list
 
 
-def Initialise_model(num_profiles):
-    '''
-    The model is ready to be initialised
-    '''
+def initialise_inputs(j=None, fname=None, num_profiles=None):
+    """Loads the provided input file and prompt the user for number of profiles if not defined
 
-    if num_profiles is None:
-        num_profiles = int(input("please indicate the number of profiles to be generated: ")) #asks the user how many profiles (i.e. code runs) he wants
-    print('Please wait...') 
-    Profile = [] #creates an empty list to store the results of each code run, i.e. each stochastically generated profile
-    
-    return (Profile, num_profiles)
-    
-def Initialise_inputs(j=None, fname=None):
-    Year_behaviour = yearly_pattern()
+    Parameters
+    ----------
+    j: int
+        the index of a .py input file that is in the format input_files/input_file_j.py
+    fname: path
+        path to a .xlsx input file
+        if provided, overrides the loading of input_files/input_file_j.py
+    num_profiles: int
+        the number of different usecase profiles which need to be generated
+
+    Returns
+    -------
+
+    """
+
+    year_behaviour = yearly_pattern()
     user_list = user_defined_inputs(j, fname)
-    
-    # Calibration parameters
-    '''
-    Calibration parameters. These can be changed in case the user has some real data against which the model can be calibrated
-    They regulate the probabilities defining the largeness of the peak window and the probability of coincident switch-on within the peak window
-    '''
-    peak_enlarg = 0.15 #percentage random enlargement or reduction of peak time range length
-    mu_peak = 0.5 #median value of gaussian distribution [0,1] by which the number of coincident switch_ons is randomly selected
-    s_peak = 0.5 #standard deviation (as percentage of the median value) of the gaussian distribution [0,1] above mentioned
-    op_factor = 0.5 #off-peak coincidence calculation parameter
-
-    return (peak_enlarg, mu_peak, s_peak, op_factor, Year_behaviour, user_list)
 
+    peak_enlarge = 0.15  # percentage random enlargement or reduction of peak time range length, corresponds to \delta_{peak} in [1], p.7
+
+    if num_profiles is None:
+        # asks the user how many profiles (i.e. code runs) they want
+        num_profiles = int(
+            input("please indicate the number of profiles to be generated: ")
+        )
+        print("Please wait...")
+
+    return (peak_enlarge, year_behaviour, user_list, num_profiles)
diff --git a/ramp/core/stochastic_process.py b/ramp/core/stochastic_process.py
index ffb18c56..eb262391 100644
--- a/ramp/core/stochastic_process.py
+++ b/ramp/core/stochastic_process.py
@@ -2,293 +2,88 @@
 
 #%% Import required libraries
 import numpy as np
-import numpy.ma as ma
 import random 
 import math
-from ramp.core.initialise import Initialise_model, Initialise_inputs
+from ramp.core.initialise import initialise_inputs
 
 #%% Core model stochastic script
 
 
-
-
-def Stochastic_Process(j=None, fname=None, num_profiles=None):
-    Profile, num_profiles = Initialise_model(num_profiles)
-    peak_enlarg, mu_peak, s_peak, op_factor, Year_behaviour, User_list = Initialise_inputs(j, fname)
-    '''
-    Calculation of the peak time range, which is used to discriminate between off-peak and on-peak coincident switch-on probability
-    Calculates first the overall Peak Window (taking into account all User classes). 
+def calc_peak_time_range(user_list, peak_enlarge=0.15):
+    """
+    Calculate the peak time range, which is used to discriminate between off-peak and on-peak coincident switch-on probability
+    Calculate first the overall Peak Window (taking into account all User classes).
+    The peak time range corresponds to `peak time frame` variable in eq. (1) of [1]
     The peak window is just a time window in which coincident switch-on of multiple appliances assumes a higher probability than off-peak
     Within the peak window, a random peak time is calculated and then enlarged into a peak_time_range following again a random procedure
-    '''
-    windows_curve = np.zeros(1440) #creates an empty daily profile
-    Tot_curve = np.zeros(1440) #creates another empty daily profile
-    for Us in User_list:
-        App_count = 0
-        for App in Us.App_list:
-            #Calculate windows curve, i.e. the theoretical maximum curve that can be obtained, for each app, by switching-on always all the 'n' apps altogether in any time-step of the functioning windows
-            single_wcurve = Us.App_list[App_count].daily_use*np.mean(Us.App_list[App_count].power)*Us.App_list[App_count].number #this computes the curve for the specific App
-            windows_curve = np.vstack([windows_curve, single_wcurve]) #this stacks the specific App curve in an overall curve comprising all the Apps within a User class
-            App_count += 1
-        Us.windows_curve = windows_curve #after having iterated for all the Apps within a User class, saves the overall User class theoretical maximum curve
-        Us.windows_curve = np.transpose(np.sum(Us.windows_curve, axis = 0))*Us.num_users
-        Tot_curve = Tot_curve + Us.windows_curve #adds the User's theoretical max profile to the total theoretical max comprising all classes
-    peak_window = np.transpose(np.argwhere(Tot_curve == np.amax(Tot_curve))) #Find the peak window within the theoretical max profile
-    peak_time = round(random.normalvariate(round(np.average(peak_window)),1/3*(peak_window[0,-1]-peak_window[0,0]))) #Within the peak_window, randomly calculate the peak_time using a gaussian distribution
-    peak_time_range = np.arange((peak_time-round(math.fabs(peak_time-(random.gauss(peak_time,(peak_enlarg*peak_time)))))),(peak_time+round(math.fabs(peak_time-random.gauss(peak_time,(peak_enlarg*peak_time)))))) #the peak_time is randomly enlarged based on the calibration parameter peak_enlarg
-    
-    '''
-    The core stochastic process starts here. For each profile requested by the software user, 
-    each Appliance instance within each User instance is separately and stochastically generated
-    '''
-    for prof_i in range(num_profiles): #the whole code is repeated for each profile that needs to be generated
-        Tot_Classes = np.zeros(1440) #initialise an empty daily profile that will be filled with the sum of the hourly profiles of each User instance
-        for Us in User_list: #iterates for each User instance (i.e. for each user class)
-            Us.load = np.zeros(1440) #initialise empty load for User instance
-            for i in range(Us.num_users): #iterates for every single user within a User class. Each single user has its own separate randomisation
-                if Us.user_preference == 0:
-                    rand_daily_pref = 0
-                    pass
-                else:
-                    rand_daily_pref = random.randint(1,Us.user_preference)
-                for App in Us.App_list: #iterates for all the App types in the given User class
-                    #initialises variables for the cycle
-                    tot_time = 0
-                    App.daily_use = np.zeros(1440)
-                    if random.uniform(0,1) > App.occasional_use: #evaluates if occasional use happens or not
-                        continue
-                    else:
-                        pass
-                    
-                    if App.pref_index == 0:
-                        pass
-                    else:
-                        if rand_daily_pref == App.pref_index: #evaluates if daily preference coincides with the randomised daily preference number
-                            pass
-                        else:
-                            continue
-                    if App.wd_we_type == Year_behaviour[prof_i] or App.wd_we_type == 2 : #checks if the app is allowed in the given yearly behaviour pattern
-                        pass
-                    else:
-                        continue
 
-                    #recalculate windows start and ending times randomly, based on the inputs
-                    rand_window_1 = np.array([int(random.uniform((App.window_1[0]-App.random_var_1),(App.window_1[0]+App.random_var_1))),int(random.uniform((App.window_1[1]-App.random_var_1),(App.window_1[1]+App.random_var_1)))])
-                    if rand_window_1[0] < 0:
-                        rand_window_1[0] = 0
-                    if rand_window_1[1] > 1440:
-                        rand_window_1[1] = 1440
-    
-                    rand_window_2 = np.array([int(random.uniform((App.window_2[0]-App.random_var_2),(App.window_2[0]+App.random_var_2))),int(random.uniform((App.window_2[1]-App.random_var_2),(App.window_2[1]+App.random_var_2)))])
-                    if rand_window_2[0] < 0:
-                        rand_window_2[0] = 0
-                    if rand_window_2[1] > 1440:
-                        rand_window_2[1] = 1440
-                            
-                    rand_window_3 = np.array([int(random.uniform((App.window_3[0]-App.random_var_3),(App.window_3[0]+App.random_var_3))),int(random.uniform((App.window_3[1]-App.random_var_3),(App.window_3[1]+App.random_var_3)))])
-                    if rand_window_3[0] < 0:
-                        rand_window_3[0] = 0
-                    if rand_window_3[1] > 1440:
-                        rand_window_3[1] = 1440
-    
-                    #redefines functioning windows based on the previous randomisation of the boundaries
-                    if App.flat == 'yes': #if the app is "flat" the code stops right after filling the newly created windows without applying any further stochasticity
-                        App.daily_use[rand_window_1[0]:rand_window_1[1]] = np.full(np.diff(rand_window_1),App.power[prof_i]*App.number)
-                        App.daily_use[rand_window_2[0]:rand_window_2[1]] = np.full(np.diff(rand_window_2),App.power[prof_i]*App.number)
-                        App.daily_use[rand_window_3[0]:rand_window_3[1]] = np.full(np.diff(rand_window_3),App.power[prof_i]*App.number)
-                        Us.load = Us.load + App.daily_use
-                        continue
-                    else: #otherwise, for "non-flat" apps it puts a mask on the newly defined windows and continues    
-                        App.daily_use[rand_window_1[0]:rand_window_1[1]] = np.full(np.diff(rand_window_1),0.001)
-                        App.daily_use[rand_window_2[0]:rand_window_2[1]] = np.full(np.diff(rand_window_2),0.001)
-                        App.daily_use[rand_window_3[0]:rand_window_3[1]] = np.full(np.diff(rand_window_3),0.001)
-                    App.daily_use_masked = np.zeros_like(ma.masked_not_equal(App.daily_use,0.001))
-                    
-                    #random variability is applied to the total functioning time and to the duration of the duty cycles, if they have been specified
-                    random_var_t = random.uniform((1-App.time_fraction_random_variability),(1+App.time_fraction_random_variability))
-                    if App.fixed_cycle == 1:
-                        App.p_11 = App.p_11*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_12 = App.p_12*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        random_cycle1 = np.concatenate(((np.ones(int(App.t_11*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_11),(np.ones(int(App.t_12*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_12))) #randomise also the fixed cycle
-                        random_cycle2 = random_cycle1
-                        random_cycle3 = random_cycle1
-                    elif App.fixed_cycle == 2:
-                        App.p_11 = App.p_11*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_12 = App.p_12*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_21 = App.p_21*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_22 = App.p_22*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        random_cycle1 = np.concatenate(((np.ones(int(App.t_11*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_11),(np.ones(int(App.t_12*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_12))) #randomise also the fixed cycle
-                        random_cycle2 = np.concatenate(((np.ones(int(App.t_21*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_21),(np.ones(int(App.t_22*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_22))) #randomise also the fixed cycle
-                        random_cycle3 = random_cycle1
-                    elif App.fixed_cycle == 3:
-                        App.p_11 = App.p_11*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_12 = App.p_12*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_21 = App.p_21*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_22 = App.p_22*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_31 = App.p_31*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        App.p_32 = App.p_32*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var))) #randomly variates the power of thermal apps, otherwise variability is 0
-                        random_cycle1 = random.choice([np.concatenate(((np.ones(int(App.t_11*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_11),(np.ones(int(App.t_12*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_12))),np.concatenate(((np.ones(int(App.t_12*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_12),(np.ones(int(App.t_11*(random.uniform((1+App.r_c1),(1-App.r_c1)))))*App.p_11)))]) #randomise also the fixed cycle
-                        random_cycle2 = random.choice([np.concatenate(((np.ones(int(App.t_21*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_21),(np.ones(int(App.t_22*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_22))),np.concatenate(((np.ones(int(App.t_22*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_22),(np.ones(int(App.t_21*(random.uniform((1+App.r_c2),(1-App.r_c2)))))*App.p_21)))])                    
-                        random_cycle3 = random.choice([np.concatenate(((np.ones(int(App.t_31*(random.uniform((1+App.r_c3),(1-App.r_c3)))))*App.p_31),(np.ones(int(App.t_32*(random.uniform((1+App.r_c3),(1-App.r_c3)))))*App.p_32))),np.concatenate(((np.ones(int(App.t_32*(random.uniform((1+App.r_c3),(1-App.r_c3)))))*App.p_32),(np.ones(int(App.t_31*(random.uniform((1+App.r_c3),(1-App.r_c3)))))*App.p_31)))])#this is to avoid that all cycles are sincronous                      
-                    else:
-                        pass
-                    rand_time = round(random.uniform(App.func_time,int(App.func_time*random_var_t)))
-                    #control to check that the total randomised time of use does not exceed the total space available in the windows
-                    if rand_time > 0.99*(np.diff(rand_window_1)+np.diff(rand_window_2)+np.diff(rand_window_3)):
-                        rand_time = int(0.99*(np.diff(rand_window_1)+np.diff(rand_window_2)+np.diff(rand_window_3)))
-                    max_free_spot = rand_time #free spots are used to detect if there's still space for switch_ons. Before calculating actual free spots, the max free spot is set equal to the entire randomised func_time
-                           
-                    while tot_time <= rand_time: #this is the key cycle, which runs for each App until the switch_ons and their duration equals the randomised total time of use of the App
-                            #check how many windows to consider
-                            if App.num_windows == 1:
-                                switch_on = int(random.choice([random.uniform(rand_window_1[0],(rand_window_1[1]))]))
-                            elif App.num_windows == 2:
-                                switch_on = int(random.choice(np.concatenate((np.arange(rand_window_1[0],rand_window_1[1]),np.arange(rand_window_2[0],rand_window_2[1])),axis=0)))
-                            else: 
-                                switch_on = int(random.choice(np.concatenate((np.arange(rand_window_1[0],
-                                                        rand_window_1[1]),
-                                              np.arange(rand_window_2[0],
-                                                        rand_window_2[1]),
-                                              np.arange(rand_window_3[0],
-                                                        rand_window_3[1]),
-                                              ), axis=0)))                            
-                            #Identifies a random switch on time within the available functioning windows
-                            if App.daily_use[switch_on] == 0.001: #control to check if the app is not already on at the randomly selected switch-on time
-                                if switch_on in range(rand_window_1[0],rand_window_1[1]):
-                                    if np.any(App.daily_use[switch_on:rand_window_1[1]]!=0.001): #control to check if there are any other switch on times after the current one    
-                                        next_switch = [switch_on + k[0] for k in np.where(App.daily_use[switch_on:]!=0.001)] #identifies the position of next switch on time and sets it as a limit for the duration of the current switch on
-                                        if (next_switch[0] - switch_on) >= App.func_cycle and max_free_spot >= App.func_cycle:
-                                            upper_limit = min((next_switch[0]-switch_on),min(rand_time,rand_window_1[1]-switch_on))
-                                        elif (next_switch[0] - switch_on) < App.func_cycle and max_free_spot >= App.func_cycle: #if next switch_on event does not allow for a minimum functioning cycle without overlapping, but there are other larger free spots, the cycle tries again from the beginning
-                                            continue
-                                        else:
-                                            upper_limit = next_switch[0]-switch_on #if there are no other options to reach the total time of use, empty spaces are filled without minimum cycle restrictions until reaching the limit                                              
-                                    else:
-                                        upper_limit = min(rand_time,rand_window_1[1]-switch_on) #if there are no other switch-on events after the current one, the upper duration limit is set this way
-                                    
-                                    if upper_limit >= App.func_cycle: #if the upper limit is higher than minimum functioning time, an array of indexes is created to be later put in the profile
-                                        indexes = np.arange(switch_on,switch_on+(int(random.uniform(App.func_cycle,upper_limit)))) #a random duration is chosen between the upper limit and the minimum cycle
-                                    else:
-                                        indexes = np.arange(switch_on,switch_on+upper_limit) #this is the case in which empty spaces need to be filled without constraints to reach the total time goal
-                                        
-                                elif switch_on in range(rand_window_2[0],rand_window_2[1]): #if random switch_on happens in windows2, same code as above is repeated for windows2
-                                    if np.any(App.daily_use[switch_on:rand_window_2[1]]!=0.001):
-                                        next_switch = [switch_on + k[0] for k in np.where(App.daily_use[switch_on:]!=0.001)]
-                                        if (next_switch[0] - switch_on) >= App.func_cycle and max_free_spot >= App.func_cycle:
-                                            upper_limit = min((next_switch[0]-switch_on),min(rand_time,rand_window_2[1]-switch_on))
-                                        elif (next_switch[0] - switch_on) < App.func_cycle and max_free_spot >= App.func_cycle:
-                                            continue
-                                        else:
-                                            upper_limit = next_switch[0]-switch_on
-                                    
-                                    else:
-                                        upper_limit = min(rand_time,rand_window_2[1]-switch_on)
-                                    
-                                    if upper_limit >= App.func_cycle:
-                                        indexes = np.arange(switch_on,switch_on+(int(random.uniform(App.func_cycle,upper_limit))))
-                                    else:    
-                                        indexes = np.arange(switch_on,switch_on+upper_limit)
-                                        
-                                else: #if switch_on is not in window1 nor in window2, it shall be in window3. Same code is repreated
-                                    if np.any(App.daily_use[switch_on:rand_window_3[1]]!=0.001):
-                                        next_switch = [switch_on + k[0] for k in np.where(App.daily_use[switch_on:]!=0.001)]
-                                        if (next_switch[0] - switch_on) >= App.func_cycle and max_free_spot >= App.func_cycle:
-                                            upper_limit = min((next_switch[0]-switch_on),min(rand_time,rand_window_3[1]-switch_on))
-                                        elif (next_switch[0] - switch_on) < App.func_cycle and max_free_spot >= App.func_cycle:
-                                            continue
-                                        else:
-                                            upper_limit = next_switch[0]-switch_on
-                                    
-                                    else:
-                                        upper_limit = min(rand_time,rand_window_3[1]-switch_on)
-                                    
-                                    if upper_limit >= App.func_cycle:
-                                        indexes = np.arange(switch_on,switch_on+(int(random.uniform(App.func_cycle,upper_limit))))
-                                    else:    
-                                        indexes = np.arange(switch_on,switch_on+upper_limit)
-                                        
-                                tot_time = tot_time + indexes.size #the count of total time is updated with the size of the indexes array
-                                
-                                if tot_time > rand_time: #control to check when the total functioning time is reached. It will be typically overcome, so a correction is applied to avoid this
-                                    indexes_adj = indexes[:-(tot_time-rand_time)] #correctes indexes size to avoid overcoming total time
-                                    if np.in1d(peak_time_range,indexes_adj).any() and App.fixed == 'no': #check if indexes are in peak window and if the coincident behaviour is locked by the "fixed" attribute
-                                        coincidence = min(App.number,max(1,math.ceil(random.gauss((App.number*mu_peak+0.5),(s_peak*App.number*mu_peak))))) #calculates coincident behaviour within the peak time range
-                                    elif np.in1d(peak_time_range,indexes_adj).any()== False and App.fixed == 'no': #check if indexes are off-peak and if coincident behaviour is locked or not
-                                        Prob = random.uniform(0,(App.number-op_factor)/App.number) #calculates probability of coincident switch_ons off-peak
-                                        array = np.arange(0,App.number)/App.number
-                                        try:
-                                            on_number = np.max(np.where(Prob>=array))+1
-                                        except ValueError:
-                                            on_number = 1 
-                                        coincidence = on_number #randomly selects how many apps are on at the same time for each app type based on the above probabilistic algorithm
-                                    else:
-                                        coincidence = App.number #this is the case when App.fixed is activated. All 'n' apps of an App instance are switched_on altogether
-                                    if App.fixed_cycle > 0: #evaluates if the app has some duty cycles to be considered
-                                        if indexes_adj.size > 0:
-                                            evaluate = round(np.mean(indexes_adj)) #calculates the mean time position of the current switch_on event, to later select the proper duty cycle
-                                        else:
-                                            evaluate = 0 
-                                        #based on the evaluate value, selects the proper duty cycle and puts the corresponding power values in the indexes range
-                                        if evaluate in range(App.cw11[0],App.cw11[1]) or evaluate in range(App.cw12[0],App.cw12[1]):
-                                            np.put(App.daily_use,indexes_adj,(random_cycle1*coincidence))
-                                            np.put(App.daily_use_masked,indexes_adj,(random_cycle1*coincidence),mode='clip')
-                                        elif evaluate in range(App.cw21[0],App.cw21[1]) or evaluate in range(App.cw22[0],App.cw22[1]):
-                                            np.put(App.daily_use,indexes_adj,(random_cycle2*coincidence))
-                                            np.put(App.daily_use_masked,indexes_adj,(random_cycle2*coincidence),mode='clip')
-                                        else:
-                                            np.put(App.daily_use,indexes_adj,(random_cycle3*coincidence))
-                                            np.put(App.daily_use_masked,indexes_adj,(random_cycle3*coincidence),mode='clip')
-                                    else: #if no duty cycles are specififed, a regular switch_on event is modelled
-                                        np.put(App.daily_use,indexes_adj,(App.power[prof_i]*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var)))*coincidence)) #randomises also the App Power if thermal_p_var is on
-                                        np.put(App.daily_use_masked,indexes_adj,(App.power[prof_i]*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var)))*coincidence),mode='clip')
-                                    App.daily_use_masked = np.zeros_like(ma.masked_greater_equal(App.daily_use_masked,0.001)) #updates the mask excluding the current switch_on event to identify the free_spots for the next iteration
-                                    tot_time = (tot_time - indexes.size) + indexes_adj.size #updates the total time correcting the previous value
-                                    break #exit cycle and go to next App
-                                else: #if the tot_time has not yet exceeded the App total functioning time, the cycle does the same without applying corrections to indexes size
-                                    if np.in1d(peak_time_range,indexes).any() and App.fixed == 'no':
-                                        coincidence = min(App.number,max(1,math.ceil(random.gauss((App.number*mu_peak+0.5),(s_peak*App.number*mu_peak)))))
-                                    elif np.in1d(peak_time_range,indexes).any() == False and App.fixed == 'no':
-                                        Prob = random.uniform(0,(App.number-op_factor)/App.number)
-                                        array = np.arange(0,App.number)/App.number
-                                        try:
-                                            on_number = np.max(np.where(Prob>=array))+1
-                                        except ValueError:
-                                            on_number = 1
-                                        coincidence = on_number
-                                    else:
-                                        coincidence = App.number
-                                    if App.fixed_cycle > 0:
-                                        if indexes.size > 0:
-                                            evaluate = round(np.mean(indexes))
-                                        else:
-                                            evaluate = 0
-                                        if evaluate in range(App.cw11[0],App.cw11[1]) or evaluate in range(App.cw12[0],App.cw12[1]):
-                                            np.put(App.daily_use,indexes,(random_cycle1*coincidence))
-                                            np.put(App.daily_use_masked,indexes,(random_cycle1*coincidence),mode='clip')
-                                        elif evaluate in range(App.cw21[0],App.cw21[1]) or evaluate in range(App.cw22[0],App.cw22[1]):
-                                            np.put(App.daily_use,indexes,(random_cycle2*coincidence))
-                                            np.put(App.daily_use_masked,indexes,(random_cycle2*coincidence),mode='clip')
-                                        else:
-                                            np.put(App.daily_use,indexes,(random_cycle3*coincidence))
-                                            np.put(App.daily_use_masked,indexes,(random_cycle3*coincidence),mode='clip')
-                                    else:
-                                        np.put(App.daily_use,indexes,(App.power[prof_i]*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var)))*coincidence))
-                                        np.put(App.daily_use_masked,indexes,(App.power[prof_i]*(random.uniform((1-App.thermal_p_var),(1+App.thermal_p_var)))*coincidence),mode='clip')
-                                    App.daily_use_masked = np.zeros_like(ma.masked_greater_equal(App.daily_use_masked,0.001))
-                                    tot_time = tot_time #no correction applied to previously calculated value
-                                                    
-                                free_spots = [] #calculate how many free spots remain for further switch_ons
-                                try:
-                                    for j in ma.notmasked_contiguous(App.daily_use_masked):
-                                        free_spots.append(j.stop-j.start)
-                                except TypeError:
-                                    free_spots = [0]
-                                max_free_spot = max(free_spots) 
-    
-                            else:
-                                continue #if the random switch_on falls somewhere where the App has been already turned on, tries again from beginning of the while cycle
-                    Us.load = Us.load + App.daily_use #adds the App profile to the User load
-            Tot_Classes = Tot_Classes + Us.load #adds the User load to the total load of all User classes
-        Profile.append(Tot_Classes) #appends the total load to the list that will contain all the generated profiles
-        print('Profile',prof_i+1,'/',num_profiles,'completed') #screen update about progress of computation
-    return(Profile)
+    Parameters
+    ----------
+    user_list: list
+        list containing all the user types
+    peak_enlarge: float
+        percentage random enlargement or reduction of peak time range length
+        corresponds to \delta_{peak} in [1], p.7
+
+    Notes
+    -----
+    [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+        Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+        Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+
+    Returns
+    -------
+    peak time range: numpy array
+    """
+
+    tot_max_profile = np.zeros(1440)  # creates an empty daily profile
+    # Aggregate each User's theoretical max profile to the total theoretical max
+    for user in user_list:
+        tot_max_profile = tot_max_profile + user.maximum_profile
+    # Find the peak window within the theoretical max profile
+    peak_window = np.squeeze(np.argwhere(tot_max_profile == np.amax(tot_max_profile)))
+    # Within the peak_window, randomly calculate the peak_time using a gaussian distribution
+    peak_time = round(random.normalvariate(
+        mu=round(np.average(peak_window)),
+        sigma=1 / 3 * (peak_window[-1] - peak_window[0])
+    ))
+    rand_peak_enlarge = round(math.fabs(peak_time - random.gauss(mu=peak_time, sigma=peak_enlarge * peak_time)))
+    # The peak_time is randomly enlarged based on the calibration parameter peak_enlarge
+    return np.arange(peak_time - rand_peak_enlarge, peak_time + rand_peak_enlarge)
+
+
+def stochastic_process(j=None, fname=None, num_profiles=None):
+    """Generate num_profiles load profile for the usecase
+
+        Covers steps 1. and 2. of the algorithm described in [1], p.6-7
+
+    Notes
+    -----
+    [1] F. Lombardi, S. Balderrama, S. Quoilin, E. Colombo,
+        Generating high-resolution multi-energy load profiles for remote areas with an open-source stochastic model,
+        Energy, 2019, https://doi.org/10.1016/j.energy.2019.04.097.
+    """
+    # creates an empty list to store the results of each code run, i.e. each stochastically generated profile
+    profiles = []
+
+    peak_enlarge, year_behaviour, user_list, num_profiles = initialise_inputs(j, fname, num_profiles)
+
+    # Calculation of the peak time range, which is used to discriminate between off-peak
+    # and on-peak coincident switch-on probability, corresponds to step 1. of [1], p.6
+    peak_time_range = calc_peak_time_range(user_list, peak_enlarge)
+
+    for prof_i in range(num_profiles):
+        # initialise an empty daily profile (or profile load)
+        # that will be filled with the sum of the daily profiles of each User instance
+        usecase_load = np.zeros(1440)
+        # for each User instance generate a load profile, iterating through all user of this instance and
+        # all appliances they own, corresponds to step 2. of [1], p.7
+        for user in user_list:
+            user.generate_aggregated_load_profile(prof_i, peak_time_range, year_behaviour)
+            # aggregate the user load to the usecase load
+            usecase_load = usecase_load + user.load
+        profiles.append(usecase_load)
+        # screen update about progress of computation
+        print('Profile', prof_i+1, '/', num_profiles, 'completed')
+    return(profiles)
+
diff --git a/ramp/core/utils.py b/ramp/core/utils.py
index 480b8ae6..4e232a6b 100644
--- a/ramp/core/utils.py
+++ b/ramp/core/utils.py
@@ -1,4 +1,6 @@
 import json
+import random
+import numpy as np
 import pandas as pd
 from openpyxl import load_workbook
 from openpyxl.worksheet.cell_range import CellRange
@@ -72,7 +74,6 @@ def read_input_file(filename):
                                 f"The provided power timeseries of the appliance '{appliance_name}' of user '{user_name}' in '{filename}' does not contain 365 values as expected\n{POSSIBLE_FORMATS}"
                             )
                         )
-
                     ts = v
                 except json.JSONDecodeError:
                     raise (
@@ -83,3 +84,82 @@ def read_input_file(filename):
             df.loc[i, "power"] = ts
             df.loc[i, "p_series"] = True
     return df
+
+
+def random_variation(var, norm=1):
+    """Pick a random variable within a uniform distribution of range [1-var, 1+var]
+
+    Parameters
+    ----------
+    var: float
+        sets the range of the uniform distribution around one
+    norm: float
+        multiplication factor of the random variable, default = 1
+
+    Returns
+    -------
+    random number close to norm
+    """
+    return norm * random.uniform((1 - var), (1 + var))
+
+
+def duty_cycle(var, t1, p1, t2, p2):
+    """Assign a two period duty cycle
+
+    concatenate an array where values equal p1 for a time (t1 +- random variation)
+    followed by values equal to p2 for a time (t2 +- random variation)
+
+    Parameters
+    ----------
+    var: float
+        sets the range of the uniform distribution around t1 and t2
+    t1: int
+        time interval of the first part of the duty cycle in minutes
+    p1: float
+        power of the first part of the duty cycle in Watt
+    t2: int
+        time interval of the second part of the duty cycle in minutes
+    p2: int
+        power of the second part of the duty cycle in Watt
+
+    Returns
+    -------
+    Power during each timestep of the duty cycle where p1 is repeated (t1 +- random variation) times and p2 is repeated (t2 +- random variation) times.
+    The duty cycle is implicitly sampled every minutes (which is the unit for t1 and t2)
+    """
+    return np.concatenate(
+        (
+            np.ones(int(random_variation(var=-var, norm=t1))) * p1,
+            np.ones(int(random_variation(var=-var, norm=t2))) * p2,
+        )
+    )
+
+
+def random_choice(var, t1, p1, t2, p2):
+    """Chooses one of two duty cycles randomly
+
+    The choice is between a normal duty cycle and a reversed duty cycle (where t1 is swapped with t2 and p1 with p2)
+
+    Parameters
+    ----------
+    var: float
+        sets the range of the uniform distribution around t1 and t2
+    t1: int
+        time interval of the first part of the duty cycle in minutes
+    p1: float
+        power of the first part of the duty cycle in Watt
+    t2: int
+        time interval of the second part of the duty cycle in minutes
+    p2: int
+        power of the second part of the duty cycle in Watt
+
+    Returns
+    -------
+    A duty cycle, see function duty_cycle
+    """
+    return random.choice(
+        [
+            duty_cycle(var, t1=t1, p1=p1, t2=t2, p2=p2),
+            duty_cycle(var, t1=t2, p1=p2, t2=t1, p2=p1),
+        ]
+    )
diff --git a/ramp/input_files/input_file_4.py b/ramp/input_files/input_file_4.py
new file mode 100644
index 00000000..d7462daa
--- /dev/null
+++ b/ramp/input_files/input_file_4.py
@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+
+#%% Definition of the inputs
+'''
+Input data definition 
+'''
+
+
+from ramp.core.core import User, np
+
+#number_of_app = int(np.random.normal(6,1))
+#power_of_app = int(np.random.normal(7,2))
+User_list = []
+
+'''
+This example input file represents an whole village-scale community,
+adapted from the data used for the Journal publication. It should provide a 
+complete guidance to most of the possibilities ensured by RAMP for inputs definition,
+including specific modular duty cycles and cooking cycles. 
+For examples related to "thermal loads", see the "input_file_2".
+'''
+
+#Create new user classes
+HI = User("high income",50)
+User_list.append(HI)
+
+#High-Income
+HI_indoor_bulb = HI.Appliance(HI,n=6, P=7, w=2, t = 120,r_t = 0.2, c = 10)
+HI_indoor_bulb.windows([1170,1440],[0,30],0.35)
+
+HI_Freezer2 = HI.Appliance(HI,1,200,1,1440,0,30,'yes',3)
+HI_Freezer2.windows([0,1440],[0,0])
+HI_Freezer2.specific_cycle_1(200,20,5,10)
+HI_Freezer2.specific_cycle_2(200,15,5,15)
+HI_Freezer2.specific_cycle_3(200,10,5,20)
+HI_Freezer2.cycle_behaviour([480,1200],[0,0],[300,479],[0,0],[0,299],[1201,1440])
+
+HI_Mixer = HI.Appliance(HI,1,50,3,30,0.1,1,occasional_use = 0.33)
+HI_Mixer.windows([420,480],[660,750],0.35,[1140,1200])
+
+
diff --git a/ramp/ramp_convert_old_input_files.py b/ramp/ramp_convert_old_input_files.py
index d8553db9..e641784d 100644
--- a/ramp/ramp_convert_old_input_files.py
+++ b/ramp/ramp_convert_old_input_files.py
@@ -71,9 +71,9 @@ def convert_old_user_input_file(
         fname = f"ramp.{fname}"
     file_module = importlib.import_module(fname)
 
-    User_list = file_module.User_list
+    user_list = file_module.User_list
 
-    UseCase(users=User_list).save(output_fname)
+    UseCase(users=user_list).save(output_fname)
 
 
 if __name__ == "__main__":
diff --git a/ramp/ramp_run.py b/ramp/ramp_run.py
index 37b0c5bf..e979e3e2 100644
--- a/ramp/ramp_run.py
+++ b/ramp/ramp_run.py
@@ -27,16 +27,16 @@
 sys.path.append('../')
 
 try:
-    from .core.stochastic_process import Stochastic_Process
+    from .core.stochastic_process import stochastic_process
     from .post_process import post_process as pp
 except ImportError:
-    from core.stochastic_process import Stochastic_Process
+    from core.stochastic_process import stochastic_process
     from post_process import post_process as pp
 
 
 def run_usecase(j=None, fname=None, num_profiles=None):
     # Calls the stochastic process and saves the result in a list of stochastic profiles
-    Profiles_list = Stochastic_Process(j=j, fname=fname, num_profiles=num_profiles)
+    Profiles_list = stochastic_process(j=j, fname=fname, num_profiles=num_profiles)
 
     # Post-processes the results and generates plots
     Profiles_avg, Profiles_list_kW, Profiles_series = pp.Profile_formatting(Profiles_list)
@@ -55,7 +55,3 @@ def run_usecase(j=None, fname=None, num_profiles=None):
 
     for i, j in enumerate(input_files_to_run):
         run_usecase(j=j)
-
-
-
-
diff --git a/tests/test_input_file_conversion.py b/tests/test_input_file_conversion.py
index 0c04a2e3..f20f4407 100644
--- a/tests/test_input_file_conversion.py
+++ b/tests/test_input_file_conversion.py
@@ -3,20 +3,16 @@
 import numpy as np
 
 from ramp.core.core import User, Appliance
-from ramp.core.initialise import Initialise_inputs
+from ramp.core.initialise import initialise_inputs
+e
 from ramp.ramp_convert_old_input_files import convert_old_user_input_file
 
 
 def load_usecase(j=None, fname=None):
-    (
-        peak_enlarg,
-        mu_peak,
-        s_peak,
-        op_factor,
-        Year_behaviour,
-        User_list,
-    ) = Initialise_inputs(j, fname)
-    return User_list
+    peak_enlarge, year_behaviour, user_list, num_profiles = initialise_inputs(
+        j, fname, num_profiles=1
+    )
+    return user_list
 
 
 class TestConversion: