def test_kde_over_time():
    results = load_results(r'./../data/eng_trans_100.cPickle', zipped=False)
    
#    kde_over_time(results, log=False)
#    kde_over_time(results, log=True)
    kde_over_time(results, group_by='policy', grouping_specifiers=['no policy', 'adaptive policy'])
    plt.show()

 def test_get_rf_feature_scores(self):
     fn = r'../data/1000 flu cases no policy.tar.gz'
     results = load_results(fn)
     
     def classify(data):
         #get the output for deceased population
         result = data['deceased population region 1']
         
         #make an empty array of length equal to number of cases 
         classes =  np.zeros(result.shape[0])
         
         #if deceased population is higher then 1.000.000 people, classify as 1 
         classes[result[:, -1] > 1000000] = 1
         
         return classes
     
     scores, forest = fs.get_rf_feature_scores(results, classify, 
                                               random_state=10)
     
     self.assertEqual(len(scores), len(results[0].dtype.fields))
     self.assertTrue(isinstance(forest, RandomForestClassifier))
     
     ooi = 'nr deaths'
     results[1][ooi] = results[1]['deceased population region 1'][:,-1]
     scores, forest = fs.get_rf_feature_scores(results, ooi, 
                                               random_state=10)
     
     self.assertEqual(len(scores), len(results[0].dtype.fields))
     self.assertTrue(isinstance(forest, RandomForestRegressor))

def test_pairs_density():
    results = load_results(r'..\data\eng_trans_100.cPickle', zipped=False)
#    pairs_density(results)
#    pairs_density(results, colormap='binary')

    pairs_density(results, group_by='policy', grouping_specifiers=['no policy'])
    plt.show()

def test_envelopes3d_group_by():
    results = expWorkbench.load_results(r'1000 flu cases.cPickle')

    envelopes3d_group_by(results, 
                         outcome='infected fraction R1', 
                         groupBy="normal interregional contact rate",
                         logSpace=True)

    def test_get_univariate_feature_scores(self):
        fn = r'../data/1000 flu cases no policy.tar.gz'
        results = load_results(fn)
        
        def classify(data):
            #get the output for deceased population
            result = data['deceased population region 1']
            
            #make an empty array of length equal to number of cases 
            classes =  np.zeros(result.shape[0])
            
            #if deceased population is higher then 1.000.000 people, classify as 1 
            classes[result[:, -1] > 1000000] = 1
            
            return classes
        
        # f classify
        scores = fs.get_univariate_feature_scores(results, classify)
        self.assertEqual(len(scores), len(results[0].dtype.fields))

        # chi2
        scores = fs.get_univariate_feature_scores(results, classify, 
                                                  score_func='chi2')
        self.assertEqual(len(scores), len(results[0].dtype.fields))
        
        # f regression
        ooi = 'nr deaths'
        results[1][ooi] = results[1]['deceased population region 1'][:,-1]
        scores = fs.get_univariate_feature_scores(results, ooi)
        self.assertEqual(len(scores), len(results[0].dtype.fields))

 def test_prepare_outcomes(self):
     fn = r'../data/1000 flu cases no policy.tar.gz'
     results = load_results(fn)
     
     # string type correct
     ooi = 'nr deaths'
     results[1][ooi] = results[1]['deceased population region 1'][:,-1]
     y, categorical = fs._prepare_outcomes(results[1], ooi)
     
     self.assertFalse(categorical)
     self.assertTrue(len(y.shape)==1)
     
     # string type not correct --> KeyError
     with self.assertRaises(KeyError):
         fs._prepare_outcomes(results[1], "non existing key")
     
     # classify function correct
     def classify(data):
         result = data['deceased population region 1']
         classes =  np.zeros(result.shape[0])
         classes[result[:, -1] > 1000000] = 1
         return classes
     
     y, categorical = fs._prepare_outcomes(results[1], classify)
     
     self.assertTrue(categorical)
     self.assertTrue(len(y.shape)==1)
     
     # neither string nor classify function --> TypeError
     with self.assertRaises(TypeError):
         fs._prepare_outcomes(results[1], 1)

def test_pairs_lines():
    results = load_results(r'..\data\eng_trans_100.cPickle', zipped=False)    
    pairs_lines(results)
#    set_fig_to_bw(pairs_lines(results)[0])
    
    pairs_lines(results, group_by='policy')
#    set_fig_to_bw(pairs_lines(results, group_by='policy')[0])
    plt.show()

