void
timbre::add_tracks(
        musly_track** tracks,
        musly_trackid* trackids,
        int length,
        bool generate_ids) {
    int num_new;
    if (generate_ids) {
        idpool.generate_ids(trackids, length);
        num_new = length;
    }
    else {
        num_new = idpool.add_ids(trackids, length);
    }

    Eigen::VectorXf sim(mp.get_normtracks()->size());
    mp.append_normfacts(num_new);
    int pos = idpool.get_size() - length;
    for (int i = 0; i < length; i++) {
        similarity_raw(tracks[i], mp.get_normtracks()->data(),
                mp.get_normtracks()->size(), sim.data());

        mp.set_normfacts(pos + i, sim);
    }
}

int main (int argc, char* argv[]) {
	// Initialize options
	int i;
	get_ops(NULL);
	for (i=1; i<argc; ++i) {
		get_ops(argv[i]);
	}

	// Seed RNG
	if (ops.seedFlag) {
		srand(ops.seed);
	} else {
		srand(time(NULL));
	}

	// Read in nuclei
	if (get_nuclei(ops.nuclei_filename)==false) {
		exit(EXIT_FAILURE);
	}

	// Read in reactions
	if (get_reactions(ops.reaclib_filename)==false) {
		exit(EXIT_FAILURE);
	}

	sim();

	free(nuclei);
	free(abun);
	free(reactions);
	return 0;
}

int main(int argc, char **argv)
{
    PseudoRandomSeed();
  
    dInitODE();

//create simulator and set robot position at (0, -8)  
    Simulator sim(0, -8);    

//create motion planner
    MotionPlanner mp(&sim);    

//create graphics
    Graphics graphics(&mp);
    
//create obstacles and terrain
    CreateWorld1(&sim);
    
//if argument specified, set as maximum motion planning time
    if(argc > 1)
	graphics.SetMotionPlannerMaxTime(atof(argv[1]));    

//print help messages
    graphics.HandleEventOnHelp();

//enter event loop
    graphics.MainLoop();
    
    return 0;    
}

int main(void)
{
	Rossler<double> ross(0.398,2.0,4.0);
	RungeKutta4<double> integrator(1e-2);
//	Euler<double> integrator(1e-2);
	Simulation<double> sim(ross,integrator);

	double ti = 0.0;
	double tf = 200.0;
	long npts = 10000;

	std::ofstream datfile("out.dat", std::ios::out | std::ios::trunc);

	ross.x(0) = 0.0;
	ross.x(1) = 0.0;
	ross.x(2) = 0.0;

	if (datfile)
	{
		//sim.run(std::cout, ti, tf);
		//sim.run(datfile, npts);
		sim.run(datfile, ti, tf);
		
		datfile.close();
	}
	else
	{
		std::cerr << "Erreur à l'ouverture du fichier !" << std::endl;
	}
	
	return 0;

}

TEST(TestASMu3D, TransferGaussPtVars)
{
  SIM3D sim(1);
  sim.opt.discretization = ASM::LRSpline;
  sim.createDefaultModel();
  ASMu3D* pch = static_cast<ASMu3D*>(sim.getPatch(1));
  size_t id[3];
  const double* xi = GaussQuadrature::getCoord(3);
  RealArray oldAr(3*3*3), newAr;
  for (size_t idx = 0; idx < 3; ++idx) {
    SIM3D simNew(1);
    simNew.opt.discretization = ASM::LRSpline;
    simNew.createDefaultModel();
    ASMu3D* pchNew = static_cast<ASMu3D*>(simNew.getPatch(1));
    pchNew->uniformRefine(idx, 1);

    for (id[2] = 0; id[2] < 3; ++id[2])
      for (id[1] = 0; id[1] < 3; ++id[1])
        for (id[0] = 0; id[0] < 3; ++id[0])
          oldAr[id[0]+(id[1]+id[2]*3)*3] = (1.0 + xi[id[idx]]) / 2.0;

    pchNew->transferGaussPtVars(pch->getVolume(), oldAr, newAr, 3);
    size_t k = 0;
    for (size_t iEl = 0; iEl < 2; ++iEl)
      for (id[2] = 0; id[2] < 3; ++id[2])
        for (id[1] = 0; id[1] < 3; ++id[1])
          for (id[0] = 0; id[0] < 3; ++id[0], ++k)
            EXPECT_FLOAT_EQ(newAr[k], 0.5*iEl + 0.5*(xi[id[idx]] + 1.0) / 2.0);
  }
}

