FVis2.py

#
# This program contains classes for managing and visualising data 
# produced by fluid simulation programs.
#
# Version: 1.1.7
#
# Last modified 20.05.2017 by Lars Frogner
#
import numpy as np
import matplotlib as mpl
mpl.use('TkAgg')                # Ensure that the Tkinter backend is used for generating figures
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size
import matplotlib.animation
import os, sys, time, datetime, struct


' *** For Windows users: Specify path to ffmpeg (required for saving animations) *** '

#plt.rcParams['animation.ffmpeg_path'] = r'<My path to ffmpeg folder>\ffmpeg\bin\ffmpeg'


class FluidVisualiser:

    'Class containing methods for visualising data generated by a fluid simulation.'

    # Public methods

    def __init__(self, printInfo=True, fontsize=11, blit=True):

        '''
        Constructor.
        '''

        self.printInfo = printInfo
        self.fontsize = float(fontsize)
        self.use_blit = blit
        self.hasSaved = False

        mpl.rcParams.update({'font.size': self.fontsize})

    def save_data(self, sim_time, update_func, rho=None, u=None, w=None, e=None, P=None, T=None, sim_fps=1, useDblPrec=False, sim_params=None, appendMode=False, folder='auto'):

        '''
        Advances the simulation by a given amount of time 
        and saves the relevant data to binary files.
        '''

        # *** Make sure that all the input is valid ***

        arr_names = ['rho', 'u', 'w', 'e', 'P', 'T']
        has_arr = []
        input_arrs = {}

        for arr, name in zip([rho, u, w, e, P, T], arr_names):

            has_arr.append(not (arr is None))
            if not (arr is None): input_arrs[name] = arr

        if not True in has_arr:

            print 'FluidVisualiser: No arrays to write. Aborting ...'
            sys.exit()

        arr_dim = len(input_arrs.values()[0].shape) # Number of dimensions for simulation
        arr_shape = input_arrs.values()[0].shape    # Shape of the data blocks to write

        # Loop through array dict and make sure that the provided arrays have the right type and shape

        for name, arr in input_arrs.iteritems():

            if not isinstance(arr, np.ndarray): raise TypeError('Input array \"%s\" has type %s. Must be a numpy array.' % (name, type(arr)))

        for arr in input_arrs.values():

            if len(arr.shape) != arr_dim: raise ValueError('Dimensions of input arrays don\'t match.')

        for arr in input_arrs.values():

            if arr.shape != arr_shape: raise ValueError('Shapes of input arrays don\'t match.')

        # Check that the other arguments have the right type
        try:

            sim_time = float(sim_time)
            sim_fps = float(sim_fps)
            folder = str(folder)

        except ValueError:

            raise ValueError('Invalid type for simulation time, fps or folder name.')

        if not hasattr(update_func, '__call__'): raise TypeError('Input argument \"update_func\" must be a function or method.')

        if not (sim_params is None):

            if not isinstance(sim_params, dict):

                raise TypeError('Input argument \"sim_params\" must be a dictionary.')

            for name in sim_params:

                if not isinstance(name, str): raise TypeError('Keys in \"sim_params\" must be strings.')


        if self.printInfo: print '\nFluidVisualiser: Simulating %g time units (%d frames will be %s binary files) ...\n' % (sim_time, sim_fps*sim_time, 'added to existing' if appendMode else 'written to new')

        # Get current working directory
        cwd = os.getcwd()

        # Initialize dictionary for writer instances
        arr_writers = {}

        if appendMode:

            # Make sure that the specified folder is present

            if folder == 'auto':

                raise ValueError('FluidVisualiser: Cannot use automatic folder naming in append mode. Please specify a folder.')

            if not os.path.isdir(folder):

                raise ValueError('FluidVisualiser: Cannot find the specified folder to append data in.')

            # Check that all the necessary data is present and compatible with the inputted data

            existing_arrs = self.get_last_data(folder)[0]

            if existing_arrs.viewkeys() != input_arrs.viewkeys():

                raise ValueError('FluidVisualiser: The set of inputted arrays doesn\'t correspond to the set of arrays in the folder.')

            if existing_arrs.values()[0].shape != arr_shape:

                raise ValueError('FluidVisualiser: The arrays to append don\'t have the same shape as the existing arrays.')

            # Move to the folder to append data in
            newdirpath = os.path.join(cwd, self.folder)
            os.chdir(newdirpath)

            # Create instances of the _BinAppender class for writing the data

            for name in input_arrs:

                arr_writers[name] = _BinAppender(name + '.pas', arr_shape)

            time_writer = _BinAppender('time.pas', (2,))

            t_init = self.t_list[-1]

        else:

            # Create output folder

            # Create a string of the current date and time (or use a user specified name)
            now = datetime.datetime.now()
            self.folder = folder if folder != 'auto' else ('FVis_output_' + now.strftime('%Y-%m-%d_%H-%M'))

            # Use the string to create a unique folder and make it the current working directory

            self.newfolder_count = 0

            def create_folder():

                try:

                    os.mkdir(self.folder)

                except OSError as e:

                    if self.newfolder_count < 5:

                        self.folder += '_new'
                        self.newfolder_count += 1
                        create_folder()

                    else:

                        raise OSError(e)

            # Keep adding _new to the name until it is unique
            create_folder()

            newdirpath = os.path.join(cwd, self.folder)
            os.chdir(newdirpath)

            # Create instances of the _BinWriter class for writing the data

            # Data type to write (only relevant when not appending)
            dtype = 'f8' if useDblPrec else 'f4'

            for name in input_arrs:

                arr_writers[name] = _BinWriter(name + '.pas', arr_shape, dtype)

            time_writer = _BinWriter('time.pas', (2,), dtype)

            t_init = 0

            # Write the initial data to the files
            for name in input_arrs:

                arr_writers[name].write_block(input_arrs[name])

            time_writer.write_block([t_init, 0])


        # *** Write data ***

        # Approximate time between each frame to save
        t_skip = 1.0/sim_fps
        total_time = t_init
        elapsed_time = 0

        # Start measuring elapsed time
        start_time = time.time()

        wasInterrupted = False

        while elapsed_time < sim_time:

            try:

                # Run an appropriate number of simulation steps

                t = 0
                skips = 0

                while t < t_skip:

                    t += update_func()
                    skips += 1

                dt_avg = t/skips

                total_time += t
                elapsed_time = total_time - t_init

                # Write the current data to the files

                for name in input_arrs:

                    arr_writers[name].write_block(input_arrs[name])

                time_writer.write_block([total_time, dt_avg])

                if wasInterrupted: break

                # Print progress
                if self.printInfo: FluidVisualiser.__print_progress(elapsed_time, sim_time, start_time)

            # If aborted by user, exit the loop and finish up the writing
            except KeyboardInterrupt:

                if self.printInfo: print '\n\nFluidVisualiser: Writing interrupted by user. Aborting ...'

                wasInterrupted = True

            except TypeError:

                raise TypeError('Input argument \"update_func\" must return the time step length (and nothing else).')

        if self.printInfo: print

        if appendMode:

            # Close files
            for name in input_arrs:

                arr_writers[name].end_append(printReport=self.printInfo)

            time_writer.end_append(printReport=self.printInfo)

        else:

            # Close files
            for name in input_arrs:

                arr_writers[name].end_write(printReport=self.printInfo)

            time_writer.end_write(printReport=self.printInfo)

            # Write simulation parameters to a text file

            info = ''

            # Write which arrays were given

            for i in xrange(len(arr_names)):

                info += 'has_%s = %d\n' % (arr_names[i], has_arr[i])

            # Write parameters

            def convert_if_num(val):

                if isinstance(val, float) or isinstance(val, int):

                    newval = '%g' % val

                else:

                    newval = str(val)

                return newval

            if not (sim_params is None):

                for name in sorted(sim_params):

                    info += name + ' = ' + convert_if_num(sim_params[name]) + '\n'

            f = open('info.txt', 'w')
            f.write(info[:-1])
            f.close()

        # Return to original working directory
        os.chdir(cwd)

        if self.printInfo: print '\nFluidVisualiser: Data stored in folder \"%s\".' % self.folder

        # Indicate that data was saved in this instance
        self.hasSaved = True

    def animate_1D(self, quantity, folder='default', extent=[0, 1], anim_fps='auto', showDeviations=True, showParams=True, height=6, aspect=1.3, title='', save=False, anim_time='auto', video_fps=30, video_name='auto', snapshots=None, units={}):

        'Creates an animation of the time evolution of a 1D simulation.'

        # *** Make sure that the inputted values are valid ***

        quantity = str(quantity)
        folder = str(folder)
        if anim_fps != 'auto': anim_fps = float(anim_fps)
        height = float(height)
        if aspect != 'auto': aspect = float(aspect)
        title = str(title)
        if anim_time != 'auto': anim_time = float(anim_time)
        video_fps = float(video_fps)
        video_name = str(video_name)
        extent = [float(extent[0]), float(extent[1])]

        save_snaps = snapshots is not None
        if save_snaps:
            snapshots = sorted([float(t) for t in list(snapshots)])
        else:
            snapshots = []

        units = {k: str(units[k]) for k in units}
            
        # Set valid folder name
        self.__set_foldername(folder)


        # *** Prepare figure ***

        extent_unit = r'm' if not 'L' in units else units['L']
        time_unit = r's' if not 't' in units else units['t']

        if time_unit == 's':
            time_displayer = lambda t: FluidVisualiser.__s_to_hms(t)
        else:
            time_displayer = lambda t: '%.3g %s' % (t, time_unit)

        # Get information required for displaying the given quantity
        self.__get_init_data(extent=extent)
        q, name, default_unit = self.__get_quantity_info(quantity)
        unit = default_unit if not quantity in units else units[quantity]
        min_val, max_val = self.__get_optimal_scaling(q, includeAll=True)

        if self.printInfo: print 'FluidVisualiser: Preparing figure ...'

        # Define figure and axis
        fig = plt.figure(figsize=(height*aspect, height))
        ax = fig.add_subplot(111)

        # Create plot for fluid
        line, = ax.plot(self.l, q(), animated=self.use_blit)
        if self.use_blit and not (save or save_snaps): ax.cla()

        # Add figure info
        textbox1, textbox2, textbox3 = self.__prepare_text(ax)
        
        # Set axis limits and labels
        ax.set_xlim(extent[0], extent[1])

        ax.set_ylim(min_val, max_val)
        ax.set_xlabel(('Vertical' if self.direction == 'z' else 'Horizontal') + ' distance' 
                      + ('' if extent_unit == '' else ' [%s]' % extent_unit))
        ax.set_ylabel(name[0].upper() + name[1:] + ('' if unit == '' else ' [%s]' % unit))
        ax.set_title(title)


        # *** Define function for updating animation ***

        if anim_fps == 'auto': anim_fps = float(self.Nt)/(self.t_list[-1] - self.t_list[0])
        t_skip = 1.0/(anim_fps) # Time between each frame to render
        N_frames = int(np.floor((self.t_list[-1] if anim_time == 'auto' else anim_time)*anim_fps)) # Total number of frames

        # Initial values of total mass and energy density
        if self.has_arr['rho']: rho_tot0 = np.sum(self.arrs['rho'])
        if self.has_arr['e']: e_tot0 = np.sum(self.arrs['e'])

        def update(i):

            # Run an appropriate number of simulation steps
            dt_avg = self.__step_time(t_skip)

            # Update fluid image data
            line.set_ydata(q())

            # Update text

            textbox1.set_text('Time: %s\ndt = %.2g %s' % (time_displayer(self.t), dt_avg, time_unit))

            if showDeviations:

                devi_text = ''
                if self.has_arr['rho']: devi_text += 'Mass deviation: %.3f %%' % (100*(np.sum(self.arrs['rho'].astype('f8'))/rho_tot0 - 1.0))
                if self.has_arr['e']: devi_text += '\nEnergy deviation: %.3f %%' % (100*(np.sum(self.arrs['e'].astype('f8'))/e_tot0 - 1.0))

                textbox2.set_text(devi_text)

            if showParams: textbox3.set_text(self.param_text)

            if save_snaps and len(snapshots) > 0 and self.t >= snapshots[0]:

                snapname = '%s_snapshot_t%g.png' % (self.folder, self.t)
                if self.printInfo: print 'FluidVisualiser: Saving snapshot as \"%s\" ...' % snapname

                fig.savefig(snapname, bbox_inches='tight')
                snapshots.pop(0)

            # Print progress if generating a movie
            if save: FluidVisualiser.__print_progress(i, N_frames, self.t0)

            # Return figure content
            return line, textbox1, textbox2, textbox3

        self.__run_animation(fig, update, save, N_frames, video_fps, video_name, show=(not save_snaps))

    def animate_2D(self, quantity, matrixLike=True, folder='default', extent=[0, 1, 0, 1], anim_fps='auto', showDeviations=True, showParams=True, showQuiver=True, quiverscale=1, N_arrows=20, interpolation='none', cmap='jet', height=6, aspect='equal', cbar_aspect=0.02, title='', save=False, anim_time='auto', video_fps=30, video_name='auto', backgrounds=None, snapshots=None, units={}):

        'Creates an animation of the time evolution of a 2D simulation.'

        # *** Make sure that the inputted values are valid ***

        quantity = str(quantity)
        folder = str(folder)
        if anim_fps != 'auto': anim_fps = float(anim_fps)
        quiverscale = float(quiverscale)
        N_arrows = int(N_arrows)
        interpolation = str(interpolation)
        cmap = str(cmap)
        height = float(height)
        if aspect != 'equal': aspect = float(aspect)
        cbar_aspect = float(cbar_aspect)
        title = str(title)
        if anim_time != 'auto': anim_time = float(anim_time)
        video_fps = float(video_fps)
        video_name = str(video_name)
        extent = [float(extent[0]), float(extent[1]), float(extent[2]), float(extent[3])]

        if not (backgrounds is None):
            if not isinstance(backgrounds, dict): raise TypeError
        else:
            backgrounds = False

        save_snaps = snapshots is not None
        if save_snaps:
            snapshots = sorted([float(t) for t in list(snapshots)])
        else:
            snapshots = []

        units = {k: str(units[k]) for k in units}

        # Set valid folder name
        self.__set_foldername(folder)

        # *** Prepare figure ***

        extent_unit_x = r'm' if not 'Lx' in units else units['Lx']
        extent_unit_z = r'm' if not 'Lz' in units else units['Lz']
        speed_unit_x = r'm$/$s' if not 'u' in units else units['u']
        speed_unit_z = r'm$/$s' if not 'w' in units else units['w']
        time_unit = r's' if not 't' in units else units['t']

        if time_unit == 's':
            time_displayer = lambda t: FluidVisualiser.__s_to_hms(t)
        else:
            time_displayer = lambda t: '%.3g %s' % (t, time_unit)

        # Get information required for displaying the given quantity
        self.__get_init_data(extent=extent, matrixLike=matrixLike)
        q, name, default_unit = self.__get_quantity_info(quantity, bg=backgrounds)
        unit = default_unit if not quantity in units else units[quantity]
        min_val, max_val = self.__get_optimal_scaling(q)

        if self.printInfo: print 'FluidVisualiser: Preparing figure ...'

        aspect = self.equal_aspect if aspect == 'equal' else float(aspect)

        # Define figure and axis
        fig = plt.figure(figsize=(height*aspect, height))
        ax = fig.add_subplot(111)

        # Create image for fluid
        img = ax.imshow(q() if matrixLike else q().T, origin='lower', norm=plt.Normalize(min_val, max_val),
                             cmap=plt.get_cmap(cmap), interpolation=interpolation, 
                             extent=extent, aspect='auto', animated=True)

        # Create quiver plot
        if extent_unit_x != extent_unit_z or speed_unit_x != speed_unit_z:
            arrowscale_unit = ''
        else:
            arrowscale_unit = '(%s)/%s' % (speed_unit_x, extent_unit_x)
        if not (self.has_arr['u'] and self.has_arr['w']): showQuiver = False
        step_q, quiver = self.__get_quiver(ax, N_arrows, quiverscale, showQuiver, arrowscale_unit)

        # Add figure info
        textbox1, textbox2, textbox3 = self.__prepare_text(ax)
        
        # Set axis limits and labels
        ax.set_xlim(extent[0], extent[1])
        ax.set_ylim(extent[2], extent[3])
        ax.set_xlabel('Horizontal distance' + ('' if extent_unit_x == '' else ' [%s]' % extent_unit_x))
        ax.set_ylabel('Vertical distance' + ('' if extent_unit_z == '' else ' [%s]' % extent_unit_z))
        ax.set_title(title, y=1.05)

        # Display colorbar
        width = axes_size.AxesY(ax, aspect=cbar_aspect)
        pad = axes_size.Fraction(0.4, width)
        divider = make_axes_locatable(ax)
        cax = divider.append_axes('right', size=width, pad=pad)
        fig.colorbar(img, cax=cax, label=(name[0].upper() + name[1:] + ('' if unit == '' else ' [%s]' % unit)))

        # Remove unnecessary figure space
        if self.Nx >= self.Nz: fig.tight_layout(pad=(1 if save else 5))


        # *** Define function for updating animation ***
        
        if anim_fps == 'auto': anim_fps = float(self.Nt)/self.t_list[-1]
        t_skip = 1.0/anim_fps    # Time between each frame to render
        N_frames = int(np.floor((self.t_list[-1] if anim_time == 'auto' else anim_time)*anim_fps)) # Total number of frames

        # Initial values of total mass and energy density
        if self.has_arr['rho']: rho_tot0 = np.sum(self.arrs['rho'])
        if self.has_arr['e']: e_tot0 = np.sum(self.arrs['e'])

        def update(i):

            # Run an appropriate number of simulation steps
            dt_avg = self.__step_time(t_skip)

            # Update fluid image data
            img.set_array(q() if matrixLike else q().T)

            # Update quiver plot
            if showQuiver:

                if matrixLike:

                    quiver.set_UVC(self.arrs['u'][::step_q, ::step_q], self.arrs['w'][::step_q, ::step_q])

                else:

                    quiver.set_UVC(self.arrs['u'][::step_q, ::step_q].T, self.arrs['w'][::step_q, ::step_q].T)

            # Update text

            textbox1.set_text('Time: %s\ndt = %.2g %s' % (time_displayer(self.t), dt_avg, time_unit))

            if showDeviations:

                devi_text = ''
                if self.has_arr['rho']: devi_text += 'Mass deviation: %.3f %%' % (100*(np.sum(self.arrs['rho'].astype('f8'))/rho_tot0 - 1.0))
                if self.has_arr['e']: devi_text += '\nEnergy deviation: %.3f %%' % (100*(np.sum(self.arrs['e'].astype('f8'))/e_tot0 - 1.0))

                textbox2.set_text(devi_text)

            if showParams: textbox3.set_text(self.param_text)

            if save_snaps and len(snapshots) > 0 and self.t >= snapshots[0]:

                snapname = '%s_snapshot_t%g.png' % (self.folder, self.t)
                if self.printInfo: print 'FluidVisualiser: Saving snapshot as \"%s\" ...' % snapname

                fig.savefig(snapname, bbox_inches='tight')
                snapshots.pop(0)

            # Print progress if generating a movie
            if save: FluidVisualiser.__print_progress(i, N_frames, self.t0)

            # Return figure content
            return img, quiver, textbox1, textbox2, textbox3

        self.__run_animation(fig, update, save, N_frames, video_fps, video_name, show=(not save_snaps))

    def animate_energyflux(self, folder='default', extent=[0, 1, 0, 1], anim_fps='auto', showParams=True, height=7, aspect=1.1, title='', save=False, anim_time='auto', video_fps=30, video_name='auto', snapshots=None, units={}):

        'Creates an animation of the horizontally averaged vertical energy flux.'

        # *** Make sure that the inputted values are valid ***

        folder = str(folder)
        if anim_fps != 'auto': anim_fps = float(anim_fps)
        height = float(height)
        if aspect != 'auto': aspect = float(aspect)
        title = str(title)
        if anim_time != 'auto': anim_time = float(anim_time)
        video_fps = float(video_fps)
        video_name = str(video_name)
        extent = [float(extent[0]), float(extent[1]), float(extent[2]), float(extent[3])]

        save_snaps = snapshots is not None
        if save_snaps:
            snapshots = sorted([float(t) for t in list(snapshots)])
        else:
            snapshots = []

        units = {k: str(units[k]) for k in units}

        # Set valid folder name
        self.__set_foldername(folder)

        # *** Prepare figure ***

        extent_unit_z = r'm' if not 'Lz' in units else units['Lz']
        time_unit = r's' if not 't' in units else units['t']

        if time_unit == 's':
            time_displayer = lambda t: FluidVisualiser.__s_to_hms(t)
        else:
            time_displayer = lambda t: '%.3g %s' % (t, time_unit)


        # *** Prepare figure ***

        # Get information required for displaying the flux
        self.__get_init_data(extent=extent)

        quantity = 'ew'
        q, name, default_unit = self.__get_quantity_info(quantity) # Get data for vertical energy flux
        q2 = lambda: np.mean(q(), axis=1)                          # Define new function returning horizontally averaged flux

        unit = default_unit if not quantity in units else units[quantity]

        min_val, max_val = self.__get_optimal_scaling(q2, includeAll=True)

        if self.printInfo: print 'FluidVisualiser: Preparing figure ...'

        # Define figure and axis
        fig = plt.figure(figsize=(height*aspect, height))
        ax = fig.add_subplot(111)

        # Create plot for fluid
        line, = ax.plot(q2(), np.linspace(extent[2], extent[3], self.Nz))

        # Add figure info
        textbox1, textbox2, textbox3 = self.__prepare_text(ax, showExtent=False)
        
        # Set axis limits and labels
        ax.set_xlim(min_val, max_val)
        ax.set_ylim(extent[2], extent[3])
        ax.set_xlabel(name[0].upper() + name[1:] + ('' if unit == '' else ' [%s]' % unit))
        ax.set_ylabel('Vertical distance' + ('' if extent_unit_z == '' else ' [%s]' % extent_unit_z))
        ax.set_title(title)


        # *** Define function for updating animation ***

        if anim_fps == 'auto': anim_fps = float(self.Nt)/(self.t_list[-1] - self.t_list[0])
        t_skip = 1.0/(anim_fps)      # Time between each frame to render
        N_frames = int(np.floor((self.t_list[-1] if anim_time == 'auto' else anim_time)*anim_fps)) # Total number of frames

        def update(i):

            # Run an appropriate number of simulation steps
            dt_avg = self.__step_time(t_skip)

            # Update fluid image data
            line.set_xdata(q2())

            # Update text

            textbox1.set_text('Time: %s\ndt = %.2g %s' % (time_displayer(self.t), dt_avg, time_unit))
            textbox2.set_text('Average: %.2g %s' % (np.mean(q2()), unit))

            if showParams: textbox3.set_text(self.param_text)

            if save_snaps and len(snapshots) > 0 and self.t >= snapshots[0]:

                snapname = '%s_snapshot_t%g.png' % (self.folder, self.t)
                if self.printInfo: print 'FluidVisualiser: Saving snapshot as \"%s\" ...' % snapname

                fig.savefig(snapname, bbox_inches='tight')
                snapshots.pop(0)

            # Print progress if generating a movie
            if save: FluidVisualiser.__print_progress(i, N_frames, self.t0)

            # Return figure content
            return line, textbox1, textbox2, textbox3

        self.__run_animation(fig, update, save, N_frames, video_fps, video_name, show=(not save_snaps))

    def plot_avg(self, quantity, folder='default', measure_time='auto', title='', showTrendline=False, relative=False, units={}):

        'Plots the time evolution of the average of a given quantity.'

        # *** Make sure that the inputted values are valid ***

        quantity = str(quantity)
        folder = str(folder)
        if measure_time != 'auto': measure_time = float(measure_time)
        units = {k: str(units[k]) for k in units}

        # Set valid folder name
        self.__set_foldername(folder)


        # *** Simulate and measure values ***

        time_unit = 's' if not 't' in units else units['t']

        self.__get_init_data()

        # Get function returning the quantity
        q, name, default_unit = self.__get_quantity_info(quantity)
        unit = default_unit if not quantity in units else units[quantity]

        # Lists for storing data points
        t_list = [self.t]
        q_list = [np.sum(q())]

        measure_time = self.t_list[-1] if measure_time == 'auto' else measure_time

        # Approximate amount of time between each data point
        t_step = measure_time/500.0

        if self.printInfo: print 'FluidVisualiser: Running simulation for %g s and measuring average %s ...' % (measure_time, name)

        t0 = time.time()

        # Loop until the given time has elapsed in the simulation
        while self.t - t_list[0] < measure_time:

            # Print progress
            FluidVisualiser.__print_progress(self.t - t_list[0], measure_time, t0)

            # Advance the simulation by the time given by t_step
            self.__step_time(t_step)

            # Stop if end of files has been reached
            if self.t < t_list[-1]:

                if (measure_time - t_list[-1])/measure_time > 0.05: print '\n\nFluidVisualiser: End of files reached. Stopping measurements.'
                break

            # Store time and quantity data
            t_list.append(self.t)
            q_list.append(np.sum(q()))

        # Close files
        self.__close_read_files()

        if self.printInfo: print '\n\nFluidVisualiser: Plotting time evolution of average %s ...' % name


        # *** Create and display plot ***

        if relative:
            data = (np.array(q_list) - np.mean(q_list))/np.mean(q_list) if q_list[0] == 0 else (np.array(q_list) - q_list[0])/q_list[0]
        else:
            data = np.array(q_list)

        x_c = np.sum(t_list)/len(t_list)
        y_c = np.sum(data)/len(data)
        a = (y_c - data[0])/x_c

        trendline = a*np.array(t_list) + data[0]

        plt.plot(t_list, data, '-')

        if showTrendline:
            plt.plot(t_list, trendline, '-', label='Trendline')
            plt.legend(loc='best')

        plt.title(title)
        plt.xlabel('Time' + ('' if time_unit == '' else ' [%s]' % time_unit))
        plt.ylabel('Average %s%s%s' % (name, '' if not relative else (' (relative to %s value)' % ('mean' if q_list[0] == 0 else 'initial')), '' if unit == '' else ' [%s]' % unit))
        plt.show()

    def delete_current_data(self):

        'Deletes the data saved by save_data in this instance.'

        if not self.hasSaved:

            print 'FluidVisualiser: Nothing to delete: no data was saved by this instance.'

        else:

            yn = raw_input('\nFluidVisualiser: Are you sure you want to delete the folder \"%s\" and all it\'s content? [Y/n] ' % (self.folder))

            while not yn.lower() in ['y', 'n']:

                print 'FluidVisualiser: Answer must be \"y\" (yes) or \"n\" (no).'
                yn = raw_input('\nFluidVisualiser: Are you sure you want to delete the folder \"%s\" and all it\'s content? [Y/n] ' % (self.folder))

            if yn.lower() == 'y':

                # Move into folder
                cwd = os.getcwd()
                newdirpath = os.path.join(cwd, self.folder)
                os.chdir(newdirpath)

                def try_remove(name):

                    try:
                        os.remove(name)

                    except OSError:

                        pass

                # Delete content
                try_remove('rho.pas')
                try_remove('u.pas')
                try_remove('w.pas')
                try_remove('e.pas')
                try_remove('P.pas')
                try_remove('T.pas')
                os.remove('time.pas')
                os.remove('info.txt')

                # Move back and delete folder
                os.chdir(cwd)
                os.rmdir(self.folder)

                print 'FluidVisualiser: All data in \"%s\" was deleted.' % (self.folder)

            else:

                print 'FluidVisualiser: No data deleted.'

    def get_last_data(self, folder):

        '''
        Reads the last arrays in the files of the given folder
        and returns them in a dictionary.
        '''

        # Set valid folder name
        self.__set_foldername(folder)

        # Set attributes required for reading the data
        self.__get_init_data()

        # Read the last arrays in the files
        self.__get_data_blocks(self.Nt - 1)

        # Close files
        self.__close_read_files()

        def adaptive_convert(val):

            if val == 'True':

                newval = True

            elif val == 'False':

                newval = False

            else:

                try:

                    newval = float(val)

                except ValueError:

                    newval = val

            return newval

        # Create a dictionary of simulation parameters
        out_params = {param_line.split(' = ')[0]:adaptive_convert(param_line.split(' = ')[1]) for param_line in self.param_text.split('\n')}

        # Return data and parameters
        return self.arrs, out_params

    # Private methods

    def __run_animation(self, fig, update, save, N_frames, video_fps, video_name, show=True):

        'Create an appropriate animation instance and run the animation procedure.'

        if save and show:

            if self.printInfo: print 'FluidVisualiser: Generating animation ...\n'

            if video_name == 'auto':

                video_name = self.folder

            while os.path.exists(video_name + '.mp4'):

                video_name += '_new'

            video_name += '.mp4'

            animation = matplotlib.animation.FuncAnimation(fig, update, blit=False, frames=N_frames)

            self.t0 = time.time()

            animation.save(video_name, writer=matplotlib.animation.FFMpegWriter(fps=video_fps, bitrate=3200, extra_args=['-vcodec', 'libx264']))

            if self.printInfo: print '\n\nFluidVisualiser: Animation saved as \"%s\".' % video_name

        elif show:

            if self.printInfo: print 'FluidVisualiser: Playing animation ...'

            animation = matplotlib.animation.FuncAnimation(fig, update, blit=self.use_blit)
            plt.show()

        else:

            for i in xrange(N_frames):
                update(i)

        # Close files
        self.__close_read_files()

    def __get_init_data(self, extent=False, matrixLike=True):

        '''
        Read the info file and headers in the saved folder and define 
        attributes required for reading the data.
        '''

        # *** Move to the folder containing the binary files to read ***

        cwd = os.getcwd()
        newdirpath = os.path.join(cwd, self.folder)
        os.chdir(newdirpath)


        # *** Get simulation parameters ***

        arr_names = ['rho', 'u', 'w', 'e', 'P', 'T']

        f = open('info.txt', 'r')
        self.lines = f.readlines()
        f.close()

        self.has_arr = {arr_names[i]:bool(int(self.lines[i].split(' = ')[1])) for i in xrange(6)}
        self.param_text = ''.join(self.lines[6:])


        # *** Create _BinReader instances for reading the data ***

        self.arr_readers = {}
        self.arrs = {}

        for name in arr_names:

            if self.has_arr[name]:

                self.arr_readers[name] = _BinReader(name + '.pas')
                self.arrs[name] = None


        # *** Read the time data ***

        time_reader = _BinReader('time.pas')
        time_dat = time_reader.read_all()
        time_reader.end_read()

        self.t_list = time_dat[:, 0]
        self.dt_avg_list = time_dat[1:, 1]
        self.Nt = len(self.t_list)

        # Set initial conditons
        self.__set_initial_conditions()

        self.dim = len(self.arrs.values()[0].shape)

        # Move back to previous working directory
        os.chdir(cwd)


        # *** Set numerical quantities ***

        if extent:

            if self.dim != len(extent)/2:

                raise ValueError('Can\'t make a %dD animation for %dD data.' %(len(extent)/2, self.dim))

            if self.dim == 2:

                self.Nx = self.arrs.values()[0].shape[1 if matrixLike else 0]    # Number of horizontal grid points
                self.Nz = self.arrs.values()[0].shape[0 if matrixLike else 1]    # Number of vertical grid points

                self.equal_aspect = float(self.Nx)/self.Nz

                # Arrays of grid points
                self.x, self.z = np.meshgrid(np.linspace(extent[0], extent[1], self.Nx), np.linspace(extent[2], extent[3], self.Nz))

            elif self.dim == 1:

                self.N = self.arrs.values()[0].shape[0]

                self.direction = 'z' if self.has_arr['w'] else 'x'

                self.l = np.linspace(extent[0], extent[1], self.N)

            else:

                raise ValueError('Invalid dimension of data files.')

    def __get_data_blocks(self, i):

        'Gets data block nr. i for each of the relevant quantities.'

        for name in self.arrs:

            self.arrs[name] = self.arr_readers[name].read_block(i)

        self.t = self.t_list[i]

    def __set_initial_conditions(self):

        'Reads the initial quantity values.'

        self.t_i = 0
        self.__get_data_blocks(self.t_i)

    def __close_read_files(self):

        'Close binary files that has been read from.'

        for reader in self.arr_readers.values():

            reader.end_read()

    def __step_time(self, min_t_skip):

        '''
        Step until a given amount of time has passed. The minimum possible 
        step time is the time between each frame in the binary files.
        '''

        # Calculate the number of steps that can be taken before the given step time is exceeded
        if min_t_skip <= self.t_list[self.t_i + 1] - self.t:

            skips = 1

        else:
            skips = np.sum(min_t_skip >= self.t_list[(self.t_i + 1):] - self.t)

        # Calculate the average time step length
        new_dt_avg = np.mean(self.dt_avg_list[self.t_i:(self.t_i + skips)])

        # Advance time
        self.t_i += skips
        if self.t_i >= self.Nt-1: self.t_i = 0

        self.__get_data_blocks(self.t_i)

        return new_dt_avg

    def __get_quantity_info(self, quantity, bg=False):

        '''
        Takes a string naming a quantity and returns a function returning 
        the quantity, a string with the unit of the quantity and values for 
        the highest and lowest value of the quantity to include in the 
        visualisation.
        '''

        if quantity == 'rho':

            # Specify name of the quantity
            name = 'mass density'

            # Specify unit for the quantity
            unit = r'kg$/$m$^3$'

            # Raise error if the required data isn't available
            if not self.has_arr['rho']: raise ValueError('No visualisation data available for %s.' % name)

            # Define function returning the relevant quantity
            q = lambda: self.arrs['rho']

        elif quantity == 'drho':

            name = 'mass density contrast'
            unit = r''

            if not self.has_arr['rho']: raise ValueError('No visualisation data available for %s.' % name)

            self.rho_bg = self.arrs['rho'].copy() if not (bg and 'rho' in bg) else bg['rho']

            q = lambda: (self.arrs['rho'] - self.rho_bg)/self.rho_bg

        elif quantity == 'u':

            name = 'horizontal velocity'
            unit = r'm$/$s'

            if not self.has_arr['u']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['u']

        elif quantity == 'w':

            name = 'vertical velocity'
            unit = r'm$/$s'

            if not self.has_arr['w']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['w']

        elif quantity == 'e':

            name = 'internal energy density'
            unit = r'J$/$m$^3$'

            if not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['e']

        elif quantity == 'de':

            name = 'internal energy density contrast'
            unit = r''

            if not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            self.e_bg = self.arrs['e'].copy() if not (bg and 'e' in bg) else bg['e']

            q = lambda: (self.arrs['e'] - self.e_bg)/self.e_bg

        elif quantity == 'es':

            name = 'specific internal energy'
            unit = r'J$/$kg'

            if not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['e']/self.arrs['rho']

        elif quantity == 'P':

            name = 'pressure'
            unit = r'Pa'

            if not self.has_arr['P']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['P']

        elif quantity == 'dP':

            name = 'pressure contrast'
            unit = r''

            if not self.has_arr['P']: raise ValueError('No visualisation data available for %s.' % name)

            self.P_bg = self.arrs['P'].copy() if not (bg and 'P' in bg) else bg['P']

            q = lambda: (self.arrs['P'] - self.P_bg)/self.P_bg

        elif quantity == 'T':

            name = 'temperature'
            unit = r'K'

            if not self.has_arr['T']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['T']

        elif quantity == 'dT':

            name = 'temperature contrast'
            unit = r''

            if not self.has_arr['T']: raise ValueError('No visualisation data available for %s.' % name)

            self.T_bg = self.arrs['T'].copy() if not (bg and 'T' in bg) else bg['T']

            q = lambda: (self.arrs['T'] - self.T_bg)/self.T_bg

        elif quantity == 'v':

            name = 'speed'
            unit = r'm$/$s'

            if not self.has_arr['u']:

                q = lambda: np.abs(self.arrs['w'])

            elif not self.has_arr['w']:

                q = lambda: np.abs(self.arrs['u'])

            else:

                q = lambda: np.sqrt(self.arrs['u']**2 + self.arrs['w']**2)

        elif quantity == 'ru':

            name = 'horizontal momentum density'
            unit = r'kg$/$sm$^2$'

            if not self.has_arr['u']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['rho']*self.arrs['u']

        elif quantity == 'rw':

            name = 'vertical momentum density'
            unit = r'kg$/$sm$^2$'

            if not self.has_arr['w']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['rho']*self.arrs['w']

        elif quantity == 'rv':

            name = 'momentum density'
            unit = r'kg$/$sm$^2$'

            if not self.has_arr['u']:

                q = lambda: self.arrs['rho']*np.abs(self.arrs['w'])

            elif not self.has_arr['w']:

                q = lambda: self.arrs['rho']*np.abs(self.arrs['u'])

            else:

                q = lambda: self.arrs['rho']*np.sqrt(self.arrs['u']**2 + self.arrs['w']**2)

        elif quantity == 'eu':

            name = 'horizontal energy flux'
            unit = r'W$/$m$^2$'

            if not self.has_arr['u'] or not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['e']*self.arrs['u']

        elif quantity == 'ew':

            name = 'vertical energy flux'
            unit = r'W$/$m$^2$'

            if not self.has_arr['w'] or not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            q = lambda: self.arrs['e']*self.arrs['w']

        elif quantity == 'ev':

            name = 'energy flux'
            unit = r'W$/$m$^2$'

            if not self.has_arr['e']: raise ValueError('No visualisation data available for %s.' % name)

            if not self.has_arr['u']:

                q = lambda: self.arrs['e']*np.sqrt(self.arrs['w'])

            elif not self.has_arr['w']:

                q = lambda: self.arrs['e']*np.sqrt(self.arrs['u'])

            else:

                q = lambda: self.arrs['e']*np.sqrt(self.arrs['u']**2 + self.arrs['w']**2)

        else:

            raise ValueError('Invalid quantity.')

        if self.printInfo: print '\nFluidVisualiser: Preparing visualisation data for %s ...' % name

        return q, name, unit

    def __get_optimal_scaling(self, q, includeAll=False):

        '''
        Advances the simulation a given time and examines 
        the resulting values of a given quantity to determine the 
        optimal limit values to use in the animation. The simulation 
        is then reset.
        '''

        if self.printInfo: print 'FluidVisualiser: Finding optimal scaling ...'

        # Define lists for the extremal quantity values
        q_max = [np.max(q())]
        q_min = [np.min(q())]

        skips = int(self.Nt)/100

        if skips == 0: skips = 1

        # Advance the simulation and store extremal values
        for i in xrange(1, self.Nt, skips):

            self.__get_data_blocks(i)

            q_max.append(np.max(q()))
            q_min.append(np.min(q()))

        # Reset simulation
        self.__set_initial_conditions()

        if includeAll:

            min_val = np.min(q_min)
            max_val = np.max(q_max)

        else:

            min_val = np.median(q_min)
            max_val = np.median(q_max)

        return min_val, max_val

    def __get_quiver(self, ax, N_arrows, quiverscale, showQuiver, arrowscale_unit):

        'Creates a quiver object for visualising the velocity field.'

        if showQuiver:

            # Find smallest axis size
            min_N = np.min([self.Nx, self.Nz])

            # Make sure smallest number of arrows isn't larger than smallest axis size
            if N_arrows > min_N:

                N_arrows = min_N

            # Calculate number of grid points to skip in order to get the right number of arrows
            step_q = int(min_N/N_arrows)

            # Set up reduced grid
            x_q = self.x[::step_q, ::step_q]
            z_q = self.z[::step_q, ::step_q]

            # Set arrow scale so that an arrow length of 5 reduced gridpoints 
            # corresponds to a given speed based on the sound speed and a modifier.

            # If the mean sound speed can be estimated, do that (assuming SI-units and gamma = 5/3)
            if (self.has_arr['P'] or self.has_arr['e']) and self.has_arr['rho']:

                cs_approx = np.mean(np.sqrt((1.67*self.arrs['P'] if self.has_arr['P'] else 1.11*self.arrs['e'])/self.arrs['rho']))

            # Otherwise use an arbitrary value
            else:

                cs_approx = 17000.0

            arrscale = 0.2*N_arrows*cs_approx/quiverscale

            # Create quiver plot
            quiver = ax.quiver(x_q, z_q, x_q, z_q, units='height', scale=arrscale, width=0.003, color='k')

            if arrowscale_unit != '':

                # Add text with arrow scale
                arrowtext = 'Arrow scale = %g %s' % (arrscale/(self.z[-1, 0] - self.z[0, 0]), arrowscale_unit)

                ax.text(0.995, 1.005, arrowtext, color='k', fontsize=self.fontsize, horizontalalignment='right', 
                                                 verticalalignment='bottom', transform=ax.transAxes)

        else:

            # Create a dummy quiver plot that doesn't show in the animation

            step_q = 0
            quiver = ax.quiver([-10], [10], [-1], [1])

        return step_q, quiver

    def __prepare_text(self, ax, showExtent=True):

        'Sets up textboxes for updateable info.'

        # Create textboxes for displaying updateable information
        textbox1 = ax.text(0.01, 0.97, '', color='k', fontsize=self.fontsize, horizontalalignment='left', 
                                           verticalalignment='top', transform=ax.transAxes)

        textbox2 = ax.text(0.99, 0.97, '', color='k', fontsize=self.fontsize, horizontalalignment='right', 
                                          verticalalignment='top', transform=ax.transAxes)

        textbox3 = ax.text(0.01, 0.03, '', color='k', fontsize=self.fontsize, horizontalalignment='left', 
                                           verticalalignment='bottom', transform=ax.transAxes)

        # Display number of data points above the animation

        if showExtent:

            if self.dim == 2:

                ax.text(0.005, 1.005, 'Nx = %d, Nz = %d' % (self.Nx, self.Nz), color='k', fontsize=self.fontsize, horizontalalignment='left', 
                                                                               verticalalignment='bottom', transform=ax.transAxes)

            elif self.dim == 1:

                ax.text(0.005, 1.005, 'N = %d' % (self.N), color='k', fontsize=self.fontsize, horizontalalignment='left', 
                                                           verticalalignment='bottom', transform=ax.transAxes)

            else:

                raise ValueError('Invalid dimension.')

        return textbox1, textbox2, textbox3

    def __set_foldername(self, folder):

        '''
        Set a valid folder name, either the one saved 
        in this instance or the provided one.
        '''

        if folder == 'default':

            if not self.hasSaved:

                raise ValueError('Name of the folder to read from must be specified unless save_data was used with this instance.')

        else:

            self.folder = folder

    @staticmethod
    def __s_to_hms(t_s):

        'Convert a time in seconds to a string showing hours, minutes and seconds.'

        hours = t_s/3600.0
        minutes = 60*(hours - int(hours))
        seconds = 60*(minutes - int(minutes))

        return '%02d:%02d:%02d' % (hours, minutes, seconds)

    @staticmethod
    def __print_progress(i, N, t0):

        'Prints progress info.'

        ETA = 0.0
        progress = i*100.0/N
        elapsed = time.time() - t0
        if progress >= 5: ETA = (100/progress - 1)*elapsed

        sys.stdout.write('Progress: %.1f%% | Elapsed: %s%s \r' % (progress, FluidVisualiser.__s_to_hms(elapsed), '' if progress < 5 else (' | ETA: %s' % FluidVisualiser.__s_to_hms(ETA))))
        sys.stdout.flush()

