def run_experiment():
    pattern = re.compile("lda_([0-9]+)\.pb")
    data_dir = "data"
    files = [
        (re.search(pattern, f).group(1), join(data_dir, f))
        for f in listdir(data_dir)
        if isfile(join(data_dir, f)) and re.match(pattern, f)
    ]
    cmd_str = "peircebayes {} -n 100 -m lda -t -s {}"
    cmd_str2 = "peircebayes {} -n 100 -m lda -t -s {} -a cgs"
    np.random.seed(1234)
    start = time.time()
    for i, f in files:
        print i
        # sample 10 times
        for j, seed in enumerate(np.random.choice(5000, 10, replace=False) + 1):
            call_cmd(cmd_str.format(f, seed))
            phi = np.load("/tmp/peircebayes/avg_samples.npz")["arr_1"]
            np.savez(join(data_dir, "phi_{}_{}".format(i, j)), **{"phi": phi})
            call_cmd("cp /tmp/peircebayes/lls.npz data/lls_{}_{}.npz".format(i, j))
            call_cmd(cmd_str2.format(f, seed))
            call_cmd("cp /tmp/peircebayes/lls.npz data/lls_cgs_{}_{}.npz".format(i, j))
    end = time.time()
    with open("data/time_pb", "w") as f:
        f.write(str(end - start))
    cmd_str_r = "Rscript run_lda.R"
    start = time.time()
    call_cmd(cmd_str_r)
    end = time.time()
    with open("data/time_r", "w") as f:
        f.write(str(end - start))

 def save(self, npz_file):
     '''Serialize to object to an npz file.'''
     np.savez(npz_file,
              chrom=self.chrom,
              sample_call_rate_full=self.sample.call_rate_full,
              sample_concordance=self.sample.concordance,
              sample_call_rate_partial=self.sample.call_rate_partial,
              sample_samples=self.sample.samples,
              
              snp_call_rate_imputed_full=self.snp.call_rate_imputed_full,
              snp_concordance_imputed_het=self.snp.concordance_imputed_het,
              snp_call_rate_imputed_partial=self.snp.call_rate_imputed_partial,
              snp_info=self.snp.info,
              snp_call_rate_training=self.snp.call_rate_training,
              snp_snps=self.snp.snps,
              snp_concordance_imputed=self.snp.concordance_imputed,
              snp_x=self.snp.x,
              snp_maf=self.snp.maf,
               
              maf_call_rate_imputed=self.maf.call_rate_imputed,
              maf_concordance_imputed_het=self.maf.concordance_imputed_het,
              maf_call_rate_training=self.maf.call_rate_training,
              maf_maf=self.maf.maf,
              maf_concordance_imputed=self.maf.concordance_imputed,
              # sample_index=self.sample_index,
              # pedigree=np.array([self.pedigree])
              )

 def save(self):
     numpy.savez('timing.npz',
                 train=self.train_timing,
                 valid=self.valid_timing,
                 test=self.test_timing,
                 k=self.k)
     self.model.save()

def recordDVM(filename='voltdata.npz',sun=False,moon=False,recordLength=np.inf,verbose=True):
    ra = 0
    dec = 0
    raArr = np.ndarray(0)
    decArr = np.ndarray(0)
    lstArr = np.ndarray(0)
    jdArr = np.ndarray(0)
    voltArr = np.ndarray(0)
    
    startTime = time.time()

    while np.less(time.time()-startTime,recordLength):
        if sun:
            raDec = sunPos()
            ra = raDec[0]
            dec = raDec[1]
        startSamp = time.time()
        currVolt = getDVMData()
        currLST = getLST()
        currJulDay = getJulDay()
        raArr = np.append(raArr,ra)
        decArr = np.append(decArr,ra)
        voltArr = np.append(voltArr,currVolt)
        lstArr = np.append(lstArr,currLST)
        jdArr = np.append(jdArr,currJulDay)

        if verbose:
            print 'Measuring voltage: ' + str(currVolt) + ' (LST: ' + str(currLST) +'  ' + time.asctime() + ')'
        
        np.savez(filename,ra=raArr,dec=decArr,jd=jdArr,lst=lstArr,volts=voltArr)
        sys.stdout.flush()
        time.sleep(np.max([0,1.0-(time.time()-startSamp)]))