TEST(TestInitialConditions, Parse)
{
  SIM2D sim({4}, false);
  EXPECT_TRUE(sim.read("src/SIM/Test/refdata/input.xinp"));

  // Recognize both comp and component attributes and correct priority
  // Boundary conditions
  for (int i = 1; i < 5; i++)
    ASSERT_FLOAT_EQ((float)i,(*sim.getSclFunc(i))(Vec3()));
  // Initial conditions
  ASSERT_TRUE(sim.getICs().begin() != sim.getICs().end());
  const std::vector<SIMdependency::ICInfo>& bar = sim.getICs().begin()->second;
  for (std::vector<SIMdependency::ICInfo>::const_iterator info = bar.begin();
       info != bar.end(); info++)
    switch (info->component) {
      case 1:
        ASSERT_EQ(info->function, "1");
        break;
      case 2:
        ASSERT_EQ(info->function, "2");
        break;
      case 3:
        ASSERT_EQ(info->function, "3");
        break;
      case 4:
        ASSERT_EQ(info->function, "4");
        break;
      default:
        ASSERT_TRUE(false);
      }
}

TEST_P(TestASMu3D, BoundaryNodes)
{
  if (GetParam() == 0 || GetParam() > 6)
    return;

  SIM3D sim(1);
  sim.opt.discretization = ASM::LRSpline;
  ASSERT_TRUE(sim.read("src/ASM/LR/Test/refdata/boundary_nodes_3d.xinp"));
  ASSERT_TRUE(sim.createFEMmodel());

  std::stringstream str;
  str << "Face" << GetParam();
  int bcode = sim.getUniquePropertyCode(str.str(),0);
  std::vector<int> vec;
  sim.getBoundaryNodes(bcode,vec);
  ASSERT_EQ(vec.size(), 16U);
  auto it = vec.begin();
  for (int i = 0; i < 4; ++i)
    for (int j = 0; j < 4; ++j)
      if (GetParam() == 1)
        EXPECT_EQ(*it++, 1+4*(4*i+j));
      else if (GetParam() == 2)
        EXPECT_EQ(*it++, 2+4*(4*i+j));
      else if (GetParam() == 3)
        EXPECT_EQ(*it++, 16*i+j+1);
      else if (GetParam() == 4)
        EXPECT_EQ(*it++, 5+16*i+j);
      else if (GetParam() == 5)
        EXPECT_EQ(*it++, 4*i+j+1);
      else if (GetParam() == 6)
        EXPECT_EQ(*it++, 17+4*i+j);
}

TEST(TestDomainDecomposition, Setup)
{
  SIM2D sim(1);
  sim.read("src/ASM/Test/refdata/DomainDecomposition_2D_1P.xinp");
  sim.preprocess();

  DomainDecomposition dd;
  // TODO: Remove after integration
  dd.setup(sim.getProcessAdm(), sim);
  const SAM* sam = sim.getSAM();

  ASSERT_EQ(dd.getMinEq(), 1);
  ASSERT_EQ(dd.getMaxEq(), sam->getNoEquations());
  ASSERT_EQ(dd.getMinNode(), 1);
  ASSERT_EQ(dd.getMaxNode(), sam->getNoNodes());
  ASSERT_EQ(dd.getMinDOF(), 1);
  ASSERT_EQ(dd.getMaxDOF(), sam->getNoDOFs());
  for (int i = 1; i <= sam->getNoEquations(); ++i)
    ASSERT_EQ(dd.getGlobalEq(i), i);
  ASSERT_EQ(dd.getGlobalEq(sam->getNoEquations()+1), 0);

  // disabled until integration in the SIM class is added
//  ASSERT_EQ(dd.getPatchOwner(1), 0);
//  ASSERT_EQ(dd.getPatchOwner(2), -1);

  ASSERT_TRUE(dd.getMLGEQ().empty());
  ASSERT_TRUE(dd.getMLGN().empty());
}

