digital_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/digital_alloy_optimizer_1d.h++>
#include <libmodelxx/cartesian_grid.h++>
#include <libmodelxx/randomness.h++>
#include <libkdotp/gcc_hints.h>

#define DIGITAL_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 {

      digital_alloy_optimizer_1d
      ::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,
           long unsigned int rng_seed,
           GList* assignments,
           struct ::kdotp::libmodelxx::postprocessing::function_exp *assignment_expr,
           GHashTable *global_symbols
           )
  : mat1(_mat1),
    mat0(_mat0),
    current_assignments(NULL),
    material_is_optimum(NULL),
    genstruct(s),
    call_set_point_callback(_set_point),
    blocksize(_blocksize),
    points_in_block(NULL),
    point_start(NULL),
    N_blocks(_N_blocks),
    previous_value(0),
    current_block(NULL),
    goal(_goal),
    seed(rng_seed)
      {
  cartesian_grid<material*, 1, 1>* matgrid = (cartesian_grid<material*, 1, 1>*)(genstruct->material_grid);
  /*First, we need to find the extents of the region.*/
  bool started=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.
   */
  for(unsigned X=0; X<matgrid->L[0]; X++) {
    if((**(matgrid->get_point_r(&X))) == *marker_material) {
      if(!started) {
        start_X = X;
        started=1;
      }
    }else if(started) {
      /*We've stopped*/
      stop_X = X;
      started=0;
    }
  }
  if(started) {
    /*We started but hit the end without switching out of the material. Odd, but solvable.*/
    stop_X = matgrid->L[0];
  }
  /*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 short unsigned[N_blocks];
  /*Previous assignments is NULL in order to know that we have no previous assignment or value.*/
  material_is_optimum = new short unsigned[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;
  if(points_remaining != 0) {
    fprintf(stderr, "N_blocks=%u\n", N_blocks);
    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(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 0, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d.c::digital_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::digital_alloy_optimizer_1d.c::digital_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.*/
  if((assignments != NULL) || (assignment_expr != NULL)) {
    /*NOTE that assignments takes precedence over assignment_expr*/
    short unsigned* state = new short unsigned[N_blocks];
    if(assignments != NULL) {
      if(assignments->next == NULL) {
        /*Assign each block based on the single user-supplied value.*/
        if(*(double*)(assignments->data) < 0.5) {
    for(unsigned i=0; i<N_blocks; i++) state[i] = 0;
        }else{
    for(unsigned i=0; i<N_blocks; i++) state[i] = 1;
        }
      }else{
        /*Assign each block based on the user's specification.*/
        GList* ll=assignments;
        for(unsigned i=0; i<N_blocks; i++) {
    state[i] = (*(double*)(ll->data) < 0.5)?0:1;
    ll = ll->next;
    if(ll == NULL) {
      *err = g_error_new(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 5, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d.c::digital_alloy_optimizer_1d(constructor): ERROR while assigning materials according to user input: there were not enough values for each block (Got %u values but had %u blocks!\n", i+1, N_blocks);
      return;
    }
        }
      }
    }else{
      /*Evaluate the expression at each point in the blocks.*/
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const start_x = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const mid_x = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const end_x = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const start_X = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const mid_X = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      struct ::kdotp::libmodelxx::postprocessing::symbol_table_entry_const end_X = {::kdotp::libmodelxx::postprocessing::CONSTANT_SYMBOL, ::kdotp::libmodelxx::postprocessing::SCALAR, NULL};
      /**TODO: need to finish symbol assignments, make the local symbol table, and then do the loop.*/
    }
    initialize(state, privdat);
  }else{
    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 digital_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 digital_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);
      }

      digital_alloy_optimizer_1d::~digital_alloy_optimizer_1d() {
  if(mat1 != NULL) delete mat1;
  if(mat0 != NULL) delete 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 digital_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;
  *current_block = get_random_unoptimized_block();
  toggle_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 digital_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();
    /**Set the previous value.
     *\note if structure is now optimum, we have the final value in previous_value.
     */
    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.*/
    toggle_block(*current_block, privdat);
    fprintf(stderr, "REVERT BLOCK %u: previous_value=%g (value=%g)\n", *current_block, previous_value, value);
    /*Note that we don't set the previous value to this current value;
     *it's not an improvement!
     */
  }
  /*We know that the block that just got changed is optimized.*/
  set_block_as_optimum(*current_block);
  /*Check to see if they're all optimized.*/
  if(N_optimum == N_blocks) {
    fprintf(stderr, "ALPHA IS LOCAL OPTIMUM: N_optimum=%u N_blocks=%u\n", N_optimum, N_blocks);
    return 1;
  }else{
    fprintf(stderr, "ALPHA IS (NOT) LOCAL OPTIMUM: N_optimum=%u N_blocks=%u\n", N_optimum, N_blocks);
  }
  /*Get the randomly-selected new block.*/
  *current_block = get_random_unoptimized_block();
  /*Set the block.*/
  toggle_block(*current_block, privdat);
  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 digital_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
       *\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* digital_alloy_optimizer_1d::do_optimization(double (*get_value)(void* privdat, unsigned assignment, long unsigned iteration), void* get_value_privdat, void* set_point_privdat) {
  long unsigned iteration=0;
  double value = get_value(get_value_privdat, 0, iteration);
  if(unlikely(isnan(value))) {
    return g_error_new(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_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(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_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(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_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(DIGITAL_ALLOY_OPTIMIZER_1D_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_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)));
  return NULL;
      }

    }
  }
}