Constraint::Constraint()
: Value(0.0),
  Type(None),
  AlignmentType(Undef),
  Name(""),
  First(GeoUndef),
  FirstPos(none),
  Second(GeoUndef),
  SecondPos(none),
  Third(GeoUndef),
  ThirdPos(none),
  LabelDistance(10.f),
  LabelPosition(0.f),
  isDriving(true),
  InternalAlignmentIndex(-1),
  isInVirtualSpace(false)
{
    // Initialize a random number generator, to avoid Valgrind false positives.
    static boost::mt19937 ran;
    static bool seeded = false;

    if (!seeded) {
        ran.seed(QDateTime::currentMSecsSinceEpoch() & 0xffffffff);
        seeded = true;
    }
    static boost::uuids::basic_random_generator<boost::mt19937> gen(&ran);

    tag = gen();
}

void UPGMpp::getRandomAssignation(CGraph &graph,
                                  map<size_t,size_t> &assignation,
                                  TInferenceOptions &options )
{
    static boost::mt19937 rng1;
    static bool firstExecution = true;

    if ( firstExecution )
    {
        rng1.seed(std::time(0));
        firstExecution = false;
    }

    vector<CNodePtr> &nodes = graph.getNodes();

    for ( size_t node = 0; node < nodes.size(); node++ )
    {
        // TODO: Check if a node has a fixed value and consider it

        size_t N_classes = nodes[node]->getType()->getNumberOfClasses();

        boost::uniform_int<> generator(0,N_classes-1);
        int state = generator(rng1);

        assignation[nodes[node]->getID()] = state;
    }

    /*map<size_t,size_t>::iterator it;
    for ( it = assignation.begin(); it != assignation.end(); it++ )
        cout << "[" << it->first << "] " << it->second << endl;*/
}

	void readSeed(boost::program_options::variables_map& variableMap, boost::mt19937& randomSource)
	{
		if(variableMap.count("seed") > 0)
		{
			randomSource.seed(variableMap["seed"].as<int>());
		}
	}

float normalvariate(float mean, float sigma)
{
    static bool seeded = false;

    if (!seeded)
    {
        // seed generator with #seconds since 1970
        rng.seed(static_cast<unsigned> (std::time(0)));
        seeded = true;
    }

    // select desired probability distribution
    boost::normal_distribution<float> norm_dist(mean, sigma);

    // bind random number generator to distribution, forming a function
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<float> >  normal_sampler(rng, norm_dist);

    // sample from the distribution
    return normal_sampler();
}