def saveEnsemble(ensemble, filename=None, **kwargs):
    """Save *ensemble* model data as :file:`filename.ens.npz`.  If *filename*
    is ``None``, title of the *ensemble* will be used as the filename, after
    white spaces in the title are replaced with underscores.  Extension is
    :file:`.ens.npz`. Upon successful completion of saving, filename is
    returned. This function makes use of :func:`numpy.savez` function."""

    if not isinstance(ensemble, Ensemble):
        raise TypeError('invalid type for ensemble, {0}'
                        .format(type(ensemble)))
    if len(ensemble) == 0:
        raise ValueError('ensemble instance does not contain data')

    dict_ = ensemble.__dict__
    attr_list = ['_title', '_confs', '_weights', '_coords']
    if isinstance(ensemble, PDBEnsemble):
        attr_list.append('_labels')
        attr_list.append('_trans')
    if filename is None:
        filename = ensemble.getTitle().replace(' ', '_')
    attr_dict = {}
    for attr in attr_list:
        value = dict_[attr]
        if value is not None:
            attr_dict[attr] = value
            
    attr_dict['_atoms'] = np.array([dict_['_atoms'], 0])
    filename += '.ens.npz'
    ostream = openFile(filename, 'wb', **kwargs)
    np.savez(ostream, **attr_dict)
    ostream.close()
    return filename

def test_npzfile_dict():
    s = StringIO.StringIO()
    x = np.zeros((3, 3))
    y = np.zeros((3, 3))

    np.savez(s, x=x, y=y)
    s.seek(0)

    z = np.load(s)

    assert 'x' in z
    assert 'y' in z
    assert 'x' in z.keys()
    assert 'y' in z.keys()

    for f, a in z.iteritems():
        assert f in ['x', 'y']
        assert_equal(a.shape, (3, 3))

    assert len(z.items()) == 2

    for f in z:
        assert f in ['x', 'y']

    assert 'x' in list(z.iterkeys())

def save(destination, train, valid, test, vocab):
    np.savez(destination,
             vocab=np.array(vocab),
             train=train,
             valid=valid,
             test=test,
             vocab_size=len(vocab))

def get_flat_distribution(y):
    start = time.time()
    logging.debug(y.shape)
    rows,cols = y.shape
    x_min = 0 
    x_max = 2500
    bin_width = 100
    n_bins = int((x_max-x_min)/bin_width)
    lst = [bin(i,i+bin_width) for i in np.arange(x_min,x_max,bin_width)]
    new_array = np.zeros((1,2))
    new_array_list = []

    for row in range(rows):
        if row == 0:
            new_array_list.append(y[row])
            [lst[i].in_bin(y[0][0][0]) for i in range(n_bins)]
        else:
            for i in range(n_bins):
                if(lst[i].in_bin(y[row][0][0]) and not lst[i].full):
                    new_array_list.append(y[row])
            if(row%1000000 == 0):
                logging.debug(row)
    stop = time.time()
    logging.debug("time elapsed running through rows"+str(stop - start))
    new_array = np.vstack(new_array_list)
    stop = time.time()
    logging.debug("time elapsed stacking"+str(stop - start))
    logging.debug("new array shape:")
    logging.debug(new_array.shape)
    
    rows,cols = new_array.shape
    [lst[i].print_object() for i in range(n_bins)]
    np.savez(file_name+"flat", new_array)
    return new_array

