# -*- coding: utf-8 -*-
"""
INITIAL COMMENTS FROM PATRIZIO (may not be up to date)
Created on Mon Mar  4 11:22:27 2019

@author: mangan36
In this version of the framework:
 - Subfunctions save to main variables
 - Thermal network has been updated (from 3-rd order max to static, via a Nyquist limit analysis)
     - Thermal predictor reduction more accurate
 - Wind model based on the work of Imre is considered. To used LUT-based model:
     - Uncomment lines that run this model in "system_solver_"
     - Comment lines that run the model of Imre in "system_solver_"
 - Net Heat is calculated from Basic Component "ideal" power (meaning without considering Rs losses). To change this:
     - Uncomment lines that activate the inclusion of Rs in "system_solver_"
 - The time vector is assumed to have constant sampling. However, the code can be easily adapted to run with a variable sampling time:
     - In "system_solver_", the last input of "Temperature_Prediction_3rd_order_Simplified", namely "Ts_DT", must be adapted at runtime to be equal to the current timestep (time_vector[sim_index]-time_vector[sim_index-1])
     - In "main", we create the predictor using the function "temperature_predictor_variable_order", that needs as input the sampling time "Ts_DT".
         - In order to properly select the thermal network order (and fulfill with Nyquist theory), the minimum input timestep must be used as input for "temperature_predictor_variable_order"
         - Small difference in sampling time do not create issues. However, if timestep changed significantly (e.g. from 1 sec up to 5 min), some errors and/or accuracy losses might happen

N.B. TO RUN PARALLEL, first of all clusters must be started in the anaconda prompt running the command
        ipcluster start -n 4 (where 4 means the number of cores we want to use)
    Also, parallel computation speed has not been fully tested, and it is probably limited by the fact that subfunctions save to main variables
"""

"""
Issues and Possible solutions:
    A) Simulation is slower than in Matlab, mainly because of:
        1. thermal coefficient run-time calculation: "lambdify" creates lambda functions, that are time-consuming to be solved:
            - coefficients should be create as standard functions (as done in matlab by means of the command "matlabFunction")
        2. interpolation is done with "interp1d". It allows extrapolation, needed to avoid big errors, but it is slow:
            - numpy interp is much faster, but do not allow extrapolation. An optimized interpolating function might be programmed
        3. code of wind model of Imre can probably be optimized for higher computation speed
    B) When compared with Spectre:
        1. If wind is not included:
            - Spectre and Matlab/Python results are extremely similar, both in terms of dynamics (no delay or significant difference in behavior) and absolute values:
                - Modeling radiation loss as variable resistance does not induce error
                - The thermal network is solved properly
        2. If wind is included (as LUT-based thermal resistors, since this is the only model implemented in Spectre):
            - Second resolution: results are quite similar
            - Minute resolution or above: some fast dynamics due to wind are obviously lost. It is important to understand:
                - if this spectre-simulated results reflect the wind effect realistically
                - if a different approach (e.g. considering wind as a "dependant source" instead of a variable resistance) leads to more realistic results (or at least to a lower error due to non-simulated dynamics):
                     - Validation with measurements instead of comparison with Spectre is needed for this!
                     - The simulation approach can be adapted accordingly to any new development of the wind cooling model
    C) Wind model of Imre needs to be fixed: sometimes it gives the error "invalid value encountered in double_scalars"
    D) Parallel processing would be probably much faster if sub-function do not save to main variable (each writing slow down computation)?
"""

import os
import multiprocess
from multiprocess import Pool

import numpy as np
import time
from scipy.interpolate import interp1d

import matplotlib.pyplot as plt
import scipy.io

from core.climatic_data import *
from core.module import *
from core.energy_yield_res import *
from core.wind_model import wind_thermal_res_front_and_back
from core.interpolate import interpolate_row_by_row
import numba
import platform

class EnergyYieldSimParams:
    """
    This class encompasses parameters that control the energy yield
    computation.
    """

    def __init__(self, nbr_of_cores=1, nbr_points_per_IV_curve=100,
    steps_accuracy=3, first_day_to_simulate=1, min_irradiance=0, debug_mode=False,
    debug_dir="debug"):
        # How many point do we want to have per IV curve
        self.nbr_points_per_IV_curve = nbr_points_per_IV_curve
        # A minimum  value of 2 is recommended
        self.steps_accuracy = steps_accuracy
        self.first_day_to_simulate = first_day_to_simulate
        self.nbr_of_cores = nbr_of_cores
        #Will impact the timestamps at which simulation start (i.e. there is no
        #point to start simulation too soon before sunrise until too far after
        #sunset)
        self.min_irradiance = min_irradiance
        self.debug = DebugEnergyYield(debug_mode, debug_dir)