def test_envelopes3d():
    results = expWorkbench.load_results(r"1000 flu cases.cPickle")
    exp, res = results

    logical = exp["policy"] == "adaptive policy"
    new_exp = exp[logical][0:100]
    new_res = {}
    for key, value in res.items():
        new_res[key] = value[logical][0:100, :]

    envelopes3d((new_exp, new_res), "infected fraction R1", logSpace=True)

def test_pairs_scatter():
    results = load_results(r'..\data\eng_trans_100.cPickle', zipped=False)    
    
    pairs_scatter(results)
#    set_fig_to_bw(pairs_scatter(results)[0])
    
    pairs_scatter(results, group_by='policy',
                  grouping_specifiers='basic policy', legend=False)
#    set_fig_to_bw(pairs_scatter(results, group_by='policy')[0])
    
    pairs_scatter(results, group_by='policy', 
                  grouping_specifiers=['no policy', 'adaptive policy'])
#    set_fig_to_bw(pairs_scatter(results, group_by='policy', legend=False)[0])
    plt.show()

def test_multiple_densities():
    results = load_results(r'..\data\eng_trans_100.cPickle', zipped=False)
    
    
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  points_in_time = [2000])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  points_in_time = [2000, 2100])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  points_in_time = [2000, 2020, 2100])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  points_in_time = [2000, 2020, 2040, 2060])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  points_in_time = [2020, 2040, 2060, 2080, 2100])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  grouping_specifiers="no policy",
#                  points_in_time = [2000, 2020, 2040, 2060, 2080, 2100],
#                  plot_type=ENV_LIN,
#                  experiments_to_show=[1,2,10])
#    multiple_densities(results, 
#                  outcome_to_show="total fraction new technologies", 
#                  group_by="policy", 
#                  grouping_specifiers="no policy",
#                  points_in_time = [2000, 2020, 2040, 2060, 2080, 2100],
#                  plot_type=ENVELOPE,
#                  experiments_to_show=[1,2,10])
    multiple_densities(results, 
#                  group_by="policy", 
    #              grouping_specifiers="no policy",
                  points_in_time = [2040, 2045, 2050, 2060,2070,2080],
                  plot_type=ENVELOPE,
                  density=KDE,
                  log=False
    #              experiments_to_show=[np.arange(0, 100, 20)]
                  )
    
    plt.show()

def test_group_results():
    results = load_results(r'./../data/eng_trans_100.cPickle', zipped=False)
    experiments, outcomes = results
    
    # test indices
    grouping_specifiers = {'set1':np.arange(0,11),
                           'set2':np.arange(11,25),
                           'set3':np.arange(25,experiments.shape[0])}
    groups = group_results(experiments, outcomes, 
                           group_by='index', 
                           grouping_specifiers=grouping_specifiers)
    total_data = 0
    for value in groups.values():
        total_data += value[0].shape[0]
    print experiments.shape[0], total_data
    
    # test continuous parameter type
    array = experiments['average planning and construction period T1']
    grouping_specifiers = make_continuous_grouping_specifiers(array, nr_of_groups=5)
    groups = group_results(experiments, outcomes, 
                           group_by='average planning and construction period T1', 
                           grouping_specifiers=grouping_specifiers) 
    total_data = 0
    for value in groups.values():
        total_data += value[0].shape[0]
    print experiments.shape[0], total_data   
    
    # test integer type
    array = experiments['seed PR T1']
    grouping_specifiers = make_continuous_grouping_specifiers(array, nr_of_groups=10)
    groups = group_results(experiments, outcomes, 
                           group_by='seed PR T1', 
                           grouping_specifiers=grouping_specifiers) 
    total_data = 0
    for value in groups.values():
        total_data += value[0].shape[0]
    print experiments.shape[0], total_data   

    
    # test categorical type
    grouping_specifiers = set(experiments["policy"])
    groups = group_results(experiments, outcomes, 
                       group_by='policy', 
                       grouping_specifiers=grouping_specifiers)
    total_data = 0
    for value in groups.values():
        total_data += value[0].shape[0]
    print experiments.shape[0], total_data   

from expWorkbench import load_results
from analysis.prim import perform_prim, write_prim_to_stdout
from analysis.prim import show_boxes_individually


def classify(data):

    result = data["total fraction new technologies"]
    classes = np.zeros(result.shape[0])
    classes[result[:, -1] > 0.8] = 1

    return classes


if __name__ == "__main__":

    results = load_results(r"CESUN_optimized_1000_new.cPickle")
    experiments, results = results
    logicalIndex = experiments["policy"] == "Optimized Adaptive Policy"
    newExperiments = experiments[logicalIndex]
    newResults = {}
    for key, value in results.items():
        newResults[key] = value[logicalIndex]
    results = (newExperiments, newResults)

    boxes = perform_prim(results, "total fraction new technologies", threshold=0.6, threshold_type=-1)

    write_prim_to_stdout(boxes)
    show_boxes_individually(boxes, results)
    plt.show()

