def main():
    """Parse command line inputs, load up files, then call run() and save() to do
    the real work"""
    parser.add_argument('output_dir', default='PDBs')
    args = parser.parse_args()

    # load...
    # project
    project = Project.load_from(args.project)

    # assignments
    try:
        assignments = io.loadh(args.assignments, 'arr_0')
    except KeyError:
        assignments = io.loadh(args.assignments, 'Data')

    # states
    if -1 in args.states:
        states = np.unique(assigments[np.where(assignments != -1)])
        logger.info('Yanking from all %d states', len(states))
    else:
        # ensure that the states are sorted, and that they're unique -- you
        # can only request each state once
        states = np.unique(args.states)
        logger.info("Yanking from the following states: %s", states)

    # extract the conformations using np.random for the randomness
    confs_by_state = project.get_random_confs_from_states(
        assignments, states=states, num_confs=args.conformations_per_state,
        replacement=args.replacement)

    # save the conformations to disk, in the requested style
    save(confs_by_state=confs_by_state, states=states, style=args.style,
         format=args.format, outdir=args.output_dir)

def test_overwrite_1():
    fid, fn = tempfile.mkstemp()
    try:
        a = np.arange(10)
        b = a + 1
        io.saveh(fn, a=a)
        io.saveh(fn, b=b)
        eq(io.loadh(fn, 'a'), a)
        eq(io.loadh(fn, 'b'), b)
    except:
        raise
    finally:
        if os.path.exists(fn):
            os.close(fid)
            os.unlink(fn)

def main():
    args, atom_indices, project, project_root = parse_cmdline()

    # load all of the data from disk
    xyzlist, sampled_frames = load_trajs(project, os.path.dirname(args.project_yaml),
                                       atom_indices, args.stride, args.fraction)
    assignments = io.loadh(args.assignments, 'arr_0')
    # pick only the assignments that had their xyz data loaded
    assignments = np.concatenate([assignments[i, sampled_frames[i]] for i in range(len(sampled_frames))])

    # make sure we didn't mess up the subsampling and get nonsense data
    assert not np.any(assignments < 0), 'assignments negative? stride/sampling messed up probs. did you use a different strid than you clustered with?'
    #assert np.all(np.unique(assignments) == np.arange(np.max(assignments)+1)), 'assignments dont go from 0 to max. did you use a different strid than you clustered with?'

    n_real_atoms = len(atom_indices)
    n_padded_atoms = xyzlist.shape[2]
    assert n_padded_atoms >= n_real_atoms

    pairwise = calculate_pairwise_rmsd(xyzlist, n_real_atoms)

    print 'computing silhouette...'
    score = silhouette_score(pairwise, assignments, metric='precomputed')
    print 'silhouette score: %f' % score

    path = os.path.join(args.output, 'silhouette.dat')
    print 'saving results to flat text file (append): %s...' % path
    if not os.path.exists(args.output):
        os.makedirs(args.output)

    with open(path, 'a') as f:
        f.write('%f\n' % score)

    def load(cls, filename):
        """
        Load a previously saved CCA object

        Parameters
        ----------
        filename : str
            filename to load data from

        Returns
        -------
        cca_object : CCA
            loaded cca_object
        """
            
        filehandler = io.loadh(filename)

        regularizer = pickle.loads(filehandler['regularizer'][0])
        eta = filehandler['eta'][0]

        cca_object = cls(regularization=regularizer, regularization_strength=eta)

        if 'sol' in filehandler.keys():
            cca_object.v = filehandler['sol']
            cca_object._has_solution = True

        return cca_object

def entry_point():
    args = parser.parse_args()

    arglib.die_if_path_exists(args.output)

    try:
        assignments = io.loadh(args.assignments, 'arr_0')
        distances = io.loadh(args.distances, 'arr_0')
    except KeyError:
        assignments = io.loadh(args.assignments, 'Data')
        distances = io.loadh(args.distances, 'Data')

    trimmed = run(assignments, distances, args.rmsd_cutoff)

    io.saveh(args.output, trimmed)
    logger.info('Saved output to %s', args.output)

