%%javascript
require(['base/js/utils'],
function(utils) {
    utils.load_extensions('IPython-notebook-extensions-3.x/usability/comment-uncomment');
    utils.load_extensions('IPython-notebook-extensions-3.x/usability/dragdrop/main');
});

"""
Explanation: Note:
This is an 'advanced user' notebook. If you don't know what is going on here, please start with the following example notebooks:

Setting up a configuration file
Using pypmj - the mie2D-project


Description
This notebook deals with convergence testing, which is a common task and an important step in every numerical simulation technique. In the FEM-method, one typically needs to make sure that the mesh is accurate enough to represent to system for the target wavelength. Another significant parameter is the polynomial degree, or FEM-degree.
In JCMsuite, there are myriads of possibilities to controle the mesh. You can learn more about these techniques in the GeoTutorial of JCMsuite. For the FEM-degree, JCMsuite provides an adaptive algorithm to automatically chose different FEM-degrees on different elements of the mesh. You can specify minimum and maximum FEM-degrees for this procedure. The target precision, which controls the choice of the FEM-degrees, can be set using the PrecisionFieldEnergy parameter in the Refinement section of the project.
In pypmj, a class called ConvergenceTest is provided, which aims on automating the task of a convergence testing in the most general way. We will use the mie2D_extended project to showcase the approach. The meshing will be controlled using side length constraints of the involved domains.
Preparations
Notebook extensions
End of explanation
"""


%matplotlib inline
import pypmj as jpy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
try:
    import seaborn as sns
    sns.set_context('notebook')
except ImportError:
    print 'You do not seem to have `seaborn` yet. You may check it out!'

"""
Explanation: Imports
End of explanation
"""


cmap = jpy._config.get('Preferences', 'colormap')

"""
Explanation: For plotting, we load the prefered colormap from the configuration.
End of explanation
"""


project = jpy.JCMProject('scattering/mie/mie2D_extended')

"""
Explanation: Simulation
Preparing and configuring the convergence test
We start by creating a JCMProject-instance. Make sure you configured the project catalog well, which is configured in the section Data under key projects in the configuration file, or specify the absolute path to the mie2D_extend project here.
End of explanation
"""


keys_test = {'constants' :{}, 
             'parameters': {'fem_degree_max': np.arange(2,6),
                            'precision_field_energy':np.array([1.e-2,1.e-4])},
             'geometry': {'radius':0.3,
                          'slc_domain': np.linspace(0.1,0.4,10),
                          'slc_circle': np.linspace(0.05,0.3,10),
                          'refine_all_circle':4}}

keys_ref = {'constants' :{}, 
            'parameters': {'fem_degree_max': 6,
                           'precision_field_energy': 1.e-9},
            'geometry': {'radius': 0.3,
                         'slc_domain': 0.1,
                         'slc_circle': 0.05,
                         'refine_all_circle': 6}}

"""
Explanation: We now specify the keys for our convergence test. In contrast two a normal SimulationSet, we need two sets of keys, one for the so-called reference simulation and one for the test simulations. The reference simulation keys must induce a single simulation and its results will be used to calculate relative deviations, so we treat them as "analytical" results. Consequently, the accuracy of the reference simulation should be much higher than any of the test simulations. The ConvergenceTest-class initializes two simulation sets behind the scenes, one for the reference simulation and one for the test simulations. They will be placed in subfolders of the storage folder.
End of explanation
"""


ctest = jpy.ConvergenceTest(project, keys_test, keys_ref, duplicate_path_levels=0,
                            storage_folder='convergence_test',
                            storage_base=os.path.abspath('tmp_storage_folder'))

"""
Explanation: Now, the ConvergenceTest can be initialized. We will store the results in a subfolder of the current working directory for now, ignoring the storage base set in the configuration file.
End of explanation
"""


ctest.make_simulation_schedule()

"""
Explanation: The syntax for preparing and running the ConvergenceTest is prefereably identical to the one for a SimulationSet. We make a simulation schedule.
End of explanation
"""


def read_scs(pp):
    results = {} #must be a dict
    results['SCS'] = pp[0]['ElectromagneticFieldEnergyFlux'][0][0].real
    return results

"""
Explanation: Running the convergence test
To run the convergence test, we need to define the processing function, which is in this case identical to the one that was used in the mie2D example. It simply reads the scattering cross section from the FluxIntegration post process.
End of explanation
"""


ctest.use_only_resources('localhost')
jpy.resources['localhost'].set_m_n(4,1)
print jpy.resources['localhost']