class DebugEnergyYield:

    def __init__(self, debug_mode=False, debug_dir="debug"):

        self.debug_mode = debug_mode
        if self.debug_mode:
            self.debug_dir = debug_dir

            if not os.path.exists(self.debug_dir):
                try:
                    os.makedirs(self.debug_dir)
                except OSError as exc:
                    if exc.errno != errno.EEXIST:
                        raise

            self.default_name_debug_file_count = 0

    def write_debug_mat_file(self, names, values, file_name=None):

        assert(self.debug_mode)

        if file_name is None:
            file_name = "proc_id_"+str(os.getpid())+"_"+\
            str(self.default_name_debug_file_count)
            self.default_name_debug_file_count = \
            self.default_name_debug_file_count +1

        debug_file_path = self.debug_dir+"/"+file_name+".mat"
        dict_to_save = dict(zip(names, values))
        scipy.io.savemat(debug_file_path, dict_to_save)

def radiation_thermal_res(T_surface, T_amb, surface_amb_T_diff, cell_area,
horiz_surface_temp_diff, accuracy_dep_idx, day_idx, sim_idx, emissivity):

    sigma = Constants.STEFAN_BOLTZMANN
    cell_area_msquared = cell_area/10**4

    if type(horiz_surface_temp_diff) == np.ndarray:
        horiz_surface_temp_diff = horiz_surface_temp_diff[day_idx, sim_idx]

    horiz_surface_amb_T_diff = T_amb[day_idx][sim_idx]-horiz_surface_temp_diff

    return surface_amb_T_diff/(cell_area_msquared*sigma*emissivity*( \
    T_surface[:, accuracy_dep_idx]**4-horiz_surface_amb_T_diff**4))

def radiation_thermal_res_front_and_back(T_surface_front, T_surface_back, T_amb,
glass_amb_delta_T, backsheet_amb_delta_T, cell, sky_temp_offset,
ground_temp_offset, accuracy_dep_idx, day_idx, sim_idx):

    rad_therm_res_front = radiation_thermal_res(T_surface_front, T_amb,
    glass_amb_delta_T, cell.area, sky_temp_offset, accuracy_dep_idx, day_idx,
    sim_idx, cell.front_glass_emissivity)
    rad_therm_res_back = radiation_thermal_res(T_surface_back, T_amb,
    backsheet_amb_delta_T, cell.area, ground_temp_offset, accuracy_dep_idx,
    day_idx, sim_idx, cell.back_sheet_emissivity)

    return (rad_therm_res_front, rad_therm_res_back)

def thermal_res_front_and_back(T_amb, T_glass, T_backsheet, wind_speed_front,
wind_speed_back, wind_direction, day_idx, accuracy_dep_idx,
sim_idx, string_of_modules, sky_temp_offset, ground_temp_offset):

    cell = string_of_modules.module.cell

    glass_amb_delta_T = T_glass[:, accuracy_dep_idx]-T_amb[day_idx, sim_idx]
    backsheet_amb_delta_T = T_backsheet[:, accuracy_dep_idx]- \
    T_amb[day_idx, sim_idx]

    R_wind_front, R_wind_back = wind_thermal_res_front_and_back(
    glass_amb_delta_T, backsheet_amb_delta_T, wind_speed_front, wind_speed_back,
    wind_direction, day_idx, sim_idx, string_of_modules)

    R_rad_front, R_rad_back = \
    radiation_thermal_res_front_and_back(T_glass, T_backsheet, T_amb,
    glass_amb_delta_T, backsheet_amb_delta_T, cell, sky_temp_offset,
    ground_temp_offset, accuracy_dep_idx, day_idx, sim_idx)

    def tot_thermal_res(R_wind, R_rad):

        R = np.zeros(R_wind.shape)

        R_wind_gt_zero = R_wind > 0
        R_wind_le_zero = np.logical_not(R_wind_gt_zero)
        R_rad_gt_zero = R_rad > 0
        R_rad_le_zero = np.logical_not(R_rad_gt_zero)

        R_wind_rad_gt_zero = np.logical_and(R_wind_gt_zero, R_rad_gt_zero)
        R_wind_rad_le_zero = np.logical_and(R_wind_le_zero, R_rad_le_zero)
        only_R_wind_gt_zero = np.logical_and(R_wind_gt_zero, R_rad_le_zero)
        only_R_rad_gt_zero = np.logical_and(R_wind_le_zero, R_rad_gt_zero)

        R[R_wind_rad_gt_zero] = (R_wind*R_rad/(R_wind+R_rad))[R_wind_rad_gt_zero]
        R[only_R_wind_gt_zero] = R_wind[only_R_wind_gt_zero]
        R[only_R_rad_gt_zero] = R_rad[only_R_rad_gt_zero]
        R[R_wind_rad_le_zero] = 10**6

        return R

    R_front_num = tot_thermal_res(R_wind_front, R_rad_front)
    R_back_num = tot_thermal_res(R_wind_back, R_rad_back)

    return (R_front_num, R_back_num)

