continuous_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/>.

NOTES:
* "block" is a set of points which are assigned a material and treated as a unit.
* "signed" vs "unsigned" blocks: because we can make a step in either a positive
     or negative direction, we have two options for every block.  Thus, we have
     2*N_blocks options to check, and the optimized-block tracking array must therefore
     be 2*N_blocks.  Ultimately, however, both positive and negative steps map to the same
     block for material assignments, which is the "unsigned" block.
*/

#ifndef ___LIBKDOTP_CONTINUOUS_ALLOY_OPTIMIZER_1D_HPP___
#define ___LIBKDOTP_CONTINUOUS_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++>
#include <libmodelxx/function_specification.h++>

namespace kdotp {
  namespace libmetacalc {
    namespace optimizer {

      class continuous_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 (ignored if _N_blocks > 0)
   *\param err GError information if something failed. Untouched if successful.
   *\param minx minimum x to be allowed
   *\param maxx maximum x to be allowed
   *\param _N_blocks Number of blocks to allocate, or <=0 to use the specified blocksize.
   *\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.
   */
  continuous_alloy_optimizer_1d(
           kdotp::libmodelxx::structure::material* marker_material, 
           const char* _mat1,
           const char* _mat0,
           kdotp::libmodelxx::structure::generic_structure* s, 
           set_point_callback _set_point,
           unsigned _blocksize, 
           enum goal _goal, 
           GError** err, 
           void* privdat, 
           struct ::kdotp::libmodelxx::structure::matdb *_mdb,
           double _minx=0.0,
           double _maxx=1.0,
           long unsigned int rng_seed=0,
           double starting_stepsize=0.1,
           unsigned _N_blocks=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*/
  char* mat1;
        char* mat0;
  /*This stores the current material assignments.*/
  double *current_assignments;
  /*This stores the set of bits which state that the material at that location is known
   *optimum for this particular point
   *\note this is 2*N_blocks: one for each way we could step (+stepsize, -stepsize)
   *\note array is ordered as this: +block0, +block1, +block2, ... , -block0, -block1, -block2 ...
   */
  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
   *\note that this is in the range [0, 2*N_blocks-1] to preserve the sign information.
   */
  unsigned* current_block;
  enum goal goal;
  /**The GSL random number generator information.*/
  gsl_rng* rng;
  unsigned long int seed;
  /**This is our current stepsize.*/
  double stepsize;
  /**previous block x, so we can reset the block to its previous state*/
  double prev_x;
  struct kdotp::libmodelxx::structure::matdb *mdb;
  /*Min and max x values*/
  double minx;
  double maxx;

  /*Internal methods*/
  inline void reset_optimum_array() {
    for(unsigned i=0; i<2*N_blocks; i++) material_is_optimum[i]=0;
    N_optimum=0;
  }
  inline void randomize_array(double* array, unsigned int N) {
    for(unsigned i=0; i<N; i++) {
      /**\note that gsl_rng_uniform returns [0,1) but that's good enough.
       */
      array[i]=minx+gsl_rng_uniform(rng)*(maxx-minx);
    }
  }
  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;
  }
  /**gets a random unoptimized block.
   *\param sign pointer to a double to store either + or -1, depending on the block.
   *\note that double pointed to by sign will be either +1 or -1; mul by stepsize to get the amount.
   *\return the block which is returned is [0, 2*N_blocks-1]
   */
  inline unsigned get_random_unoptimized_block(double* sign) {
    unsigned check = (unsigned)floor(gsl_rng_uniform(rng)*(2*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.
     */
    /**\note we first walk through the first half, then the second.*/
    *sign = 1.0;
    for(unsigned block=0; block<N_blocks; block++) {
      if(!material_is_optimum[block]) {
        if(j++ == check) {
    return block;
        }
      }
    }
    /**we're now in the negative blocks.*/
    *sign = -1.0;
    for(unsigned block=0; block<N_blocks; block++) {
      if(!material_is_optimum[block]) {
        if(j++ == check) {
    return N_blocks+block;
        }
      }
    }
    return 2*N_blocks;
  }
  /**Steps a block by the specified step
   *\param signed block
   *\param privdat private data for the set_point callback.
   *\return true if stepping succeeded; false if stepping failed.
   */
  inline bool step_block(unsigned block, void* privdat) {
    double sign = -1;
    if(block < N_blocks) sign=1;
    if(current_assignments[block%N_blocks]+sign*stepsize > maxx) return false;
    if(current_assignments[block%N_blocks]+sign*stepsize < minx) return false;
    fprintf(stderr, "STEPPING: step_block: sign=%g stepsize=%g block[%u->%u]=%g->", sign, stepsize, block, block%N_blocks, current_assignments[block%N_blocks]);
    current_assignments[block%N_blocks] += sign*stepsize;
    fprintf(stderr, "%g\n", current_assignments[block%N_blocks]);
    set_block(block%N_blocks, privdat);
    return true;
  }
  /**Gets the x value of a block
   *\param signed block
   */
  inline double get_x(unsigned block) {
    return current_assignments[block%N_blocks];
  }
  /**Sets a block to the specified x value
   *\param signed block
   */
  inline void set_block(unsigned block, void* privdat, double x) {
    current_assignments[block%N_blocks] = x;
    set_block(block%N_blocks, privdat);
  }
  /**\note block is an unsigned block*/
  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++) {
	    ::kdotp::libmodelxx::structure::material** m;
      if((m=matgrid->get_point_r(&final_point)) != NULL) delete *m;
	    ::kdotp::libmodelxx::structure::material* mat = new ::kdotp::libmodelxx::structure::material(mdb, mat1, mat0, current_assignments[block]);
      matgrid->set_point(&final_point, mat);
      if(call_set_point_callback != NULL)
        call_set_point_callback(final_point, mat, 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, double step_fraction=0.33, double min_stepsize=0.001);
  /**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*/
  ~continuous_alloy_optimizer_1d();
      private:
  continuous_alloy_optimizer_1d() {}
  bool operator==(const continuous_alloy_optimizer_1d& o) {return false;}
  continuous_alloy_optimizer_1d operator=(const continuous_alloy_optimizer_1d& o) {return *this;}
  continuous_alloy_optimizer_1d(const continuous_alloy_optimizer_1d& o){}
      };



    }
  }
}


#endif