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 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 smooth(self, path='', parameters='dummy', span=0. ):
        """ Smooths SPECFEM2D kernels by convolving them with a Gaussian
        """
        from seisflows.tools.array import meshsmooth, stack

        #assert parameters == self.parameters

        # implementing nproc > 1 would be straightforward, but a bit tedious
        #assert self.mesh.nproc == 1

        kernels = self.load(path, suffix='_kernel')
        if not span:
            return kernels

        # set up grid
        _,x = loadbypar(PATH.MODEL_INIT, ['x'], 0)
        _,z = loadbypar(PATH.MODEL_INIT, ['z'], 0)
        mesh = stack(x[0], z[0])

        for key in self.parameters:
            kernels[key] = [meshsmooth(kernels[key][0], mesh, span)]

        unix.rm(path + '_nosmooth')
        unix.mv(path, path + '_nosmooth')
        self.save(path, kernels, suffix='_kernel')

 def forward(self):
     """ Calls SPECFEM3D_GLOBE forward solver
     """
     solvertools.setpar('SIMULATION_TYPE', '1')
     solvertools.setpar('SAVE_FORWARD', '.true.')
     self.mpirun('bin/xspecfem3D')
     unix.mv(self.data_wildcard, 'traces/syn')

    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 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 save(self, path, v, backup=None):
        if backup:
            src = path +'/'+ 'gradient'
            dst = path +'/'+ 'gradient_'+backup
            unix.mv(src, dst)

        solver.save(path +'/'+ 'gradient',
                    solver.split(v),
                    suffix='_kernel')

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

        # work around conflicting name conventions
        self.rename_kernels()

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

 def clean(self):
     isready = self.solver_status()
     if isready:
         unix.rm(PATH.GRAD)
         unix.mv(PATH.FUNC, PATH.GRAD)
         unix.mkdir(PATH.FUNC)
         unix.rm(PATH.SOLVER)
         unix.mv(PATH.SOLVER+'_best', PATH.SOLVER)
     else:
         super(thrifty_inversion, self).clean()

    def clean(self):
        # can forward simulations from line search be carried over?
        self.update_status()

        if self.status==1:
            unix.rm(PATH.GRAD)
            unix.mv(PATH.FUNC, PATH.GRAD)
            unix.mkdir(PATH.FUNC)
        else:
            super(thrifty_inversion, self).clean()

    def forward(self, path='traces/syn'):
        """ Calls SPECFEM3D_GLOBE forward solver
        """
        solvertools.setpar('SIMULATION_TYPE', '1')
        solvertools.setpar('SAVE_FORWARD', '.true.')
        call_solver(system.mpiexec(), 'bin/xspecfem3D')

        if PAR.FORMAT in ['ASCII', 'ascii']:
            src = glob('OUTPUT_FILES/*.sem.ascii')
            dst = path
            unix.mv(src, dst)

    def update_aggregate_gradient(self):
        """ Update stored aggregate gradient.
        """
        names = self.check_source_names()
        subset = [names[isrc] for isrc in self._source_subset]
        print subset

        for source_name in subset:
            src = glob(join(PATH.GRAD, 'kernels', source_name, '*'))
            dst = join(PATH.GRAD_AGG, source_name)
            unix.mv(src, dst)

    def iterate_search(self):
        super(thrifty_inversion, self).iterate_search()

        isdone = optimize.isdone
        isready = self.solver_status()

        # to avoid redundant forward simulation, save solver files associated
        # with 'best' trial model
        if isready and isdone:
            unix.rm(PATH.SOLVER+'_best')
            unix.mv(PATH.SOLVER, PATH.SOLVER+'_best')

    def generate_data(self, **model_kwargs):
        """ Generates data
        """
        self.generate_mesh(**model_kwargs)

        unix.cd(self.getpath)
        setpar('SIMULATION_TYPE', '1')
        setpar('SAVE_FORWARD', '.true.')
        self.call('bin/xspecfem3D')

        unix.mv(self.data_wildcard, 'traces/obs')
        self.export_traces(PATH.OUTPUT, 'traces/obs')

    def forward(self, path='traces/syn'):
        """ Calls SPECFEM3D forward solver
        """
        setpar('SIMULATION_TYPE', '1')
        setpar('SAVE_FORWARD', '.true.')
        call_solver(system.mpiexec(), 'bin/xgenerate_databases')
        call_solver(system.mpiexec(), 'bin/xspecfem3D')

        if PAR.FORMAT in ['SU', 'su']:
            src = glob('OUTPUT_FILES/*_d?_SU')
            dst = path
            unix.mv(src, dst)

    def forward(self, path='traces/syn'):
        """ Calls SPECFEM2D forward solver
        """
        setpar('SIMULATION_TYPE', '1')
        setpar('SAVE_FORWARD', '.true.')

        call_solver(system.mpiexec(), 'bin/xmeshfem2D')
        call_solver(system.mpiexec(), 'bin/xspecfem2D')

        if PAR.FORMAT in ['SU', 'su']:
            filenames = glob('OUTPUT_FILES/*.su')
            unix.mv(filenames, path)

    def eval_func(self, path='', export_traces=False):
        """ Evaluates misfit function by carrying out forward simulation and
            comparing observations and synthetics.
        """
        unix.cd(self.getpath)
        self.import_model(path)

        self.forward()
        unix.mv(self.data_wildcard, 'traces/syn')
        preprocess.prepare_eval_grad(self.getpath)

        self.export_residuals(path)
        if export_traces:
            self.export_traces(path, prefix='traces/syn')

    def clip(self, path='', thresh=1.):
        """ Clips SPECFEM2D kernels
        """
        parts = self.load(path)
        if thresh >= 1.:
            return parts

        for key in self.parameters:
            # scale to [-1,1]
            minval = parts[key][0].min()
            maxval = parts[key][0].max()
            np.clip(parts[key][0], thresh*minval, thresh*maxval, out=parts[key][0])
        unix.mv(path, path + '_noclip')
        self.save(path, parts)

    def apply_hess(self, path=''):
        """ Computes action of Hessian on a given model vector.
        """
        # a gradient evaluation must have already been carried out
        unix.cd(self.getpath)
        unix.mkdir('traces/lcg')

        self.import_model(path)
        self.forward()
        unix.mv(self.data_wildcard, 'traces/lcg')
        preprocess.prepare_apply_hess(self.getpath)

        self.adjoint()
        self.export_kernels(path)

    def finalize_search(self):
        """ Cleans working directory and writes updated model
        """
        unix.cd(PATH.OPTIMIZE)

        m = self.load("m_new")
        g = self.load("g_new")
        p = self.load("p_new")
        s = loadtxt("s_new")

        x = self.step_lens()
        f = self.func_vals()

        # clean working directory
        unix.rm("alpha")
        unix.rm("m_try")
        unix.rm("f_try")

        if self.iter > 1:
            unix.rm("m_old")
            unix.rm("f_old")
            unix.rm("g_old")
            unix.rm("p_old")
            unix.rm("s_old")

        unix.mv("m_new", "m_old")
        unix.mv("f_new", "f_old")
        unix.mv("g_new", "g_old")
        unix.mv("p_new", "p_old")
        unix.mv("s_new", "s_old")

        # write updated model
        alpha = x[f.argmin()]
        savetxt("alpha", alpha)
        self.save("m_new", m + alpha * p)
        savetxt("f_new", f.min())

        # append latest output
        self.writer("factor", -self.dot(g, g) ** -0.5 * (f[1] - f[0]) / (x[1] - x[0]))
        self.writer("gradient_norm_L1", np.linalg.norm(g, 1))
        self.writer("gradient_norm_L2", np.linalg.norm(g, 2))
        self.writer("misfit", f[0])
        self.writer("restarted", self.restarted)
        self.writer("slope", (f[1] - f[0]) / (x[1] - x[0]))
        self.writer("step_count", self.step_count)
        self.writer("step_length", x[f.argmin()])
        self.writer("theta", 180.0 * np.pi ** -1 * angle(p, -g))

        self.stepwriter.newline()