MeshKit  1.0
IAWeights.cpp
Go to the documentation of this file.
00001 // IAWeights.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 // Weights
00005 
00006 #include "meshkit/IAWeights.hpp"
00007 
00008 #include <math.h>
00009 #include <limits>
00010 #include <algorithm>
00011 #include <assert.h>
00012 #include <cstdio>
00013 
00014 namespace MeshKit 
00015 {
00016 
00017 
00018 // generate integers from 0..n-1
00019 struct c_unique {
00020   int current;
00021   c_unique() {current=0;}
00022   int operator()() {return current++;}
00023 } UniqueNumber;
00024 
00025 
00026 class FabsWeightComparer
00027 {
00028 public: 
00029   std::vector<double> *w;
00030   FabsWeightComparer(std::vector<double> *weightvec) {w = weightvec;}
00031   bool operator() (const int lhs, const int rhs)
00032   {
00033     return fabs((*w)[lhs]) < fabs((*w)[rhs]);
00034   }
00035 };
00036 
00037 IAWeights::IAWeights() : std::vector<double>(), 
00038 debugging(true)  
00039 {}
00040 
00041 IAWeights::~IAWeights() 
00042 {}
00043   
00044   // static
00045 double IAWeights::rand_excluded_middle()
00046 {
00047   // generate double in [-1,-0.5] U [.5,1]
00048   double d = ((double) rand() / RAND_MAX) - 0.5;
00049   if (d<0.)
00050     d -= 0.5;
00051   else
00052     d += 0.5;
00053   assert( d >= -1. );
00054   assert( d <= 1. );
00055   assert( d <= -0.5 || d >= 0.5 );
00056   return d;
00057 }
00058 
00059 void IAWeights::uniquify(const double lo, const double hi)
00060 {
00061   assert( hi >= lo );
00062   assert( lo >= 0. );
00063 
00064   // find min an max of input
00065   double fabs_min_weight = std::numeric_limits<double>::max();
00066   double fabs_max_weight = 0.;
00067   for (unsigned int i = 0; i < size(); ++i)
00068   {
00069     const double w = (*this)[i];
00070     const double fabsw = fabs(w);
00071     if (fabsw < fabs_min_weight)
00072       fabs_min_weight = fabsw;
00073     if (fabsw > fabs_max_weight)
00074       fabs_max_weight = fabsw;
00075   }
00076   
00077   // relative range of input and output 
00078   const double input_range = fabs_max_weight - fabs_min_weight; 
00079   assert( input_range >= 0.);
00080   const double output_range = hi - lo;
00081   assert( output_range >= 0.);
00082   
00083   // scale the weightvec so | max | is 1.e4, and min is 1
00084   // the range should be well below the ipopt solver tolerance, which is 1.0e-7
00085   // "typically" the raw weights are between 1.e-4 and 10
00086   if (fabs_max_weight < 1.)
00087     fabs_max_weight = 1.;
00088   //was const double s = 1.e4 / fabs_max_weight;
00089   double s = output_range / input_range; // could be nan, so limit in next line
00090   if ( s > 1.e8 )
00091     s = 1.e8;
00092   for (unsigned int i = 0; i < size(); ++i)
00093   {
00094     const double fabsw = lo + ( (fabs((*this)[i]) - fabs_min_weight) * s );
00095     (*this)[i] = (*this)[i] > 0. ? fabsw : -fabsw;
00096 
00097     /* was
00098     weightvec[i] *= s;
00099     if (fabs(weightvec[i]) < 1.)
00100       if (weightvec[i] < 0.)
00101         weightvec[i] = -1.;
00102       else
00103         weightvec[i] = 1.;
00104      */
00105   }
00106   
00107   // uniquify the weights. 
00108   // ensure a random minimum ratio between consecutive weights
00109   // we'd really like no weight to be the sum of other weights, but randomization should catch most of the cases.
00110   // We randomize rather than make a deterministict fraction.
00111   
00112   // get the indices of the sorted order of the weights
00113   std::vector<int> sorted_fabs_weights(size());
00114   std::generate(sorted_fabs_weights.begin(), sorted_fabs_weights.end(), UniqueNumber );
00115   std::sort( sorted_fabs_weights.begin(), sorted_fabs_weights.end(), FabsWeightComparer(this) );
00116   
00117   
00118   srand(9384757); 
00119   double prior_fw = 0.;
00120   for (unsigned int i = 0; i < size(); ++i)
00121   {
00122     // detect consecutive identical weights and modify the later one
00123     const int j = sorted_fabs_weights[i]; // index of weight in weights
00124     double w = (*this)[j];
00125     double fw = fabs(w);
00126     if (fw - prior_fw < lo * 0.01) // relative tolerance
00127     {
00128       const double eps = lo * (0.01 + 0.02 * ((double) rand() / RAND_MAX));
00129       fw = prior_fw + eps; // use prior_fw rather than fw to ensure a min gap, uniform in [0.01, 0.03]
00130       (*this)[j] = w = (w<0.) ? -fw : fw;
00131     }
00132     prior_fw = fw;
00133     //printf("%d: w_%d %10.4f\n", i, j, weightvec[j]); 
00134   }
00135   
00136   // scale again if max was exceeded
00137   if (prior_fw > hi)
00138   {
00139     // assume minimum is lo, ignoring the random bit we added
00140     const double ss = output_range / ( prior_fw - lo );
00141     for (unsigned int i = 0; i < size(); ++i)
00142     {
00143       double w = (fabs((*this)[i]) - lo) * ss + lo;
00144       assert( w >= 0. );
00145       // with roundoff, could be slightly above hi, force it
00146       if ( w > hi )
00147         w = hi;
00148       if ( w < lo )
00149         w = lo;
00150       (*this)[i] = ((*this)[i] < 0.) ? -w : w; 
00151       assert( w <= hi );
00152       assert( w >= lo );
00153     }
00154   }
00155   
00156   if (0) // debugging
00157   {
00158     printf("unique weights with fabs in [%e, %e]\n", lo, hi);
00159     for (unsigned int i = 0; i < size(); ++i)
00160     {
00161       const int j = sorted_fabs_weights[i]; // index of weight in weightvec
00162       const double w = (*this)[j];
00163       printf("%d: w_%d %10.4f\n", i, j, w); 
00164       assert( fabs(w) <= hi );
00165       assert( fabs(w) >= lo );
00166     }
00167   }
00168   // exit(1); //zzyk
00169 }
00170 
00171 // debug
00172 void IAWeights::print() const
00173 {
00174   printf("weights:\n");
00175   for (unsigned int i = 0; i < size(); ++i)
00176   {
00177     const double w = (*this)[i];
00178     printf("w[%u] = %f\n", i, w);
00179   }
00180   printf("\n");
00181 
00182 }
00183 
00184 } // namespace MeshKit
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines