digital_alloy_optimizer_1d_mpi.c++

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

*/

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

#ifndef HAVE_MPI
#error 'digital_alloy_optimizer_1d_mpi should ONLY be built if MPI support is included, but you do not have MPI support!'
#endif

#define DIGITAL_ALLOY_OPTIMIZER_1D_MPI_GERROR_DOMAIN 74

#define EPSILON 1e-9

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

namespace kdotp {
  namespace libmetacalc {
    namespace optimizer {

      digital_alloy_optimizer_1d_mpi
      ::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,
           long unsigned int rng_seed,
           GList* assignments,
           struct ::kdotp::libmodelxx::postprocessing::function_exp *assignment_expr,
           double temperature,
           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),
    N_nodes(_N_nodes),
          node(_node),
    T(temperature),
    node_assignments(NULL),
    results(NULL),
    all_results(NULL)
      {
  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) {
    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_MPI_GERROR_DOMAIN, 0, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d_mpi.c::digital_alloy_optimizer_1d_mpi(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*/
    if(node == 0) {
      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_mpi.c::digital_alloy_optimizer_1d_mpi(constructor): Failed to get random seed:\n\t");
        seed=0;
        MPI_Bcast(&seed, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
        return;
      }
      fprintf(stderr, "rng_seed was zero, so got new one (seed=%lu rng_seed=%lu)\n", seed, rng_seed);
      MPI_Bcast(&seed, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
    }else{
      int result = MPI_Bcast(&seed, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
      if((result != 0) || (seed == 0)) {
        g_propagate_prefixed_error(err, e, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d_mpi.c::digital_alloy_optimizer_1d_mpi(constructor): Failed to get random seed (either MPI error or seed(%lu)==0):\n\t", seed);
        return;
      }
    }
  }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_MPI_GERROR_DOMAIN, 5, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d_mpi.c::digital_alloy_optimizer_1d_mpi(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);
  }
  /*Set up the nodes' data*/
  node_assignments = new int[N_nodes+1];
  node_assignments[0] = 0;
  for(unsigned n=1; n<N_nodes+1; ++n) {
    node_assignments[n]=node_assignments[n-1]+N_blocks/N_nodes;
    /*we need <= because we're setting the upper limit for the *previous* node!*/
    if(n <= N_blocks%N_nodes) node_assignments[n]++;
  }
  if(node == 0) {
    fprintf(stderr, "Assignments for %u blocks across %u nodes (%u/%u=%u; %u%%%u=%u): \n", N_blocks, N_nodes, N_blocks, N_nodes, N_blocks/N_nodes, N_blocks, N_nodes, N_blocks%N_nodes);
    for(unsigned i=0; i<N_nodes; ++i) {
      fprintf(stderr, "node %u: %u\n", i, node_assignments[i]);
    }
    fprintf(stderr, "node %u: %u\n", N_nodes, node_assignments[N_nodes]);
  }
  results = new double[node_assignments[node+1]-node_assignments[node]];
  all_results = new double[N_blocks];
  /*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_mpi::initialize(void* privdat) {
  if(current_block != NULL) {
    delete current_block;
    current_block = NULL;
  }
  randomize_array(current_assignments, N_blocks);
  set_material_from_current_assignments(privdat);
      }

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

      digital_alloy_optimizer_1d_mpi::~digital_alloy_optimizer_1d_mpi() {
  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;
  if(results != NULL) delete[] results;
  if(node_assignments != NULL) delete[] node_assignments;
  if(all_results != NULL) delete[] all_results;
      }

      int digital_alloy_optimizer_1d_mpi::do_step(void* privdat) {
  int converged = 0;
  int* offsets = new int[N_nodes];
  int* receive_counts = new int[N_nodes];
  fprintf(stderr, "assignment for %u: %u (%u-%u)\n", node, node_assignments[node+1]-node_assignments[node], node_assignments[node], node_assignments[node+1]);
  int offset = 0;
  for(unsigned i=0; i<N_nodes; ++i) {
    offsets[i] = offset;
    receive_counts[i] = node_assignments[i+1]-node_assignments[i];
    offset += receive_counts[i];
  }
  int result = MPI_Gatherv(results, node_assignments[node+1]-node_assignments[node], MPI_DOUBLE, all_results, receive_counts, offsets, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  unsigned selected_mod = 0;
  if(node == 0) {
    fprintf(stderr, "offsets\treceive_counts\n");
    for(unsigned i=0; i<N_nodes; ++i) {
      fprintf(stderr, "%d\t%d\n", offsets[i], receive_counts[i]);
    }
    for(unsigned i=0; i<N_blocks; ++i){
      fprintf(stderr, "%u %g\n", i, all_results[i]);
    }
  }
  delete[] offsets;
  delete[] receive_counts;
  if(result != 0) {
    fprintf(stderr, "Node %u: Got result %d from MPI_Gatherv!\n", node, result);
    /*If we got an MPI error there, all the nodes should've gotten the error*/
    return 5;
  }
  /*Check for convergence*/
  if(node == 0) {
    /*The mod which was selected*/
    /*A mod was selected*/
    int mod_was_selected=0;
    double sum = 0;
    /*These are in rank order, so they're in the right order; they represent the change at each block*/
    int N_improved=0;
    for(unsigned i=0; i<N_blocks; ++i) {
      double dist = dist_to_goal(all_results[i]);
      /*Improvements will be >= 0, regardless of the goal*/
      /*\note we put the n_improved incrementing up here because we need to check improvement independent of the sum (because of the annealing variant).*/
      if(dist > 0) {
        ++N_improved;
      }
      if(T != 0) {
        dist = exp(dist/T);
      }
      fprintf(stderr, "%u: sum=%g dist=%g ", i, sum, dist);
      if(dist > 0) {
        sum += dist;
        /*If this was the only one, this tells us which one it was*/
        selected_mod=i;
      }
      fprintf(stderr, " -> sum=%g\n", sum);
    }
    if(N_improved == 0) {
      converged=3;
    }else{
      if(N_improved == 1) {
        mod_was_selected=1;
        /*selected_mod was set in the loop above*/
      }else{
        double rand_selection = gsl_rng_uniform(rng)*sum;
        sum=0;
        for(unsigned i=0; i<N_blocks; ++i) {
    double dist = dist_to_goal(all_results[i]);
    if(T != 0) {
      dist = exp(dist/T);
    }
    fprintf(stderr, "%u: sum=%g dist=%g rand=%g ", i, sum, dist, rand_selection);
    if(dist > 0) {
      if((rand_selection >= sum) && (rand_selection <= (sum+dist))) {
        selected_mod = i;
        mod_was_selected=1;
        fprintf(stderr, "-> SELECTED\n");
        break;
      }
      sum += dist;
    }
    fprintf(stderr, " -> sum=%g\n", sum);
        }
        if(!mod_was_selected) {
    fprintf(stderr, "ERROR: no modification was selected!\n");
    converged=2;
        }
      }
    }
    /*Previous_value only matters on the root node; the others just calc what they're told.*/
    fprintf(stderr, "converged=%d mod_was_selected=%d selected_mod=%d previous_value=%g", converged, mod_was_selected, selected_mod, previous_value);
    previous_value = all_results[selected_mod];
    fprintf(stderr, "->%g\n", previous_value);
  }
  result = MPI_Bcast(&converged, 1, MPI_INT, 0, MPI_COMM_WORLD);
  if(result != 0) return result;
  fprintf(stderr, "node %u: converged=%d\n", node, converged);
  if(converged) {
    return converged;
  }
  /*If we've not yet converged, we need to send the modification to the nodes*/
  result = MPI_Bcast(&selected_mod, 1, MPI_INT, 0, MPI_COMM_WORLD);
  if(result != 0) {
    fprintf(stderr, "node %u: error %d bcasting the selected mod\n", node, result);
    return result;
  }
  toggle_block(selected_mod, privdat);
  return 0;
      }

      /*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 anneal_start starting teperature for a simulated annealing program
       *\param anneal_end ending temperature for a simulated annealing program
       *\param n_anneal_steps number of annealing steps per annealing cycle
       *\param n_anneal_cyles number of annealing cycles to perform. Must be at least 1 to use a simulated annealing program
       *\param num_calculations_per_annealing_step
       *\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_mpi::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) {
  long unsigned iteration=0;
  if(node == 0) {
    double value = get_value(get_value_privdat, 0, iteration);
    if(unlikely(isnan(value))) {
      return g_error_new(DIGITAL_ALLOY_OPTIMIZER_1D_MPI_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d_mpi.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_MPI_GERROR_DOMAIN, 101, "kdotp::libkdotp::optimizer::digital_alloy_optimizer_1d_mpi.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");
    }
    previous_value = value;
  }
  if(n_anneal_cycles == 0) {
    fprintf(stderr, "(iter %lu) assignment for node %u: %u (%u-%u)\n", iteration, node, node_assignments[node+1]-node_assignments[node], node_assignments[node], node_assignments[node+1]);
    do {
      ++iteration;
      for(int i=0; i<(node_assignments[node+1]-node_assignments[node]); ++i) {
        toggle_block(node_assignments[node]+i, set_point_privdat);
        results[i] = get_value(get_value_privdat, i, iteration);
        toggle_block(node_assignments[node]+i, set_point_privdat);
        fprintf(stderr, "(iter %lu) node %u: assignment %u: %g\n", iteration, node, i, results[i]);
      }
    }while(!do_step(set_point_privdat));
  }else{
    for(unsigned anneal_cycle=0; anneal_cycle < n_anneal_cycles; ++anneal_cycle) {
      T = anneal_start;
      /*We iterate over a cycle-total calculation step; temperature is set based on the annealing step*/
      for(unsigned calc_step=0; calc_step < n_anneal_steps*num_calculations_per_annealing_step; ++calc_step) {
        ++iteration;
        /*Parentheses around division ensures it's done as integer, leading to rounding down*/
        T = anneal_start + (calc_step/num_calculations_per_annealing_step)*(anneal_end-anneal_start)/(n_anneal_steps-1);
        for(int i=0; i<=(node_assignments[node+1]-node_assignments[node]); ++i) {
    toggle_block(node_assignments[node]+i, set_point_privdat);
    results[i] = get_value(get_value_privdat, i, iteration);
    toggle_block(node_assignments[node]+i, set_point_privdat);
        }
        /*If we converge before we hit the end of an annealing cycle, terminate the cycle early.*/
        if(do_step(set_point_privdat)) break;
      }
    }
  }
  return NULL;
      }
    }
  }
}