#@numba.jit
def energy_yield_for_one_day(climatic_data, string_of_modules, index_min,
index_max, steps_accuracy, nbr_cells_in_basic_comp, sky_temp_offset,
ground_temp_offset, time_step, thermal_predictor, Heat_Irr_mat, Iph_mat,
Vdiode_BC, R_in, R_out, debug, day_index):

    irradiance = climatic_data.global_irradiance
    T_amb = climatic_data.ambient_temperature
    wind_direction = climatic_data.wind_direction
    wind_speed_back = climatic_data.wind_speed_back
    wind_speed_front = climatic_data.wind_speed_front

    nbr_of_basic_components = irradiance.shape[1]
    nbr_of_timestamps_per_day = irradiance.shape[2]

    cell = string_of_modules.module.cell
    cell_elec = cell.electrical_prop

    I_BC_mat = np.zeros((nbr_of_basic_components, len(Vdiode_BC)))
    Pmpp_system = np.zeros(nbr_of_timestamps_per_day)
    Vmpp_system = np.zeros(nbr_of_timestamps_per_day)
    Pmpp_BC = np.zeros((nbr_of_basic_components, nbr_of_timestamps_per_day))

    T_amb_at_day = T_amb[day_index]
    Tcell_mat = np.reshape(T_amb_at_day, (1, nbr_of_timestamps_per_day))
    Tcell_mat = np.repeat(Tcell_mat, nbr_of_basic_components, axis=0)
    Tglass_mat = Tcell_mat.copy()
    Tbacksheet_mat = Tcell_mat.copy()
    Net_Heat = Heat_Irr_mat[day_index].copy()

    debug_names = []
    debug_values = []

    for sim_index in range(int(index_min[day_index]),int(index_max[day_index]+1)):

        for index_accuracy in range(steps_accuracy):
            if index_accuracy == 0:
                accuracy_dep_index = sim_index-1
            else:
                accuracy_dep_index = sim_index

            R_front_num, R_back_num = thermal_res_front_and_back(T_amb,
            Tglass_mat, Tbacksheet_mat, wind_speed_front, wind_speed_back,
            wind_direction, day_index, accuracy_dep_index, sim_index,
            string_of_modules, sky_temp_offset, ground_temp_offset)

            Ppv_Pred_for_NetHeat = Pmpp_BC[:, accuracy_dep_index]/nbr_cells_in_basic_comp

            temp_prediction_3rd_order_with_back_rad_loss(
            Ppv_Pred_for_NetHeat,R_front_num,R_back_num,day_index,
            sim_index,time_step, thermal_predictor, Heat_Irr_mat, Net_Heat,
            T_amb, Tcell_mat, Tglass_mat, Tbacksheet_mat)

            I_BC_mat = BC_IV((Iph_mat[day_index, :, sim_index])[:, np.newaxis],
            (Tcell_mat[:, sim_index])[:, np.newaxis], cell_elec,
            nbr_cells_in_basic_comp, Vdiode_BC)

            system_obj(R_in, R_out, Vdiode_BC, I_BC_mat, day_index, sim_index,
            Pmpp_system, Vmpp_system, Pmpp_BC, nbr_of_basic_components, debug,
            debug_names, debug_values)

            # Uncomment next line to include Rs in Net Heat calculation
            #Pmpp_BC[day_index, :, sim_index] = Pmpp_BC[day_index, :, sim_index]\
            #- Rs*nbr_cells_in_basic_comp*(Pmpp_system[day_index][sim_index]/Vmpp_system[day_index][sim_index])**2
            Net_Heat[:, sim_index] =  Heat_Irr_mat[day_index,: , sim_index] - \
            Pmpp_BC[:, sim_index]/nbr_cells_in_basic_comp

            if debug.debug_mode:
                debug_names.extend(["index_accuracy_"+str(day_index)+"_"+str(sim_index),
                "R_front_num_"+str(day_index)+"_"+str(sim_index),
                "R_back_num_"+str(day_index)+"_"+str(sim_index),
                "Ppv_Pred_for_NetHeat_"+str(day_index)+"_"+str(sim_index),
                "I_BC_mat_"+str(day_index)+"_"+str(sim_index)])
                debug_values.extend([index_accuracy, R_front_num,
                R_back_num, Ppv_Pred_for_NetHeat, I_BC_mat])

    if debug.debug_mode:
        debug_names.extend(["Tglass_mat_"+str(day_index), "Tbacksheet_mat_"+str(day_index)])
        debug_values.extend([Tglass_mat, Tbacksheet_mat])
        file_name = "sim_loop_"+str(day_index)
        debug.write_debug_mat_file(debug_names, debug_values, file_name)
    return (Pmpp_system, Vmpp_system, Pmpp_BC, Tcell_mat, Net_Heat)

    #print('day #%s done' %(day_index+1))

