""" Runs  Sensitivity Analysis on a finished MCE project.

"""

# Copyright (c) 2014 MTA Centre for Ecological Research
# Distributed under the GNU GPL v3. For full terms see the file LICENSE.

import numpy
import datetime
import ReportGenerator
import FileManager
import bbgc_db
import os
import wsconf

class BbgcSA:

    def __init__(self, input_arr, output_arr):
        self.in_arr = numpy.array(input_arr).astype('float')
        self.out_arr = numpy.array(output_arr)
        # Diff from mean
        in_var_means = numpy.mean(self.in_arr, axis=0)
        vecsub = numpy.vectorize(lambda a, mn: a - mn)
        self.diff_arr = vecsub(self.in_arr[:,[0]], in_var_means[0])
        for i in range(1,len(in_var_means)):
            self.diff_arr = numpy.concatenate((self.diff_arr, 
                        vecsub(self.in_arr[:,[i]], in_var_means[i])), axis=1)
        #variance of diff of input variables
        self.diff_variance_arr = numpy.var(self.diff_arr, axis=0)
        
    def get_sens_coeff_arr(self, col_index):
        # Lin. regression
        #print('Diff arr. dimensions: ' + str(self.diff_arr.shape))
        #print(self.diff_arr)
        #print('Out dim: ' + str(self.out_arr.shape))
        #print(self.out_arr)
        x = numpy.linalg.lstsq(self.diff_arr, self.out_arr[:,[col_index]].T[0])[0]
        #print('X')
        #print(x)
        #Overall variance
        vec_variance_mult = numpy.vectorize(lambda a, b: numpy.power(a, 2) * b)
        sigma_all = numpy.sum(vec_variance_mult(x, self.diff_variance_arr))
        # sensitivity coefficients for variables
        vec_sens_coeff = numpy.vectorize(lambda a, b, c: (numpy.power(a, 2) *b) / c * 100)
        sens_coeff_arr = vec_sens_coeff(x, self.diff_variance_arr, sigma_all)
        #print(sens_coeff_arr)
        return sens_coeff_arr

#-----------------------------------------------------------------------
def create_sa_report(var_names, var_values, output_var_ids, 
                     p_id, sa_pid, mce_data, sa_data, bdb):
    """ Create a PDF report from Sensitivity Analysis results.

    """

    p_id = str(p_id)
    sa_pid = str(sa_pid)
    report_fn = 'SA_report_' + p_id +'_'+ sa_pid + '.pdf'
    title = ['Biome-BGC SA', 'Sensitivity Analysis Report', 'BioVeL 2011 - 2014']
    rg = ReportGenerator.ReportGenerator(title,
             report_fn)
    #rg.add_text(title)
    rg.add_spacer()
    rg.add_section_title('Project summary')
    rg.add_text('Project title: ' + sa_data.project_data['title'])
    pr_url =  wsconf.drupal_url + 'node/' + sa_pid
    rg.add_url_line(pr_url, pr_url, 'SA project ID: ' + sa_pid + ' [', ']')
    rg.add_spacer()
    rg.add_text(sa_data.description)
    rg.add_spacer()

    # Get label - in vars, out vars
    #-------------------------------------
    in_var_labels = []
    for name in var_names:
        tmp_name = name if name.startswith('field_ds') else 'field_ds_mce_'+name+'1'
        in_var_labels.append(bdb.get_fields_label(tmp_name))
    #print(in_var_labels)
    out_var_labels = []
    for oid in output_var_ids:
        #print('Output id: ' + str(oid) + '\n')
        tmp_dict = bdb.get_parameter_data(oid)
        try:
            out_var_labels.append(tmp_dict['field_ov_name'] + ', ' + tmp_dict['title'])
        except KeyError:
            out_var_labels.append('(' + str(oid) + ')')
            
    # Selected Output variables summary
    #--------------------------------------
    rg.add_text('Selected annual outputs:')
    for outvar_label in out_var_labels:
        rg.add_text(' - ' + outvar_label)

    # Other summary
    #--------------------------------------
    rg.add_spacer()
    start_dt = sa_data.project_data['field_pr_start_time']
    now_dt = datetime.datetime.now()
    rg.add_text('Generated: ' + now_dt.strftime('%Y-%m-%d %H:%M'))
    delta_t = now_dt - start_dt
    hours, remainder = divmod(delta_t.seconds, 3600)
    minutes, seconds = divmod(remainder, 60)
    rg.add_text('Duration: %s:%s:%s' % (hours, minutes, seconds))
    rg.add_spacer()
    rg.add_text('This analysis is based on the outputs of the following Biome-BGC Monte Carlo Experiment (MCE):')
    rg.add_spacer()
    pr_url =  wsconf.drupal_url + 'node/' + p_id
    rg.add_url_line(pr_url, pr_url, 'MCE project ID: ' + p_id + ' [', ']')

    rg.add_text(mce_data.project_data['title'])
    rg.add_text('Model type: ' + str(mce_data.model_data['title']))
    rg.add_text('MCE ini ID: ' + str(mce_data.project_data['field_pr_dsin_mceini']))
    rg.add_text('Repetition: ' + str(mce_data.repeat_count))
    rg.add_text('Randomized input variables:')
    for rand_in_name in mce_data.mcefield_list.keys():
        rg.add_text(' - ' + rand_in_name + ' [' + mce_data.mcefield_list[rand_in_name]['minimum'] + ' - '
                    + mce_data.mcefield_list[rand_in_name]['maximum'] + '] ' 
                    + mce_data.mcefield_list[rand_in_name]['dimension'])
    rg.new_page()

    # Chart per output variable
    #-------------------------------------
    for var_n in range(len(output_var_ids)):
        rg.add_section_title(out_var_labels[var_n])
        table = []
        for i in range(0, len(in_var_labels)):
            table.append([str(i+1),in_var_labels[i],
                          var_values[var_n][i]])
        rg.add_table(table)
        rg.new_page()
        #nums = [str(x+1) for x in range(len(var_values[var_n]))]
        ebar = ReportGenerator.BarChart('Sensitivity coefficients of input variables on\n Output: ' 
                                        + out_var_labels[var_n], 'Input variables','%',
                                        [0] + var_values[var_n])
        image_name = 'sens_coeffs' + str(var_n) + '.png'
        ebar.create_image(image_name)
        rg.add_figure(image_name,
            'Sensitivity coefficients of randomised input variables (numbers are from the table above)')
        rg.new_page()
    rg.save_report()
    FileManager.copy_file2proj_dir(report_fn, p_id)
    bdb.insert_result_file(report_fn, p_id, 'report', sa_pid)
    os.system(wsconf.drupal_cache_clear)
    

