def read_trj(trj_file):
    """
    return: Simulation
    parameters:
        trj_file: string | name of trj file
    """
    simulation = Simulation()
    with open(trj_file, 'r') as trj:

        while True:

            line = trj.readline()
            if not line:
                break
            lattice = Lattice()
            lattice.set_a(np.array(line.split(), dtype=float))
            lattice.set_b(np.array(trj.readline().split(), dtype=float))
            lattice.set_c(np.array(trj.readline().split(), dtype=float))

            configuration = Configuration(lattice=lattice)

            atom_types = trj.readline().split()
            atom_counts = np.array(trj.readline().split(), dtype=int)
            natom = np.sum(atom_counts)

            for i in xrange(natom):
                atom_record = trj.readline().split()
                atom_name = atom_record[0]
                atom_position = np.array(atom_record[1:], dtype=float)
                configuration.insert_atom(Atom(atom_name, atom_position))

            simulation.insert_configuration(configuration)

    return simulation

def do_main(graph, teams, model):
  print "\n\nGraph:", graph, "  Teams:", teams

  # Create the adjacency list for the graph.
  adj_list = create_adj_list(graph)
  # Read in the node selection for each team.
  team_nodes = read_nodes(graph, adj_list.keys(), teams)

  # Run the simulation and output the run to file.
  simulation = Simulation(model, team_nodes, adj_list)
  (output, results) = simulation.run()
  output_filename = graph + "-" + str(time.time()) + ".txt"
  output_file = open(OUTPUT_FOLDER + output_filename, "w")
  output_file.write(str(json.dumps(output)))
  output_file.close()

  # Get the final results of teams to their nodes and update their points in
  # the database.
  scores = update_points(results)
  db.test.runs.insert({ \
    "teams": teams, \
    "scores": scores, \
    "graph": graph, \
    "file": output_filename \
  })

    def __init__(self,**kwargs):
        if not 'sims' in kwargs:
            self.error('sims must be provided')
        #end if
        sims = kwargs['sims']
        self.sims = sims
        del kwargs['sims']
        files = set()
        for sim in sims:
            files = files | sim.files
        #end for
        kwargs['files'] = files

        inputs = []
        filenames = []
        for sim in sims:
            inputs.append(sim.input)
            filenames.append(sim.infile)
        #end for
        kwargs['input'] = BundledQmcpackInput(inputs=inputs,filenames=filenames)

        Simulation.__init__(self,**kwargs)
        deps = []
        for sim in sims:
            for dep in sim.dependencies:
                deps.append((dep.sim,'other'))
            #end for
        #end for
        self.depends(*deps)

def explore(fpath):
    _, ext = splitext(fpath)
    ftype = 'data' if ext in ('.h5', '.hdf5') else 'simulation'
    print("Using %s file: '%s'" % (ftype, fpath))
    if ftype == 'data':
        globals_def, entities = entities_from_h5(fpath)
        data_source = H5Data(None, fpath)
        h5in, _, globals_data = data_source.load(globals_def, entities)
        h5out = None
        simulation = Simulation(globals_def, None, None, None, None,
                                entities.values(), None)
        period, entity_name = None, None
    else:
        simulation = Simulation.from_yaml(fpath)
        h5in, h5out, globals_data = simulation.load()
        period = simulation.start_period + simulation.periods - 1
        entity_name = simulation.default_entity
    entities = simulation.entities_map
    if entity_name is None and len(entities) == 1:
        entity_name = entities.keys()[0]
    if period is None and entity_name is not None:
        entity = entities[entity_name]
        period = max(entity.output_index.keys())
    eval_ctx = EvaluationContext(simulation, entities, globals_data, period,
                                 entity_name)
    try:
        c = Console(eval_ctx)
        c.run()
    finally:
        h5in.close()
        if h5out is not None:
            h5out.close()

    def do_the_walk(self):
        sim_param = self.sim_parameters

        # current insertion genes: (new genes)
        current_insertion_gene = sim_param.num_genes + 1

        # start genome:
        current_genome = Genome.identity(sim_param.num_genes, sim_param.num_chr, name="G_0")

        # add copy number information to track orthologous/paralogous, when duplications are present:
        for chromosome in current_genome.chromosomes:
            chromosome.copy_number = [1] * len(chromosome.gene_order)
        current_copy_number = current_genome.gene_count()

        # do some pre-dups if necessary:
        if sim_param.pre_duplications > 0:
            for i in range(sim_param.pre_duplications):
                Simulation.apply_random_segmental_duplication(current_genome, range(1, param.duplication_length+1), current_copy_number)

        self.genomes.append(current_genome)

        for key in RandomWalk.number_of_events:
            self.events[key].append(0)

        # walk:
        for step in range(self.length):
            # apply random event on current;
            current_genome = current_genome.clone("G_%d" % (step+1))
            n_rearrangements, n_insertions, n_deletions, n_duplications, current_insertion_gene = \
                Simulation.apply_random_events(sim_param, current_genome, self.steps, current_insertion_gene, current_copy_number)
            for key, value in zip(RandomWalk.number_of_events,
                                  [self.steps, n_rearrangements, n_insertions, n_deletions, n_duplications]):
                self.events[key].append(value)
            self.genomes.append(current_genome)

def main(argv):
    R = 500

    try:
        opts, args = getopt.getopt(argv, "hr:", ["help", "reps="])
    except getopt.GetoptError:
                sys.exit(2)
    for opt, arg in opts:
        if opt in ('-h', '--help'):
            usage()
            sys.exit(2)
        if opt in ('-r', '-replications'):
            R = int(arg)

    # Instantiate ruler ideology parameters
    rid0 = float(np.random.beta(2, 2, 1))
    rid1 = float(np.random.beta(2, 2, 1))
    params0 = {'ideology': 1, 'quality': 0, 'seniority': 0}
    params1 = {'ideology': 1, 'quality': 0, 'seniority': 0}
    leonidasr = Ruler(rid0, params0)
    spartar = Army(3, 3, 4, 30, [2, 4], leonidasr)
    leonidasl = Ruler(rid1, params1)
    spartal = Army(3, 3, 4, 30, [4, 2], leonidasl)

    print('Replication: ruler0-params {}, \
                        ruler1-params {}'.
          format(rid0,
                 rid1))
    for oo in [True]:
        # print 'Inits: {}, Ordered: {}'.format(params, oo)
        population = Population().population
        sargs = {'R': R, 'method': 'none'}
        simp = Simulation(spartar, spartal, population, sargs)
        simp.run()
        simp.write()

def sim_from_file(filename):
    """
    The purpose of this function is to:

    + open a file containing a pickled dictionary of input values to a simulation,

    + initialize the objects which the corresponding `py_qcode.Simulation` takes as input,
    
    + run the simulation, and 

    + save the results to a file of the same name as the input, with a different extension.  
    """
    #Obsolete, scavenging code for pickle-independent implementation
    #(you can't pickle functions).
    with open(filename,'r') as phil:
        sim_dict = pkl.load(phil)
    sim = Simulation(**sim_dict)
    sim.run()

    split_name = filename.split('.')
    try:
        file_prefix, file_ext = split_name
    except ValueError:
        raise ValueError('Filenames are assumed to be of the form'+\
        ' "prefix.ext".')

    output_name = '.'.join([file_prefix, 'out'])

    sim.save(output_name)

def test_spinstell_spinstell_ampa():
    netdata = TraubFullNetData()
    sim = Simulation('spinstell_spinstell_synapse')
    spinstell_index = netdata.celltype.index('SpinyStellate')
    pre = SpinyStellate(SpinyStellate.prototype, sim.model.path + '/SpinyStellate1')
    spinstell = SpinyStellate(SpinyStellate.prototype, sim.model.path + '/SpinyStellate2')
    precomp = pre.comp[SpinyStellate.presyn]
    postcomp = spinstell.comp[5] # 5 is among the allowed post synaptic compartments in spiny stellate cell
    synchan = precomp.makeSynapse(postcomp, 
                                  name='ampa_from_SPINSTELL', 
                                  classname='SynChan', 
                                  Ek=0.0,
                                  Gbar=netdata.g_ampa_baseline[spinstell_index][spinstell_index],
                                  tau1=netdata.tau_ampa[spinstell_index][spinstell_index],
                                  tau2=netdata.tau_ampa[spinstell_index][spinstell_index],
                                  delay = synapse.SYNAPTIC_DELAY_DEFAULT
                                  )
    stim = pre.soma.insertPulseGen('stimulus', sim.model, firstLevel=1e-9, firstDelay=200e-3, firstWidth=2e-3)
    pre_soma_tab = pre.soma.insertRecorder('stim', 'Vm', sim.data)
    ss_soma_tab = spinstell.soma.insertRecorder('post_soma', 'Vm', sim.data)
    ss_dend_tab = postcomp.insertRecorder('post_dend', 'Vm', sim.data)
    sim.schedule()
    sim.run(1.0)
    pylab.plot(numpy.linspace(0, 1.0, len(pre_soma_tab)), pre_soma_tab, label='pre_soma')
    pylab.plot(numpy.linspace(0, 1.0, len(ss_soma_tab)), ss_soma_tab, label='ss_soma')
    pylab.plot(numpy.linspace(0, 1.0, len(ss_dend_tab)), ss_dend_tab, label='ss_dend')
    pylab.legend()
    pylab.show()