def system_obj(Rin, Rout, voltage, current, day_index, sim_index, Pmpp_system,
Vmpp_system, Pmpp_BC, nbr_of_basic_components, debug, debug_names, debug_values):

    # Series connection of BCs' I-V curves
    voltages_at_common_currents, common_currents = IV_series_connection(
    voltage, current)
    voltage_tot = np.sum(voltages_at_common_currents, axis=0)

    # Adding total series resistance
    voltage_tot_real = IV_real(Rin, Rout, voltage_tot, common_currents)
    Pmpp = voltage_tot_real*common_currents
    index_max = np.argmax(Pmpp) # find MPP index
    Pmpp_system[sim_index] = Pmpp[index_max] # Calculate MPP Power
    Vmpp_system[sim_index] = voltage_tot_real[index_max] # Calculate MPP voltage

    #find operating voltage
    Vmax_BC = voltages_at_common_currents[:, index_max]
    Imax_BC = common_currents[index_max]
    Pmpp_BC[:, sim_index] = Vmax_BC*Imax_BC

    if debug.debug_mode:
        debug_names.extend(
        ["voltages_at_common_currents_"+str(day_index)+"_"+str(sim_index),
        "common_currents_"+str(day_index)+"_"+str(sim_index),
        "voltage_tot_"+str(day_index)+"_"+str(sim_index),
        "voltage_tot_real_"+str(day_index)+"_"+str(sim_index),
        "Pmpp_"+str(day_index)+"_"+str(sim_index),
        "Pmpp_max_"+str(day_index)+"_"+str(sim_index),
        "Vmpp_max_"+str(day_index)+"_"+str(sim_index),
        "Vmax_BC_"+str(day_index)+"_"+str(sim_index),
        "Imax_BC_"+str(day_index)+"_"+str(sim_index),
        "Pmpp_BC_"+str(day_index)+"_"+str(sim_index)])
        debug_values.extend([voltages_at_common_currents, common_currents,
        voltage_tot, voltage_tot_real, Pmpp, Pmpp_system[sim_index],
        Vmpp_system[sim_index], Vmax_BC, Imax_BC, Pmpp_BC[:, sim_index]])


def comp_Vdiode_BC(cell_elec, sim_params, Iph_STC, nbr_cells_in_basic_comp,
nbr_basic_components):
    """
    Basic Component Diode Voltage for I-V curve calculation
    """

    #Maximum voltage at which Basic Component's IV curve is evaluated (without Rs) [V]
    Voc_BC_max = cell_elec.V_max_open_circuit*nbr_cells_in_basic_comp
    nbr_points_per_IV_curve = sim_params.nbr_points_per_IV_curve
    Rsh = cell_elec.R_shunt
    Rs = cell_elec.R_series
    Vdiode_BC = np.append(-(Rsh+Rs)*nbr_cells_in_basic_comp*2.5*Iph_STC,
    np.linspace(0,Voc_BC_max+Rs*nbr_cells_in_basic_comp*Iph_STC, nbr_points_per_IV_curve))
    Vdiode_BC = np.reshape(Vdiode_BC, (1, Vdiode_BC.size))
    Vdiode_BC = np.repeat(Vdiode_BC, nbr_basic_components, axis=0)
    return Vdiode_BC

def temp_prediction_3rd_order_with_back_rad_loss(Ppv_Pred_for_NetHeat,
R_front_num,R_back_num,day_index, sim_index,Ts, thermal_predictor, Heat_Irr_mat,
Net_Heat, Tamb_mat, Tcell_mat, Tglass_mat, Tbacksheet_mat):

    thermal_coefs = np.array(list(map(
    lambda therm_coef_fun: therm_coef_fun(R_front_num, R_back_num,Ts),
    thermal_predictor)))
    T_cell_coefs = thermal_coefs[0:4]
    T_front_coefs = thermal_coefs[4:8]
    T_back_coefs = thermal_coefs[8:12]
    denom_coefs = thermal_coefs[12:15]

    prev_sim_idx = sim_index-np.arange(1,4)
    heat_irr_minus_ppv_pred_for_net_heat = Heat_Irr_mat[day_index,:,sim_index]-\
    Ppv_Pred_for_NetHeat

    Tcell_mat[:, sim_index] = predict_temperature(T_cell_coefs, denom_coefs,
    heat_irr_minus_ppv_pred_for_net_heat, Net_Heat, Tamb_mat, Tcell_mat,
    day_index, sim_index, prev_sim_idx)
    Tglass_mat[:, sim_index] = predict_temperature(T_front_coefs, denom_coefs,
    heat_irr_minus_ppv_pred_for_net_heat, Net_Heat, Tamb_mat, Tglass_mat,
    day_index, sim_index, prev_sim_idx)
    Tbacksheet_mat[:, sim_index] = predict_temperature(T_back_coefs, denom_coefs,
    heat_irr_minus_ppv_pred_for_net_heat, Net_Heat, Tamb_mat, Tbacksheet_mat,
    day_index, sim_index, prev_sim_idx)

def predict_temperature(coefs, coefs_denom, heat_irr_minus_ppv_pred_for_net_heat,
Net_Heat, Tamb_mat, prev_temp, day_index, sim_index, prev_sim_idx):
    delta_T_amb_prev = Tamb_mat[day_index, prev_sim_idx]-prev_temp[:, prev_sim_idx]
    res = coefs[0]*heat_irr_minus_ppv_pred_for_net_heat+Tamb_mat[day_index, sim_index]
    for i in range(0,3):
        res = res+coefs[i+1]*Net_Heat[:, prev_sim_idx[i]]+coefs_denom[i]*delta_T_amb_prev[:,i]
    return res