def write_potential(N=2.5, pphw=20, amplitude=1.0, sigmax=1e-1, sigmay=1e-1,
                    L=100., W=1.0, x_R0=0.05, y_R0=0.4, loop_type='Bell',
                    init_phase=0.0, shape='RAP', plot=True,
                    plot_dimensions=False, direction='right',
                    boundary_only=False, with_boundary=False, boundary_phase=0.0,
                    theta=0.0, smearing=False, verbose=True, linearized=False):

    p = Potential(N=N, pphw=pphw, amplitude=amplitude, sigmax=sigmax,
                  sigmay=sigmay, x_R0=x_R0, y_R0=y_R0, init_phase=init_phase,
                  shape=shape, L=L, W=W, loop_type=loop_type,
                  direction=direction, boundary_only=boundary_only,
                  with_boundary=with_boundary, theta=theta,
                  verbose=verbose, linearized=linearized)

    if not boundary_only:
        imag, imag_vector = p.imag, p.imag_vector
        real, real_vector = p.real, p.real_vector
    X, Y = p.X, p.Y

    if not boundary_only:
        if plot:
            import matplotlib.pyplot as plt
            if plot_dimensions:
                plt.figure(figsize=(L, W))
            plt.pcolormesh(X, Y, imag, cmap='RdBu_r')
            plt.savefig("imag.png")
            plt.pcolormesh(X, Y, real, cmap='RdBu_r')
            plt.savefig("real.png")

        np.savetxt("potential_imag.dat", zip(xrange(len(imag_vector)), imag_vector),
                   fmt=["%i", "%.12f"])
        np.savetxt("potential_real.dat", zip(xrange(len(real_vector)), real_vector),
                   fmt=["%i", "%.12f"])
        if shape != 'science':
            np.savez("potential_imag_xy.npz", X=X, Y=Y, P=imag_vector,
                     X_nodes=p.xnodes, Y_nodes=p.ynodes,
                     sigmax=sigmax, sigmay=sigmay)

    if shape == 'RAP':
        xi_lower, xi_upper = p.WG.get_boundary(theta=theta, smearing=smearing,
                                               boundary_phase=boundary_phase)
        # set last element to 0 (xi_lower) or W (xi_upper)
        print "WARNING: end of boundary not set zero!"
        # xi_lower[-1] = 0.0
        # xi_upper[-1] = W
        np.savetxt("upper.boundary", zip(xrange(p.nx), xi_upper))
        np.savetxt("lower.boundary", zip(xrange(p.nx), xi_lower))
        eps, delta = p.WG.get_cycle_parameters()
        np.savetxt("boundary.eps_delta", zip(eps, delta))
    if shape == 'RAP_TQD':
        eps_prime, delta_prime, theta_prime = p.WG.get_quantum_driving_parameters()
        xi_lower, xi_upper = p.WG.get_boundary(eps=eps_prime, delta=delta_prime,
                                               theta=theta_prime,
                                               smearing=smearing)
        # set last element to 0 (xi_lower) or W (xi_upper)
        xi_lower[-1] = 0.0
        xi_upper[-1] = W
        np.savetxt("upper.boundary", zip(xrange(p.nx), xi_upper))
        np.savetxt("lower.boundary", zip(xrange(p.nx), xi_lower))
        np.savetxt("boundary.eps_delta_theta", zip(eps_prime, delta_prime, theta_prime))

    def __call__(self, u, x, t, n):
        # Save solution u to a file using numpy.savez
        if self.filename is not None:
            name = 'u%04d' % n  # array name
            kwargs = {name: u}
            fname = '.' + self.filename + '_' + name + '.dat'
            np.savez(fname, **kwargs)
            self.t.append(t[n])  # store corresponding time value
            if n == 0:           # save x once
                np.savez('.' + self.filename + '_x.dat', x=x)

        # Animate
        if n % self.skip_frame != 0:
            return
        # Plot u and mark medium x=x_L and x=x_R
        x_L, x_R = self.medium
        umin, umax = self.yaxis
        title = 'Nx=%d' % (x.size-1)
        if self.title:
            title = self.title + ' ' + title
        if self.backend is None:
            # native matplotlib animation
            if n == 0:
                self.plt.ion()
                self.lines = self.plt.plot(
                    x, u, 'r-',
                    [x_L, x_L], [umin, umax], 'k--',
                    [x_R, x_R], [umin, umax], 'k--')
                self.plt.axis([x[0], x[-1],
                               self.yaxis[0], self.yaxis[1]])
                self.plt.xlabel('x')
                self.plt.ylabel('u')
                self.plt.title(title)
                self.plt.text(0.75, 1.0, 'C=0.25')
                self.plt.text(0.32, 1.0, 'C=1')
                self.plt.legend(['t=%.3f' % t[n]])
            else:
                # Update new solution
                self.lines[0].set_ydata(u)
                self.plt.legend(['t=%.3f' % t[n]])
                self.plt.draw()
        else:
            # scitools.easyviz animation
            self.plt.plot(x, u, 'r-',
                          [x_L, x_L], [umin, umax], 'k--',
                          [x_R, x_R], [umin, umax], 'k--',
                          xlabel='x', ylabel='u',
                          axis=[x[0], x[-1],
                                self.yaxis[0], self.yaxis[1]],
                          title=title,
                          show=self.screen_movie)
        # pause
        if t[n] == 0:
            time.sleep(2)  # let initial condition stay 2 s
        else:
            if self.pause is None:
                pause = 0.2 if u.size < 100 else 0
            time.sleep(pause)

        self.plt.savefig('frame_%04d.png' % (n))