def square_toric_code_sim(size, error_rate, n_trials, filename):
    """
    This function is square in more than one sense; it does everything
    the most vanilla way possible, and it uses a square grid to define 
    the torus. You put in an integer size, an error rate and a number
    of trials to execute, and it produces a pickled dict containing 
    the input to a simulation object in a file.
    """
    
    sim_lattice = SquareLattice((size,size))
    sim_dual_lattice = SquareLattice((size,size), is_dual=True)
    sim_model = depolarizing_model(error_rate)
    sim_code = toric_code(sim_lattice, sim_dual_lattice)
    sim_decoder = mwpm_decoder(sim_lattice, sim_dual_lattice)
    sim_log_ops = toric_log_ops((size,size))

    sim_keys = ['lattice', 'dual_lattice', 'error_model', 'code', 
                            'decoder', 'logical_operators', 'n_trials']

    sim_values = [sim_lattice, sim_dual_lattice, sim_model, sim_code, 
                                sim_decoder, sim_log_ops, n_trials]
    
    sim_dict = dict(zip(sim_keys, sim_values))

    sim = Simulation(**sim_dict)
    sim.run()
    sim.save(filename + '.sim')

def basic_test():
	test_geom = Geometry(1,1,[Rod(0,0,Dielectric(11.8), \
	occupancy_radius(0.3,1))])
	sim = Simulation('basic_test', test_geom)
	sim.run_simulation()
	sim.post_process()
	return True

    def test_single_cell(cls):
        """Simulates a single superficial pyramidal regular spiking
        cell and plots the Vm and [Ca2+]"""

        config.LOGGER.info("/**************************************************************************")
        config.LOGGER.info(" *")
        config.LOGGER.info(" * Simulating a single cell: %s" % (cls.__name__))
        config.LOGGER.info(" *")
        config.LOGGER.info(" **************************************************************************/")
        sim = Simulation(cls.__name__)
        mycell = SupPyrRS(SupPyrRS.prototype, sim.model.path + "/SupPyrRS")
        config.LOGGER.info('Created cell: %s' % (mycell.path))
        vm_table = mycell.comp[SupPyrRS.presyn].insertRecorder('Vm_suppyrrs', 'Vm', sim.data)
        pulsegen = mycell.soma.insertPulseGen('pulsegen', sim.model, firstLevel=3e-10, firstDelay=50e-3, firstWidth=50e-3)

        sim.schedule()
        if mycell.has_cycle():
            config.LOGGER.warning("WARNING!! CYCLE PRESENT IN CICRUIT.")
        t1 = datetime.now()
        sim.run(200e-3)
        t2 = datetime.now()
        delta = t2 - t1
        if config.has_pylab:
            mus_vm = config.pylab.array(vm_table) * 1e3
            mus_t = linspace(0, sim.simtime * 1e3, len(mus_vm))
            try:
                nrn_vm = config.pylab.loadtxt('../nrn/mydata/Vm_deepLTS.plot')
                nrn_t = nrn_vm[:, 0]
                nrn_vm = nrn_vm[:, 1]
                config.pylab.plot(nrn_t, nrn_vm, 'y-', label='nrn vm')
            except IOError:
                print 'NEURON Data not available.'
            config.pylab.plot(mus_t, mus_vm, 'g-.', label='mus vm')
            config.pylab.legend()
            config.pylab.show()