# Private classes (used for input and output)

class _BinWriter:

    '''
    Class for writing array data of arbitrary dimensionality to a binary file.
    '''

    def __init__(self, filename, block_shape, dtype_dbytes):

        '''
        Opens binary file and writes a header. The header consists of the following entries:
    
        block_dim:                             Number of dimensions of the data block arrays
        block_num:                             Total number of data blocks in the body
        block_shape[0], ... , block_shape[-1]: Number of elements along each dimension of a data block
        dtype:                                 Character indicating the data type of the elements in a data block
        dbytes:                                Number of bytes used by an element in a data block
        packtype:                              Character indicating the way the arrays have been flattened
                                               ("F": Fortran style, "C": C++ style)

        All header entries are 32 bit integers except dtype and packtype, which are 8 bit characters.
        '''

        self.block_shape = np.asarray(block_shape, dtype='i4') # Set data block shape
        self.dtype_dbytes = dtype_dbytes                       # Set element data type and precision

        self.block_num = 0                          # Set counter to zero
        self.write_time = 0                         # Set total writing time to zero

        # Make sure that the correct file extension is used

        fn_ext = '' if not '.' in filename else filename.split('.')[1]
        fn_noext = filename.split('.')[0]

        if len(fn_ext) > 0:

            if fn_ext != 'pas':

                print '_BinWriter warning (\"%s\"): Invalid file extension. Using \".pas\" instead.' \
                      % filename

        self.filename = fn_noext + '.pas'

        # Check that the input string specifying the data type and precision is complete

        if len(dtype_dbytes) < 2:

            print '_BinWriter error (\"%s\"): Invalid input \"%s\" for data type/precision. ' % (self.filename, dtype_dbytes) \
                  + 'Must be a string consisting of a letter (describing the type) followed ' \
                  + 'by a number (describing the precision), e.g. "i4" for the 4 byte signed integer type.'
            sys.exit()

        try:

            self.dtype = dtype_dbytes[0]
            self.dbytes = int(dtype_dbytes[1:])

        except ValueError:

            print '_BinWriter error (\"%s\"): Invalid input \"%s\" for data type/precision. ' % (self.filename, dtype_dbytes) \
                  + 'The part following \"%s\" must be a number.' % (dtype_dbytes[0])
            sys.exit()

        # Make sure that the specified data type and precision is valid

        if not self.dtype in ['S', 'i', 'f', 'u', 'c']:

            print '_BinWriter error (\"%s\"): Invalid data type (\"%s\") . ' % (self.filename, self.dtype) \
                  + 'Valid types are \"S\" (character), \"i\" (signed integer), \"f\" (float), \"u\" (unsigned integer) and \"c\" (complex).'
            sys.exit()

        if self.dtype == 'S' and self.dbytes != 1:

            print '_BinWriter error (\"%s\"): Invalid data precision (%d bytes) for character type. ' % (self.filename, self.dbytes) \
                  + 'The only valid precision is 1 byte.'
            sys.exit()

        elif self.dtype == 'i' and not self.dbytes in [1, 2, 4, 8]:

            print '_BinWriter error (\"%s\"): Invalid data precision (%d bytes) for signed integer type. ' % (self.filename, self.dbytes) \
                  + 'Valid precisions are 1, 2, 4 and 8 bytes.'
            sys.exit()

        elif self.dtype == 'f' and not self.dbytes in [2, 4, 8]:

            print '_BinWriter error (\"%s\"): Invalid data precision (%d bytes) for float type. ' % (self.filename, self.dbytes) \
                  + 'Valid precisions are 2, 4 and 8 bytes.'
            sys.exit()

        elif self.dtype == 'u' and not self.dbytes in [1, 2, 4, 8]:

            print '_BinWriter error (\"%s\"): Invalid data precision (%d bytes) for unsigned integer type. ' % (self.filename, self.dbytes) \
                  + 'Valid precisions are 1, 2, 4 and 8 bytes.'
            sys.exit()

        elif self.dtype == 'c' and not self.dbytes in [8, 16]:

            print '_BinWriter error (\"%s\"): Invalid data precision (%d bytes) for complex type. ' % (self.filename, self.dbytes) \
                  + 'Valid precisions are 8 and 16 bytes.'
            sys.exit()

        # Block array dimension
        self.block_dim = len(self.block_shape)

        # Start measuring time
        self.start_time = time.time()

        # Open binary file for writing output
        self.f = open(self.filename, 'wb')

        # Write first header entry
        self.f.write(struct.pack('i4', self.block_dim))

        # Skip second header entry (not yet determined)
        self.f.seek(8, os.SEEK_SET)

        # Write remaining header entries
        self.block_shape.tofile(self.f, sep='')
        self.f.write(struct.pack('c1', self.dtype))
        self.f.write(struct.pack('i4', self.dbytes))
        self.f.write(struct.pack('c1', 'C'))

    def write_block(self, block):

        'Adds a data block to the end of the file.'

        # Get initial time
        t0 = time.time()

        # Write data block
        np.asarray(block, dtype=self.dtype_dbytes).tofile(self.f, sep='')

        # Add elapsed time to total time
        self.write_time += time.time() - t0

        # Increase counter
        self.block_num += 1

    def print_header_info(self):

        print '\n*** Header info for \"%s\" ***' % (self.filename)
        print 'Data block dimension:: %dD' % (self.block_dim)
        print 'Total number of data blocks: %d' % (self.block_num)
        print 'Data block shape: ' + _Common.get_shape_string(self.block_shape, useParan=True)
        print 'Data type: ' + _Common.dtype_descript[self.dtype]
        print 'Number of bytes per element: %d' % (self.dbytes)
        print 'Pack type: C++/Python style'

        print 'I live in _BinWriter!'

    def end_write(self, printReport=False):

        'Adds the second header number, closes the file and prints writing info.'

        # Go to second header number
        self.f.seek(4, os.SEEK_SET)

        # Write number of blocks written
        self.f.write(struct.pack('i4', self.block_num))

        # Close file
        self.f.close()

        # Stop measuring time
        self.stop_time = time.time()

        if printReport:

            # Print title
            print '\n*** Writing report for \"%s\" ***' % (self.filename)

            # Print block info
            print 'Data blocks: %s %s (%d bit precision)' % (_Common.get_shape_string(self.block_shape), 
                                                             _Common.dtype_descript[self.dtype], 
                                                             8*self.dbytes)

            # Print number of blocks written
            print 'Number of blocks written: %d' % self.block_num

            # Print total amount of data written
            tot_size = self.dbytes*self.block_num*np.prod(self.block_shape.astype('i8'))/1024000.0
            print 'Amount of data written: %g MB' % tot_size

            if self.write_time != 0:

                # Print writing time
                print 'Total writing time: %g s' % (self.write_time)

                # Print writing speed
                print 'Average writing speed: %.3g MB/s' % (tot_size/self.write_time)

            # Print amount of time file was open for
            print 'File open for: %.3g s' % (self.stop_time - self.start_time)