def test_read_stride(get_fn):
    # Read a binpos with stride=3
    fn_binpos = get_fn('frame0.binpos')
    with BINPOSTrajectoryFile(fn_binpos) as f:
        xyz = f.read(stride=3)
    xyz2 = io.loadh(get_fn('frame0.binpos.h5'), 'xyz')
    assert eq(xyz, xyz2[::3])

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('assignments', default='Macro4/MacroAssignments.h5', help='Path to an assignments file. (default=Macro4/MacroAssignments.h5)')
    parser.add_argument('--project', default='ProjectInfo.yaml', help='Path to ProjectInfo.yaml file. (default=ProjectInfo.yaml)')
    args = parser.parse_args()

    project = Project.load_from(args.project)
    t = reduce(operator.add, (project.load_traj(i) for i in range(project.n_trajs)))

    phi_angles = md.compute_dihedrals(t, [PHI_INDICES]) * 180.0 / np.pi
    psi_angles = md.compute_dihedrals(t, [PSI_INDICES]) * 180.0 / np.pi
    state_index = np.hstack(io.loadh(args.assignments)['arr_0'])

    for i in np.unique(state_index):
        pp.plot(phi_angles[np.where(state_index == i)],
                psi_angles[np.where(state_index == i)],
                'x', label='State %d' % i)


    pp.title("Alanine Dipeptide Macrostates")
    pp.xlabel(r"$\phi$")
    pp.ylabel(r"$\psi$")
    annotate()

    pp.legend(loc=1, labelspacing=0.075, prop={'size': 8.0}, scatterpoints=1,
              markerscale=0.5, numpoints=1)
    pp.xlim([-180, 180])
    pp.ylim([-180, 180])
    pp.show()

    def load(cls, tica_fn):
        """
        load a tICA solution to use in projecting data.

        Parameters:
        -----------
        tica_fn : str
            filename pointing to tICA solutions

        """
        # the only variables we need to save are the two matrices
        # and the eigenvectors / values as well as the lag time

        logger.warn("NOTE: You can only use the tICA solution, you will "
                    "not be able to continue adding data")
        f = io.loadh(tica_fn)

        metric = cPickle.loads(f["metric_string"][0])

        tica_obj = cls(f['lag'][0], prep_metric=metric)
        # lag entry is an array... with a single item

        tica_obj.timelag_corr_mat = f['timelag_corr_mat']
        tica_obj.cov_mat = f['cov_mat']

        tica_obj.vals = f['vals']
        tica_obj.vecs = f['vecs']

        tica_obj._sort()

        return tica_obj

def reduce_shells(system, num_shells):
    """
    Combine shells into final maps for average intensity, <I>/sigma(I), and number of pixels.
    Also compute mean of these parameters per resolution shell; return as shell_stats.
    """
    
    map_shape = (len(system['bins']['h']), len(system['bins']['k']), len(system['bins']['l']))
    map_keys = ["I", "I_sigI", "n_pixels"]

    file_glob = glob.glob(system['map_path'] + "temp/grid_rshell*.h5")
    filelist = sorted(file_glob, key = lambda name: int(name.split('rshell')[-1].split('_')[0]))
    assert len(filelist) == num_shells

    hkl = np.array(list(itertools.product(system['bins']['h'], system['bins']['k'], system['bins']['l'])))
    hkl_res = map_utils.compute_resolution(system['space_group'], system['cell'], hkl)
    hkl_res = hkl_res.reshape(map_shape)
    
    combined_maps, shell_stats = dict(), dict()
    shell_stats['resolution'] = np.zeros(num_shells)
    
    for key in map_keys:
        print "on key %s" %key
        combined_maps[key] = np.zeros(map_shape)
        shell_stats[key] = np.zeros(num_shells)

        for shell in range(len(filelist)):
            data = io.loadh(filelist[shell], key)
            combined_maps[key] += data
            shell_stats[key][shell] = np.median(data[data>0])

            if key == "I":
                shell_stats['resolution'][shell] = np.median(hkl_res[data>0])
            
    return combined_maps, shell_stats

def test_read_stride_2(get_fn):
    # Read a binpos with stride=3 when n_frames is supplied (different code path)
    fn_binpos = get_fn('frame0.binpos')
    with BINPOSTrajectoryFile(fn_binpos) as f:
        xyz = f.read(n_frames=1000, stride=3)
    xyz2 = io.loadh(get_fn('frame0.binpos.h5'), 'xyz')
    assert eq(xyz, xyz2[::3])

def test_read_atom_indices(get_fn):
    "Read a binpos with atom_indices as a list"
    fn_binpos = get_fn('frame0.binpos')
    with BINPOSTrajectoryFile(fn_binpos) as f:
        xyz = f.read(atom_indices=[0, 1, 2])
    xyz2 = io.loadh(get_fn('frame0.binpos.h5'), 'xyz')
    assert eq(xyz, xyz2[:, [0, 1, 2], :])