"""
Explanation: Running a convergence test is nothing else than running two simulation sets, that's why you can simply pass the same arguments to the run-method. But there are two more arguments which can be specified. If save_run is True, it runs both simulation sets using the utility-function run_simusets_in_save_mode. More interesting is the run_ref_with_max_cores argument, which is set to 'AUTO' as default. As the reference simulation is a single simulation which typically needs massive computation power, it makes sense to use as many cores as possible for it. The default behavior is to detect the resource with the most available cores automatically and to use only this machine with a multiplicity of 1. For the number of threads, the product of the currently configured multiplicity and number of threads is used.
Let's make that clear. Let's configure the localhost as our only resource and set its multiplicity to 4 and the number of threads to 1:
End of explanation
"""


ctest.run(N=200, processing_func=read_scs, wdir_mode='zip')

"""
Explanation: There are now 4 cores available and there is no other machine available, as we restricted the resources to localhost. The run-method will now use this machine with a multiplicity of 1 and 4 threads for the reference simulation, which is the most efficient choice. For the test simulations, the normal configuration as specified above will be used.
We now run the convergence test, 200 simulations at a time, while zipping the working directories.
End of explanation
"""


analysis = ctest.analyze_convergence_results('SCS')
analysis.head()

"""
Explanation: Analyzing the results
Analyzing and interpreting a convergence test is a task which requires some experience and maybe sometimes also intuition. So the name of this section and of the analyze_convergence_results-method may be misleading in some way. The method just automates the most common step in this procedure: it calculates the relative deviation of the test data from the reference data for given result columns.
We only calculated the scattering cross section ('SCS'). In other cases, we could pass a list of result keys, i.e. column names, for the dev_columns-argument of the analyze_convergence_results-method. The sort_by-argument specifies the column of which the deviation data is used to sort the results in the end. If it is None, the first one is used.
End of explanation
"""


fig = plt.figure(figsize=(14,4))
ax=plt.gca()
analysis.plot.scatter(ax=ax, x='deviation_SCS', y='AccumulatedCPUTime', 
                      c='fem_degree_max', s=50, logx=True,
                      cmap=cmap)
plt.autoscale(tight=True)

# Save for later:
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()

plt.show()

"""
Explanation: The returned data frame (which is also stored in the class attribute analyzed_data) contains a column called 'deviation_SCS' and is sorted by this column in ascending order. So the first row represents the simulation with the smallest deviation.
Example of an extended analysis and interpretation of the data
In general, the task is to find the simulation parameters which satisfy our requirements on accuracy, while causing the smallest computational cost, e.g. CPU time. As the computational costs are automatically recorded by pypmj, we have full access to that data and can plot, for example, the distribution of CPU time as a function of the SCS deviation. We color the scatter dots by the used FEM degree.
End of explanation
"""


from mpl_toolkits.axes_grid1 import ImageGrid
sns.set(style="ticks")

analysis['slc_average'] = (analysis['slc_circle'] + analysis['slc_domain'])/2.
max_fems = np.sort(np.unique(analysis['fem_degree_max']))
precisions = np.sort(np.unique(analysis['precision_field_energy']))[::-1]

nrows = len(precisions)
ncols = len(max_fems)

fig = plt.figure(figsize=(15,8))
grid = ImageGrid(fig, 111, # similar to subplot(111)
                 nrows_ncols = (nrows, ncols), # NxM grid of axes
                 axes_pad = (0.25,0.5), # pad between axes in inch.
                 share_all=True,
                 label_mode = 'L',
                 cbar_location = 'right',
                 cbar_mode = 'single',
                 cbar_pad=0.2,
                 cbar_size = '5%',
                 aspect=None)

vmin, vmax = (np.min(analysis['slc_average']), np.max(analysis['slc_average']))

for index, ax in enumerate(grid):
    row, col = divmod(index,ncols)
    prec = precisions[row]
    fem = max_fems[col]
    
    data_ = analysis[analysis['fem_degree_max'] == fem]
    data_ = data_[data_['precision_field_energy'] == prec]
    scat = ax.scatter(x=data_['deviation_SCS'], y=data_['AccumulatedCPUTime'], 
                      c=data_['slc_average'], s=50, cmap=cmap, vmin=vmin, 
                      vmax=vmax)
    ax.semilogx(True)
    ax.set_xlim((xmin,xmax))
    ax.set_ylim((ymin,ymax))
    ax.set_ylabel('AccumulatedCPUTime')
    ax.set_xlabel('deviation_SCS')
    ax.set_title('max FEM={}, precision={}'.format(fem, prec), y=1.03)