def packageMergedSpec():
    dataDir = getPackageDir('SIMS_SKYBRIGHTNESS_DATA')
    outDir = os.path.join(dataDir, 'ESO_Spectra/MergedSpec')

    # A large number of the background components only depend on Airmass, so we can merge those together
    npzs = ['LowerAtm/Spectra.npz',
            'ScatteredStarLight/scatteredStarLight.npz',
            'UpperAtm/Spectra.npz']
    files = [os.path.join(dataDir, 'ESO_Spectra', npz) for npz in npzs]


    temp = np.load(files[0])

    wave = temp['wave'].copy()
    spec = temp['spec'].copy()
    spec['spectra'] = spec['spectra']*0.
    spec['mags'] = spec['mags']*0.

    for filename in files:
        restored = np.load(filename)
        spec['spectra'] += restored['spec']['spectra']
        try:
            flux = 10.**(-0.4*(restored['spec']['mags']-np.log10(3631.)))
        except:
            import pdb ; pdb.set_trace()
        flux[np.where(restored['spec']['mags'] == 0.)] = 0.
        spec['mags'] += flux

    spec['mags'] = -2.5*np.log10(spec['mags'])+np.log10(3631.)

    np.savez(os.path.join(outDir,'mergedSpec.npz'), spec=spec, wave=wave, filterWave=temp['filterWave'])

    def save(self, path):
        savedir = smartutils.create_folder(pjoin(path, type(self).__name__))
        smartutils.save_dict_to_json_file(pjoin(savedir, "hyperparams.json"), self.hyperparameters)

        params = {param.name: param.get_value() for param in self.parameters}
        assert len(self.parameters) == len(params)  # Implies names are all unique.
        np.savez(pjoin(savedir, "params.npz"), **params)

def convert_npys_to_npzs(npy_files, arr_key, output_dir):
    """Create a number of NPY files.

    Parameters
    ----------
    npy_files = list of str
        Paths to the created set of files.

    arr_key : str
        Name to write the array under in the npz archive.

    output_dir : str
        Path under which to write data.

    Returns
    -------
    npz_files : list of str
        Newly created NPZ files.
    """
    npz_files = []
    for fpath in npy_files:
        data = {arr_key: np.load(fpath)}
        npz_path = os.path.join(output_dir, "{}.npz".format(filebase(fpath)))
        np.savez(npz_path, **data)
        npz_files.append(npz_path)

    return npz_files

    def save(self, file):
        """
        Saves data from a CorpusSent object as an `npz` file.
        
        :param file: Designates the file to which to save data. See
            `numpy.savez` for further details.
        :type file: str-like or file-like object
            
        :returns: None

        :See Also: :class: Corpus, :meth: Corpus.save, :meth: numpy.savez
        """
	
	print 'Saving corpus as', file
        arrays_out = dict()
        arrays_out['corpus'] = self.corpus
        arrays_out['words'] = self.words
        arrays_out['sentences'] = self.sentences
        arrays_out['context_types'] = np.asarray(self.context_types)

        for i,t in enumerate(self.context_data):
            key = 'context_data_' + self.context_types[i]
            arrays_out[key] = t

        np.savez(file, **arrays_out)

	def save(self, directory, suffix=""):
		fname = "mlp_" + ("_".join([str(l) for l in self.__layer_sizes])) + suffix + ".npz"
		path = os.path.join(directory, fname)
		with open(path, "w") as f:
			mats_to_save = [np.array(self.__layer_sizes)] + self.__W + self.__b
			np.savez(f, *mats_to_save)
		return path

