C++ Experiment::microcircuit方法代码示例

int main(int argc, const char* argv[])
{
    const URI blue_config_filename(argv[1]);
    try {
    	// LOADING CURCUIT
        Experiment experiment;
        experiment.open(blue_config_filename);
        Microcircuit & microcircuit = experiment.microcircuit();
        const Targets & targets = experiment.targets();
        const Cell_Target target = targets.cell_target("Column");
        microcircuit.load(target, NEURONS | MORPHOLOGIES);

        // PRELOAD the Trees amd Neuron Morphology Mapping
        ISpatialIndex *trees[MORPHOLOGIES_COUNT];
        ISpatialIndex *neurons[NEURONS_COUNT];
        global_transformer *transforms[NEURONS_COUNT];

        int cm=0;
        int cn=0;
        string baseName = "";
        Morphologies & myMorphologies = microcircuit.morphologies();
        Neurons & myNeurons = microcircuit.neurons();

        cout << "PreLoading Mappings \n";
        Morphologies::iterator myMorphologiesEnd = myMorphologies.end();
          for (Morphologies::iterator m = myMorphologies.begin(); m != myMorphologiesEnd; ++m)
          {
        	  baseName = m->label();
              IStorageManager* diskfile = StorageManager::loadDiskStorageManager(baseName);
              trees[cm] = RTree::loadRTree(*diskfile, 1);

              std::cout << "Checking R-tree structure... ";
              if (!trees[cm]->isIndexValid()) std::cerr << "R-tree internal checks failed!\n"; else std::cout << "OK\n";
              IStatistics * tree_stats;
              trees[cm]->getStatistics (&tree_stats);
              cout << *tree_stats;

              Neurons::iterator myNeuronsEnd = myNeurons.end();
              for (Neurons::iterator n = myNeurons.begin(); n != myNeuronsEnd; ++n)
                  {
            	  if (strcmp(n->morphology().label().c_str(),m->label().c_str())==0)
            		  {
            		  transforms[cn] = n->global_transform().inverse();
            		  neurons[cn] = trees[cm];
            		  }
            	  cn++;
            	  if (cn>=NEURONS_COUNT) break;
                  }
              cn=0;cm++;
          }
          /*
          // SINGLE QUERY
          cout << "Start Querying \n";
          range_query_visitor visitor;
		  SpatialIndex::Region query;

		  micron_vector plow,phigh;
		  plow[0] = 0; plow[1] = 0; plow[2] = 0;
		  phigh[0]=90; phigh[1]=90; phigh[2]=90;

          for (int i=0;i<NEURONS_COUNT;i++)
          {
    		  get_transformed_cube_mbr(plow,phigh,*transforms[i],&query);
    		  neurons[i]->intersectsWithQuery(query,visitor);
    		  visitor.inc_neuron();
          }
          visitor.print_stats();
*/

          // PERFORMANCE EVALUATION RANDOM RANGE QUERIES
          {
          cout << "Start Range Query Analysis \n";
          range_query_visitor visitor;
  		  double plow[3],phigh[3];
  		  plow[0]  = 98.2538; plow[1]  = 1005.14; plow[2]  = 92.9046;
  		  phigh[0] = 452.301; phigh[1] = 1385.56; phigh[2] = 456.743;
  	 	  //plow[0]  = -1698.38; plow[1] = -1065.03; plow[2] = -1724.75;
  	 	 // phigh[0] = 2248.56;  phigh[1]= 1894.28;  phigh[2]= 2276.71;
          SpatialIndex::Region query_region = SpatialIndex::Region(plow,phigh,3);
  		  boost::mt11213b generator (42u);
  		  const double x[3] = {0.0, 0.0, 0.0};
  		  SpatialIndex::Point rnd_point1 (x, query_region.m_dimension);
  		  SpatialIndex::Point rnd_point2 (x, query_region.m_dimension);
  		  SpatialIndex::Region query;

  		  for (int j=0;j<QUERIES_FOR_ANALYSIS;j++)
  		  {
  		     for (size_t i = 0; i < query_region.m_dimension; i++)
  	 	     {
  	 	          boost::uniform_real<> uni_dist (query_region.m_pLow[i],query_region.m_pHigh[i]);
  	 	          boost::variate_generator<boost::mt11213b &,boost::uniform_real<> > uni(generator, uni_dist);
  	 	          rnd_point1.m_pCoords[i] = uni();

  	 	          boost::uniform_real<> uni_dist1 (rnd_point1.m_pCoords[i],query_region.m_pHigh[i]);
  	 	          boost::variate_generator<boost::mt11213b &,boost::uniform_real<> > uni1(generator, uni_dist1);
  	 	          rnd_point2.m_pCoords[i] = uni1();
  	 	     }
  			  visitor.new_query();
  			  visitor.reset_neuron();
  			  for (int i=0;i<NEURONS_COUNT;i++)
//.........这里部分代码省略.........

int main(int argc, const char* argv[])
{
	/*
	 *  Build 48 Index with Links
	 */


	// Load Circuit
    Experiment experiment;
    experiment.open(blue_config_filename);
    Microcircuit & microcircuit = experiment.microcircuit();
    const Targets & targets = experiment.targets();
    const Cell_Target target = targets.cell_target("Column");
    microcircuit.load(target, NEURONS | MORPHOLOGIES);

    //Make Neuron Rtrees
    ISpatialIndex *neuronTrees[MORPHOLOGIES_COUNT];
    string *morphologyLabels[MORPHOLOGIES_COUNT];

    int cm=0;
    Morphologies & myMorphologies = microcircuit.morphologies();
    Morphologies::iterator myMorphologiesEnd = myMorphologies.end();
      for (Morphologies::iterator i = myMorphologies.begin(); i != myMorphologiesEnd; ++i)
      {
    	  morphologyLabels[cm] = i->label();
    	  neuronTrees[cm] = RTree::createNewRTree (createNewMemoryStorageManager(), 0.7, 127, 127, 3,RTree::RV_RSTAR,indexIdentifier);
    	  cm++;
      }

    Neurons & myNeurons = microcircuit.neurons();
    Neurons::iterator myNeuronsEnd = myNeurons.end();
    for (Neurons::iterator i = myNeurons.begin(); i != myNeuronsEnd; ++i)
    {
    	cm=0;
    	for (cm=0;cm<MORPHOLOGIES_COUNT;cm++)
    		if (strcmp(i->morphology().label(),morphologyLabels[cm])==0) break;

    	Transform_3D<Micron> trafo = i->global_transform();
    	Sections mySections = i->morphology().all_sections();
    	Sections::iterator mySectionsEnd = mySections.end();
    	for (Sections::iterator s = mySections.begin(); s != mySectionsEnd; ++s)
    	    {
    		 Segments segments = s->segments();
    		 Segments::const_iterator segments_end = segments.end();
    		 for (Segments::const_iterator j = segments.begin(); j != segments_end ; ++j)
    		     {
     			 vect plow, phigh;
     			 get_segment_mbr (*j, trafo, &plow, &phigh);
     			 SpatialIndex::Region mbr = SpatialIndex::Region(plow.data(),phigh.data(),3);

     			 std::stringstream strStream;
     			 strStream << i->gid() <<"-"<< s->id()<< "-" << j->id();
     			 neuronTrees[cm]->insertData (strStream.str().length(), (byte*)(strStream.str().c_str()), mbr, segmentid);
    		     }
    	    }
    }

    // Make Morphology Rtrees
    Morphologies & myMorphologies = microcircuit.morphologies();
    Morphologies::iterator myMorphologiesEnd = myMorphologies.end();
      for (Morphologies::iterator i = myMorphologies.begin(); i != myMorphologiesEnd; ++i)
      {
      	cout << "Indexing Morphology: " << i->label();
      	string baseName = i->label();
        IStorageManager* diskfile = StorageManager::createNewDiskStorageManager(baseName, 4096);
        ISpatialIndex *tree = RTree::createNewRTree (*diskfile, 0.7, 127, 127, 3,RTree::RV_RSTAR,indexIdentifier);
        indexIdentifier++; segmentid=0;

      	Sections mySections = i->all_sections();
      	Sections::iterator mySectionsEnd = mySections.end();
      	for (Sections::iterator s = mySections.begin(); s != mySectionsEnd; ++s)
      	{
            Segments segments = s->segments();
      		Segments::const_iterator segments_end = segments.end();
      		for (Segments::const_iterator j = segments.begin(); j != segments_end ; ++j)
      		{

      			Box<bbp::Micron> Mbr = AABBCylinder::calculateAABBForCylinder(j->begin().center(),
      								  j->begin().radius(),j->end().center(),j->begin().radius());

      			vect plow, phigh;

      			plow[0] = Mbr.center().x() - Mbr.dimensions().x() / 2;
      			phigh[0] = Mbr.center().x() + Mbr.dimensions().x() / 2;
      			plow[1] = Mbr.center().y() - Mbr.dimensions().y() / 2;
      			phigh[1] = Mbr.center().y() + Mbr.dimensions().y() / 2;
      			plow[2] = Mbr.center().z() - Mbr.dimensions().z() / 2;
      			phigh[2] = Mbr.center().z() + Mbr.dimensions().z() / 2;

				SpatialIndex::Region mbr = SpatialIndex::Region(plow.data(),phigh.data(),3);

      			std::stringstream strStream;
      			strStream << s->id()<< "-" << j->id();

      			tree->insertData (strStream.str().length(), (byte*)(strStream.str().c_str()), mbr, segmentid);
      			segmentid++;
      		}
      	}
      	cout << ".. Total Segments: " << segmentid << "\n";
        tree->~ISpatialIndex();
//.........这里部分代码省略.........

/*!
    \todo Clarify that spatio_temporal_average specifies milliseconds from 
    simulation start, not from simulation time window that is reported in 
    the file. also produce an error message if interval not in file. Should 
    actually probably go in the Stream_Playback_Control then.
*/
void spatiotemporal_average (URI            blue_config_filename, 
							 Filename		volume_compartment_mapping_filename,
							 Filename		output_filename_prefix,
							 Millisecond	begin_time,
							 Millisecond	end_time,
							 Millisecond	average_time_window,
                             bool           portable = false)
{
	Count frames_per_second(10);

    // experiment
    Experiment  experiment;
    experiment.open(blue_config_filename);
	Cell_Target circuit_target = experiment.cell_target("Column");

    std::cout << "Microcircuit: loading target Column" << std::endl;
    Count j = 0;

    for (Cell_Target::const_iterator i = circuit_target.begin();
        i != circuit_target.end(); ++i)
    {
        Neuron_Ptr new_neuron(new Neuron(experiment.microcircuit(), *i , j));
	    experiment.microcircuit().neurons().insert(new_neuron);
        ++j;
    }
	
	Volume_Temporal_Average<Millivolt, Micron, Millisecond>	
		temporal_average(Count (average_time_window * frames_per_second));

    Filter_Data<Microcircuit> microcircuit(experiment.microcircuit_ptr());
    
    // pipeline setup
    Pipeline    pipeline;
        
    std::cout << "Loading compartment report..." << std::endl;
    Reports_Specification::const_iterator result = 
        experiment.reports().find("allCompartments");
    if (result == experiment.reports().end())
    {
        std::cerr << "No \"allCompartments\" report found" << std::endl;
        exit(1);
    }
    Compartment_Report_Stream_Reader<Millivolt> binary_report_reader(* result, 2);
    binary_report_reader.cell_target(circuit_target);
    Compartment_Report_Printer<Millivolt>   printer;
    Compartment_Voltages_to_Microcircuit    report_link;
    Compartment_Voltage_Voxel_Average voltage_volume_average(volume_compartment_mapping_filename, false);
    Volume_Voltage_to_Byte voltage_to_byte;

	binary_report_reader.time_window(begin_time, end_time);
	binary_report_reader.speed(1.0);
	voltage_to_byte.dynamic_range(-65.0, -20.0);

	std::cout << "Pipeline 1: Setup spatio-temporal average." << std::endl;
	pipeline = binary_report_reader
		  + printer
		  + report_link
		  + microcircuit
		  + voltage_volume_average
		  + temporal_average;

	std::cout << "Pipeline 1: Started processing spatio-temporal average." 
		<< std::endl;
	pipeline.start();
	pipeline.wait_until_finished();
	pipeline.stop();
	std::cout << "Pipeline 1: Temporal average finished." << std::endl;


    std::cout << "Pipeline 2: Setup NRRD volume file writer." << std::endl;

    Pipeline    pipeline2;
    Filter_Reader<Volume<Millivolt, Micron, Millisecond> >   volume_buffer_reader;
    volume_buffer_reader.input_ptr(temporal_average.output_ptr());
	Volume_NRRD_File_Stream_Writer nrrd_file_writer(output_filename_prefix);

    pipeline2 = volume_buffer_reader 
                + voltage_to_byte 
                + nrrd_file_writer;

    std::cout << "Pipeline 2: Started writing volumes to disk." << std::endl;
	pipeline2.start();
	pipeline2.wait_until_finished();
	std::cout << "Pipeline 2: Written volumes to disk." << std::endl;
	
}