def smooth(self, input_path='', output_path='', parameters=[], span=0.):
        """ Smooths kernels by convolving them with a Gaussian.  Wrapper over 
            xsmooth_sem utility.
        """
        if not exists(input_path):
            raise Exception

        if not exists(output_path):
            unix.mkdir(output_path)

        # apply smoothing operator
        unix.cd(self.cwd)
        for name in parameters or self.parameters:
            print ' smoothing', name
            call_solver(
                system.mpiexec(),
                PATH.SPECFEM_BIN +'/'+ 'xsmooth_sem '
                + str(span) + ' '
                + str(span) + ' '
                + name + '_kernel' + ' '
                + input_path + '/ '
                + output_path + '/ ',
                output='/dev/null')

        print ''

        # rename output files
        files = glob(output_path+'/*')
        unix.rename('_smooth', '', files)

    def _launch(self, classname, funcname, hosts='all'):
        unix.mkdir(PATH.SYSTEM)

        with open(PATH.SYSTEM+'/'+'job_id', 'w') as f:
            subprocess.call('sbatch '
                + '%s ' % PAR.SLURMARGS
                + '--job-name=%s ' % PAR.TITLE
                + '--nodes=%d ' % math.ceil(PAR.NPROC/float(PAR.NODESIZE))
                + '--ntasks-per-node=%d ' % PAR.NODESIZE
                + '--ntasks=%d ' % PAR.NPROC
                + '--time=%d ' % PAR.STEPTIME
                + self._launch_args(hosts)
                + findpath('seisflows.system') +'/'+ 'wrappers/run '
                + PATH.OUTPUT + ' '
                + classname + ' '
                + funcname + ' ',
                shell=1,
                stdout=f)

        # retrieve job ids
        with open(PATH.SYSTEM+'/'+'job_id', 'r') as f:
            line = f.readline()
            job = line.split()[-1].strip()
        if hosts == 'all' and PAR.NTASK > 1:
            return [job+'_'+str(ii) for ii in range(PAR.NTASK)]
        else:
            return [job]

    def apply_hessian(self, m, dm, h):
        """ Computes the action of the Hessian on a given vector through
          solver calls
        """
        system = sys.modules['seisflows_system']
        solver = sys.modules['seisflows_solver']
        postprocess = sys.modules['seisflows_postprocess']

        self.save('m_lcg', m + h*dm)

        solver.save(solver.split(m + h*dm), 
                PATH.HESS+'/'+'model')

        system.run('optimize', 'apply_hess',
                path=PATH.HESS)

        postprocess.write_gradient(
                path=PATH.HESS)

        self.save('g_lcg', solver.merge(solver.load(
                PATH.HESS+'/'+'gradient', suffix='_kernel')))

        # uncomment for debugging
        #if True:
        #    unix.rm(PATH.HESS+'_debug')
        #    unix.mv(PATH.HESS, PATH.HESS+'_debug')
        #    unix.mkdir(PATH.HESS)

        unix.rm(PATH.HESS)
        unix.mkdir(PATH.HESS)

        return self.hessian_product(h)

    def save(self, dict, path, parameters=['vp','vs','rho'],
             prefix='', suffix=''):
        """ 
          Saves SPECFEM2D/3D models or kernels

          :input dict :: model stored as a dictionary or Container
          :input path :: directory to which model is written
          :input parameters :: list of material parameters to be written
          :input prefix :: optional filename prefix
          :input suffix :: optional filename suffix, eg '_kernel'
        """
        unix.mkdir(path)

        # fill in any missing parameters
        missing_keys = diff(parameters, dict.keys())
        for iproc in range(self.mesh_properties.nproc):
            for key in missing_keys:
                dict[key] += self.io.read_slice(
                    PATH.MODEL_INIT, prefix+key+suffix, iproc)

        # write slices to disk
        for iproc in range(self.mesh_properties.nproc):
            for key in parameters:
                self.io.write_slice(
                    dict[key][iproc], path, prefix+key+suffix, iproc)

    def combine(self, input_path='', output_path='', parameters=[]):
        """ Sums individual source contributions. Wrapper over xcombine_sem
            utility.
        """
        if not exists(input_path):
            raise Exception

        if not exists(output_path):
            unix.mkdir(output_path)

        unix.cd(self.cwd)

        names = self.check_source_names()
        subset = [names[isrc] for isrc in self._source_subset]

        with open('kernel_paths', 'w') as f:
            f.writelines([join(input_path, dir)+'\n' for dir in subset])

        # SAGA component - include contributions from reference gradient
        remainder = list(set(self._source_names) - set(subset))

        with open('kernel_paths', 'a') as f:
            f.writelines([join(PATH.GRAD_AGG, dir)+'\n' for dir in remainder])

        for name in parameters or self.parameters:
            call_solver(
                system.mpiexec(),
                PATH.SPECFEM_BIN +'/'+ 'xcombine_sem '
                + name + '_kernel' + ' '
                + 'kernel_paths' + ' '
                + output_path)

    def submit(self, workflow):
        """Submits job
        """
        unix.mkdir(PATH.OUTPUT)
        unix.cd(PATH.OUTPUT)

        # save current state
        self.checkpoint()

        # construct resource list
        nodes = int(PAR.NTASK / PAR.NODESIZE)
        cores = PAR.NTASK % PAR.NODESIZE
        hours = int(PAR.WALLTIME / 60)
        minutes = PAR.WALLTIME % 60
        resources = 'walltime=%02d:%02d:00'%(hours, minutes)
        if nodes == 0:
            resources += ',mem=%dgb,nodes=1:ppn=%d'%(PAR.MEMORY, cores)
        elif cores == 0:
            resources += ',mem=%dgb,nodes=%d:ppn=%d'%(PAR.MEMORY, nodes, PAR.NODESIZE)
        else:
            resources += ',mem=%dgb,nodes=%d:ppn=%d+1:ppn=%d'%(PAR.MEMORY, nodes, PAR.NODESIZE, cores)

        # construct arguments list
        unix.run('qsub '
                + '-N %s '%PAR.TITLE
                + '-o %s '%(PATH.SUBMIT +'/'+ 'output.log')
                + '-l %s '%resources
                + '-j %s '%'oe'
                + findpath('system') +'/'+ 'wrappers/submit '
                + '-F %s '%PATH.OUTPUT)

 def save(self, path):
     """ Save objects to disk for later reference
     """
     fullpath = self.fullpath(path)
     unix.mkdir(fullpath)
     for key in self.objects:
         saveobj(fullpath +'/'+ key+'.p', sys.modules[key])

    def submit(self, workflow):
        """ Submits workflow
        """
        unix.mkdir(PATH.OUTPUT)
        unix.cd(PATH.OUTPUT)
        unix.mkdir(PATH.SUBMIT+'/'+'output.pbs')

        self.checkpoint()

        hours = PAR.WALLTIME/60
        minutes = PAR.WALLTIME%60
        walltime = 'walltime=%02d:%02d:00 ' % (hours, minutes)

        ncpus = PAR.NODESIZE
        mpiprocs = PAR.NODESIZE

        # prepare qsub arguments
        unix.run( 'qsub '
                + '%s ' % PAR.PBSARGS
                + '-l select=1:ncpus=%d:mpiprocs=%d ' % (ncpus, mpiprocs)
                + '-l %s ' % walltime
                + '-N %s ' % PAR.TITLE
                + '-j %s '%'oe'
                + '-o %s ' % (PATH.SUBMIT+'/'+'output.log')
                + '-V '
                + ' -- ' + findpath('seisflows.system') +'/'+ 'wrappers/submit '
                + PATH.OUTPUT)

 def save(self, path):
     """ Save objects to disk for later reference
     """
     unix.mkdir(_full(path))
     for name in self.names:
         fullfile = join(_full(path), name+'.p')
         saveobj(fullfile, sys.modules[name])

    def launch(self, classname, funcname, hosts='all'):
        unix.mkdir(PATH.SYSTEM)

        # submit job
        with open(PATH.SYSTEM+'/'+'job_id', 'w') as f:
            subprocess.call('bsub '
                + '%s ' % PAR.LSFARGS
                + '-n %d ' % PAR.NPROC 
                + '-R "span[ptile=%d]" ' % PAR.NODESIZE
                + '-W %d:00 ' % PAR.STEPTIME
                + '-J "%s' %PAR.TITLE
                + self.launch_args(hosts)
                + findpath('seisflows.system') +'/'+ 'wrapper/run '
                + PATH.OUTPUT + ' '
                + classname + ' '
                + funcname + ' ',
                shell=True,
                stdout=f)

        # retrieve job ids
        with open(PATH.SYSTEM+'/'+'job_id', 'r') as f:
            # reads one entire line from the file
            line = f.readline()
            job_buf = line.split()[1].strip()
            job = job_buf[1:-1]
        if hosts == 'all' and PAR.NSRC > 1:
            nn = range(1,PAR.NSRC+1)
            #return [job+'_'+str(ii) for ii in nn]
            return [job+'['+str(ii)+']' for ii in nn]
        else:
            return [job]

    def export_residuals(self, path):
        # hack deals with problems on parallel filesystem
        unix.mkdir(join(path, 'residuals'), noexit=True)

        src = join(self.getpath, 'residuals')
        dst = join(path, 'residuals', basename(self.getpath))
        unix.mv(src, dst)

    def export_traces(self, path, prefix='traces/obs'):
        # hack deals with problems on parallel filesystem
        unix.mkdir(join(path, 'traces'), noexit=True)

        src = join(self.getpath, prefix)
        dst = join(path, 'traces', basename(self.getpath))
        unix.cp(src, dst)

    def clip(self, path='', parameters=[], minval=-np.inf, maxval=np.inf):
        """ Clips kernels by convolving them with a Gaussian.  Wrapper over 
            xclip_sem utility.
        """
        assert exists(path)
        assert len(parameters) > 0

        unix.cd(self.getpath)
        for name in self.parameters:
            self.mpirun(
                PATH.SPECFEM_BIN +'/'+ 'xclip_sem '
                + str(minval) + ' '
                + str(maxval) + ' '
                + name + '_kernel' + ' '
                + path + '/ '
                + path + '/ ')

        # move input files
        src = path
        dst = path + '_noclip'
        unix.mkdir(dst)
        for name in self.parameters:
            unix.mv(glob(src+'/*'+name+'.bin'), dst)

        # rename output files
        unix.rename('_clip', '', glob(src+'/*'))

    def export_kernels(self, path):
        unix.cd(self.kernel_databases)

        # work around conflicting name conventions
        files = []
        files += glob('*proc??????_alpha_kernel.bin')
        files += glob('*proc??????_alpha[hv]_kernel.bin')
        files += glob('*proc??????_reg1_alpha_kernel.bin')
        files += glob('*proc??????_reg1_alpha[hv]_kernel.bin')
        unix.rename('alpha', 'vp', files)

        files = []
        files += glob('*proc??????_beta_kernel.bin') 
        files += glob('*proc??????_beta[hv]_kernel.bin')
        files += glob('*proc??????_reg1_beta_kernel.bin')
        files += glob('*proc??????_reg1_beta[hv]_kernel.bin')
        unix.rename('beta', 'vs', files)

        # hack to deal with problems on parallel filesystem
        unix.mkdir(join(path, 'kernels'), noexit=True)

        unix.mkdir(join(path, 'kernels', basename(self.getpath)))
        src = join(glob('*_kernel.bin'))
        dst = join(path, 'kernels', basename(self.getpath))
        unix.mv(src, dst)

    def smooth(self, path='', parameters=[], span=0.):
        """ Smooths kernels by convolving them with a Gaussian.  Wrapper over 
            xsmooth_sem utility.
        """
        assert exists(path)
        assert len(parameters) > 0

        # apply smoothing operator
        unix.cd(self.getpath)
        for name in parameters:
            print ' smoothing', name
            self.mpirun(
                PATH.SPECFEM_BIN +'/'+ 'xsmooth_sem '
                + str(span) + ' '
                + str(span) + ' '
                + name + '_kernel' + ' '
                + path + '/ '
                + path + '/ ',
                output=self.getpath+'/'+'OUTPUT_FILES/output_smooth_sem.txt')

        print ''

        # move input files
        src = path
        dst = path + '_nosmooth'
        unix.mkdir(dst)
        for name in self.parameters:
            unix.mv(glob(src+'/*'+name+'.bin'), dst)

        # rename output files
        unix.rename('_smooth', '', glob(src+'/*'))

    def _launch(self, classname, funcname, hosts='all'):
        unix.mkdir(PATH.SYSTEM)

        # submit job
        with open(PATH.SYSTEM+'/'+'job_id', 'w') as f:
            hours = PAR.WALLTIME/60
            minutes = PAR.WALLTIME%60
            resources = 'walltime=%02d:%02d:00 '%(hours, minutes)
            args = ('/opt/pbs/12.1.1.131502/bin/qsub '
                + '-l select=1:ncpus=32:mpiprocs=32 '
                + '-l %s '%resources
                + '-q standard '
                + '-A ERDCH38424KSC '
                + '-J 0-%s ' % PAR.NTASK
                + '-N %s ' % PAR.TITLE
                + ' -- ' + findpath('system') +'/'+ 'wrappers/runpbsdsh '
                + PATH.OUTPUT + ' '
                + classname + ' '
                + funcname + ' ')

            subprocess.call(args, shell=1, stdout=f)

        # retrieve job ids
        with open(PATH.SYSTEM+'/'+'job_id', 'r') as f:
            line = f.readline()
            job = line.split()[-1].strip()
        if hosts == 'all' and PAR.NTASK > 1:
            nn = range(PAR.NTASK)
            return [job+'_'+str(ii) for ii in nn]
        else:
            return [job]

    def setup(self):
        """ Sets up nonlinear optimization machinery
        """
        unix.mkdir(PATH.OPTIMIZE)

        # prepare output writers
        self.writer = Writer(path=PATH.OUTPUT)

        self.stepwriter = StepWriter(path=PATH.SUBMIT)

        # prepare algorithm machinery
        if PAR.SCHEME in ["NLCG"]:
            self.NLCG = NLCG(path=PATH.OPTIMIZE, maxiter=PAR.NLCGMAX, thresh=PAR.NLCGTHRESH, precond=self.precond)

        elif PAR.SCHEME in ["LBFGS"]:
            self.LBFGS = LBFGS(
                path=PATH.OPTIMIZE,
                memory=PAR.LBFGSMEM,
                maxiter=PAR.LBFGSMAX,
                thresh=PAR.LBFGSTHRESH,
                precond=self.precond,
            )

        # write initial model
        if exists(PATH.MODEL_INIT):
            src = PATH.MODEL_INIT
            dst = join(PATH.OPTIMIZE, "m_new")
            savenpy(dst, solver.merge(solver.load(src)))

    def submit(self, workflow):
        """Submits job
        """
        unix.mkdir(PATH.OUTPUT)
        unix.cd(PATH.OUTPUT)

        # save current state
        self.checkpoint()

        # construct resource list
        nodes = PAR.NTASK/PAR.NODESIZE
        cores = PAR.NTASK%PAR.NODESIZE
        hours = PAR.WALLTIME/60
        minutes = PAR.WALLTIME%60
        resources = 'walltime=%02d:%02d:00 '%(hours, minutes)
        if nodes == 0:
            resources += ',nodes=1:ppn=%d'%cores
        elif cores == 0:
            resources += ',nodes=%d:ppn=16'%nodes
        else:
            resources += ',nodes=%d:ppn=16+1:ppn=%d'%(nodes, cores)

        # construct arguments list
        unix.run('qsub '
                + '-N %s '%PAR.TITLE
                + '-o %s '%(PATH.SUBMIT +'/'+ 'output.log')
                + '-l %s '%resources
                + '-j %s '%'oe'
                + findpath('system') +'/'+ 'pbs/wrapper_qsub '
                + PATH.OUTPUT,
                shell=True)

 def checkpoint(self, path, classname, method, args, kwargs):
     """ Writes information to disk so tasks can be executed remotely
     """
     argspath = join(path, 'kwargs')
     argsfile = join(argspath, classname+'_'+method+'.p')
     unix.mkdir(argspath)
     saveobj(argsfile, kwargs)
     save()

    def write_gradient(self, path):
        """ Reads kernels and writes gradient of objective function
        """
        if not exists(path):
            unix.mkdir(path)

        self.combine_kernels(path, solver.parameters)
        self.process_kernels(path, solver.parameters)