def explore(fpath):
    _, ext = splitext(fpath)
    ftype = 'data' if ext in ('.h5', '.hdf5') else 'simulation'
    print("Using {} file: '{}'".format(ftype, fpath))
    if ftype == 'data':
        globals_def, entities = entities_from_h5(fpath)
        simulation = Simulation(globals_def, None, None, None, None,
                                entities.values(), 'h5', fpath, None)
        period, entity_name = None, None
    else:
        simulation = Simulation.from_yaml(fpath)
        # use output as input
        simulation.data_source = H5Source(simulation.data_sink.output_path)
        period = simulation.start_period + simulation.periods - 1
        entity_name = simulation.default_entity
    dataset = simulation.load()
    data_source = simulation.data_source
    data_source.as_fake_output(dataset, simulation.entities_map)
    data_sink = simulation.data_sink
    entities = simulation.entities_map
    if entity_name is None and len(entities) == 1:
        entity_name = entities.keys()[0]
    if period is None and entity_name is not None:
        entity = entities[entity_name]
        period = max(entity.output_index.keys())
    eval_ctx = EvaluationContext(simulation, entities, dataset['globals'],
                                 period, entity_name)
    try:
        c = Console(eval_ctx)
        c.run()
    finally:
        data_source.close()
        if data_sink is not None:
            data_sink.close()

 def skip_test_handling_page_table_full(self):
     Simulation._load_memory_accesses_from_file = stub_empty_method
     Simulation._load_replacement_algorithm = stub_empty_method
     s = Simulation('dummy path', 'dummy alg', 4)
     a = MemoryAccess('00000000', 'R', 0)
     e = PageTableFullException()
     s._handle_page_fault(e, a)
     assert s.page_fault_counter == 1

 def run_simulations(self):
     """Run NUM_SIMULATIONS simulations"""
     self.results = []
     append = self.results.append
     for _ in xrange(self.num_simulations):
         simulation = Simulation(self.num_simulations, self.attrition, self.iterations_per_simulation, 
             self.promotion_bias, self.num_positions_at_level, self.bias_favors_this_gender)
         simulation.run()
         append(simulation.get_result())

 def testPBAndDIncreaseWithLambda(self):
     print "======= TEST: PB and D increase when lambda increases ======="
     lambdas = [10, 20, 40, 80, 100]
     for lamb in lambdas:
         print "Lambda = " + str(lamb)
         sim = Simulation(0, 3, 1000, 20, 5, 30, "poisson", lamb, 0, 0)
         sim.simulate()
         print ""
     print "============================================================="

 def __init__(self,**sim_args):
     has_group_atoms = 'group_atoms' in sim_args
     group_atoms = True
     if has_group_atoms:
         group_atoms = sim_args['group_atoms']
         del sim_args['group_atoms']
     #end if
     Simulation.__init__(self,**sim_args)
     if group_atoms and isinstance(self.system,PhysicalSystem):
         self.system.structure.group_atoms()

 def __init__(self,**sim_args):
     has_group_atoms = 'group_atoms' in sim_args
     group_atoms = False
     if has_group_atoms:
         group_atoms = sim_args['group_atoms']
         del sim_args['group_atoms']
     #end if
     Simulation.__init__(self,**sim_args)
     if group_atoms and isinstance(self.system,PhysicalSystem):
         self.warn('requested grouping by atomic species, but pwscf does not group atoms anymore!')

class Billiard:
	def new_game(self):
		self.scores = [0, 0]
		self.sim.add_ball(0, 0.5, 0.5, 0.0, 0.0)

		#build balls triangle
		for i in range(6):
			for j in range(i):
				self.sim.add_ball((i*(i-1))/2+j+1, 1.3+i*0.06, 0.5+j*0.06-i*0.03, 0.0, 0.0)
		
	def __init__(self):
		self.gui = Interface()
		self.gui.start()
		self.sim = Simulation()
		clock = pygame.time.Clock()
		self.gui.current_player = 0

		while not self.gui.done:
			current_player = self.gui.current_player

			#start new game if requested
			if self.gui.new_game_request:
				self.gui.new_game_request = False
				self.new_game()
			self.gui.balls = {}

			#has current player changed?
			if not self.gui.stable and self.sim.is_stable():
				current_player = (current_player+1)%2
				self.gui.current_player = current_player
			self.gui.stable = self.sim.is_stable()

			#update ball positions
			for label, ball in self.sim.balls.iteritems():
				self.gui.balls[label] = ball.pos

			#read shot command from interface and execute them
			if len(self.gui.shots) != 0:
				(angle, power) = self.gui.shots.pop()
				v = Conf.VMAX*power
				angle = (angle/180.0)*math.pi
				self.sim.balls[0].x_velocity = -v*math.sin(angle)/Conf.FPS
				self.sim.balls[0].y_velocity = -v*math.cos(angle)/Conf.FPS
			
			#check if player hit any pockets and update score
			res = self.sim.next_iter()
			if 0 in [p[0] for p in res]:
				self.sim.add_ball(0, 0.5, 0.5, 0.0, 0.0)
				self.scores[current_player] -= 1
			for ball, pocket in res:
				if ball != 0:
					self.scores[current_player] += 1
			self.gui.scores = self.scores

			clock.tick(Conf.FPS)

def test_simulation_can_have_a_room_and_roombas(room):
    roomba = Roomba()
    sim = Simulation(room=room, roombas=[roomba])


    assert sim.room is room
    assert sim.roombas == [roomba]
    sim.step()

    assert roomba.position != start_position
    assert roomba.angle != 90

def bs2d_test():
	test_geom = wheel(1,1,5,0.3,0,Dielectric('si')\
	,priority='Occupancy')
	draw_geometry(test_geom,'test_2d')
	test_kspace = KSpace(2,x_res=50,y_res=50)
	sim = Simulation('test_2d',test_geom,test_kspace,numbands=5, \
	resolution=64)
	sim.runSimulation()
	sim.postProcess()
	draw_bandstructure('test_2d',test_kspace,5,filled=False)
	return True