digital_alloy_optimizer_1d_mpi.h++

Go to the documentation of this file.

/*
 *Optimizes a material according to some number.


    Copyright (C) 2011 Joseph Pingenot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#ifndef ___LIBKDOTP_DIGITAL_ALLOY_OPTIMIZER_MPI_1D_HPP___
#define ___LIBKDOTP_DIGITAL_ALLOY_OPTIMIZER_MPI_1D_HPP___

#include <math.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <libmodelxx/material.h++>
#include <libmodelxx/generic_structure.h++>

namespace kdotp {
  namespace libmetacalc {
    namespace optimizer {

      class digital_alloy_optimizer_1d_mpi {
      public:
  typedef void (*set_point_callback)(unsigned X, const kdotp::libmodelxx::structure::material* mat, void* privdat);
  enum goal {
    MAXIMIZE,
    MAXIMIZE_ABS,
    MINIMIZE,
    MAKE_ZERO
  };

  /**constructor and initializer in one
   *\param marker_material material used to mark the region to be optimized
   *\param _mat1 material to be used if the block is zero
   *\param _mat0 material to be used if the block is one
   *\param s generic_structure to be optimized
   *\param _set_point callback to use when setting a point, or NULL if you don't want to do anything
   *\param _blocksize number of points per block
   *\param _goal the goal of the optimization
   *\param err GError information if something failed. Untouched if successful.
   *\param privdat private data for the set_points callback (used when initializing)
   *\param N_nodes number of nodes being executed on
   *\param node node number for this program.
   *\param rng_seed seed for the random number generator (default is 0, i.e. no seed)
   *\param assignements list of initial material assignments for each block (or for all blocks)
   *\param assignement_expr function to evaluate to get the initial material assignments.
   *\param global_symbols global symbol table
   *\param temperature if specified or non-zero, use simulated annealing at the specified temperature.  If the temperature is never set, the steepest-descent-related Monte Carlo algorithm will be used.
   *\note Starting and stopping points are fetched from the generic structure.
   *\note that the assignments are done on the block level,
   *i.e. in blocks of points.  This is to improve efficiency and
   *accuracy, since structures can only be grown so finely, but
   *having a grid spacing which is a multiple of the minimum
   *feature size is desirable for high-accuracy calculations.
   *\note that any incommensurate number of points will be split
   *amongst the points as evenly as possible.
   *\note that this modifies the generic_structure, removing the marker material and replacing it with materials from the optimizer state.
   *\note that this constructor calls _set_point, since it has to initialize the material.
   */
  digital_alloy_optimizer_1d_mpi(
           kdotp::libmodelxx::structure::material* marker_material, 
           kdotp::libmodelxx::structure::material* _mat1, 
           kdotp::libmodelxx::structure::material* _mat0, 
           kdotp::libmodelxx::structure::generic_structure* s, 
           set_point_callback _set_point,
           unsigned _blocksize, 
           enum goal _goal, 
           GError** err, 
           void* privdat, 
           int _N_nodes,
           int _node,
           unsigned _N_blocks=0,
           long unsigned int rng_seed=0,
           GList* assignments=NULL,
           struct ::kdotp::libmodelxx::postprocessing::function_exp *assignment_expr=NULL,
           double temperature=0,
           GHashTable *global_symbols=NULL
           );
  
  /*The two materials we can be using at each site*/
  kdotp::libmodelxx::structure::material* mat1;
  kdotp::libmodelxx::structure::material* mat0;
  /*This stores the current material assignments.*/
  short unsigned *current_assignments;
  /*This stores the set of bits which state that the material at that location is known optimum for this particular point*/
  short unsigned *material_is_optimum;
  /*Pointer to the generic structure we're to modify.
   *\note Please note that although it's perfectly doable to just set the genstruct
   *  and then use it to create the Hamiltonian over and over and over again, that it's
   *  much much more efficient to use the callback to also set the Hamiltonian directly.
   *\note To manipulate the Hamiltonian directly, you'll need to set the material yourself
   *  for now.  I may later add an interface through the generic hamltonian class to accomplish
   *  this.
   */
  kdotp::libmodelxx::structure::generic_structure* genstruct;
  /*The callback to use when setting a point*/
  set_point_callback call_set_point_callback;
  /*Starting and ending points. Derived in the constructor from
    the generic_structure.  These points are *in*clusive. */
  unsigned start_X;
  /*This one is exclusive, i.e. this is the first point after the optimizer region.*/
  unsigned stop_X;
  /*Size of a block, in points*/
  unsigned blocksize;
  /*Stores the number of points in each block*/
  unsigned* points_in_block;
  /*This stores the point that the start of the block corresponds to*/
  unsigned* point_start;
  /**the number of blocks*/
  unsigned N_blocks;
  /**Previous value we're to track, to compare how we're changing.
   *\note this is a *pointer*; if it's NULL, the previous calculation was the first one, so all of the blocks had been changed.
   */
  double previous_value;
  /**Which block was last changed.
   *\note this is a *pointer*; if it's NULL, the previous calculation was the first one, so all of the blocks had been changed.
   *\note this is the block which was just tested
   */
  unsigned* current_block;
  enum goal goal;
  /**The GSL random number generator information.*/
  gsl_rng* rng;
  unsigned long int seed;
        unsigned N_nodes;
  unsigned node;
  /*Temperature for simulated annealing*/
  double T;

  /*array containing the first block assigned to each node.  There are N_nodes+1 values (the last is 
    the off the end of the array!*/
  int* node_assignments;
  /*The results of this node's calculation*/
  double* results;
  /*The results of modifying the requested sites.  Only on the root node; there are N_blocks nodes.
   */
  double* all_results;

  inline void randomize_array(short unsigned* array, unsigned int N) {
    for(unsigned i=0; i<N; i++) {
      /**\note that gsl_rng_uniform returns [0,1), so 0.5 will be considered in the upper half,
       *  so that both options will be equally probable
       *\note Upper half of [0,1) => 1; lower => 0.
       */
      if(gsl_rng_uniform(rng) >= 0.5) {
        array[i]=1;
      }else{
        array[i]=0;
      }
    }
  }
  inline bool is_improved(double value) {
    switch(goal) {
    case MAXIMIZE:
      return (value > previous_value);
    case MAXIMIZE_ABS:
      return (fabs(value) > fabs(previous_value));
    case MINIMIZE:
      return (value < previous_value);
    case MAKE_ZERO:
      return (fabs(value) < fabs(previous_value));
    }
    return false;
  }
  inline double dist_to_goal(double value) {
    switch(goal) {
    case MAXIMIZE:
      return value - previous_value;
    case MAXIMIZE_ABS:
      return fabs(value) - fabs(previous_value);
    case MINIMIZE:
      return previous_value - value;
    case MAKE_ZERO:
      return fabs(previous_value) - fabs(previous_value);
    }
    return false;
  }
  inline void toggle_block(unsigned block, void* privdat) {
    if(current_assignments[block] == 0) {
      current_assignments[block] = 1;
    }else{
      current_assignments[block] = 0;
    }
    set_block(block, privdat);
  }
  inline void set_block(unsigned block, void* privdat) {
	  ::kdotp::libmodelxx::grid::cartesian_grid< ::kdotp::libmodelxx::structure::material*, 1, 1>* matgrid
      = (::kdotp::libmodelxx::grid::cartesian_grid< ::kdotp::libmodelxx::structure::material*, 1, 1>*)
      (genstruct->material_grid);
    unsigned final_point = point_start[block];
    for(unsigned pt=0; pt<points_in_block[block]; pt++) {
      matgrid->set_point(&final_point, (current_assignments[block]==0)?&mat0:&mat1);
      if(call_set_point_callback != NULL)
        call_set_point_callback(final_point, (current_assignments[block]==0)?mat0:mat1, privdat);
      final_point++;
    }
  }

  inline void set_material_from_current_assignments(void* privdat) {
    for(unsigned block=0; block<N_blocks; block++) {
      set_block(block, privdat);
    }
  }
  /**Does the first step in the algorithm.
   *\param privdat pointer to data for the set_point callback.
   *\return 0 if we've not converged, or 1 if we've converged.
   */
  int do_step(void* privdat);

  /*This performs the calculation without stepping through manually.
   *At each step, get_value is called with privdat
   *\param get_value callback to use after setting the materials.
   *\param privdat private data for the callback
   *\return NULL if we converged; GError* if there was an error.
   *\note that if NaN is returned, we stop and return an error
   *\note this should eventually deal with GErrors for better error reporting.
   */
  GError*  do_optimization(double (*get_value)(void* privdat, unsigned assignment, long unsigned iteration), void* get_value_privdat, void* set_point_privdat, double anneal_start, double anneal_end, unsigned n_anneal_steps, unsigned n_anneal_cycles, unsigned num_calculations_per_annealing_step);
  /**Resets the state to get ready to do an optimization.
   *\note that this randomizes the current assignments array, so this is what
   *  you'll want to call to start from another initial condition
   */
  void initialize(void* privdat);
  /**Initializes to a known state.
   *\param state the state to initialize to (array of N_blocks short unsigned ints)
   */
  void initialize(const short unsigned* state, void* privdat);
  /*this should be virtual if we had any real inheritance, but we don't*/
  ~digital_alloy_optimizer_1d_mpi();
      private:
  digital_alloy_optimizer_1d_mpi() {}
  bool operator==(const digital_alloy_optimizer_1d_mpi& o) {return false;}
  digital_alloy_optimizer_1d_mpi operator=(const digital_alloy_optimizer_1d_mpi& o) {return *this;}
  digital_alloy_optimizer_1d_mpi(const digital_alloy_optimizer_1d_mpi& o){}
      };



    }
  }
}


#endif