#-----------------------------------------------------------------------
def run_sensitivity_analysis(input_arr, output_arr, randin_names, mce_pid, sa_pid):
    """ Run Sensitivity Analysis. """
    
    # az output arr egy 3D tomb (futas * evek * outputok), amibol
    # az evek szerinti atlagolassal 2D-s tombot kell kesziteni.
    bdb = bbgc_db.DB()
    
    out_mean_arr = []
    #print(numpy.array(output_arr).shape)
    for run_matrix in output_arr:
        out_mean_arr.append(numpy.mean(numpy.array(run_matrix),axis=0))
    out_mean_arr = numpy.array(out_mean_arr)
    in_arr = numpy.array(input_arr)
    #print(in_arr.shape)
    
    # Output Variable Names
    #-----------------------
    mce_data = bbgc_db.MceIniData(mce_pid, bdb)
    annout_dict = mce_data.normal_data['field_out_ann_vars']
    mce_annout_list = [annout_dict[k] for k in sorted(annout_dict)]
    sa_data = bbgc_db.ProjectSAData(sa_pid, bdb)
    selected_sa_outputs = sa_data.project_data['field_out_ann_vars'].values()
    csv_fn = 'sensitivity.csv'
    f = open(csv_fn,'w')
    #print(randin_names)
    #print(wsconf.csv_separator.join(randin_names))
    ret_str = 'output_id' + wsconf.csv_separator + \
              wsconf.csv_separator.join(randin_names) + '\n' 
    # for n in randin_names:
    #     ret_str += n + wsconf.csv_separator
    # ret_str += '\n'
    sa_obj = BbgcSA(in_arr, out_mean_arr)
    #coeff_matr = []
    var_values = []
    selected_outvar_names = []
    for col_index in range(0, len(out_mean_arr[0])):
        if mce_annout_list[col_index] in selected_sa_outputs:
            coeffs = sa_obj.get_sens_coeff_arr(col_index)
            ret_str += str(mce_annout_list[col_index]) + wsconf.csv_separator
            selected_outvar_names.append(mce_annout_list[col_index])
            tmp_values = []
            for c in coeffs:
                tmp_values.append("{0:0.2f}".format(c))
                ret_str += str(c) + wsconf.csv_separator
            var_values.append(tmp_values)
            ret_str += '\n'
    f.write(ret_str)
    f.close()
    FileManager.copy_file2proj_dir(csv_fn, str(mce_pid))
    bdb.insert_result_file(csv_fn, mce_pid, 'result', sa_pid)
    create_sa_report(randin_names, var_values, selected_outvar_names,
                     mce_pid, sa_pid, mce_data, sa_data, bdb)
    
    bdb.close()
    return(ret_str)

# a create_report meghivasa:
# create_report(wsconf.input_dir + pdir, pdir, stat_dict)
# A stat_dict a fentiben a project_dir-ben levo pickled fajl
# mcpr = ResultEvaluator.MCEProjectResults(pdir, pid, stat_dict)
# run_sensitivity_analysis(mcpr.randin_list, mcpr.annout_list, mcpr.randin_names)