ema_logging.log_to_stderr(level=ema_logging.INFO)

def classify(data):
    #get the output for deceased population
    result = data['deceased population region 1']
    
    #make an empty array of length equal to number of cases 
    classes =  np.zeros(result.shape[0])
    
    #if deceased population is higher then 1.000.000 people, classify as 1 
    classes[result[:, -1] > 1000000] = 1
    
    return classes

#load data
results = load_results(r'./data/1000 flu cases.bz2')
experiments, results = results

#extract results for 1 policy
logical = experiments['policy'] == 'no policy'
new_experiments = experiments[ logical ]
new_results = {}
for key, value in results.items():
    new_results[key] = value[logical]

results = (new_experiments, new_results)

#perform prim on modified results tuple

prim = prim.Prim(results, classify, threshold=0.8, threshold_type=1)
box_1 = prim.find_box()

def perform_loop_knockout():    
    unique_edges = [['In Goods', 'lost'],
                    ['loss unprofitable extraction capacity', 'decommissioning extraction capacity'],
                    ['production', 'In Goods'],
                    ['production', 'lost'],
                    ['production', 'Supply'],
                    ['Real Annual Demand', 'substitution losses'],
                    ['Real Annual Demand', 'price elasticity of demand losses'],
                    ['Real Annual Demand', 'desired extraction capacity'],
                    ['Real Annual Demand', 'economic demand growth'],
                    ['average recycling cost', 'relative market price'],
                    ['recycling fraction', 'lost'],
                    ['commissioning recycling capacity', 'Recycling Capacity Under Construction'],
                    ['maximum amount recyclable', 'recycling fraction'],
                    ['profitability recycling', 'planned recycling capacity'],
                    ['relative market price', 'price elasticity of demand losses'],
                    ['constrained desired recycling capacity', 'gap between desired and constrained recycling capacity'],
                    ['profitability extraction', 'planned extraction capacity'],
                    ['commissioning extraction capacity', 'Extraction Capacity Under Construction'],
                    ['desired recycling', 'gap between desired and constrained recycling capacity'],
                    ['Installed Recycling Capacity', 'decommissioning recycling capacity'],
                    ['Installed Recycling Capacity', 'loss unprofitable recycling capacity'],
                    ['average extraction costs', 'profitability extraction'],
                    ['average extraction costs', 'relative attractiveness recycling']]
    
    unique_cons_edges = [['recycling', 'recycling'],
                           ['recycling', 'supply demand ratio'],
                           ['decommissioning recycling capacity', 'recycling fraction'],
                           ['returns to scale', 'relative attractiveness recycling'],
                           ['shortage price effect', 'relative price last year'],
                           ['shortage price effect', 'profitability extraction'],
                           ['loss unprofitable extraction capacity', 'loss unprofitable extraction capacity'],
                           ['production', 'recycling fraction'],
                           ['production', 'constrained desired recycling capacity'],
                           ['production', 'new cumulatively recycled'],
                           ['production', 'effective fraction recycled of supplied'],
                           ['loss unprofitable recycling capacity', 'recycling fraction'],
                           ['average recycling cost', 'loss unprofitable recycling capacity'],
                           ['recycling fraction', 'new cumulatively recycled'],
                           ['substitution losses', 'supply demand ratio'],
                           ['Installed Extraction Capacity', 'Extraction Capacity Under Construction'],
                           ['Installed Extraction Capacity', 'commissioning extraction capacity'],
                           ['Installed Recycling Capacity', 'Recycling Capacity Under Construction'],
                           ['Installed Recycling Capacity', 'commissioning recycling capacity'],
                           ['average extraction costs', 'profitability extraction']]
    
#    CONSTRUCTING THE ENSEMBLE AND SAVING THE RESULTS
    ema_logging.log_to_stderr(ema_logging.INFO)
    results = load_results(r'base.cPickle')

#    GETTING OUT THOSE BEHAVIOURS AND EXPERIMENT SETTINGS
#    Indices of a number of examples, these will be looked at.
    runs = [526,781,911,988,10,780,740,943,573,991]
    VOI = 'relative market price'
    
    results_of_interest = experiment_settings(results,runs,VOI)
    cases_of_interest = experiments_to_cases(results_of_interest[0])
    behaviour_int = results_of_interest[1][VOI]
    
#    CONSTRUCTING INTERVALS OF ATOMIC BEHAVIOUR PATTERNS
    ints = intervals(behaviour_int,False)