class _BinReader:

    '''
    Class for reading binary files.
    '''

    def __init__(self, filename):

        '''
        Opens file and reads header.
        '''

        # Make sure that the correct file extension is used

        fn_ext = '' if not '.' in filename else filename.split('.')[1]
        fn_noext = filename.split('.')[0]

        if len(fn_ext) > 0:

            if fn_ext != 'pas':

                print '_BinReader error (\"%s\"): Invalid file extension. Can only read \".pas\" files.' \
                      % filename
                sys.exit()

        else:

            print '_BinReader warning (\"%s\"): No file extension supplied. Assuming extension is \".pas\".' % filename

        self.filename = fn_noext + '.pas'

        # Start measuring time
        self.start_time = time.time()

        # Open binary file
        self.f = open(self.filename, 'rb')

        # Read array dimensions from header
        self.block_dim = int(np.fromfile(self.f, count=1, dtype='i4'))

        # Move read indicator to next header entry
        self.f.seek(4, os.SEEK_SET)

        # Read array shape numbers
        self.body_shape = np.fromfile(self.f, count=(self.block_dim + 1), dtype='i4')

        # Move read indicator to next header entry
        self.f.seek(4*(self.block_dim + 2), os.SEEK_SET)

        # Read element data type
        self.dtype = str(np.fromfile(self.f, count=1, dtype='a1')[0])

        # Move read indicator to next header entry
        self.f.seek(4*(self.block_dim + 2) + 1, os.SEEK_SET)

        # Read element precision (number of bytes for each element)
        self.dbytes = int(np.fromfile(self.f, count=1, dtype='i4'))

        # Move read indicator to final header entry
        self.f.seek(4*(self.block_dim + 3) + 1, os.SEEK_SET)

        # Read pack type
        self.packtype = str(np.fromfile(self.f, count=1, dtype='a1')[0])

        # Print error and abort if the data type or precision in the file cannot be read

        if not self.dtype in ['S', 'i', 'f', 'u', 'c']:

            print '_BinReader error (\"%s\"): Cannot read the data type of this file (\"%s\") . ' % (self.filename, self.dtype) \
                  + 'Can only read \"S\" (character), \"i\" (signed integer), \"f\" (float), \"u\" (unsigned integer) and \"c\" (complex).'
            self.f.close()
            sys.exit()

        if self.dtype == 'S' and self.dbytes != 1:

            print '_BinReader error (\"%s\"): Cannot read the data precision of this file (%d byte). ' % (self.filename, self.dbytes) \
                  + 'Can only read 1 byte precision for character type.'
            self.f.close()
            sys.exit()

        elif self.dtype == 'i' and not self.dbytes in [1, 2, 4, 8]:

            print '_BinReader error (\"%s\"): Cannot read the data precision of this file (%d byte). ' % (self.filename, self.dbytes) \
                  + 'Can only read 1, 2, 4 and 8 byte precision for signed integer type.'
            self.f.close()
            sys.exit()

        elif self.dtype == 'f' and not self.dbytes in [2, 4, 8]:

            print '_BinReader error (\"%s\"): Cannot read the data precision of this file (%d byte). ' % (self.filename, self.dbytes) \
                  + 'Can only read 2, 4 and 8 byte precision for float type.'
            self.f.close()
            sys.exit()

        elif self.dtype == 'u' and not self.dbytes in [1, 2, 4, 8]:

            print '_BinReader error (\"%s\"): Cannot read the data precision of this file (%d byte). ' % (self.filename, self.dbytes) \
                  + 'Can only read 1, 2, 4 and 8 byte precision for unsigned integer type.'
            self.f.close()
            sys.exit()

        elif self.dtype == 'c' and not self.dbytes in [8, 16]:

            print '_BinReader error (\"%s\"): Cannot read the data precision of this file (%d byte). ' % (self.filename, self.dbytes) \
                  + 'Can only read 8 and 16 byte precision for complex type.'
            self.f.close()
            sys.exit()

        # Print warning if the stated pack type is unknown

        if not self.packtype.upper() in ['C', 'F']:

            print '_BinReader warning (\"%s\"): Pack type \"%s\" not recognized. Assuming C++/Python style flattening.' \
                  % self.packtype

        self.dtype_dbytes = self.dtype + ('%d' % self.dbytes)   # Data type and precision in body array
        self.start = 4*(self.block_dim + 3) + 2                 # Position of first value in body
        self.block_shape = self.body_shape[1:]                  # Shape of the data blocks
        self.block_num = self.body_shape[0]                     # Total number of data blocks
        self.block_len = np.prod(self.block_shape.astype('i8')) # Number of elements in each data block
        self.block_size = self.dbytes*self.block_len            # Number of bytes in each data block

        self.read_count = 0
        self.read_time = 0

    def read_all(self):

        '''
        Reads the entire body array and stores the array in an attribute.
        '''

        # Move read indicator to beginning of body
        self.f.seek(self.start, os.SEEK_SET)

        try:

            # Get initial time
            t0 = time.time()

            # Create array containing all the body data
            data = np.fromfile(self.f, dtype=self.dtype_dbytes)

            # Reshape array
            if self.block_num == 1:

                data = data.reshape(self.block_shape, order=self.packtype)

            else:

                data = data.reshape(self.body_shape, order=self.packtype)

        except IOError:

            print '_BinReader error (\"%s\"): File is too large to read all at once. Read in blocks instead.' \
                   % (self.filename)
            self.f.close()
            sys.exit()

        # Add elapsed time to total time
        self.read_time += time.time() - t0

        # Set number of blocks read
        self.read_count = self.block_num

        return data

    def read_block(self, i):

        '''
        Reads block nr. i of the body array and stores the array in an attribute.
        '''

        # Move read indicator to beginning of data block i
        self.f.seek(self.start + i*self.block_size, os.SEEK_SET)

        # Get initial time
        t0 = time.time()

        # Create array containing data in block i
        data = np.fromfile(self.f, count=self.block_len, dtype=self.dtype_dbytes).reshape(self.block_shape, order=self.packtype)

        # Add elapsed time to total time
        self.read_time += time.time() - t0

        # Increase counter
        self.read_count += 1

        return data

    def print_header_info(self):

        print '\n*** Header info for \"%s\" ***' % (self.filename)
        print 'Data block dimension: %dD' % (self.block_dim)
        print 'Total number of data blocks: %d' % (self.block_num)
        print 'Data block shape: ' + _Common.get_shape_string(self.block_shape, useParan=True)
        print 'Data type: ' + _Common.dtype_descript[self.dtype]
        print 'Number of bytes per element: %d' % (self.dbytes)
        print 'Pack type: ' + ('Fortran style' if self.packtype.upper() == 'F' else 'C++/Python style')

        print 'I live in _BinReader!'

    def end_read(self, printReport=False):

        '''
        Closes the file and prints reading info.
        '''

        # Close file
        self.f.close()

        # Stop measuring time
        self.stop_time = time.time()

        if printReport:

            # Print title
            print '\n*** Reading report for \"%s\" ***' % (self.filename)

            # Print block info
            print 'Data blocks: %s %s (%d bit precision)' % (_Common.get_shape_string(self.block_shape), 
                                                             _Common.dtype_descript[self.dtype], 
                                                             8*self.dbytes)

            # Print number of blocks read
            print 'Number of blocks read: %d/%d' % (self.read_count, self.block_num)

            # Print amount of data read
            read_size = self.read_count*self.block_size/1024000.0
            tot_size = self.block_num*self.block_size/1024000.0
            print 'Amount of data read: %g/%g MB' % (read_size, tot_size)

            if self.read_time != 0:

                # Print reading time
                print 'Total reading time: %g s' % (self.read_time)

                # Print reading speed
                print 'Average reading speed: %g MB/s' % (read_size/self.read_time)

            # Print amount of time file was open for
            print 'File open for: %.3g s' % (self.stop_time - self.start_time)