def sky_ground_temp_offsets(string_of_modules):

    tilt_angle = string_of_modules.tilt_angle
    cell = string_of_modules.module.cell
    #Sky temperature offset values, weighted by solid angle of sky.
    #Ground is assumed to be at ambient temperature
    offset_front = cell.sky_temp_horiz_surface_temp_diff*\
    (1+np.cos(np.deg2rad(tilt_angle)))/2
    offset_back = cell.sky_temp_horiz_surface_temp_diff*\
    (1+np.cos(np.deg2rad(180-tilt_angle)))/2
    return (offset_front, offset_back)

def dispatch_work_among_proc(climatic_data, sim_params, timestamp_index_min,
timestamp_index_max, nbr_cells_in_basic_comp, string_of_modules,
real_sim_nbr_bc_ratio, nbr_sim_modules, debug):
    #dill.settings['recurse'] = True

    irradiance = climatic_data.global_irradiance
    time_vector = climatic_data.time_vector

    assert(time_vector.shape[0] > 1)

    steps_accuracy = sim_params.steps_accuracy
    first_day_to_simulate = sim_params.first_day_to_simulate
    nbr_days_to_simulate = irradiance.shape[0]
    nbr_basic_components = irradiance.shape[1]
    nbr_timestamps_per_day = irradiance.shape[2]

    module = string_of_modules.module
    cell = module.cell
    cell_elec = cell.electrical_prop

    Iph_STC = cell_elec.I_photo_generated
    Vdiode_BC = comp_Vdiode_BC(cell_elec, sim_params, Iph_STC,
    nbr_cells_in_basic_comp, nbr_basic_components)

    #Total heat as a function of irradiation
    Heat_STC = cell.heat_generated_at_standard_cond
    Heat_Irr_mat = Heat_STC*irradiance/1000
    #Iph as a function of irradiation (no thermal dependency assumed)
    Iph_mat = Iph_STC*irradiance/1000
    sky_temp_offset, ground_temp_offset = sky_ground_temp_offsets(
    string_of_modules)

    # Sampling time for temperature prediction = sampling time of input
    time_step = time_vector[1]-time_vector[0]
    #thermal_predictor = temperature_predictor_variable_order(time_step,
    #cell.thermal_prop)
    thermal_predictor = \
    cell.thermal_prop.temperature_predictor_variable_order(time_step)

    day_indices = list(np.arange(first_day_to_simulate-1,
    first_day_to_simulate+nbr_days_to_simulate-1))

    nbr_cores =  sim_params.nbr_of_cores
    if nbr_cores is None:
        nbr_cores = min(multiprocess.cpu_count(), len(day_indices))

    #print("Spawning "+str(nbr_cores)+" processes for energy yield computation.")

    R_cabling_tot = string_of_modules.R_cabling_per_module*nbr_sim_modules

    Rs = cell_elec.R_series
    R_in, R_out = Rs*nbr_cells_in_basic_comp*nbr_basic_components/2+R_cabling_tot

    #Useful when multiprocessing done at a "upper level", because
    #daemonic processes are not allowed to have children (even one.)
    if nbr_cores == 1:
        results = list(map(functools.partial(
        energy_yield_for_one_day, climatic_data, string_of_modules,
        timestamp_index_min, timestamp_index_max, steps_accuracy,
        nbr_cells_in_basic_comp, sky_temp_offset, ground_temp_offset, time_step,
        thermal_predictor, Heat_Irr_mat, Iph_mat, Vdiode_BC, R_in, R_out, debug),
        day_indices))
    else:
        proc_pool = Pool(nbr_cores)
        results = list(proc_pool.map(functools.partial(
        energy_yield_for_one_day, climatic_data, string_of_modules,
        timestamp_index_min, timestamp_index_max, steps_accuracy,
        nbr_cells_in_basic_comp, sky_temp_offset, ground_temp_offset, time_step,
        thermal_predictor, Heat_Irr_mat, Iph_mat, Vdiode_BC, R_in, R_out, debug),
        day_indices))
        proc_pool.close()

    Pmpp_system = np.zeros((nbr_days_to_simulate, nbr_timestamps_per_day))
    Vmpp_system = np.zeros((nbr_days_to_simulate, nbr_timestamps_per_day))
    Pmpp_BC = np.zeros(irradiance.shape)
    T_cell = np.zeros(irradiance.shape)
    Net_Heat = np.zeros(irradiance.shape)

    for day_index in day_indices:
        day_res = results[day_index]
        Pmpp_system[day_index] = day_res[0]
        Vmpp_system[day_index] = day_res[1]
        Pmpp_BC[day_index] = day_res[2]
        T_cell[day_index] = day_res[3]
        Net_Heat[day_index] = day_res[4]

    Pmpp_system = Pmpp_system*real_sim_nbr_bc_ratio
    Vmpp_system = Vmpp_system*real_sim_nbr_bc_ratio

    wpeak_watt = comp_system_wpeak(cell, Vdiode_BC,
    climatic_data.basic_component, nbr_cells_in_basic_comp,
    nbr_basic_components, real_sim_nbr_bc_ratio)

    if debug.debug_mode:
        debug.write_debug_mat_file(["nbr_days_to_simulate", "nbr_basic_components",
        "nbr_timestamps_per_day", "Iph_STC", "Vdiode_BC", "Heat_STC", "Heat_Irr_mat",
        "Iph_mat", "sky_temp_offset", "ground_temp_offset", "time_step",
        "day_indices", "R_cabling_tot", "Rs", "R_in", "R_out", "wpeak_watt"],
        [nbr_days_to_simulate, nbr_basic_components, nbr_timestamps_per_day,
        Iph_STC, Vdiode_BC, Heat_STC, Heat_Irr_mat, Iph_mat, sky_temp_offset,
        ground_temp_offset, time_step, day_indices, R_cabling_tot, Rs, R_in,
        R_out, wpeak_watt])

    Pmpp_system = Pmpp_system*2 if isinstance(cell, HalfCellTsec) else Pmpp_system

    return (Pmpp_system, Vmpp_system, Pmpp_BC, T_cell, wpeak_watt, Net_Heat,
    Heat_Irr_mat)