def getData(gps, dt, target_name, channel_names, chunk_size):
    print 'Acquiring Data...\n'

    target = getChannel('trend', target_name, gps, dt)
    bad = []
    chunks = [np.asarray(range(len(channel_names))[i:i + chunk_size])
              for i in xrange(0, len(channel_names), chunk_size)]

    for i in xrange(len(chunks)):
        gc.collect()
        good = []
        aux = np.empty((dt, len(chunks[i])), object)
        for j in xrange(len(chunks[i])):
            aux[:, j] = getChannel('trend', channel_names[chunks[i][j]], gps,
                                   dt).data
            if aux[:, j][0] != np.mean(aux[:, j]):
                good.append(j)
            else:
                print 'Constant Channel: ' + channel_names[chunks[i][j]]
                bad.append(i * chunk_size + j)
            if (i * len(chunks[i]) + j) % 10 == 0:
                print '\nLoaded ' + str(int(((i * chunk_size + j) /
                                        float(len(channel_names))) * 100))
                                        + '% : ' + str(i * chunk_size + j)
                                        + ' of ' + str(len(channel_names))
                                        + '\n'

        aux = aux[:, good]
        channel_names[chunks[i]] = channel_names[chunks[i][good]]
        np.savez('nonna_storage/data_chunk_' + str(i), aux=aux,
                 channel_names=channel_names[chunks[i]])

    def assertCubeDataAlmostEqual(self, cube, reference_filename, *args, **kwargs):
        reference_path = self.get_result_path(reference_filename)
        if self._check_reference_file(reference_path):
            kwargs.setdefault("err_msg", "Reference file %s" % reference_path)

            result = np.load(reference_path)
            if isinstance(result, np.lib.npyio.NpzFile):
                self.assertIsInstance(cube.data, ma.MaskedArray, "Cube data was not a masked array.")
                # Avoid comparing any non-initialised array data.
                data = cube.data.filled()
                np.testing.assert_array_almost_equal(data, result["data"], *args, **kwargs)
                np.testing.assert_array_equal(cube.data.mask, result["mask"])
            else:
                np.testing.assert_array_almost_equal(cube.data, result, *args, **kwargs)
        else:
            self._ensure_folder(reference_path)
            logger.warning("Creating result file: %s", reference_path)
            if isinstance(cube.data, ma.MaskedArray):
                # Avoid recording any non-initialised array data.
                data = cube.data.filled()
                with open(reference_path, "wb") as reference_file:
                    np.savez(reference_file, data=data, mask=cube.data.mask)
            else:
                with open(reference_path, "wb") as reference_file:
                    np.save(reference_file, cube.data)

def curvature(outfile, subject, smooth=20, **kwargs):
    curvs = []
    for pts, polys in db.get_surf(subject, "fiducial"):
        surf = polyutils.Surface(pts, polys)
        curv = surf.smooth(surf.mean_curvature(), smooth)
        curvs.append(curv)
    np.savez(outfile, left=curvs[0], right=curvs[1])

    def mutation_histogram(self, records, mut, filename):
        """Records is a pd.Dataframe."""
        if os.path.exists(filename + '.npz'):
            logging.critical(filename + '.npz esists.')
            return filename
        if records.shape[0] < self.min_seqs:
            return ''

        igs = [IgRecord(x.to_dict()) for _, x in records.iterrows()]
        igsimilarity_learn = copy.deepcopy(self.igsimilarity)
        igsimilarity_learn.correct = self.correction
        igsimilarity_learn.rm_duplicates = True
        if not self.correction:
            igsimilarity_learn.tol = 1000
        else:
            igsimilarity_learn.correct_by = self.correction

        sim_func = igsimilarity_learn.pairwise
        logging.info("Computing %s", filename)
        dnearest = parallel_distance.dnearest_intra_padding(
            igs, sim_func, filt=lambda x: x > 0, func=max)

        if not os.path.exists(filename.split('/')[0]):
            os.makedirs(filename.split('/')[0])
        np.savez(filename, X=dnearest, mut=mut)

        # Plot distance distribution
        title = "Similarities for {:.3f}-{:.3f}%" \
                .format(np.min(records['MUT']), np.max(records['MUT']))
        plot_hist(dnearest, self.bins, title, filename)
        return filename

def tissots_indicatrix(outfile, sub, radius=10, spacing=50, maxfails=100): 
    tissots = []
    allcenters = []
    for hem in ["lh", "rh"]:
        fidpts, fidpolys = db.get_surf(sub, "fiducial", hem)
        #G = make_surface_graph(fidtri)
        surf = polyutils.Surface(fidpts, fidpolys)
        nvert = fidpts.shape[0]
        tissot_array = np.zeros((nvert,))

        centers = [np.random.randint(nvert)]
        cdists = [surf.geodesic_distance(centers)]
        while True:
            ## Find possible vertices
            mcdist = np.vstack(cdists).min(0)
            possverts = np.nonzero(mcdist > spacing)[0]
            #possverts = np.nonzero(surf.geodesic_distance(centers) > spacing)[0]
            if not len(possverts):
                break
            ## Pick random vertex
            centervert = possverts[np.random.randint(len(possverts))]
            centers.append(centervert)
            print("Adding vertex %d.." % centervert)
            dists = surf.geodesic_distance([centervert])
            cdists.append(dists)

            ## Find appropriate set of vertices
            selverts = dists < radius
            tissot_array[selverts] = 1

        tissots.append(tissot_array)
        allcenters.append(np.array(centers))

    np.savez(outfile, left=tissots[0], right=tissots[1], centers=allcenters)