cax = grid.cbar_axes[0]
cbar = plt.colorbar(scat, cax=cax, label='slc_average')
plt.show()

"""
Explanation: We observe that a higher maximum FEM-degree leads to a better accuracy. However, the relative deviation seems to saturate already for a fem_degree_max of 4 if other parameters are chosen well. This could be due to a smaller side length constraint (SLC). We could check this by calculating an average SLC and use this for the coloring. But we still have far to many parameters to make sense of the data by just looking at a single plot.
In a more advanced step, we can try to divide our data into subsets which have a constant maximum FEM degree and a constant precision field energy. Each of these subsets will be plotted seperately using the plotting technique above, this time coloring by the average SLC.
End of explanation
"""


candidates = analysis[(analysis['fem_degree_max']==4) & 
                      (analysis['precision_field_energy']==0.01) &
                      (analysis['deviation_SCS']<=1.e-4)]

"""
Explanation: Sure, this afforded some serious plotting and it may already be close to the limit of how many parameters we can study by this technique. However, we now see different tendencies as expected:

a larger maximum FEM degree causes smaller deviations, but it saturates if the target precision is not high enough
smaller average SLCs cause a lot of CPU time

We also spotted candidates which give a very high accuracy but small computation times, in the third column of the first row:
End of explanation
"""


relevant_cols = ['AccumulatedCPUTime', 'TotalMemory_GB', 'Unknowns',
                 'fem_degree_max', 'precision_field_energy', 'slc_circle',
                 'slc_domain', 'deviation_SCS']
optimum_1 = candidates.ix[candidates['AccumulatedCPUTime'].argmin()]
optimum_1[relevant_cols]

"""
Explanation: The candidate with the smallest accumulated CPU time has the following parameters:
End of explanation
"""


subset = analysis[analysis['deviation_SCS'] <= 0.1].copy(deep=True)
for slc in ['slc_circle', 'slc_domain']:
    subset.ix[:,slc] = np.round(subset[slc], 2)

"""
Explanation: Using seaborn to make sense of the convergence data
The seaborn package provides a lot of functionality to analyze data sets with lots of attributes, i.e. categories. This can be an enormous advantage and save much time. As you have seen above, plotting can be tedious in such cases and we even helped ourselves with the average SLC, which might not work out well in all cases.
Let's see what seaborn can do for us. We restrict our data to those simulations which have a deviation smaller than 0.1. (we only round the SLC values here for better readability in the plot below)
End of explanation
"""


sns.set(style='darkgrid', font_scale=1.2)
grid = sns.factorplot(x='fem_degree_max', y='deviation_SCS', hue='slc_domain', 
                      col='precision_field_energy', row = 'slc_circle',
                      data=subset, kind='swarm',
                      palette=cmap, margin_titles=True, size=2, aspect=2)
plt.show()

"""
Explanation: The factorplot-method can produce grids of plots by only specifying the keys that we want to use for generating the rows and the columns. Each plot will be a swarmplot, showing us the deviations grouped by the FEM degree and coloured by the SLC of the domain. We will have a column for each value of precision_field_energy and a row for each slc_circle.
End of explanation
"""


subset['performance'] = 1./(subset['deviation_SCS']*subset['AccumulatedCPUTime'])
grid = sns.factorplot(x='fem_degree_max', y='performance', hue='slc_domain', 
                      col='precision_field_energy', row = 'slc_circle',
                      data=subset, kind='swarm',
                      palette=cmap, margin_titles=True, size=2, aspect=2)
plt.show()

"""
Explanation: We could also use the inverse of the product of deviation_SCS and AccumulatedCPUTime as a performance measure, as this is the quantity we want to minimize.
End of explanation
"""


print 'The new optimum:'
optimum_2 = analysis.loc[subset['performance'].argmax()]
optimum_2[relevant_cols]

"""
Explanation: If we look for the performance maximum in this subset of our data, we find the following result:
End of explanation
"""


print 'The old optimum:'
optimum_1[relevant_cols]

"""
Explanation: If we compare to the previous result, we find that we now found a candidate with a smaller deviation at smaller CPU time. This is achieved by a much smaller value for the SLC of the circle, which seems to have the dominating effect. The domain does not need an as fine grid, which is expected because of the smaller permittivity. When looking at the second figure that was created using seaborn, this behavior can be seen very clearly.
End of explanation
"""