#@numba.njit #slows down a bit
def IV_series_connection(voltage, current):
    """
    Series connection of I-V curves.
    N.B. When it comes to interpolation, numpy.interp is faster but cannot
    extrapolate, leading to big error
    """
    nbr_of_basic_components = current.shape[0]
    if nbr_of_basic_components == 1:
        return (voltage, current.flatten())
    else:
        #Sort current by increasing order (row by row) and change positions of
        #values in voltage accordingly
        sort_indices = np.argsort(current, axis=1)
        sort_current = np.take_along_axis(current, sort_indices, axis=1)
        voltage_sort_by_current = np.take_along_axis(voltage, sort_indices,
        axis=1)
        flat_current = sort_current.flatten()
        curr = np.unique(np.sort(flat_current))
        volt = np.zeros((nbr_of_basic_components, curr.size))
        #Initialize internal vector of interpolate here, because numba fails to
        #np.zeros and because it allows to decreases costly call to malloc
        #inside interpolate_row_by_row called many times.
        bound_right_indices = np.zeros(curr.size).astype(int)

        # Old interpolation code using scipy (much slower)
        # for i in range(voltage_sort_by_current.shape[0]):
        #     itp_volt = interp1d(sort_current[i,:], voltage_sort_by_current[i,:],
        #     kind='linear', fill_value="extrapolate", copy=False)
        #     volt[i,:] = itp_volt(curr)

        interpolate_row_by_row(sort_current, voltage_sort_by_current, curr,
        bound_right_indices, volt)
        return (volt, curr)

# Parallel connection of I-V curves:
def IV_parallel_connection(Voltage,Current):
    if len(Voltage.shape)==1:
        return(Voltage, sum(Current))
    else:
        volt = np.unique(np.sort(Voltage.flatten()))[::-1] # voltage vector for interpolation
        curr = np.zeros([Voltage.shape[0],len(volt)])
        for i in range(Voltage.shape[0]):
            itp_curr = interp1d(Voltage[i][:],Current[i][:], kind='linear',
            fill_value="extrapolate", copy=False) # create a scipy.interpolate.interp1d instance
            curr[i,:] = itp_curr(volt) # filling current matrix by interpolation
            #curr[i][:] = np.interp(volt,Voltage[i][:],Current[i][:])
        return(volt, sum(curr))

# Find operating voltage given the operating current
def find_operating_voltage(common_current, voltage_sort_by_current,
sort_current):

    itp_volt = interp1d(Current, Voltage, kind='linear',
    fill_value="extrapolate", copy=False) # create a scipy.interpolate.interp1d instance
    return itp_volt(common_current)

    # interpolate_row_by_row(sort_current, voltage_sort_by_current,
    # common_current, itp_volt)
    # itp_volt = interp1d(Current, Voltage, kind='linear',
    # fill_value="extrapolate", copy=False) # create a scipy.interpolate.interp1d instance
    # return itp_volt(common_current)

    #return (np.interp(common_current,Current,Voltage))

# Find operating current given the operating voltage
def find_operating_current(common_voltage,Voltage,Current):
    itp_curr = interp1d(Voltage,Current, kind='linear',fill_value="extrapolate",
    copy=False) # create a scipy.interpolate.interp1d instance
    return (itp_curr(common_voltage))
    #return (np.interp(common_voltage,Voltage,Current))

