continuous_alloy_optimizer_1d.c++

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

*/

#include <gsl/gsl_errno.h>
#include <libmetacalc/continuous_alloy_optimizer_1d.h++>
#include <libmodelxx/cartesian_grid.h++>
#include <libmodelxx/randomness.h++>
#include <libkdotp/gcc_hints.h>
#include <libmodelxx/function_specification.h++>

#define CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN 73

#define EPSILON 1e-9

using namespace kdotp::libmodelxx::grid;
using namespace kdotp::libmodelxx::structure;

namespace kdotp {
  namespace libmetacalc {
    namespace optimizer {

      continuous_alloy_optimizer_1d
      ::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 matdb *_mdb,
           double _minx,
           double _maxx,
           long unsigned int rng_seed,
           double starting_stepsize,
           unsigned _N_blocks,
           GList* assignments,
           struct ::kdotp::libmodelxx::postprocessing::function_exp *assignment_expr,
           GHashTable *global_symbols
           )
  : mat1(NULL),
    mat0(NULL),
    current_assignments(NULL),
    material_is_optimum(NULL),
    genstruct(s),
    call_set_point_callback(_set_point),
    blocksize(_blocksize),
    points_in_block(NULL),
    point_start(NULL),
    previous_value(0),
    current_block(NULL),
    goal(_goal),
    seed(rng_seed),
    stepsize(starting_stepsize),
    mdb(_mdb),
    minx(_minx),
    maxx(_maxx)
      {
  cartesian_grid<material*, 1, 1>* matgrid = (cartesian_grid<material*, 1, 1>*)(genstruct->material_grid);
  /*Assign material strings.*/
  GString *str = g_string_new(_mat1);
  mat1=str->str;
  g_string_free(str, FALSE);
  str = g_string_new(_mat0);
  mat0=str->str;
  g_string_free(str, FALSE);
  /*Easy way to make this a GList: start off with list=NULL, create a new (start,end) struct
   *  on start, save it to the list on end.  leave the started variable as it is.
   */
  /*If marker_mat_in_grid is NULL, we've not started yet.*/
  material** marker_mat_in_grid = NULL;
  material* nullmat = NULL;
  for(unsigned X=0; X<matgrid->L[0]; X++) {
    if((**(matgrid->get_point_r(&X))) == *marker_material) {
      /*Set this point to NULL to prevent deleting this material*/
      matgrid->set_point_(&X, nullmat);
      if(marker_mat_in_grid == NULL) {
        start_X = X;
        marker_mat_in_grid = matgrid->get_point_r(&X);
      }
    }else if(marker_mat_in_grid != NULL) {
      /*We've stopped*/
      stop_X = X;
      /*Mark this region with NULL so we don't re-free the material.*/
      delete *marker_mat_in_grid;
      marker_mat_in_grid = NULL;
    }
  }
  if(marker_mat_in_grid != NULL) {
    /*We started but hit the end without switching out of the material. Odd, but solvable.*/
    stop_X = matgrid->L[0];
    /*Mark this region with NULL so we don't re-free the material.*/
    delete *marker_mat_in_grid;
  }
  /*We could probably do this with rounding, but this only happens once and works.*/
  if(_N_blocks > 0) {
    N_blocks = (unsigned)_N_blocks;
    blocksize = ((stop_X-start_X) - (stop_X-start_X)%N_blocks)/N_blocks;
  }else{
    N_blocks = ((stop_X-start_X) - (stop_X-start_X)%blocksize)/blocksize;
  }
  /*Allocate arrays*/
  points_in_block = new unsigned[N_blocks];
  point_start = new unsigned[N_blocks];
  current_assignments = new double[N_blocks];
  /*Previous assignments is NULL in order to know that we have no previous assignment or value.
   */
  material_is_optimum = new short unsigned[2*N_blocks];
  /*Map the points to the blocks*/
  for(unsigned i=0; i<N_blocks; i++) {
    points_in_block[i] = blocksize;
  }
  unsigned points_remaining = (stop_X - start_X)-blocksize*N_blocks;
  if(points_remaining != 0) {
    unsigned pos=points_remaining/N_blocks;
    for(unsigned pt=0; pt<points_remaining; pt++) {
      points_in_block[pos]++;
      pos += points_remaining/N_blocks;
      if(pos >= N_blocks) pos = pos%N_blocks;
    }
  }
  /*Put in the offsets.*/
  unsigned final_point = start_X;
  for(unsigned block=0; block<N_blocks; block++) {
    point_start[block] = final_point;
    final_point += points_in_block[block];
  }

  /*allocate and then seed the random number generator*/
  GError *e = NULL;
  gsl_rng_env_setup();
  rng = gsl_rng_alloc(gsl_rng_default);
  if(rng == NULL) {
    *err = g_error_new(CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 0, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::continuous_alloy_optimizer_1d(constructor): Failed to set up GSL random number generator.");
    return;
  }
  /*Set the seed randomly if it's not already set*/
  if(rng_seed == 0) {
    /*A seed value of 0 is special to GSL, so don't allow a seed of 0. Loop until it's not, or an error has occurred*/
    seed=0;
    while((seed == 0) && (e == NULL)) {
      seed = kdotp::libmodelxx::randomness::get_urandom_unsigned_long_int(&e);
    }
    if(e != NULL) {
      g_propagate_prefixed_error(err, e, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::continuous_alloy_optimizer_1d(constructor): Failed to get random seed:\n\t");
      return;
    }
    fprintf(stderr, "rng_seed was zero, so got new one (seed=%lu rng_seed=%lu)\n", seed, rng_seed);
  }else{
    fprintf(stderr, "rng_seed was non-zero, so keeping it (seed=%lu rng_seed=%lu)\n", seed, rng_seed);
  }
  fprintf(stderr, "\n\nseed is %lu\n\n", seed);
  /*We now seed the rng*/
  gsl_rng_set(rng, seed);
  /*Begin initialization.*/
  initialize(privdat);
  /*At this point, we should be set up.
   *The user should perform their first calculation, and call us back using one of the do_*_step funx.*/
      }

      void continuous_alloy_optimizer_1d::initialize(void* privdat) {
  if(current_block != NULL) {
    delete current_block;
    current_block = NULL;
  }
  reset_optimum_array();
  randomize_array(current_assignments, N_blocks);
  set_material_from_current_assignments(privdat);
      }

      void continuous_alloy_optimizer_1d::initialize(const short unsigned* state, void* privdat) {
  if(current_block != NULL) {
    delete current_block;
    current_block = NULL;
  }
  reset_optimum_array();
  for(unsigned i=0; i<N_blocks; i++) current_assignments[i] = state[i];
  set_material_from_current_assignments(privdat);
      }

      continuous_alloy_optimizer_1d::~continuous_alloy_optimizer_1d() {
  if(mat1 != NULL) g_free(mat1);
  if(mat0 != NULL) g_free(mat0);
  if(current_assignments != NULL) delete[] current_assignments;
  if(material_is_optimum != NULL) delete[] material_is_optimum;
  if(points_in_block != NULL) delete[] points_in_block;
  if(current_block != NULL) delete current_block;
      }

      /**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 continuous_alloy_optimizer_1d::do_first_step(double value, void* privdat) {
  /*There's no way we can yet be converged.  We don't even know if this
   *  point is optimal.  This is our first point.*/
  previous_value=value;
  current_block = new unsigned;
  double sign;
  *current_block = get_random_unoptimized_block(&sign);
  /*Save this block's current composition.*/
  prev_x = get_x(*current_block);
  step_block(*current_block, privdat);
  return 0;
      }
      /**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 continuous_alloy_optimizer_1d::do_subsequent_step(double value, void* privdat) {
  if(is_improved(value)) {
    /*This is an improved point!*/
    /*We don't know right now if any other points are optimized.*/
    reset_optimum_array();
    /*This direction (+ or -) succeeded, so the *other* direction will fail and therefore is
     *  already optimized.
     *So set that the *other* direction is optimized.
     */
    set_block_as_optimum((*current_block + N_blocks)%(2*N_blocks));
    previous_value = value;
    fprintf(stderr, "RESET ARRAY: previous_value=%g (value=%g)\n", previous_value, value);
  }else{
    /*Choosing this point worsened the structure.  Go back, set this as optimized, and choose a new point.*/
    /*Algorithm to get sign without a branch:
     *1) Get the real block
     *2) Get the offset between the signed and real block.  This will either be 0 or N_blocks
     *3) Divide by N_blocks/2 (maps to range [0, 2])
     *4) Subtract 1 (maps to range [-1, 1])
     *\note that we need -1 times this, to undo the previous step
     *\note that we now reset the block to its previous x instead of adding and then subtracting
     */
    set_block(*current_block, privdat, prev_x);
    fprintf(stderr, "REVERT BLOCK %u: previous_value=%g (value=%g)\n", *current_block, previous_value, value);
    /**We know that the block that just got tested was already optimal.
     *\note that we don't know whether the other direction (+ or -) is optimized!  */
    set_block_as_optimum(*current_block);
    /*Note that we don't set the previous value to this current value;
     *it's not an improvement!
     */
  }
  /*Check to see if they're all optimized.*/
  if(N_optimum == 2*N_blocks) {
    fprintf(stderr, "ALPHA IS LOCAL OPTIMUM: N_optimum=%u 2*N_blocks=%u\n", N_optimum, 2*N_blocks);
    return 1;
  }else{
    fprintf(stderr, "ALPHA IS (NOT) LOCAL OPTIMUM: N_optimum=%u 2*N_blocks=%u\n", N_optimum, 2*N_blocks);
  }
  /*Get the randomly-selected new block.*/
  double sign;
  bool success;
  do {
    *current_block = get_random_unoptimized_block(&sign);
    prev_x = get_x(*current_block);
    /*Set the block.*/
    success = step_block(*current_block, privdat);
    if(!success) {
      /*If we failed, this block is optimum!*/
      set_block_as_optimum(*current_block);
      /*We may have finished optimizing!*/
      if(N_optimum == 2*N_blocks) {
        fprintf(stderr, "stepping block %u failed, and ALPHA IS LOCAL OPTIMUM: N_optimum=%u 2*N_blocks=%u\n",  *current_block, N_optimum, 2*N_blocks);
        return 1;
      }else{
        fprintf(stderr, "stepping block %u failed, but ALPHA IS (NOT) LOCAL OPTIMUM: N_optimum=%u 2*N_blocks=%u\n", *current_block, N_optimum, 2*N_blocks);
      }
    }
  }while(!success);
  return 0;
      }
      /**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 continuous_alloy_optimizer_1d::do_generic_step(double value, void* privdat) {
  if(current_block == NULL) {
    return do_first_step(value, privdat);
  }
  return do_subsequent_step(value, 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
       *\param step_fraction when we're optimized at the current stepsize, we multiply the stepsize by this factor to get the new stepsize
       *\param min_stepsize minimum stepsize to iterate over.  When stepsize gets below this value, we stop.
       *\return NULL if we converged; GError* if there was an error.
       *\note that if NaN is returned, we stop and return an error
       *\note that initialize() should have been called before calling this function, either directly or implicitly in the constructor.
       */
      GError* continuous_alloy_optimizer_1d::do_optimization(double (*get_value)(void* privdat, unsigned assignment, long unsigned iteration), void* get_value_privdat, void* set_point_privdat, double step_fraction, double min_stepsize) {
  long unsigned iteration=0;
  do {
    double value = get_value(get_value_privdat, 0, iteration++);
    if(unlikely(isnan(value))) {
      return g_error_new(CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::do_optimization: Got NaN (%g) from get_value callback in the %lu%s iteration.  Iteration has been halted.\n", value, iteration, (iteration%10==1&&iteration!=11)?"st":(iteration%10==2&&iteration!=12)?"nd":(iteration%10==3&&iteration!=13)?"rd":"th");
    }else if(unlikely(isinf(value))) {
      return g_error_new(CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::do_optimization: Got inf (%g) from get_value callback in the %lu%s iteration.  Iteration has been halted.\n", value, iteration, (iteration%10==1&&iteration!=11)?"st":(iteration%10==2&&iteration!=12)?"nd":(iteration%10==3&&iteration!=13)?"rd":"th");
    }
    do_first_step(value, set_point_privdat);
    do {
      value = get_value(get_value_privdat, 0, iteration++);
      if(unlikely(isnan(value))) {
        return g_error_new(CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::do_optimization: Got NaN (%g) from get_value callback in the %lu%s iteration.  Iteration has been halted.\n", value, iteration, (iteration%10==1&&iteration!=11)?"st":(iteration%10==2&&iteration!=12)?"nd":(iteration%10==3&&iteration!=13)?"rd":"th");
      }else if(unlikely(isinf(value))) {
        return g_error_new(CONTINUOUS_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::continuous_alloy_optimizer_1d.c::do_optimization: Got inf (%g) from get_value callback in the %lu%s iteration.  Iteration has been halted.\n", value, iteration, (iteration%10==1&&iteration!=11)?"st":(iteration%10==2&&iteration!=12)?"nd":(iteration%10==3&&iteration!=13)?"rd":"th");
      }
    }while(!(do_subsequent_step(value, set_point_privdat)));
    /*We can't optimize anymore at this step size.*/
    stepsize *= step_fraction;
    /*Free the current_block so it'll be reallocated after first iteration
     *Then reset our optimized steps..
     */
    delete current_block;
    reset_optimum_array();
  }while(stepsize >= min_stepsize);
  return NULL;
      }

    }
  }
}