C++ Grid_continuous_property类代码示例

bool Unary_property_item_action::execute(GsTL_object_item* item) {

    Error_messages_handler error_messages;
    Grid_continuous_property* prop = dynamic_cast<Grid_continuous_property*>(item);
    if(prop == 0) { //May be inside a region
        GsTL_object_property_item* prop_in_group = dynamic_cast<GsTL_object_property_item*>(item);
        if(prop_in_group)
            prop = prop_in_group->property();
    }

    if(prop==0) {
        GsTLcerr << "Command " << unary_action_name_.toStdString() << " could not be performed; a property was not correctly selected";
        GsTLcerr << gstlIO::end;
        return false;
    }

    std::string grid_name = item->parent()->parent()->item_data(0).toString().toStdString();  //prop->group->grid
    std::string params = grid_name+"::"+prop->name();
    bool ok = project_->execute(unary_action_name_.toStdString(), params , &error_messages);

    if (!ok)
    {
        GsTLcerr << "Command " << unary_action_name_.toStdString() << " could not be performed: \n";
        if (!error_messages.empty())
            GsTLcerr << error_messages.errors() << "\n";
        GsTLcerr << gstlIO::end;
    }
    return ok;
}

void Cosisim::reset_indicator_properties() {
  for( unsigned int i = 0; i < indicators_.size() ; i++ ) {
    Grid_continuous_property* prop = indicators_[i];
    for( unsigned int j = 0; j < prop->size() ; j++ ) {
      if( !prop->is_harddata( j ) )
        prop->set_not_informed( j );
    }
  }
}

void Swap_to_disk_property_item_action::initialize(GsTL_object_item* item) {
    Grid_continuous_property* prop = dynamic_cast<Grid_continuous_property*>(item);
    if(prop == 0) {
        this->setVisible(false);
    }
    else if(prop->is_in_memory()) {
        this->setVisible(true);
    }
    else {
        this->setVisible(false);
    }
}

int Cosisim::execute( GsTL_project* ) {
  // Initialize the global random number generator
  Global_random_number_generator::instance()->seed( seed_ );

  // Set up a progress notifier	
  int total_steps = simul_grid_->size() * (nb_of_realizations_);
  int frequency = std::max( total_steps / 20, 1 );
  SmartPtr<Progress_notifier> progress_notifier = 
    utils::create_notifier( "Running CoSisim", 
			    total_steps, frequency );


  simul_grid_->init_random_path();

  // loop on all realizations
  for( int nreal = 0; nreal < nb_of_realizations_ ; nreal ++ ) {

    // update the progress notifier
    progress_notifier->message() << "working on realization " 
                                 << nreal+1 << gstlIO::end;
    if( !progress_notifier->notify() ) return 1;


    // Create a new property to hold the realization and tell the simulation 
    // grid to use it as the current property 
    appli_message( "Creating new realization" );
    Grid_continuous_property* prop = multireal_property_->new_realization();
    prop->set_parameters(parameters_);
    simul_grid_->select_property( prop->name() );

    typedef Geostat_grid::random_path_iterator RandomPathIterator;
    RandomPathIterator path_begin = simul_grid_->random_path_begin();
    RandomPathIterator path_end = simul_grid_->random_path_end();

    int status = 0;
    if( do_median_ik_ )
      status = median_ik( path_begin, path_end, progress_notifier.raw_ptr() );
    else
      status = full_ik( path_begin, path_end, progress_notifier.raw_ptr() );

    reset_indicator_properties();

    // check if the simulation was aborted
    if( status == 1 ) {
      clean( prop );
      return 1;
    }

  } // end loop on realizations

  remove_temporary_properties();
  return 0;
}