int
main(int argc, char *argv[])
{
    // Parse and read the command line arguments.
    if (argc != 5)
    {
        cerr << "Usage: " << PNAME
             << " <min-weight> <max-weight> <#vertices> <#tests>\n";
        exit(1);
    }

    long int min_weight = strtol(argv[1], 0, 0);
    long int max_weight = strtol(argv[2], 0, 0);
    long int n_vertices = strtol(argv[3], 0, 0);
    long int n_tests = strtol(argv[4], 0, 0);

    if (errno)
    {
        cerr << "Bad arguments\n";
        exit(1);
    }

    // Check consistency of command line arguments.
    if (n_tests < 0)
    {
        cerr << "the number of tests must be a positive integer\n";
        exit(1);
    }
    
    if (n_vertices > numeric_limits<unsigned>::digits)
    {
        cerr << "Too many vertices\n"
             << "Maximum n.o. vertices is "
             << numeric_limits<unsigned>::digits << "\n";
        exit(1);
    }

    if (n_vertices < 2)
    {
        cerr << "Too few vertices\n"
             << "must be at least two\n";
        exit(1);
    }

    if (min_weight > max_weight)
    {
        cerr << "min weight must be greater than max weight\n";
        exit(1);
    }

    if (abs(max_weight) >= numeric_limits<int>::max() / n_vertices ||
        abs(min_weight) >= numeric_limits<int>::max() / n_vertices)
    {
        cerr << "weights too large\n";
        exit(1);
    }

    // Initialize the random number generator.
    timeval tv;
    gettimeofday(&tv, 0);
    g_generator.seed(tv.tv_sec * 1000000 + tv.tv_usec);
    boost::uniform_int<int> uniform_int_dist(min_weight, max_weight);
    boost::variate_generator<mt19937, uniform_int<int> > rng_i(g_generator, uniform_int_dist);
    
    // Perform the tests.
    enum {MAX_SPAN_0, MAX_SPAN_1, MAX_SPAN_2,
          MAX_NOSPAN_0, MAX_NOSPAN_1, MAX_NOSPAN_2,
          MIN_SPAN_0, MIN_SPAN_1, MIN_SPAN_2,
          MIN_NOSPAN_0, MIN_NOSPAN_1, MIN_NOSPAN_2,
          N_CASES};
    string case_names[] = {
        "MAX_SPAN_0", "MAX_SPAN_1", "MAX_SPAN_2",
        "MAX_NOSPAN_0", "MAX_NOSPAN_1", "MAX_NOSPAN_2",
        "MIN_SPAN_0", "MIN_SPAN_1", "MIN_SPAN_2",
        "MIN_NOSPAN_0", "MIN_NOSPAN_1", "MIN_NOSPAN_2"
    };

    complete_graph      g(n_vertices);
    multi_array<int, 2> weights(extents[n_vertices][n_vertices]);
    vector<Vertex>      parent(n_vertices);
    Vertex              roots[] = {0, 1};

    vector<Edge>        branching;
    vector< vector< Edge > > all_branchings(N_CASES);
    vector<int>         ans(N_CASES);

    while (n_tests--)
    {
        // Clear all the branchings from previous iteration.
        BOOST_FOREACH (vector<Edge> &vec, all_branchings)
        {
            vec.clear();
        }

        // create new weights
        int *end = weights.origin() + n_vertices * n_vertices;
        for (int *ip = weights.origin(); ip != end; ++ip)
        {
            *ip = rng_i();
        }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
//, 0,  0,  0,  1,  1

///*--------------------------------- Graph 23 ----------------------------------*/

//    Eigen::MatrixXi g23_adj(11,11);
//, 0,  1,  1
//, 1,  0,  1
//, 1,  1,  0

///*--------------------------------- Graph 24 ----------------------------------*/

//    Eigen::MatrixXi g24_adj(11,11);

//    g24_adj <<  0,  1,  1,  0,  1,  1,  0,  0,  0,
//                1,  0,  1,  1,  1,  0,  0,  0,  0,
//                1,  1,  0,  1,  1,  0,  1,  0,  0,
//                0,  1,  1,  0,  1,  1,  1,  1,  0,
//                1,  1,  1,  1,  0,  0,  1,  0,  0,
//                1,  0,  0,  1,  0,  0,  0,  0,  0,
//                0,  0,  1,  1,  1,  0,  0,  0,  0,
//                0,  0,  0,  1,  0,  0,  0,  0,  1,
//                0,  0,  0,  0,  0,  0,  0,  1,  0;

///*--------------------------------- Graph 25 ----------------------------------*/

//    Eigen::MatrixXi g25_adj(11,11);

//    g25_adj <<  0,  1,  1,  0,  0,
//                1,  0,  1,  1,  1,
//                1,  1,  0,  0,  0,
//                0,  1,  0,  0,  1,
//                0,  1,  0,  1,  0;

    static boost::mt19937 rng;
    rng.seed(std::time(0));

    size_t N_graphsForTraining = 9;

    size_t absoluteSuccess_Greedy  = 0;
    size_t absoluteSuccess_ICM     = 0;
    size_t absoluteSuccess_AlphaExpansion = 0;
    size_t absoluteSuccess_AlphaBetaSwap  = 0;
    size_t absoluteSuccess_LBP = 0;
    size_t absoluteSuccess_TRP = 0;
    size_t absoluteSuccess_RBP = 0;
    size_t absoluteNumberOfNodes = 0;

    size_t N_repetitions = 100;

    for ( size_t rep = 0; rep < N_repetitions; rep++ )
    {

    vector<size_t> graphs_to_train;
    size_t N_graphsAdded       = 0;

    while ( N_graphsAdded < N_graphsForTraining )
    {
        boost::uniform_int<> int_generator(0,11);
        int rand = int_generator(rng);

        if ( std::find(graphs_to_train.begin(),graphs_to_train.end(),rand)
             == graphs_to_train.end() )
        {
            graphs_to_train.push_back(rand);
            N_graphsAdded++;
        }

int main(int argc, char** argv)
{
	int args = 1;
#ifdef USE_OPENCL
	if (argc < 10) {
#else
	if (argc < 7) {
#endif
		std::cout << "Not enough arguments" << std::endl;
		system("pause");
		return 1;
	}

	DIM = util::toInt(argv[args++]);
	N = util::toInt(argv[args++]);
	K = util::toInt(argv[args++]);
	ITERATIONS = util::toInt(argv[args++]);
	RUNS = util::toInt(argv[args++]);
#ifdef USE_OPENCL
	AM_LWS = util::toInt(argv[args++]);
	RP_LWS = util::toInt(argv[args++]);
	CT_LWS = util::toInt(argv[args++]);

	USE_ALL_DEVICES = util::toInt(argv[args++]);
#else
	device_count = util::toInt(argv[args++]);
#endif


	std::cout << "DIM = " << DIM << std::endl;
	std::cout << "N = " << N << std::endl;
	std::cout << "K = " << K << std::endl;
	std::cout << "ITERATIONS = " << ITERATIONS << std::endl;
	std::cout << "RUNS = " << RUNS << std::endl;
#ifdef USE_OPENCL
	std::cout << "AM_LWS = " << AM_LWS << std::endl;
	std::cout << "RP_LWS = " << RP_LWS << std::endl;
	std::cout << "CT_LWS = " << CT_LWS << std::endl;
	std::cout << "USE_ALL_DEVICES = " << USE_ALL_DEVICES << std::endl << std::endl;
#else
	std::cout << "device_count = " << device_count << std::endl << std::endl;
#endif


#ifdef _WIN32
	rng.seed();
	srand(GetTickCount());
#else
	rng.seed();
	srand(getTimeMs());
#endif

	u = boost::uniform_real<float>(0.0f, 1000000.0f);
	gen = new boost::variate_generator<boost::mt19937&, boost::uniform_real<float> >(rng, u);

#ifdef USE_OPENCL
	cl_int clError = CL_SUCCESS;
	initCL();

	for (int i = 0; i < clDevices.size(); ++i) {

		clInputBuf.push_back(cl::Buffer(clContext, CL_MEM_READ_ONLY, N * DIM * sizeof(float), NULL, &clError));
		if (clError != CL_SUCCESS) std::cout << "OpenCL Error: Could not create buffer" << std::endl;

		clCentroidBuf.push_back(cl::Buffer(clContext, CL_MEM_READ_WRITE, K * DIM * sizeof(float), NULL, &clError));
		if (clError != CL_SUCCESS) std::cout << "OpenCL Error: Could not create buffer" << std::endl;

		clMappingBuf.push_back(cl::Buffer(clContext, CL_MEM_READ_WRITE, N * sizeof(int), NULL, &clError));
		if (clError != CL_SUCCESS) std::cout << "OpenCL Error: Could not create buffer" << std::endl;

		clReductionBuf.push_back(cl::Buffer(clContext, CL_MEM_WRITE_ONLY, N * sizeof(float), NULL, &clError));
		if (clError != CL_SUCCESS) std::cout << "OpenCL Error: Could not create buffer" << std::endl;

		clClusterAssignment[i].setArgs(clInputBuf[i](), clCentroidBuf[i](), clMappingBuf[i]());
		clClusterReposition[i].setArgs(clInputBuf[i](), clMappingBuf[i](), clCentroidBuf[i]());
		clClusterReposition_k[i].setArgs(clInputBuf[i](), clMappingBuf[i](), clCentroidBuf[i]());
		//clClusterReposition_k_c[i].setArgs(clInputBuf[i](), clMappingBuf[i](), clCentroidBuf[i](), clConvergedBuf[i]());
		clComputeCost[i].setArgs(clInputBuf[i](), clCentroidBuf[i](), clMappingBuf[i](), clReductionBuf[i]());

	}

	device_count = clDevices.size();
#endif

	util::Clock clock;
	clock.reset();

	for (int i = 0; i < RUNS; ++i) {
		mapping_list.push_back(NULL);
		centroids_list.push_back(NULL);
		cost_list.push_back(0.0f);
	}

	float* source = new float[N*DIM];
	for (int i = 0; i < N*DIM; ++i)
		source[i] = gen_random_float();

	input_list.push_back(source);

	for (int i = 1; i < device_count; ++i) {
//.........这里部分代码省略.........

void UPGMpp::getSpanningTree( CGraph &graph, std::vector<size_t> &tree)
{
    // TODO: The efficiency of this method can be improved

    // Reset tree
    tree.clear();

    static boost::mt19937 rng;
    static bool firstExecution = true;

    if ( firstExecution )
    {
        rng.seed(time(0));
        firstExecution = false;
    }

    vector<CNodePtr> &v_nodes = graph.getNodes();
    vector<CNodePtr> v_nodesToExplore;
    map<size_t,size_t> nodesToExploreMap;
    size_t N_nodes = v_nodes.size();

    //cout << "Random: ";

    for ( size_t i_node = 0; i_node < N_nodes; i_node++ )
    {        
        boost::uniform_real<> real_generator(0,1);
        real_generator.reset();
        double rand = real_generator(rng);

        //cout << rand << " ";

        if ( rand > 0.25 )
        {
            v_nodesToExplore.push_back( v_nodes[i_node] );
            nodesToExploreMap[ v_nodes[i_node]->getID() ] =  v_nodesToExplore.size()-1;
        }

    }

    //cout << endl;

    //SHOW_VECTOR_NODES_ID("Nodes to explore: ", v_nodesToExplore)

    size_t N_nodesToExplore = v_nodesToExplore.size();

    if ( !N_nodesToExplore )
        return;

    bool nodeAdded = true;
    vector<bool> v_nodesAdded( N_nodesToExplore, false );

    //
    // Randomly select the root node
    //

    boost::uniform_int<> generator(0,N_nodesToExplore-1);
    int rand = generator(rng);

    //cout << "Root:" << rand << endl;

    v_nodesAdded[rand] = true;
    multimap<size_t, CEdgePtr> &edgesF = graph.getEdgesF();

    while ( nodeAdded )
    {
        nodeAdded = false;

        for ( size_t i_node = 0; i_node < N_nodesToExplore; i_node++ )
        {
            //cout << "Node " << i_node << " ID "<< v_nodesToExplore[i_node]->getID() << " Added? " << v_nodesAdded[i_node] << endl;

            // Check that the node has not been added to the tree yet
            if ( !v_nodesAdded[i_node] )
            {
                size_t nodeID = v_nodesToExplore[i_node]->getID();

                NEIGHBORS_IT neighbors = edgesF.equal_range(nodeID);

                for ( multimap<size_t,CEdgePtr>::iterator itNeigbhor = neighbors.first;
                      itNeigbhor != neighbors.second;
                      itNeigbhor++ )
                {

                    CEdgePtr edgePtr = (*itNeigbhor).second;
                    size_t neighborID;

                    if ( !edgePtr->getNodePosition( nodeID ) )
                        neighborID = edgePtr->getSecondNodeID();
                    else
                        neighborID = edgePtr->getFirstNodeID();

                    //cout << "Current node ID: " << nodeID << " neighbor: " << neighborID << endl;

                    if ( v_nodesAdded[nodesToExploreMap[neighborID]] )
                    {
                        v_nodesAdded[i_node] = true;
                        nodeAdded = true;
                    }
                }

//.........这里部分代码省略.........

std::vector<T1> generateRandomSet(const unsigned int size_,
								  const T1 min_,
								  const T1 max_,
								  const bool allowRepetition_)
{
	generator.seed(std::rand());

#if BOOST_MINOR_VERSION <= 46
	NumberType range(min_, max_);
	boost::variate_generator<Generator&, NumberType> dist(generator, range);
#else
	Generator dist(min_, max_);
#endif

	std::vector<T1> numbers;
	numbers.reserve(size_);

	std::map<T1, bool> used;
	while (numbers.size() < size_)
	{
#if BOOST_MINOR_VERSION <= 46
		T1 number = dist();
#else
		T1 number = dist(generator);
#endif

		if (allowRepetition_ || used.find(number) == used.end())
		{
			used[number] = true;
			numbers.push_back(number);
		}
	}

	return numbers;
}

	void setUp() {
		randomGen.seed(time(NULL));
		eventLoop = new DummyEventLoop();
		timerFactory = new DummyTimerFactory();
		connection = boost::make_shared<MockeryConnection>(failingPorts, true, eventLoop);
		//connection->onDataSent.connect(boost::bind(&SOCKS5BytestreamServerSessionTest::handleDataWritten, this, _1));
		//stream1 = boost::make_shared<ByteArrayReadBytestream>(createByteArray("abcdefg")));
//		connection->onDataRead.connect(boost::bind(&SOCKS5BytestreamClientSessionTest::handleDataRead, this, _1));
	}

	void setUp() {
		crypto = boost::shared_ptr<CryptoProvider>(PlatformCryptoProvider::create());
		destination = "092a44d859d19c9eed676b551ee80025903351c2";
		randomGen.seed(static_cast<unsigned int>(time(NULL)));
		eventLoop = new DummyEventLoop();
		timerFactory = new DummyTimerFactory();
		connection = boost::make_shared<MockeryConnection>(failingPorts, true, eventLoop);
		//connection->onDataSent.connect(boost::bind(&SOCKS5BytestreamServerSessionTest::handleDataWritten, this, _1));
		//stream1 = boost::make_shared<ByteArrayReadBytestream>(createByteArray("abcdefg")));
//		connection->onDataRead.connect(boost::bind(&SOCKS5BytestreamClientSessionTest::handleDataRead, this, _1));
	}

  void init ()
  {
    std::chrono::high_resolution_clock::time_point now =
      std::chrono::high_resolution_clock::now();

    std::chrono::nanoseconds time =
      std::chrono::duration_cast<std::chrono::nanoseconds>
      (now.time_since_epoch () );

    ran.seed (time.count() );

    pid = getpid();
  }

void parse_args(int argc, char* argv[]) {
  uint32_t seed = 0;
  bool has_seed = false;

  struct option long_opts[] = {
    { "help", false, nullptr, 'h' },
    { "seed", true, nullptr, 's' },
    { nullptr, 0, nullptr, 0 }
  };

  while (true) {
    optopt = 1;
    int optchar = getopt_long(argc, argv, "hs:", long_opts, nullptr);
    if (optchar == -1) {
      break;
    }

    switch (optchar) {
      case 's': {
        char *endptr;
        seed = strtol(optarg, &endptr, 0);
        if (endptr == optarg || *endptr != '\0') {
          fprintf(stderr, "invalid seed value \"%s\": must be a positive "
                  "integer\n", optarg);
          exit(1);
        }
        has_seed = true;
        break;
      }
      case 'h':
        print_usage(stdout, argv[0]);
        exit(0);
      case '?':
        exit(1);
      default:
        // Only happens if someone adds another option to the optarg string,
        // but doesn't update the switch statement to handle it.
        fprintf(stderr, "unknown option \"-%c\"\n", optchar);
        exit(1);
    }
  }

  if (!has_seed) {
    seed = time(nullptr);
  }

  printf("seed: %" PRIu32 "\n", seed);
  rng.seed(seed);
}

boost::unit_test::test_suite* init_unit_test_suite(int argc, char* argv[]) {
  THRIFT_UNUSED_VARIABLE(argc);
  THRIFT_UNUSED_VARIABLE(argv);
  uint32_t seed = static_cast<uint32_t>(time(NULL));
  printf("seed: %" PRIu32 "\n", seed);
  rng.seed(seed);

  boost::unit_test::test_suite* suite = &boost::unit_test::framework::master_test_suite();
  suite->p_name.value = "ZlibTest";

  uint32_t buf_len = 1024 * 32;
  add_tests(suite, gen_uniform_buffer(buf_len, 'a'), buf_len, "uniform");
  add_tests(suite, gen_compressible_buffer(buf_len), buf_len, "compressible");
  add_tests(suite, gen_random_buffer(buf_len), buf_len, "random");

  suite->add(BOOST_TEST_CASE(test_no_write));

  return NULL;
}

bool init_unit_test_suite() {
  uint32_t seed = static_cast<uint32_t>(time(NULL));
#ifdef HAVE_INTTYPES_H
  printf("seed: %" PRIu32 "\n", seed);
#endif
  rng.seed(seed);

  boost::unit_test::test_suite* suite = &boost::unit_test::framework::master_test_suite();
  suite->p_name.value = "ZlibTest";

  uint32_t buf_len = 1024 * 32;
  add_tests(suite, gen_uniform_buffer(buf_len, 'a'), buf_len, "uniform");
  add_tests(suite, gen_compressible_buffer(buf_len), buf_len, "compressible");
  add_tests(suite, gen_random_buffer(buf_len), buf_len, "random");

  suite->add(BOOST_TEST_CASE(test_no_write));

  return true;
}

int main (int argc, char const *argv[]) {
  // The default maze size is 20x10.  A different size may be specified on
  // the command line.
  std::size_t x = 20;
  std::size_t y = 10;

  if (argc == 3) {
    x = boost::lexical_cast<std::size_t>(argv[1]);
    y = boost::lexical_cast<std::size_t>(argv[2]);
  }

  random_generator.seed(std::time(0));
  maze m(x, y);
  random_maze(m);
  if (m.solve())
    std::cout << "Solved the maze." << std::endl;
  else
    std::cout << "The maze is not solvable." << std::endl;
  std::cout << m << std::endl;
  return 0;
}

int main() {
    rnd.seed(42);
    vertex_iterator vi,vi_end;
    Graph m_graph;

    auto start = chrono::high_resolution_clock::now();
    generate_random_graph(m_graph, 1000/*00*/, 500/*00*/, rnd); // reduced load for Coliru

    std::cout << "Generated " << num_vertices(m_graph) << " vertices and " << num_edges(m_graph) << " edges in "
        << chrono::duration_cast<chrono::milliseconds>(chrono::high_resolution_clock::now() - start).count() << "ms\n"
        << "The graph has a cycle? " << std::boolalpha << has_cycle(m_graph) << "\n"
        << "starting selective removal...\n";

    start = chrono::high_resolution_clock::now();
    size_t count = 0;
    for (boost::tie(vi, vi_end) = boost::vertices(m_graph); vi!=vi_end;)
    {
        if (m_graph[*vi].guid.part1 == 0) {
            count++;
            clear_vertex(*vi, m_graph);
            //std::cout << "." << std::flush;
#if defined(STABLE_IT)
            auto toremove = vi++;
            boost::remove_vertex(*toremove,m_graph);
#else
            boost::remove_vertex(*vi,m_graph);
            boost::tie(vi, vi_end) = boost::vertices(m_graph);
#endif
        } else 
            ++vi;
    }

    std::cout << "Done in " << chrono::duration_cast<chrono::milliseconds>(
            chrono::high_resolution_clock::now() - start).count() << "ms\n";

    std::cout << "After: " << num_vertices(m_graph) << " vertices and " << num_edges(m_graph) << " edges\n";
}

int main(int argc, char ** argv)
{
	using namespace std;
	using namespace boost;

	if(argc <= 3) {
		cout << "usage: " << argv[0] << " <grid size> <time> <runs>" << endl;
		return -1;
	}

	int system_random = open("/dev/random", O_RDONLY);
	uint32_t seed;
	read(system_random, &seed, 4);
	close(system_random);

	rng.seed(seed);

	for(int count = 0; count < atoi(argv[3]);) {
		grid_lattice grid(atoi(argv[1]));
		while(true) {
			double dt = grid.next_event();
			grid.time += dt;
			if(grid.time > atof(argv[2])) break;
			grid.flip();
			if(grid.empty()) grid.restart();
		}

		if(!grid.empty()) {
			cout << grid.size() << endl;
//			cout << grid << endl;
			count++;
		}
	}

	return 0;
}

void RandBipartiteGraph::GenerateGraph()
{
    // make n left-vertices
    for(unsigned int i = 0; i < n; i++){
        AddVariable(i);
    }

    vector<Vertex*> availableVertexList;
    availableVertexList.reserve(n);
    gen.seed(time(NULL));

    if (c%d != 0){
    // make c/d n right-vertices
        for(unsigned int i = 0; i < GetRightVerticesNumber(); i++){
            availableVertexList.push_back(static_cast<Vertex*>(AddConstraint(i)));
        }
    }

    for (unsigned int i = 0; i < n; i++){
        for(unsigned int j = 0; j < c; j++){
        	AddEdge(static_cast<Vertex*>(variables[i]), GetRandRightVertex(availableVertexList));
        }
    }
}

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
//int main(int ac, char* av[]){

  // global objects for the program:
  /*Cogaps_options_class Cogaps_options(ac, av);
  try {
    // Cogaps_options_class Cogaps_options(ac, av); // WS -- change from Sasha's
    if (Cogaps_options.to_help){
      Cogaps_options.help(cout);
      return 0;
    }
    //cout << Cogaps_options; 
  } 
  catch( const exception & e) {
    cerr <<e.what()<<endl;
    return 1;
  }
*/
  // ===========================================================================
  // Initialization of the random number generator.
  // Different seeding methods:
  // --- fixed seed 
  //std::vector<unsigned long> ve(2);
  //ve[0]=198782;ve[1]=89082;
  //boost::random::seed_seq seq(ve);
  //rng.seed(seq);
  // --- seeded with time
  rng.seed(static_cast<boost::uint32_t>(std::time(0)));
  //---------------------

  // ===========================================================================
  // Part 1) Initialization: 
  // In this section, we read in the system parameters from the paremter file
  // parameter.txt, and matrices D and S from datafile.txt. 
  // Then we initialize A and P in both their atomic domains and 
  // matrix forms.
  // ===========================================================================
	
	//Matlab pointers for inputs
    mxArray* DArrayElems;
	mxArray* SArrayElems;
	mxArray* InputColsPtr;
	mxArray* InputRowsPtr;

	//Matlab pointers for config 
	mxArray* nFactorPtr;
	mxArray* simulationIDPtr;
	mxArray* nEquilPtr;
	mxArray* nSamplePtr;
	mxArray* QAtomicPtr;
	mxArray* alphaAPtr;
	mxArray* nIterAPtr;
	mxArray* nMaxAPtr;
	mxArray* maxGibbsAPtr;
	mxArray* lambdaAScalePtr;
	mxArray* alphaPPtr;
	mxArray* nIterPPtr;
	mxArray* nMaxPPtr;
	mxArray* maxGibbsPPtr;
	mxArray* lambdaPScalePtr;
	
	//Transferring inputs from mex gateway to pointers
	DArrayElems = mxGetCell(prhs[0], 0);
	SArrayElems = mxGetCell(prhs[0], 1);
	InputColsPtr = mxGetCell(prhs[0], 2);
	InputRowsPtr = mxGetCell(prhs[0], 3);
	nFactorPtr = mxGetCell(prhs[0], 4);
	simulationIDPtr = mxGetCell(prhs[0], 5);
	nEquilPtr = mxGetCell(prhs[0], 6);
	nSamplePtr = mxGetCell(prhs[0], 7);
	QAtomicPtr = mxGetCell(prhs[0], 8);
	alphaAPtr = mxGetCell(prhs[0], 9);
	nIterAPtr = mxGetCell(prhs[0], 10);
	nMaxAPtr = mxGetCell(prhs[0], 11);
	maxGibbsAPtr = mxGetCell(prhs[0], 12);
	lambdaAScalePtr = mxGetCell(prhs[0], 13);
	alphaPPtr = mxGetCell(prhs[0], 14);
	nIterPPtr = mxGetCell(prhs[0], 15);
	nMaxPPtr = mxGetCell(prhs[0], 16);
	maxGibbsPPtr = mxGetCell(prhs[0], 17);
	lambdaPScalePtr = mxGetCell(prhs[0], 18);
	
	
	//creating native C style data types
	double* DArray;
	double* SArray;
	double* numInputRows;
	double* numInputCols;
	
	//Establish the C types for the Config Information
	double* nFactorC = mxGetPr(nFactorPtr);
	
	int ptrLength; 
	char* simulationIDC;
	ptrLength = mxGetNumberOfElements(simulationIDPtr) + 1;
	simulationIDC = (char *)(mxCalloc(ptrLength, sizeof(char)));
	mxGetString(simulationIDPtr, simulationIDC, ptrLength);
	
//.........这里部分代码省略.........