def diode_current(cell_elec, cell_temp, nbr_cells_in_basic_comp, Vdiode_BC):

    Rs = cell_elec.R_series
    Rsh = cell_elec.R_shunt
    C = cell_elec.saturation_cst
    n = cell_elec.ideality_factor
    k = Constants.ELECTRON_CHARGE_DIV_BY_BOLTZMANN_CONST
    a = cell_elec.a
    b = cell_elec.b

    return C*cell_temp**3*\
    np.exp(-k*(1.1557-(a*cell_temp**2)/(cell_temp+b))/cell_temp)* \
    (np.exp(k*Vdiode_BC/(n*nbr_cells_in_basic_comp*cell_temp))-1)

# I-V curve calculation
def BC_IV(Iph, cell_temp, cell_elec, nbr_cells_in_basic_comp, Vdiode_BC):
    Idiode = diode_current(cell_elec, cell_temp, nbr_cells_in_basic_comp, Vdiode_BC)
    return (Iph-Idiode-Vdiode_BC/(cell_elec.R_shunt*nbr_cells_in_basic_comp)) # BC current

# Adding series resistance to any I-V curve
def IV_real(Rin,Rout,Voltage,Current):
    return (Voltage-(Rin+Rout)*Current) # Voltage of the real I-V curve

## ###################################### AUXILIARY FUNCTIONS ##############################################################################################################################################################################################
# Check if a variable is not empty
def is_not_empty(any_structure):
    return any_structure

# Check if a variable is empty
def is_empty(any_structure):
    return not is_not_empty(any_structure)

def comp_nbr_cells_in_basic_comp(basic_component_type, module):
    """
    A BASIC COMPONENT IS THE COMPONENT FOR WHICH WE GENERATE THE INPUT. IF INPUT
    IS DEFINED PER CELL, CELL IS THE BASIC COMPONENT. IF INPUT IS DEFINED BY
    MODULE, MODULE IS THE BASIC COMPONENT. A BASIC COMPONENT IN THIS VERSION OF
    THE FRAMEWORK IS MADE OF num_cells SERIES-CONNECTED SOLAR CELLS
    """
    if basic_component_type == "cell":
        return 1
    elif basic_component_type == "module":
        return module.number_of_cells
    else:
        raise NotImplementedError()

def comp_system_wpeak(cell, Vdiode_BC, basic_component, nbr_cells_in_basic_comp,
nbr_basic_comp, real_sim_nbr_bc_ratio):

    cell_temp = 298.15
    cell_elec = cell.electrical_prop
    Rsh = cell_elec.R_shunt
    Iph_STC = cell_elec.I_photo_generated
    Rs = cell_elec.R_series

    Idiode = diode_current(cell_elec, cell_temp, nbr_cells_in_basic_comp, Vdiode_BC)
    I_basic_comp = (Iph_STC-Idiode-Vdiode_BC/(Rsh*nbr_cells_in_basic_comp)) # BC current

    V_basic_comp = Vdiode_BC-Rs*nbr_cells_in_basic_comp*I_basic_comp
    wpeak_watt_basic_comp = np.max(V_basic_comp*I_basic_comp)

    #if basic_component == "module":
    wpeak_watt = wpeak_watt_basic_comp*nbr_basic_comp
    if isinstance(cell, HalfCellTsec):
        wpeak_watt = wpeak_watt*2

    return wpeak_watt*real_sim_nbr_bc_ratio

def time_interval_with_irr_gt_min(sim_params, irradiances):
    """
    For each day and for each basic component find the indices of the timestamps
    that define the temporal interval during which the irradiance is superior to
    a minimum irradiance.
    """

    min_irradiance = sim_params.min_irradiance
    first_day_to_simulate = sim_params.first_day_to_simulate
    nbr_days_to_simulate = irradiances.shape[0]
    nbr_basic_components = irradiances.shape[1]

    index_min = np.zeros(nbr_days_to_simulate)
    index_max = np.zeros(nbr_days_to_simulate)
    index_min_temp = np.zeros((nbr_days_to_simulate, nbr_basic_components))
    index_max_temp = np.zeros((nbr_days_to_simulate, nbr_basic_components))

    for day_index in range(first_day_to_simulate-1,
    first_day_to_simulate + nbr_days_to_simulate-1):
        for BC_index in range(0, nbr_basic_components):

            buffer_condition = np.where(
            irradiances[day_index][BC_index][:] > min_irradiance)

            if buffer_condition[0].size > 0:
                index_min_temp[day_index][BC_index] = max([3, buffer_condition[0][0]])
                index_max_temp[day_index][BC_index] = buffer_condition[0][-1]
            else:
                index_min_temp[day_index][BC_index] = 10
                index_max_temp[day_index][BC_index] = 10

        index_min[day_index] = min(index_min_temp[day_index][:])
        index_max[day_index] = max(index_max_temp[day_index][:])

    return (index_min, index_max)