void PostKriging::initialize_operation(std::vector<Grid_continuous_property*>& props, 
		std::vector<float>& vals,Error_messages_handler* errors, 
		const Parameters_handler* parameters, std::string tag)
{
  std::string base_prop_name = parameters->value( tag+"_prop.value" );
	errors->report(base_prop_name.empty(), tag+"_prop", "No name specified" );
  std::string val_names = parameters->value( tag+"_vals.value" );
	errors->report( val_names.empty(), tag+"_vals", "No values specified" );
  if( errors->empty() ) {
	  vals = String_Op::to_numbers<float>(val_names);

	  for(std::vector< float >::iterator it = vals.begin();it != vals.end(); ++it )
	  {
		  std::ostringstream prop_name;
		  prop_name << base_prop_name <<"_tr"<< *it;
      Grid_continuous_property* prop = geostat_utils::add_property_to_grid( grid_,prop_name.str() );
      prop->set_parameters(parameters_);
		  props.push_back(prop);
  		 
	  }
  }
}

int LU_sim::execute( GsTL_project* ) {
  
  // Initialize the global random number generator
  Global_random_number_generator::instance()->seed( seed_ );
  
  
  // Set up a progress notifier	
  int total_steps = simul_grid_->size() * (nb_of_realizations_);
  int frequency = std::max( total_steps / 20, 1 );
  SmartPtr<Progress_notifier> progress_notifier = 
    utils::create_notifier( "Running LU_sim", 
			    total_steps, frequency );

  // In LU simulation, the marginal is a Gaussian cdf, 
  // with mean 0 and variance 1.
  Gaussian_cdf marginal( 0.0, 1.0 );

  //typedef Geostat_grid::random_path_iterator iterator;
  typedef Gval_iterator< TabularMapIndex > iterator;

  std::vector< int > id_unknown;
  std::vector< int > id_data;

  //harddata_grid_->select_property(harddata_property_->name());
  if( harddata_grid_ ) {
    for( int i=0; i< harddata_property_->size(); i++ ) {
      if(!harddata_property_->is_informed(i) ) continue;
      id_data.push_back( i );
    }
  }
  for( int i=0; i< simul_grid_->size(); i++ ) {
    if( simul_grid_ == harddata_grid_ && harddata_property_->is_informed(i) ) continue;
    id_unknown.push_back( i );
  }

  iterator begin_d( harddata_grid_, harddata_property_ ,0,  
			      id_data.size(), TabularMapIndex(&id_data) );
  iterator end_d( harddata_grid_, harddata_property_ ,id_data.size(),  
			      id_data.size(), TabularMapIndex(&id_data) );

  Grid_continuous_property* prop = multireal_property_->new_realization();
  
  iterator begin_u( simul_grid_, prop ,0,  
			      id_unknown.size(), TabularMapIndex(&id_unknown) );
  iterator end_u( simul_grid_, prop ,id_unknown.size(),  
			      id_unknown.size(), TabularMapIndex(&id_unknown) );

  Random_number_generator gen;
  LU_simulator<
        Covariance<Location>,
        Random_number_generator,
        Geostat_grid> 
    lu_sim( covar_,gen, simul_grid_ );
  lu_sim.initialize_matrix(begin_u, end_u, begin_d, end_d);

  // loop on all realizations
  for( int nreal = 0; nreal < nb_of_realizations_ ; nreal ++ ) {

    // update the progress notifier
    progress_notifier->message() << "working on realization " 
                                 << nreal+1 << gstlIO::end;
    if( !progress_notifier->notify() ) return 1;


    // Create a new property to hold the realization and tell the simulation 
    // grid to use it as the current property 
    appli_message( "Creating new realization" );
    //Grid_continuous_property* prop = multireal_property_->new_realization();
    if(nreal>0) prop = multireal_property_->new_realization();
    prop->set_parameters(parameters_);
    simul_grid_->select_property( prop->name() );


    // initialize the new realization with the hard data, if that was requested 
    if( property_copier_ ) {
      property_copier_->copy( harddata_grid_, harddata_property_,
                              simul_grid_, prop );
    }


    iterator begin_unk( simul_grid_, prop ,0,  
			   id_unknown.size(), TabularMapIndex(&id_unknown) );
    iterator  end_unk( simul_grid_, prop ,id_unknown.size(),  
			   id_unknown.size(), TabularMapIndex(&id_unknown) );

    int status = lu_sim(begin_unk,end_unk,marginal);

    // initialize the new realization with the hard data, if that was requested 
  //  if( property_copier_ ) {
  //    property_copier_->copy( harddata_grid_, harddata_property_,
  //                            simul_grid_, prop );
      //initializer_->assign( prop, harddata_grid_, harddata_property_->name() );
  //  }


    if( status == -1 ) {
      clean( prop );
      return 1;
    }
    // back-transform if needed
//.........这里部分代码省略.........

int Indicator_kriging::full_ik( Progress_notifier* progress_notifier ) {
  bool ok = true;
  bool order_relation_problems = false;

  // create all the properties we will populate
  std::vector< Grid_continuous_property* > simul_properties;
  for( int thres = 0; thres < thres_count_; thres++ ) {
    Grid_continuous_property* prop = multireal_property_->new_realization();
    prop->set_parameters(parameters_);
    simul_properties.push_back( prop );
  }

  std::vector<double> krig_weights;
  SK_constraints Kconstraints;
  typedef std::vector<double>::const_iterator weight_iterator;
  typedef SK_combiner< weight_iterator, Neighborhood > SKCombiner;

  // the following line could probably be omitted
  simul_grid_->select_property( simul_properties[0]->name() );

  Geostat_grid::iterator begin = simul_grid_->begin();
  Geostat_grid::iterator end = simul_grid_->end();

  // loop over the grid

  for( ; begin != end; ++begin ) {
    if( !progress_notifier->notify() ) return 1;

    if( begin->is_informed() ) continue;
    

    // for each threshold / class:
    Non_parametric_cdf<float>::p_iterator p_it = ccdf_->p_begin();
    for( int thres = 0; thres < thres_count_; thres++, ++p_it ) {
      neighborhoods_vector_[thres]->find_neighbors( *begin );
 
      if( neighborhoods_vector_[thres]->is_empty() ){
        //if we don't have any conditioning data, use the marginal
        *p_it = marginal_probs_[thres];
        continue;
      }
      
      int status = kriging_weights( krig_weights,
	                          			  *begin,
				                            *(neighborhoods_vector_[thres].raw_ptr()),
				                            covar_vector_[thres], Kconstraints );
      
      if( status != 0 ) {
        // the kriging system could not be solved, issue a warning and skip the
        // node
        ok = false;
        *p_it = marginal_probs_[thres];
        continue;
      }
      
      SKCombiner combiner( marginal_probs_[thres] );
      double estimate = combiner( krig_weights.begin(), 
		                          		krig_weights.end(),
			  	                        *(neighborhoods_vector_[thres].raw_ptr()) );
      *p_it = estimate;
    }

    // make sure the ccdf is a valid cdf
    if( !ccdf_->make_valid() ) {
      // there was a problem making the cdf a valid cdf:
      // leave the node un-estimated and set the flag so that an error will
      // be reported
      order_relation_problems = true; 
      continue;
    }

    GsTLInt node_id = begin->node_id();

    // output the ccdf probabilities to the grid properties
    p_it = ccdf_->p_begin();
    for( int thres2 = 0; thres2 < thres_count_; thres2++, ++p_it ) {
  	  simul_properties[ thres2 ]->set_value( *p_it, node_id );
    }
  }


  if( !ok )
    GsTLcerr << "The kriging system could not be solved for every node\n" 
             << gstlIO::end; 
  if( order_relation_problems ) {
    GsTLcerr << "A cdf could not be estimated for all nodes because of major "
             << "order-relation problems (all probabilities < 0 )" 
             << gstlIO::end;
  }

  return 0;
}

int Indicator_kriging::median_ik( Progress_notifier* progress_notifier ) {
  bool ok = true;
  
  // create all the properties we will populate
  std::vector< Grid_continuous_property* > simul_properties;
  for( int thres = 0; thres < thres_count_; thres++ ) {
    Grid_continuous_property* prop = multireal_property_->new_realization();
    prop->set_parameters(parameters_);
    simul_properties.push_back( prop );
  }

  std::vector<double> krig_weights;
  SK_constraints Kconstraints;
  typedef std::vector<double>::const_iterator weight_iterator;
  typedef SK_combiner< weight_iterator, Neighborhood > SKCombiner;

  // the following line could probably be omitted
  simul_grid_->select_property( simul_properties[0]->name() );

  Geostat_grid::iterator begin = simul_grid_->begin();
  Geostat_grid::iterator end = simul_grid_->end();
  
  
  for( ; begin != end; ++begin ) {
    if( !progress_notifier->notify() ) return 1;

    if( begin->is_informed() ) continue;
      
    neighborhood_->find_neighbors( *begin );
//    if( neighborhood_->is_empty() ){
    if( neighborhood_->size() < min_neigh_ ){
      //if we don't have any conditioning data, skip the node
      continue;
    }
    else {
      int status = kriging_weights( krig_weights,
				    *begin,
				    *(neighborhood_.raw_ptr()),
				    covar_, Kconstraints );
      
      if(status == 0) {
	      // the kriging system could be solved
	      // Since we're using the same covariance and the 
	      // same neighborhood for all thresholds, we can re-use the same
	      // weights for all thresholds 

       	GsTLInt node_id = begin->node_id();
      	Non_parametric_cdf<float>::p_iterator p_it = ccdf_->p_begin();

      	for( int thres = 0; thres < thres_count_; thres++, ++p_it ) {

      	  // tell the neighbors to work on the correct property
	        for( Neighborhood::iterator it = neighborhood_->begin();
	             it != neighborhood_->end(); it++ ) {
      	    it->set_property_array( hdata_properties_[ thres ] );
	        }

    	    SKCombiner combiner( marginal_probs_[thres] );
	        double estimate = combiner( krig_weights.begin(), 
				                              krig_weights.end(),
				                              *(neighborhood_.raw_ptr()) );
	        *p_it = estimate;
	      }

      	// make sure the ccdf is a valid cdf
      	ccdf_->make_valid();

      	// output the ccdf probabilities to the grid properties
      	p_it = ccdf_->p_begin();
      	for( int thres2 = 0; thres2 < thres_count_; thres2++, ++p_it ) {
      	  simul_properties[ thres2 ]->set_value( *p_it, node_id );
      	}
      }
      else {
      	// the kriging system could not be solved, issue a warning and skip the
      	// node
        ok = false;
      }
    }
  }

  if( !ok )
    GsTLcerr << "The kriging system could not be solved for every node\n" << gstlIO::end; 


  return 0;
}

/* Convert_continuous_to_categorical_property grid_name::prop1::[prop2::]
* will copy convert the property from continuous to categorical
*/
bool Convert_continuous_to_categorical_property::
init( std::string& parameters, GsTL_project* proj,
      Error_messages_handler* errors ) {

  std::vector< std::string > params =
    String_Op::decompose_string( parameters, Actions::separator,
                      				   Actions::unique );

  if( params.size() <2 ) {
    errors->report( "some parameters are missing, Needs at least two parameters (grid + property)" );
    return false;
  }

  SmartPtr<Named_interface> grid_ni =
    Root::instance()->interface( gridModels_manager + "/" + params[0] );
  Geostat_grid* grid = dynamic_cast<Geostat_grid*>( grid_ni.raw_ptr() );
  if( !grid ) {
    std::ostringstream message;
    message << "No grid called \"" << params[0] << "\" was found";
    errors->report( message.str() );
    return false;
  }

  for(int i=1; i<params.size(); ++i) {
    Grid_categorical_property* cprop = grid->categorical_property(params[i]);
    if(cprop) continue;
    Grid_continuous_property* prop = grid->property(params[i]);
    if(prop == 0) continue;
    std::set<int> cat_codes;
    std::string prop_name = prop->name()+" - categorical";
    cprop = grid->add_categorical_property(prop_name);
    while(!cprop) {
      prop_name.append("_0");
      cprop = grid->add_categorical_property(prop_name);
    }
    for(int nodeid=0; nodeid < prop->size(); ++nodeid) {
      if( prop->is_informed(nodeid)) {
        int code = static_cast<int>(prop->get_value(nodeid));
        cprop->set_value(code,nodeid);
        cat_codes.insert(code);
      }
    }
    // Check is sequential coding
    bool is_sequential_coding = false;
    std::set<int>::iterator it = cat_codes.begin();
    if( *it == 0  ) {
      std::advance(it, cat_codes.size()-1);
      if( *it == cat_codes.size()-1 ) {
        is_sequential_coding = true;
      }
      
    }

    if( !is_sequential_coding  ) {  // Need to build a categorical definition
      CategoricalPropertyDefinitionName* cat_def = 0;
      std::string catdef_name = grid->name()+"-"+prop->name();
      while(!cat_def) {
        SmartPtr<Named_interface> ni = Root::instance()->new_interface( "categoricaldefinition://"+catdef_name,categoricalDefinition_manager +"/"+catdef_name );
        cat_def = dynamic_cast<CategoricalPropertyDefinitionName*>(ni.raw_ptr());
          if(!cat_def) catdef_name.append("_0");
      }

      std::set<int>::iterator it_cat_codes = cat_codes.begin();
      for(; it_cat_codes != cat_codes.end(); ++it_cat_codes) {
        cat_def->add_category(*it_cat_codes, QString("Code %1").arg(*it_cat_codes).toStdString() );
      }
      cprop->set_category_definition( cat_def );
    }

  }

  proj->update( params[0] );
  return true;
}

int dssim::execute( GsTL_project* ) {
  
  // Initialize the global random number generator
  Global_random_number_generator::instance()->seed( seed_ );
  
  
  // Set up a progress notifier	
  int total_steps = simul_grid_->size() * (nb_of_realizations_);
  int frequency = std::max( total_steps / 20, 1 );
  SmartPtr<Progress_notifier> progress_notifier = 
    utils::create_notifier( "Running dssim", 
			    total_steps, frequency );

  // Initialize the marginal cdf for the direct sequential simulation
 // LogNormal_cdf marginal( 0.1, 1.0 );
 // LogNormal_cdf ccdf;


  // set up the cdf-estimator
  typedef First2_moments_cdf_Kestimator< Covariance<Location>,
                                   Neighborhood,
                                   geostat_utils::KrigingConstraints
                                  >    Kriging_cdf_estimator;
  Kriging_cdf_estimator cdf_estimator( covar_,
				       *Kconstraints_,
				       *combiner_ );

  // set up the sampler
  Random_number_generator gen;
  Monte_carlo_sampler_t< Random_number_generator > sampler( gen );
  
  // compute the random path
  Grid_path path(simul_grid_, simul_grid_->selected_property(), target_grid_region_);
  //simul_grid_->init_random_path();

  // loop on all realizations
  for( int nreal = 0; nreal < nb_of_realizations_ ; nreal ++ ) {

    // update the progress notifier
    progress_notifier->message() << "working on realization " 
                                 << nreal+1 << gstlIO::end;
    if( !progress_notifier->notify() ) return 1;


    // Create a new property to hold the realization and tell the simulation 
    // grid to use it as the current property 
    appli_message( "Creating new realization" );
    Grid_continuous_property* prop = multireal_property_->new_realization();
    prop->set_parameters(parameters_);
    simul_grid_->select_property( prop->name() );
    neighborhood_->select_property( prop->name() );

    // initialize the new realization with the hard data, if that was requested 
    if( property_copier_ ) {
      property_copier_->copy( harddata_grid_, harddata_property_,
                              simul_grid_, prop );
      //initializer_->assign( prop, harddata_grid_, harddata_property_->name() );
    }

    path.randomize();
    path.set_property(prop->name());

    appli_message( "Doing simulation" );
    // do the simulation
    int status = 
      sequential_simulation( path.begin(),
			     path.end(),
			     *(neighborhood_.raw_ptr()),
			     *ccdf_,
			     cdf_estimator,
			     *marginal_,
			     sampler, progress_notifier.raw_ptr()
			     );
    if( status == -1 ) {
      clean( prop );
      return 1;
    }
  }

  clean();

  return 0;
}