digital_alloy_optimizer_1d.h++

Go to the documentation of this file.

/*
 *Optimizes a material according to some number.


    Copyright (C) 2010 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_1D_HPP___
#define ___LIBKDOTP_DIGITAL_ALLOY_OPTIMIZER_1D_HPP___

#include <math.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 {
      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 a material to be used if the block is zero
   *\param b 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 err GError information if something failed. Untouched if successful.
   *\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.
   *\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(
           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, 
           unsigned _N_blocks=0,
           long unsigned int rng_seed=0,
           GList* assignments=NULL,
           struct ::kdotp::libmodelxx::postprocessing::function_exp *assignment_expr=NULL,
           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;
  /**This could be gleaned by counting material_is_optimum, but is faster.*/
  unsigned N_optimum;
  /**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;


  /*Internal methods*/
  inline void reset_optimum_array() {
    for(unsigned i=0; i<N_blocks; i++) material_is_optimum[i]=0;
    N_optimum=0;
  }
  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);
      break;
    case MAXIMIZE_ABS:
      fprintf(stderr, "maximize_abs: value=%g previous_value=%g\n", value, previous_value);
      return (fabs(value) > fabs(previous_value));
      break;
    case MINIMIZE:
      return (value < previous_value);
      break;
    case MAKE_ZERO:
      return (fabs(value) < fabs(previous_value));
      break;
    }
    return false;
  }
  inline unsigned get_random_unoptimized_block() {
    unsigned check = (unsigned)floor(gsl_rng_uniform(rng)*(N_blocks - N_optimum-1));
    unsigned j=0;
    /*we need to find the check-th un-optimized block.
     *algorithm: walk through the list incrementing j when we have an unoptimized block
     *  return the block where j becomes check.
     */
    for(unsigned block=0; block<N_blocks; block++) {
      if(!material_is_optimum[block]) {
        if(j++ == check) return block;
      }
    }
    return N_blocks;
  }
  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_block_as_optimum(unsigned block) {
    material_is_optimum[block]=1;
    N_optimum++;
  }

  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 value the value that we're tracking
   *\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_first_step(double value, void* privdat);
  /**Does the later steps in the algorithm
   *\param value the value that we're tracking
   *\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_second_step(double value, void* privdat);
  /**Does the later steps in the algorithm
   *\param value the value that we're tracking
   *\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_subsequent_step(double value, void* privdat);
  /**Does any step in the algorithm, but checks to see if we've taken the first step yet.
   *\param value the value that we're tracking
   *\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_generic_step(double value, 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);
  /**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();
      private:
  digital_alloy_optimizer_1d() {}
  bool operator==(const digital_alloy_optimizer_1d& o) {return false;}
  digital_alloy_optimizer_1d operator=(const digital_alloy_optimizer_1d& o) {return *this;}
  digital_alloy_optimizer_1d(const digital_alloy_optimizer_1d& o){}
      };



    }
  }
}


#endif