#    GETTING OUT ONLY THOSE OF MAXIMUM LENGTH PER BEHAVIOUR
    max_intervals = intervals_interest(ints)
    
#    THIS HAS TO DO WITH THE MODEL FORMULATION OF THE SWITCHES/VALUES
    double_list = [6,9,11,17,19]
    
    indCons = len(unique_edges)
#    for elem in unique_cons_edges:
#        unique_edges.append(elem)
    
    current = os.getcwd()

    for beh_no in range(0,10):
#        beh_no = 0 # Varies between 0 and 9, index style.
        interval = max_intervals[beh_no]
    
        rmp = behaviour_int[beh_no]
    #    rmp = rmp[interval[0]:interval[1]]
        x = range(0,len(rmp))
        fig = plt.figure()
        ax = fig.add_subplot(111)
    
        vensim.be_quiet()
    #    for loop_index in range(7,8):
        for loop_index in range(1,len(unique_edges)+1):
    
            if loop_index-indCons > 0:
                model_location = current + r'\Models\Consecutive\Metals EMA.vpm'
            elif loop_index == 0:
                model_location = current + r'\Models\Base\Metals EMA.vpm'
            else:
                model_location = current + r'\Models\Switches\Metals EMA.vpm'
        
            serie = run_interval(model_location,loop_index,
                                  interval,'relative market price',
                                  unique_edges,indCons,double_list,
                                  cases_of_interest[beh_no])
            
#.........这里部分代码省略.........

def test_lines3d():
    results = expWorkbench.load_results(r'eng_trans_100.cPickle')
    lines3d(results, outcomes=['installed capacity T1',
                               'installed capacity T2'])

'''
Created on 26 sep. 2011

@author: jhkwakkel
'''
import matplotlib.pyplot as plt

from expWorkbench import load_results
from analysis.graphs import lines

data = load_results(r'../../../src/analysis/1000 flu cases.cPickle')
fig = lines(data, column='fatality ratio region 1', density=False)
plt.show()

'''
Created on Sep 8, 2011

@author: gonengyucel, jhkwakkel
'''
import matplotlib.pyplot as plt

from analysis.clusterer import cluster

from expWorkbench import load_results
from expWorkbench import ema_logging

ema_logging.log_to_stderr(ema_logging.INFO)

#load the data
data = load_results(r'..\gallery\data\100 flu cases no policy.cPickle')

# specify the number of desired clusters
# note: the meaning of cValue is tied to the value for cMethod
cValue = 5

#perform cluster analysis
dist, clusteraloc, runlog, z = cluster(data=data, 
                                    outcome='deceased population region 1', 
                                    distance='gonenc', 
                                    interClusterDistance='complete', 
                                    cMethod = 'maxclust',
                                    cValue = cValue,
                                    plotDendrogram=False, 
                                    plotClusters=False, 
                                    groupPlot=False,

Created on Sep 8, 2011

.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
                gonengyucel

'''
import matplotlib.pyplot as plt

from analysis.clusterer import cluster

from expWorkbench import load_results, ema_logging

ema_logging.log_to_stderr(ema_logging.INFO)

#load the data
data = load_results(r'..\examples\100 flu cases no policy.cPickle')

# specify the number of desired clusters
# note: the meaning of cValue is tied to the value for cMethod
cValue = 5

#perform cluster analysis
dist, clusteraloc, runlog, z = cluster(data=data, 
                                    outcome='deceased population region 1', 
                                    distance='gonenc', 
                                    interClusterDistance='complete', 
                                    cMethod = 'maxclust',
                                    cValue = cValue,
                                    plotDendrogram=False, 
                                    plotClusters=False, 
                                    groupPlot=False,

'''
Created on 16 okt. 2012

@author: localadmin
'''
from expWorkbench import load_results
import numpy as np
import matplotlib.pyplot as plt

from analysis.plotting import lines
from analysis import plotting

results = load_results(r'.\data\2000 flu cases no policy.bz2')

experiments, outcomes = results

# get indices with of worst runs in terms of deaths and max fraction of 
# population that is ill at any given point in time
deaths = outcomes['deceased population region 1'][:, -1]
peak = np.max(outcomes['infected fraction R1'], axis=1)

deaths = [(deaths[i], i) for i in range(deaths.shape[0])]
deaths = sorted(deaths, reverse=True)
death_indices = [death[1] for death in deaths]

peak = [(peak [i], i) for i in range(peak.shape[0])]
peak = sorted(peak, reverse=True)
peak_indices = [element[1] for element in peak]

# combine the top 20 of both
indices = death_indices[0:20]