int main(int argc,char *argv[]) {
    if(argc != 5) {
            printf("Wrong Number of Arguments\n");
            exit(0);
    }

    //double alpha1 = argv[1]-0.0;
    double alpha1, alpha2, beta1, beta2, gamma1, gamma2, phi1, phi2;
    alpha1 = alpha2 = atof( argv[1] );
    beta1 = 0.35;
    beta2 = atof( argv[2] );
    gamma1 = gamma2 = 1.0/5.0;
    phi1 = phi2 = atof( argv[4]);
    int intro_time;
    intro_time = atoi( argv[3] );

    int num_reps = 5;

    Network net = Network("gillespie toy", Network::Undirected);
    net.populate(10000);
    net.connect_all_nodes();
    Gillespie_TwoStrain_Network_Sim sim(&net, alpha1, alpha2, gamma1, gamma2, beta1, beta2, phi1, phi2, intro_time);
    for(int i =1; i <= num_reps; i++){
        cout << "Simulation number: " << i << endl;
        sim.reset();
        sim.rand_infect(5, 1);
        sim.run_simulation(10000.0);
    }
    return 0;
}

TEST(Simplifier, JumpConstFold) {
  BCMarker dummy = BCMarker::Dummy();
  IRUnit unit(test_context);
  Simplifier sim(unit);

  // Folding JmpZero and JmpNZero.
  {
    auto tester = [&] (SSATmp* val, Opcode op) {
      auto jmp = unit.gen(op, dummy, unit.defBlock(), val);
      return sim.simplify(jmp, false);
    };

    auto resultFalseZero  = tester(unit.cns(false),  JmpZero);
    auto resultFalseNZero = tester(unit.cns(false), JmpNZero);
    auto resultTrueZero   = tester(unit.cns(true),   JmpZero);
    auto resultTrueNZero  = tester(unit.cns(true),  JmpNZero);

    EXPECT_SINGLE_OP(resultFalseNZero, Nop);
    EXPECT_SINGLE_OP(resultTrueZero, Nop);

    EXPECT_SINGLE_OP(resultFalseZero, Jmp);
    EXPECT_SINGLE_OP(resultTrueNZero, Jmp);
  }

  // Folding query jumps.
  {
    auto jmpeqTaken = unit.gen(JmpEq, dummy, unit.cns(10), unit.cns(10));
    auto result = sim.simplify(jmpeqTaken, false);
    EXPECT_SINGLE_OP(result, Jmp);

    auto jmpeqNTaken = unit.gen(JmpEq, dummy, unit.cns(10), unit.cns(400));
    result = sim.simplify(jmpeqNTaken, false);
    EXPECT_SINGLE_OP(result, Nop);
  }
}

/**
 * @brief Add a new simulation instance
 * @param name unique internal name for this instance
 * @param lattice file name
 */
void flamePrepare(const char *name, const char *lattice)
{
    if(!name) return;
    try {
        Guard G(SimGlobal.lock);

        if(SimGlobal.sims.find(name)!=SimGlobal.sims.end())
            throw std::runtime_error("Sim name already in use");

        std::auto_ptr<Config> conf;
        {
            GLPSParser parser;
            conf.reset(parser.parse_file(lattice));
        }

        std::auto_ptr<Sim> sim(new Sim(name));

        sim->machine.reset(new Machine(*conf));

        SimGlobal.sims[name] = sim.get();

        sim.release();

    }catch(std::exception& e){
        fprintf(stderr, "Error: %s\n", e.what());
    }
}

void FpuStackAllocator::pop_always(LIR_Op* op, LIR_Opr opr) {
  assert(op->fpu_pop_count() == 0, "fpu_pop_count alredy set");
  assert(tos_offset(opr) == 0, "can only pop stack top");

  op->set_fpu_pop_count(1);
  sim()->pop();
}