def readData(FN):
    """
This function will read data from a filename based on it's extension.
Inputs:
    1) FN: filename to find data
Outputs:
    2) data: numpy array containing data read in.

The function tries load and loadtxt from numpy, and throws an error if it cannot read the file.
Additionally if the filename is lh5 it will load the data as a Trajectory.Trajectory Object but return the array of coordinates.

    """
    if FN.split('.')[-1] == 'npy':
    	data = load(FN)
    elif FN.split('.')[-1] == 'lh5': # Assume this is a trajectory in msmbuilder
    	data = md.load(FN).xyz # Only return the coordinates.
    elif FN.split('.')[-1] == 'h5':
    	data = io.loadh( FN )
        try: data = data['arr_0']
        except: data = data['Data']
    	data = data[ where( data != -1 ) ]
    else:
    	try:
    		data = loadtxt(FN)
    	except:
    		print "\n\n dataIO.readData: Cannot read %s. Use numpy.save or numpy.savetxt. Exiting" % FN
    		exit()
    if data.shape == (): # If there is only one data point, then turn it into a list
    	data = array( [ data ] )
    return data

def test_read_atom_indices_slice(get_fn):
    "Read a binpos with atom_indices as a slice"
    fn_binpos = get_fn('frame0.binpos')
    with BINPOSTrajectoryFile(fn_binpos) as f:
        xyz = f.read(atom_indices=slice(0, 10, None))
    xyz2 = io.loadh(get_fn('frame0.binpos.h5'), 'xyz')
    assert eq(xyz, xyz2[:, 0:10, :])

def test_read_atomindices_1(get_fn, fn_xtc):
    iofile = io.loadh(get_fn('frame0.xtc.h5'), deferred=False)
    with XTCTrajectoryFile(fn_xtc) as f:
        xyz, time, step, box = f.read(atom_indices=[0, 1, 2])
    assert eq(xyz, iofile['xyz'][:, [0, 1, 2]])
    assert eq(step, iofile['step'])
    assert eq(box, iofile['box'])
    assert eq(time, iofile['time'])

def test_read_atomindices_2():
    iofile = io.loadh(get_fn('frame0.xtc.h5'), deferred=False)
    with XTCTrajectoryFile(fn_xtc) as f:
         xyz, time, step, box = f.read(atom_indices=slice(None, None, 2))
    yield lambda: eq(xyz, iofile['xyz'][:, ::2])
    yield lambda: eq(step, iofile['step'])
    yield lambda: eq(box, iofile['box'])
    yield lambda: eq(time, iofile['time'])

def test_read_1():
    with BINPOSTrajectoryFile(fn_binpos, chunk_size_multiplier=0.5) as f:
        xyz = f.read()
    with DCDTrajectoryFile(fn_dcd) as f:
        xyz2 = f.read()[0]
    xyz3 = io.loadh(get_fn('frame0.binpos.h5'), 'xyz')

    yield lambda: eq(xyz[1:], xyz2)
    yield lambda: eq(xyz, xyz3)

def test_read_chunk2():
    with XTCTrajectoryFile(fn_xtc, 'r', chunk_size_multiplier=1) as f:
        xyz, time, step, box = f.read()

    iofile = io.loadh(get_fn('frame0.xtc.h5'), deferred=False)
    yield lambda: eq(xyz, iofile['xyz'])
    yield lambda: eq(step, iofile['step'])
    yield lambda: eq(box, iofile['box'])
    yield lambda: eq(time, iofile['time'])

def test_read_stride_2():
    "read xtc with stride with n_frames"
    iofile = io.loadh(get_fn('frame0.xtc.h5'), deferred=False)
    with XTCTrajectoryFile(fn_xtc) as f:
         xyz, time, step, box = f.read(n_frames=1000, stride=3)
    yield lambda: eq(xyz, iofile['xyz'][::3])
    yield lambda: eq(step, iofile['step'][::3])
    yield lambda: eq(box, iofile['box'][::3])
    yield lambda: eq(time, iofile['time'][::3])

def test_read_chunk3(get_fn, fn_xtc):
    with XTCTrajectoryFile(fn_xtc, chunk_size_multiplier=2) as f:
        xyz, time, step, box = f.read(n_frames=100)

    iofile = io.loadh(get_fn('frame0.xtc.h5'), deferred=False)
    assert eq(xyz, iofile['xyz'][:100])
    assert eq(step, iofile['step'][:100])
    assert eq(box, iofile['box'][:100])
    assert eq(time, iofile['time'][:100])

def test_read_0():
    with BINPOSTrajectoryFile(fn_binpos) as f:
        xyz = f.read()
    with DCDTrajectoryFile(fn_dcd) as f:
        xyz2 = f.read()[0]
    xyz3 = io.loadh(get_fn("frame0.binpos.h5"), "xyz")

    yield lambda: eq(xyz[1:], xyz2)
    yield lambda: eq(xyz, xyz3)