class _Restrictor:

    '''
    Class for overwriting some methods from _BinWriter and _BinReader 
    to prevent them from being called from a _BinAppender instance.
    '''

    def end_write(self, printReport=False):

        raise AttributeError

    def read_all(self):

        raise AttributeError

    def read_block(self, i):

        raise AttributeError

    def end_read(self, printReport=False):

        raise AttributeError


class _BinAppender(_Restrictor, _BinWriter, _BinReader):

    '''
    Class for appending array data of arbitrary dimensionality to 
    a binary file. Data blocks can be appended by calling the 
    method "write_block", inherited from the _BinWriter class.
    '''

    def __init__(self, filename, block_shape):

        'Opens a binary file, reads the header and prepares to append data to the file.'

        # Call constructor of the writing class to get header data
        _BinReader.__init__(self, filename)

        # Close the file (which is opened in reading mode)
        self.f.close()

        # Check that the data to append has the same format as the existing data

        if (np.asarray(block_shape, dtype='i4') != self.block_shape).any():

            print '_BinAppender error (\"%s\"): Input block shape [%s] must match that of the blocks already written to the file [%s].' \
                  % (self.filename, _Common.get_shape_string(np.asarray(block_shape, dtype='i4')), self.block_shape)
            sys.exit()

        if self.packtype.upper() != 'C':

            print '_BinAppender error (\"%s\"): Pack type of file to append to must be C++/Python style.' % (self.filename)
            sys.exit()

        # Start measuring time
        self.start_time = time.time()

        # Open the file in append mode
        self.f = open(self.filename, 'ab')

        self.write_time = 0
        self.block_num_orig = self.block_num
        self.block_num = 0

    def end_append(self, printReport=False):

        'Updates the second header number, closes the file and prints writing info.'

        # Close the file (which is in append mode)
        self.f.close()

        # Stop measuring time
        self.stop_time = time.time()

        # Reopen the file in a mode that permits overwriting
        self.f = open(self.filename, 'r+b')

        # Go to second header number
        self.f.seek(4, os.SEEK_SET)

        # Update total number of blocks written
        self.f.write(struct.pack('i4', self.block_num_orig + self.block_num))

        # Close file
        self.f.close()

        if printReport:

            # Print title
            print '\n*** Writing report for \"%s\" (append mode) ***' % (self.filename)

            # Print block info
            print 'Data blocks: %s %s (%d bit precision)' % (_Common.get_shape_string(self.block_shape), 
                                                             _Common.dtype_descript[self.dtype], 
                                                             8*self.dbytes)

            # Print number of blocks appended
            print 'Number of blocks appended: %d (total is now %d)' % (self.block_num, 
                                                                       self.block_num + self.block_num_orig)

            # Print total amount of data appended
            block_size_MB = self.dbytes*np.prod(self.block_shape.astype('i8'))/1024000.0
            tot_size = self.block_num*block_size_MB
            print 'Amount of data appended: %g MB (total is now %g)' % (tot_size, 
                                                                        tot_size + self.block_num_orig*block_size_MB)

            if self.write_time != 0:

                # Print writing time
                print 'Total writing time: %g s' % (self.write_time)

                # Print writing speed
                print 'Average writing speed: %.3g MB/s' % (tot_size/self.write_time)

            # Print amount of time file was open for
            print 'File open for: %.3g s' % (self.stop_time - self.start_time)


class _Common:

    'Namespace available to all the classes in this module.'
    
    # Dictionary with type descriptions
    dtype_descript = {'S':'character', 'i':'signed integer', 'f':'float', 'u':'unsigned integer', 'c':'complex'}

    @staticmethod
    def get_shape_string(arr, useParan=False):

        'Returns a pretty string representation of an input 1D array.'

        symb = ', ' if useParan else ' x '

        arr_str = '%d' % (arr[0])

        for val in arr[1:]:
            arr_str += '%s%d' % (symb, val)

        if useParan:
            arr_str = '(' + arr_str + ')'

        return arr_str

if __name__ == '__main__':
    
    None