int main() {
    
    // Construct Network
    Network net("name", Network::Undirected);
    net.populate(10000);
    net.erdos_renyi(5);

    cout << "days_ahead" << "," << "transmissibility1" << "," << "transmissibility2" << "," << "strain1" << "," << "strain2" << endl;

    for (int days_ahead = 0; days_ahead < 20; days_ahead ++){
        for (double transmissibility1 = 0.2; transmissibility1 <=  0.5; transmissibility1 += 0.05){
            for (double transmissibility2 = 0.2; transmissibility2 <= 0.5; transmissibility2 += 0.05){
                int i = 0;
                while ( i < 100){
                    // Choose and run simulation
                    Two_Strain_Percolation_Sim sim(&net);
                    sim.set_transmissibility1(transmissibility1);
                    sim.set_transmissibility2(transmissibility2);
                    sim.rand_infect(5, 0);
                    sim.run_headstart_simulation(days_ahead, 5);

                    if(sim.epidemic_size1() > 50){
                        cout << days_ahead << "," << transmissibility1 << "," << transmissibility2 << "," << sim.epidemic_size1() << "," << sim.epidemic_size2() << endl; 
                        i++;
                    }
                    
                    sim.reset();
                }
            }
        }
    }

    return 0;
}

Else::PhysicsSimulatorPtr
Else::PhysicsSimulatorFactory::create( void )
{
    PhysicsSimulatorPtr sim( new PhysicsSimulator() );
    sim->setTimeStep( myTimeStep );
    return sim;
}

void process_enable_simulator(LLMessageSystem *msg, void **user_data)
{
    // enable the appropriate circuit for this simulator and
    // add its values into the gSimulator structure
    U64		handle;
    U32		ip_u32;
    U16		port;

    msg->getU64Fast(_PREHASH_SimulatorInfo, _PREHASH_Handle, handle);
    msg->getIPAddrFast(_PREHASH_SimulatorInfo, _PREHASH_IP, ip_u32);
    msg->getIPPortFast(_PREHASH_SimulatorInfo, _PREHASH_Port, port);

    // which simulator should we modify?
    LLHost sim(ip_u32, port);

    // Viewer trusts the simulator.
    msg->enableCircuit(sim, TRUE);
    LLWorld::getInstance()->addRegion(handle, sim);

    // give the simulator a message it can use to get ip and port
    llinfos << "simulator_enable() Enabling " << sim << " with code " << msg->getOurCircuitCode() << llendl;
    msg->newMessageFast(_PREHASH_UseCircuitCode);
    msg->nextBlockFast(_PREHASH_CircuitCode);
    msg->addU32Fast(_PREHASH_Code, msg->getOurCircuitCode());
    msg->addUUIDFast(_PREHASH_SessionID, gAgent.getSessionID());
    msg->addUUIDFast(_PREHASH_ID, gAgent.getID());
    msg->sendReliable(sim);
}

int main() {
    // Construct Network
    Network net("name", false);
    net.populate(10000);

    // Let's connect the network using an arbitrary degree distribution
    // Out of 10 nodes, none will have degrees 0 or 1, 1 will have degree 2
    // 3 with degree 3, 2 with degree 4, 1 with degree 5, and 3 with degree 3
    vector<double> degree_dist;
    double tmp_array[] = {0, 0, 1, 3, 2, 1, 3};
    degree_dist.assign(tmp_array,tmp_array+7); // 5 == the length of tmp_array

    // Now make those probabilities that sum to 1
    degree_dist = normalize_dist(degree_dist, sum(degree_dist));

    // Finally, connect up the network using that user-defined degree distribution
    // Note that poisson, exponential, and scale-free random network generators are built-in
    net.rand_connect_user(degree_dist);

    // Simulation parameters
    int infectious_period = 10; // 10 days, hours,  whatever 
    double T = 0.05; // per timestep, per outbound edge probability of transmission

    for (int i = 0; i < 10; i++){
        // Choose and run simulation
        ChainBinomial_Sim sim(&net, infectious_period, T);
        sim.rand_infect(1);
        sim.run_simulation();
        cout << sim.epidemic_size() << endl;
        sim.reset();
    }
    return 0;
}

__interrupt void TIM4_Updata_IRQHandler(void)
{
    sim();//关总中断

    
    TIM4->SR1 &=~BIT(0);//清标志位
    rim();//开总中断
}