def comp_real_sim_nbr_bc_ratio(string_of_modules, climatic_data):
    """
    Compute ratio between simulated number of basic components and
    the number of basic components in the real solar system.
    """

    basic_component = climatic_data.basic_component
    nbr_basic_comp = climatic_data.nbr_of_components

    err_msg = "energy_yield.py, real_sim_nbr_bc_ratio: Number of \
    basic components "+str(nbr_basic_comp)+" does not match \
    string_of_modules.number_of_modules "+\
    str(string_of_modules.number_of_modules)

    if basic_component == "module":
        if string_of_modules.number_of_modules == nbr_basic_comp:
            return 1
        elif nbr_basic_comp == 1:
            return string_of_modules.number_of_modules
        else:
            raise Exception(err_msg)
    elif basic_component == "cell":
        total_syst_cell_nbr = string_of_modules.number_of_modules*\
        string_of_modules.module.number_of_cells
        if nbr_basic_comp == total_syst_cell_nbr:
            return 1
        elif nbr_basic_comp == string_of_modules.module.number_of_cells:
            return string_of_modules.number_of_modules
        elif nbr_basic_comp == 1:
            return total_syst_cell_nbr
        else:
            raise Exception(err_msg)

def comp_nbr_simulated_modules(climatic_data, module):
    """
    Compute the number of simulated modules, which may not be equal to the
    number of modules of the real system (string_of_modules.number_of_modules).
    """
    if climatic_data.basic_component == "module":
        return climatic_data.nbr_of_components
    elif climatic_data.basic_component == "cell":
        return climatic_data.nbr_of_components//module.number_of_cells
    else:
        raise NotImplementedError

def validate_inputs(climatic_data, string_of_modules):
    """
    Several supported cases:
    1) We have one set of climatic data ("for 1 module"), we make the
    thermo-electrical simulation for 1 module:
    a) We want output power for 1 module (climatic_data.nbr_of_components
    == string_of_modules.number_of_modules).
    b) We want output power for several modules (climatic_data.nbr_of_components
    > string_of_modules.number_of_modules).

    2) We have one set of climatic data for N modules, we make the
    thermo-electrical simulation for N modules. climatic_data.nbr_of_components
    == string_of_modules.number_of_modules

    3) We have one set of climatic data per cell, for N cells of 1 module. That
    is N = module.number_of_cells.
    We make the thermo-electrical simulation for N cells.
    a) We want output power for 1 module. (climatic_data.nbr_of_components == N,
    == string_of_modules.module.number_of_cells,
    string_of_modules.number_of_modules == 1)
    b) We want output power for M modules (climatic_data.nbr_of_components == N,
    == string_of_modules.module.number_of_cells,
    string_of_modules.number_of_modules == M).

    4) We have one set of climatic data per cell, for N*M cells, where
    is N = module.number_of_cells and M = string_of_modules.number_of_modules.
    Then climatic_data.nbr_of_components == N*M.
    """

    assert(climatic_data.basic_component == "module" or \
    climatic_data.basic_component == "cell")

    if climatic_data.basic_component == "module":
        assert(climatic_data.nbr_of_components == 1 or \
        climatic_data.nbr_of_components == string_of_modules.number_of_modules)
    else:
        assert(climatic_data.nbr_of_components == 1 or \
        climatic_data.nbr_of_components == string_of_modules.number_of_modules*\
        string_of_modules.module.number_of_cells or \
        climatic_data.nbr_of_components == \
        string_of_modules.module.number_of_cells)


def compute_energy_yield(climatic_data, string_of_modules, output_data_path=None,
sim_params=None, debug=False):

    validate_inputs(climatic_data, string_of_modules)

    if sim_params is None:
        sim_params = EnergyYieldSimParams()

    module = string_of_modules.module
    nbr_cells_in_basic_comp = comp_nbr_cells_in_basic_comp(
    climatic_data.basic_component, module)
    debug = sim_params.debug
    timestamp_index_min, timestamp_index_max = \
    time_interval_with_irr_gt_min(sim_params, climatic_data.global_irradiance)

    real_sim_nbr_bc_ratio = comp_real_sim_nbr_bc_ratio(string_of_modules,
    climatic_data)
    nbr_sim_modules = comp_nbr_simulated_modules(climatic_data, module)

    if debug.debug_mode:
        debug.write_debug_mat_file(["nbr_cells_in_basic_comp", "timestamp_index_min",
        "timestamp_index_max", "real_sim_nbr_bc_ratio", "nbr_sim_modules"],
        [nbr_cells_in_basic_comp, timestamp_index_min,
        timestamp_index_max, real_sim_nbr_bc_ratio, nbr_sim_modules])

    time_before_sim = time.time()
    result = dispatch_work_among_proc(climatic_data, sim_params,
    timestamp_index_min, timestamp_index_max, nbr_cells_in_basic_comp,
    string_of_modules, real_sim_nbr_bc_ratio, nbr_sim_modules, debug)
    simulation_time = time.time() - time_before_sim

    Pmpp_system, Vmpp_system, Pmpp_BC, temp_BC, wpeak_watt, net_heat, \
    heat_irr = result

    return EnergyYieldRes(Pmpp_system, Vmpp_system, Pmpp_BC, temp_BC,
    net_heat, heat_irr, simulation_time, wpeak_watt, climatic_data.time_vector)