void test_sim(void){

  uint8_t instr[] = {0XAB};
  
  TEST_ASSERT_EQUAL_INT8(1, sim(instr));
  TEST_ASSERT_EQUAL_INT8(1, I0);
  TEST_ASSERT_EQUAL_INT8(1, I1);
}

int main()
{
	int m=1;
	int l=1;

	double dt=0.01;

	pendulum1d_cost::N(0,0)=0.3;
	pendulum1d_cost::N(1,1)=0.3;

	pendulum1d_cost::Q=0.1;

	pendulum1d::Cost cost=pendulum1d_cost::cost;
	pendulum1d::DCost dcost=pendulum1d_cost::dcost;

	pendulum1d::Cost costN=pendulum1d_cost::costN;
	pendulum1d::DCost dcostN=pendulum1d_cost::dcostN;

	SAC<pendulum1d>::Functions func(cost, dcost, costN, dcostN);

	SAC<pendulum1d>::Params params;

	params.alpha=-20;
	params.beta=0.6;
	
	params.Tc=0.64;
	params.Tp=2.5;
	params.Ts=0.02;
	params.Tcal=0.01;
	params.kmax=6;

	pendulum1d::System sys(m,l);

	pendulum1d::State state0(-2,7);

	std::list<pendulum1d::Params> list_pars(1000, pendulum1d::State(0, 0));
	std::list<Eigen::Matrix<double,1,1> > list_u;

	Sim<pendulum1d>::State_Sim sim(sys,dt);
	sim.init(state0);

	SAC<pendulum1d> sac(sys,dt);
	
	for(double t=0; t<30; t+=params.Ts)
	{
		sac.init(sim.get_state(), params, func, list_pars.begin(), list_pars.end(), &list_u);
		sac.sac_ctrl(list_u);

		sim.update(list_u.front());
		list_u.pop_front();

		sim.update(list_u.front());
		list_u.pop_front();
	}

	sim.save("/media/fantaosha/Documents/JHU/Summer 2015/results/result3.mat");
	return 0;
}

int
main(int argc, char* argv[])
{
	namespace po = boost::program_options;

	try {

		po::options_description desc = commonOptions();
		desc.add_options()
			("neurons,n", po::value<unsigned>()->default_value(1000), "number of neurons")
			("synapses,m", po::value<unsigned>()->default_value(1000), "number of synapses per neuron")
			("dmax,d", po::value<unsigned>()->default_value(1), "maximum excitatory delay,  where delays are uniform in range [1, dmax]")
		;

		po::variables_map vm = processOptions(argc, argv, desc);

		unsigned ncount = vm["neurons"].as<unsigned>();
		unsigned scount = vm["synapses"].as<unsigned>();
		unsigned dmax = vm["dmax"].as<unsigned>();
		unsigned duration = vm["duration"].as<unsigned>();
		unsigned stdp = vm["stdp-period"].as<unsigned>();
		unsigned verbose = vm["verbose"].as<unsigned>();
		bool runBenchmark = vm.count("benchmark") != 0;

		std::ofstream file;
		std::string filename;

		if(vm.count("output-file")) {
			filename = vm["output-file"].as<std::string>();
			file.open(filename.c_str()); // closes on destructor
		}

		std::ostream& out = filename.empty() ? std::cout : file;

		LOG(verbose, "Constructing network");
		boost::scoped_ptr<nemo::Network> net(nemo::random::construct(ncount, scount, dmax, stdp != 0));
		LOG(verbose, "Creating configuration");
		nemo::Configuration conf = configuration(vm);
		LOG(verbose, "Simulation will run on %s", conf.backendDescription());
		LOG(verbose, "Creating simulation");
		boost::scoped_ptr<nemo::Simulation> sim(nemo::simulation(*net, conf));
		LOG(verbose, "Running simulation");
		if(runBenchmark) {
			benchmark(sim.get(), ncount, scount, vm);
		} else {
			simulate(sim.get(), duration, stdp, out);
		}
		LOG(verbose, "Simulation complete");
		return 0;
	} catch(std::exception& e) {
		std::cerr << e.what() << std::endl;
		return -1;
	} catch(...) {
		std::cerr << "random: An unknown error occurred\n";
		return -1;
	}

}