Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/04/16 09:47:19 (8 years ago)
Author:
bwerth
Message:

#1087 moved project bugfixes

Location:
branches/HeuristicLab.Problems.MultiObjectiveTestFunctions/HeuristicLab.Problems.MultiObjectiveTestFunctions/3.3/Calculators
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/HeuristicLab.Problems.MultiObjectiveTestFunctions/HeuristicLab.Problems.MultiObjectiveTestFunctions/3.3/Calculators/Crowding.cs

    r13936 r13988  
    2424using System.Linq;
    2525
    26 namespace HeuristicLab.Problems.MultiObjectiveTestFunctions
    27 {
     26namespace HeuristicLab.Problems.MultiObjectiveTestFunctions {
    2827
    29     /// <summary>
    30     /// Crowding distance d(x,A) is usually defined between a point x and a set of points A
    31     /// d(x,A) is then a weighted sum over all dimensions where for each dimension the next larger and the next smaller Point to x are subtracted
    32     /// I extended the concept and defined the Crowding distance of a front A as the mean of the crowding distances of every point x in A
    33     /// C(A) = mean(d(x,A)) where x in A  and d(x,A) is not infinite
    34     /// </summary>
    35     public class Crowding
    36     {
     28  /// <summary>
     29  /// Crowding distance d(x,A) is usually defined between a point x and a set of points A
     30  /// d(x,A) is then a weighted sum over all dimensions where for each dimension the next larger and the next smaller Point to x are subtracted
     31  /// I extended the concept and defined the Crowding distance of a front A as the mean of the crowding distances of every point x in A
     32  /// C(A) = mean(d(x,A)) where x in A  and d(x,A) is not infinite
     33  /// </summary>
     34  public class Crowding {
    3735
    38         public static double Calculate(IEnumerable<double[]> front, double[,] bounds)
    39         {
    40             double sum = 0;
    41             int c = 0;
    42             double[] pointsums = new double[front.Count()];
     36    public static double Calculate(IEnumerable<double[]> front, double[,] bounds) {
     37      double sum = 0;
     38      int c = 0;
     39      double[] pointsums = new double[front.Count()];
    4340
    44             for (int dim = 0; dim < front.First().Length; dim++)
    45             {
    46                 double[] arr = front.Select(x => x[dim]).ToArray();
    47                 Array.Sort(arr);
    48                 double fmax = bounds[dim % bounds.GetLength(0), 1];
    49                 double fmin = bounds[dim % bounds.GetLength(0), 0];
    50                 int pointIdx = 0;
    51                 foreach (double[] point in front)
    52                 {
    53                     double d = 0;
    54                     int pos = Array.BinarySearch(arr, point[dim]);
    55                     if (pos != 0 && pos != arr.Count() - 1)
    56                     {
    57                         d = (arr[pos + 1] - arr[pos - 1]) / (fmax - fmin);
    58                         pointsums[pointIdx++] += d;
    59                     }
    60                 }
    61             }
     41      for (int dim = 0; dim < front.First().Length; dim++) {
     42        double[] arr = front.Select(x => x[dim]).ToArray();
     43        Array.Sort(arr);
     44        double fmax = bounds[dim % bounds.GetLength(0), 1];
     45        double fmin = bounds[dim % bounds.GetLength(0), 0];
     46        int pointIdx = 0;
     47        foreach (double[] point in front) {
     48          double d = 0;
     49          int pos = Array.BinarySearch(arr, point[dim]);
     50          if (pos != 0 && pos != arr.Count() - 1) {
     51            d = (arr[pos + 1] - arr[pos - 1]) / (fmax - fmin);
     52            pointsums[pointIdx++] += d;
     53          }
     54        }
     55      }
    6256
    63             foreach (double d in pointsums)
    64             {
    65                 if (!Double.IsInfinity(d))
    66                 {
    67                     sum += d;
    68                     c++;
    69                 }
    70             }
    71             return c == 0 ? Double.PositiveInfinity : sum / c;
     57      foreach (double d in pointsums) {
     58        if (!Double.IsInfinity(d)) {
     59          sum += d;
     60          c++;
    7261        }
     62      }
     63      return c == 0 ? Double.PositiveInfinity : sum / c;
     64    }
    7365
    74     }
     66  }
    7567}

  • branches/HeuristicLab.Problems.MultiObjectiveTestFunctions/HeuristicLab.Problems.MultiObjectiveTestFunctions/3.3/Calculators/HyperVolume.cs

    r13936 r13988  
    2424using System.Linq;
    2525
    26 namespace HeuristicLab.Problems.MultiObjectiveTestFunctions
    27 {
    28 
    29     public class Hypervolume
    30     {
    31 
    32         /// <summary>
    33         /// The Hyprevolume-metric is defined as the Hypervolume enclosed between a given reference point,
    34         /// that is fixed for every evaluation function and the evaluated front.
    35         ///
    36         /// Example:
    37         /// r is the reference Point at (1|1)  and every Point p is part of the evaluated front
    38         /// The filled Area labled HV is the 2 diensional Hypervolume enclosed by this front.
    39         ///
    40         /// (0|1)                (1|1)
    41         ///   +      +-------------r
    42         ///   |      |###### HV ###|
    43         ///   |      p------+######|
    44         ///   |             p+#####|
    45         ///   |              |#####|
    46         ///   |              p-+###|
    47         ///   |                p---+
    48         ///   |                 
    49         ///   +--------------------1
    50         /// (0|0)                (1|0)
    51         ///
    52         ///  Please note that in this example both dimensions are minimized. The reference Point need to be dominated by EVERY point in the evaluated front
    53         ///
    54         /// </summary>
    55         ///
    56         public static double Calculate(IEnumerable<double[]> points, IEnumerable<double> reference, bool[] maximization)
    57         {
    58             var front = NonDominatedSelect.removeNonReferenceDominatingVectors(points, reference.ToArray<double>(), maximization, false);
    59             if (maximization.Length == 2)
    60             { //Hypervolume analysis only with 2 objectives for now
    61                 return Hypervolume.Calculate2D(front, reference, maximization);
    62             }
    63             else if (Array.TrueForAll(maximization, x => !x))
    64             {
    65                 return Hypervolume.CalculateMD(front, reference);
    66             }
    67             else
    68             {
    69                 throw new NotImplementedException("Hypervolume calculation for more than two dimensions is supported only with minimization problems");
    70             }
    71         }
    72 
    73 
    74         private static double Calculate2D(IEnumerable<double[]> front, IEnumerable<double> reference, bool[] maximization)
    75         {
    76             List<double> list = new List<double>();
    77             foreach (double d in reference) list.Add(d);
    78             double[] refp = list.ToArray<double>();
    79             if (front == null) throw new ArgumentException("Fronts must not be null");
    80             double[][] set = front.ToArray();   //Still no Good
    81             if (set.Length == 0) throw new ArgumentException("Fronts must not be empty");
    82             if (refp.Length != set[0].Length) throw new ArgumentException("Front and referencepoint need to be of the same dimensionality");
    83             Array.Sort<double[]>(set, Utilities.getDimensionComparer(0, maximization[0]));
    84             double[] last = set[set.Length - 1];
    85             CheckConsistency(last, 0, refp, maximization);
    86             CheckConsistency(last, 1, refp, maximization);
    87 
    88             double sum = 0;
    89             for (int i = 0; i < set.Length - 1; i++)
    90             {
    91                 CheckConsistency(set[i], 1, refp, maximization);
    92                 sum += Math.Abs((set[i][0] - set[i + 1][0])) * Math.Abs((set[i][1] - refp[1]));
    93             }
    94 
    95             sum += Math.Abs(refp[0] - last[0]) * Math.Abs(refp[1] - last[1]);
    96             return sum;
    97         }
    98 
    99         private static void CheckConsistency(double[] point, int dim, double[] reference, bool[] maximization)
    100         {
    101             if (!maximization[dim] && point[dim] > reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
    102             if (maximization[dim] && point[dim] < reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
    103             if (point.Length != 2) throw new ArgumentException("Only 2-dimensional cases are supported yet");
    104         }
    105 
    106         private static double CalculateMD(IEnumerable<double[]> points, IEnumerable<double> reference)
    107         {
    108             double[] referencePoint = reference.ToArray<double>();
    109             if (referencePoint == null || referencePoint.Length < 3) throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
    110             if (!IsDominated(referencePoint, points))
    111             {
    112                 throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
    113             }
    114             int objectives = referencePoint.Length;
    115             List<double[]> lpoints = new List<double[]>();
    116             foreach (double[] p in points)
    117             {
    118                 lpoints.Add(p);
    119             }
    120             lpoints.StableSort(Utilities.getDimensionComparer(objectives - 1, false));
    121 
    122             double[] regLow = new double[objectives];
    123             for (int i = 0; i < objectives; i++)
    124             {
    125                 regLow[i] = 1E15;
    126             }
    127             foreach (double[] p in lpoints)
    128             {
    129                 for (int i = 0; i < regLow.Length; i++)
    130                 {
    131                     if (regLow[i] > p[i]) regLow[i] = p[i];
    132                 }
    133             }
    134 
    135             return Stream(regLow, referencePoint, lpoints, 0, referencePoint[objectives - 1], (int)Math.Sqrt(points.Count()), objectives);
    136         }
    137 
    138         private static bool IsDominated(double[] referencePoint, IEnumerable<double[]> points)
    139         {
    140             foreach (double[] point in points)
    141             {
    142                 for (int i = 0; i < referencePoint.Length; i++)
    143                 {
    144                     if (referencePoint[i] < point[i])
    145                     {
    146                         return false;
    147                     }
    148                 }
    149             }
    150             return true;
    151         }
    152 
    153         private static double Stream(double[] regionLow, double[] regionUp, List<double[]> points, int split, double cover, int sqrtNoPoints, int objectives)
    154         {
    155             double coverOld = cover;
    156             int coverIndex = 0;
    157             int coverIndexOld = -1;
    158             int c;
    159             double result = 0;
    160 
    161             double dMeasure = GetMeasure(regionLow, regionUp, objectives);
    162             while (cover == coverOld && coverIndex < points.Count())
    163             {
    164                 if (coverIndexOld == coverIndex) break;
    165                 coverIndexOld = coverIndex;
    166                 if (Covers(points[coverIndex], regionLow, objectives))
    167                 {
    168                     cover = points[coverIndex][objectives - 1];
    169                     result += dMeasure * (coverOld - cover);
    170                 }
    171                 else coverIndex++;
    172 
    173             }
    174 
    175             for (c = coverIndex; c > 0; c--) if (points[c - 1][objectives - 1] == cover) coverIndex--;
    176             if (coverIndex == 0) return result;
    177 
    178             bool allPiles = true;
    179             int[] piles = new int[coverIndex];
    180             for (int i = 0; i < coverIndex; i++)
    181             {
    182                 piles[i] = IsPile(points[i], regionLow, regionUp, objectives);
    183                 if (piles[i] == -1)
    184                 {
    185                     allPiles = false;
    186                     break;
    187                 }
    188             }
    189 
    190             if (allPiles)
    191             {
    192                 double[] trellis = new double[regionUp.Length];
    193                 for (int j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];
    194                 double current = 0;
    195                 double next = 0;
    196                 int i = 0;
    197                 do
    198                 {
    199                     current = points[i][objectives - 1];
    200                     do
    201                     {
    202                         if (points[i][piles[i]] < trellis[piles[i]]) trellis[piles[i]] = points[i][piles[i]];
    203                         i++;
    204                         if (i < coverIndex) next = points[i][objectives - 1];
    205                         else { next = cover; break; }
    206                     } while (next == current);
    207                     result += ComputeTrellis(regionLow, regionUp, trellis, objectives) * (next - current);
    208                 } while (next != cover);
    209             }
    210             else
    211             {
    212                 double bound = -1;
    213                 double[] boundaries = new double[coverIndex];
    214                 double[] noBoundaries = new double[coverIndex];
    215                 int boundIdx = 0;
    216                 int noBoundIdx = 0;
    217 
    218                 do
    219                 {
    220                     for (int i = 0; i < coverIndex; i++)
    221                     {
    222                         int contained = ContainesBoundary(points[i], regionLow, split);
    223                         if (contained == 0) boundaries[boundIdx++] = points[i][split];
    224                         else if (contained == 1) noBoundaries[noBoundIdx++] = points[i][split];
    225                     }
    226                     if (boundIdx > 0) bound = GetMedian(boundaries, boundIdx);
    227                     else if (noBoundIdx > sqrtNoPoints) bound = GetMedian(noBoundaries, noBoundIdx);
    228                     else split++;
    229                 } while (bound == -1.0);
    230 
    231                 List<double[]> pointsChildLow, pointsChildUp;
    232                 pointsChildLow = new List<double[]>();
    233                 pointsChildUp = new List<double[]>();
    234                 double[] regionUpC = new double[regionUp.Length];
    235                 for (int j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];
    236                 double[] regionLowC = new double[regionLow.Length];
    237                 for (int j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];
    238 
    239                 for (int i = 0; i < coverIndex; i++)
    240                 {
    241                     if (PartCovers(points[i], regionUpC, objectives)) pointsChildUp.Add(points[i]);
    242                     if (PartCovers(points[i], regionUp, objectives)) pointsChildLow.Add(points[i]);
    243                 }
    244                 //this could/should be done in Parallel
    245 
    246                 if (pointsChildUp.Count() > 0) result += Stream(regionLow, regionUpC, pointsChildUp, split, cover, sqrtNoPoints, objectives);
    247                 if (pointsChildLow.Count() > 0) result += Stream(regionLowC, regionUp, pointsChildLow, split, cover, sqrtNoPoints, objectives);
    248             }
    249             return result;
    250         }
    251 
    252         private static double GetMedian(double[] vector, int length)
    253         {
    254             if (vector.Length != length)
    255             {
    256                 double[] vec = new double[length];
    257                 Array.Copy(vector, vec, length);
    258                 vector = vec;
    259             }
    260             return vector.Median();
    261 
    262         }
    263 
    264         private static double ComputeTrellis(double[] regionLow, double[] regionUp, double[] trellis, int objectives)
    265         {
    266             bool[] bs = new bool[objectives - 1];
    267             for (int i = 0; i < bs.Length; i++) bs[i] = true;
    268 
    269             double result = 0;
    270             uint noSummands = BinarayToInt(bs);
    271             int oneCounter; double summand;
    272             for (uint i = 1; i <= noSummands; i++)
    273             {
    274                 summand = 1;
    275                 IntToBinary(i, bs);
    276                 oneCounter = 0;
    277                 for (int j = 0; j < objectives - 1; j++)
    278                 {
    279                     if (bs[j])
    280                     {
    281                         summand *= regionUp[j] - trellis[j];
    282                         oneCounter++;
    283                     }
    284                     else
    285                     {
    286                         summand *= regionUp[j] - regionLow[j];
    287                     }
    288                 }
    289                 if (oneCounter % 2 == 0) result -= summand;
    290                 else result += summand;
    291 
    292             }
    293             return result;
    294         }
    295 
    296         private static void IntToBinary(uint i, bool[] bs)
    297         {
    298             for (int j = 0; j < bs.Length; j++) bs[j] = false;
    299             uint rest = i;
    300             int idx = 0;
    301             while (rest != 0)
    302             {
    303                 bs[idx] = rest % 2 == 1;
    304                 rest = rest / 2;
    305                 idx++;
    306             }
    307 
    308         }
    309 
    310         private static uint BinarayToInt(bool[] bs)
    311         {
    312             uint result = 0;
    313             for (int i = 0; i < bs.Length; i++)
    314             {
    315                 result += bs[i] ? ((uint)1 << i) : 0;
    316             }
    317             return result;
    318         }
    319 
    320         private static int IsPile(double[] cuboid, double[] regionLow, double[] regionUp, int objectives)
    321         {
    322             int pile = cuboid.Length;
    323             for (int i = 0; i < objectives - 1; i++)
    324             {
    325                 if (cuboid[i] > regionLow[i])
    326                 {
    327                     if (pile != objectives) return 1;
    328                     pile = i;
    329                 }
    330             }
    331             return pile;
    332         }
    333 
    334         private static double GetMeasure(double[] regionLow, double[] regionUp, int objectives)
    335         {
    336             double volume = 1;
    337             for (int i = 0; i < objectives - 1; i++)
    338             {
    339                 volume *= (regionUp[i] - regionLow[i]);
    340             }
    341             return volume;
    342         }
    343 
    344         private static int ContainesBoundary(double[] cub, double[] regionLow, int split)
    345         {
    346             if (regionLow[split] >= cub[split]) return -1;
    347             else
    348             {
    349                 for (int j = 0; j < split; j++)
    350                 {
    351                     if (regionLow[j] < cub[j]) return 1;
    352                 }
    353             }
    354             return 0;
    355         }
    356 
    357         private static bool PartCovers(double[] v, double[] regionUp, int objectives)
    358         {
    359             for (int i = 0; i < objectives - 1; i++)
    360             {
    361                 if (v[i] >= regionUp[i]) return false;
    362             }
    363             return true;
    364         }
    365 
    366         private static bool Covers(double[] v, double[] regionLow, int objectives)
    367         {
    368             for (int i = 0; i < objectives - 1; i++)
    369             {
    370                 if (v[i] > regionLow[i]) return false;
    371             }
    372             return true;
    373         }
    374 
    375 
    376     }
     26namespace HeuristicLab.Problems.MultiObjectiveTestFunctions {
     27
     28  public class Hypervolume {
     29
     30    /// <summary>
     31    /// The Hyprevolume-metric is defined as the Hypervolume enclosed between a given reference point,
     32    /// that is fixed for every evaluation function and the evaluated front.
     33    ///
     34    /// Example:
     35    /// r is the reference Point at (1|1)  and every Point p is part of the evaluated front
     36    /// The filled Area labled HV is the 2 diensional Hypervolume enclosed by this front.
     37    ///
     38    /// (0|1)                (1|1)
     39    ///   +      +-------------r
     40    ///   |      |###### HV ###|
     41    ///   |      p------+######|
     42    ///   |             p+#####|
     43    ///   |              |#####|
     44    ///   |              p-+###|
     45    ///   |                p---+
     46    ///   |                 
     47    ///   +--------------------1
     48    /// (0|0)                (1|0)
     49    ///
     50    ///  Please note that in this example both dimensions are minimized. The reference Point need to be dominated by EVERY point in the evaluated front
     51    ///
     52    /// </summary>
     53    ///
     54    public static double Calculate(IEnumerable<double[]> points, IEnumerable<double> reference, bool[] maximization) {
     55      var front = NonDominatedSelect.removeNonReferenceDominatingVectors(points, reference.ToArray<double>(), maximization, false);
     56      if (maximization.Length == 2) { //Hypervolume analysis only with 2 objectives for now
     57        return Hypervolume.Calculate2D(front, reference, maximization);
     58      } else if (Array.TrueForAll(maximization, x => !x)) {
     59        return Hypervolume.CalculateMD(front, reference);
     60      } else {
     61        throw new NotImplementedException("Hypervolume calculation for more than two dimensions is supported only with minimization problems");
     62      }
     63    }
     64
     65
     66    private static double Calculate2D(IEnumerable<double[]> front, IEnumerable<double> reference, bool[] maximization) {
     67      List<double> list = new List<double>();
     68      foreach (double d in reference) list.Add(d);
     69      double[] refp = list.ToArray<double>();
     70      if (front == null) throw new ArgumentException("Fronts must not be null");
     71      double[][] set = front.ToArray();   //Still no Good
     72      if (set.Length == 0) throw new ArgumentException("Fronts must not be empty");
     73      if (refp.Length != set[0].Length) throw new ArgumentException("Front and referencepoint need to be of the same dimensionality");
     74      Array.Sort<double[]>(set, Utilities.getDimensionComparer(0, maximization[0]));
     75      double[] last = set[set.Length - 1];
     76      CheckConsistency(last, 0, refp, maximization);
     77      CheckConsistency(last, 1, refp, maximization);
     78
     79      double sum = 0;
     80      for (int i = 0; i < set.Length - 1; i++) {
     81        CheckConsistency(set[i], 1, refp, maximization);
     82        sum += Math.Abs((set[i][0] - set[i + 1][0])) * Math.Abs((set[i][1] - refp[1]));
     83      }
     84
     85      sum += Math.Abs(refp[0] - last[0]) * Math.Abs(refp[1] - last[1]);
     86      return sum;
     87    }
     88
     89    private static void CheckConsistency(double[] point, int dim, double[] reference, bool[] maximization) {
     90      if (!maximization[dim] && point[dim] > reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
     91      if (maximization[dim] && point[dim] < reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
     92      if (point.Length != 2) throw new ArgumentException("Only 2-dimensional cases are supported yet");
     93    }
     94
     95    private static double CalculateMD(IEnumerable<double[]> points, IEnumerable<double> reference) {
     96      double[] referencePoint = reference.ToArray<double>();
     97      if (referencePoint == null || referencePoint.Length < 3) throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
     98      if (!IsDominated(referencePoint, points)) {
     99        throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
     100      }
     101      int objectives = referencePoint.Length;
     102      List<double[]> lpoints = new List<double[]>();
     103      foreach (double[] p in points) {
     104        lpoints.Add(p);
     105      }
     106      lpoints.StableSort(Utilities.getDimensionComparer(objectives - 1, false));
     107
     108      double[] regLow = new double[objectives];
     109      for (int i = 0; i < objectives; i++) {
     110        regLow[i] = 1E15;
     111      }
     112      foreach (double[] p in lpoints) {
     113        for (int i = 0; i < regLow.Length; i++) {
     114          if (regLow[i] > p[i]) regLow[i] = p[i];
     115        }
     116      }
     117
     118      return Stream(regLow, referencePoint, lpoints, 0, referencePoint[objectives - 1], (int)Math.Sqrt(points.Count()), objectives);
     119    }
     120
     121    private static bool IsDominated(double[] referencePoint, IEnumerable<double[]> points) {
     122      foreach (double[] point in points) {
     123        for (int i = 0; i < referencePoint.Length; i++) {
     124          if (referencePoint[i] < point[i]) {
     125            return false;
     126          }
     127        }
     128      }
     129      return true;
     130    }
     131
     132    private static double Stream(double[] regionLow, double[] regionUp, List<double[]> points, int split, double cover, int sqrtNoPoints, int objectives) {
     133      double coverOld = cover;
     134      int coverIndex = 0;
     135      int coverIndexOld = -1;
     136      int c;
     137      double result = 0;
     138
     139      double dMeasure = GetMeasure(regionLow, regionUp, objectives);
     140      while (cover == coverOld && coverIndex < points.Count()) {
     141        if (coverIndexOld == coverIndex) break;
     142        coverIndexOld = coverIndex;
     143        if (Covers(points[coverIndex], regionLow, objectives)) {
     144          cover = points[coverIndex][objectives - 1];
     145          result += dMeasure * (coverOld - cover);
     146        } else coverIndex++;
     147
     148      }
     149
     150      for (c = coverIndex; c > 0; c--) if (points[c - 1][objectives - 1] == cover) coverIndex--;
     151      if (coverIndex == 0) return result;
     152
     153      bool allPiles = true;
     154      int[] piles = new int[coverIndex];
     155      for (int i = 0; i < coverIndex; i++) {
     156        piles[i] = IsPile(points[i], regionLow, regionUp, objectives);
     157        if (piles[i] == -1) {
     158          allPiles = false;
     159          break;
     160        }
     161      }
     162
     163      if (allPiles) {
     164        double[] trellis = new double[regionUp.Length];
     165        for (int j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];
     166        double current = 0;
     167        double next = 0;
     168        int i = 0;
     169        do {
     170          current = points[i][objectives - 1];
     171          do {
     172            if (points[i][piles[i]] < trellis[piles[i]]) trellis[piles[i]] = points[i][piles[i]];
     173            i++;
     174            if (i < coverIndex) next = points[i][objectives - 1];
     175            else { next = cover; break; }
     176          } while (next == current);
     177          result += ComputeTrellis(regionLow, regionUp, trellis, objectives) * (next - current);
     178        } while (next != cover);
     179      } else {
     180        double bound = -1;
     181        double[] boundaries = new double[coverIndex];
     182        double[] noBoundaries = new double[coverIndex];
     183        int boundIdx = 0;
     184        int noBoundIdx = 0;
     185
     186        do {
     187          for (int i = 0; i < coverIndex; i++) {
     188            int contained = ContainesBoundary(points[i], regionLow, split);
     189            if (contained == 0) boundaries[boundIdx++] = points[i][split];
     190            else if (contained == 1) noBoundaries[noBoundIdx++] = points[i][split];
     191          }
     192          if (boundIdx > 0) bound = GetMedian(boundaries, boundIdx);
     193          else if (noBoundIdx > sqrtNoPoints) bound = GetMedian(noBoundaries, noBoundIdx);
     194          else split++;
     195        } while (bound == -1.0);
     196
     197        List<double[]> pointsChildLow, pointsChildUp;
     198        pointsChildLow = new List<double[]>();
     199        pointsChildUp = new List<double[]>();
     200        double[] regionUpC = new double[regionUp.Length];
     201        for (int j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];
     202        double[] regionLowC = new double[regionLow.Length];
     203        for (int j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];
     204
     205        for (int i = 0; i < coverIndex; i++) {
     206          if (PartCovers(points[i], regionUpC, objectives)) pointsChildUp.Add(points[i]);
     207          if (PartCovers(points[i], regionUp, objectives)) pointsChildLow.Add(points[i]);
     208        }
     209        //this could/should be done in Parallel
     210
     211        if (pointsChildUp.Count() > 0) result += Stream(regionLow, regionUpC, pointsChildUp, split, cover, sqrtNoPoints, objectives);
     212        if (pointsChildLow.Count() > 0) result += Stream(regionLowC, regionUp, pointsChildLow, split, cover, sqrtNoPoints, objectives);
     213      }
     214      return result;
     215    }
     216
     217    private static double GetMedian(double[] vector, int length) {
     218      if (vector.Length != length) {
     219        double[] vec = new double[length];
     220        Array.Copy(vector, vec, length);
     221        vector = vec;
     222      }
     223      return vector.Median();
     224
     225    }
     226
     227    private static double ComputeTrellis(double[] regionLow, double[] regionUp, double[] trellis, int objectives) {
     228      bool[] bs = new bool[objectives - 1];
     229      for (int i = 0; i < bs.Length; i++) bs[i] = true;
     230
     231      double result = 0;
     232      uint noSummands = BinarayToInt(bs);
     233      int oneCounter; double summand;
     234      for (uint i = 1; i <= noSummands; i++) {
     235        summand = 1;
     236        IntToBinary(i, bs);
     237        oneCounter = 0;
     238        for (int j = 0; j < objectives - 1; j++) {
     239          if (bs[j]) {
     240            summand *= regionUp[j] - trellis[j];
     241            oneCounter++;
     242          } else {
     243            summand *= regionUp[j] - regionLow[j];
     244          }
     245        }
     246        if (oneCounter % 2 == 0) result -= summand;
     247        else result += summand;
     248
     249      }
     250      return result;
     251    }
     252
     253    private static void IntToBinary(uint i, bool[] bs) {
     254      for (int j = 0; j < bs.Length; j++) bs[j] = false;
     255      uint rest = i;
     256      int idx = 0;
     257      while (rest != 0) {
     258        bs[idx] = rest % 2 == 1;
     259        rest = rest / 2;
     260        idx++;
     261      }
     262
     263    }
     264
     265    private static uint BinarayToInt(bool[] bs) {
     266      uint result = 0;
     267      for (int i = 0; i < bs.Length; i++) {
     268        result += bs[i] ? ((uint)1 << i) : 0;
     269      }
     270      return result;
     271    }
     272
     273    private static int IsPile(double[] cuboid, double[] regionLow, double[] regionUp, int objectives) {
     274      int pile = cuboid.Length;
     275      for (int i = 0; i < objectives - 1; i++) {
     276        if (cuboid[i] > regionLow[i]) {
     277          if (pile != objectives) return 1;
     278          pile = i;
     279        }
     280      }
     281      return pile;
     282    }
     283
     284    private static double GetMeasure(double[] regionLow, double[] regionUp, int objectives) {
     285      double volume = 1;
     286      for (int i = 0; i < objectives - 1; i++) {
     287        volume *= (regionUp[i] - regionLow[i]);
     288      }
     289      return volume;
     290    }
     291
     292    private static int ContainesBoundary(double[] cub, double[] regionLow, int split) {
     293      if (regionLow[split] >= cub[split]) return -1;
     294      else {
     295        for (int j = 0; j < split; j++) {
     296          if (regionLow[j] < cub[j]) return 1;
     297        }
     298      }
     299      return 0;
     300    }
     301
     302    private static bool PartCovers(double[] v, double[] regionUp, int objectives) {
     303      for (int i = 0; i < objectives - 1; i++) {
     304        if (v[i] >= regionUp[i]) return false;
     305      }
     306      return true;
     307    }
     308
     309    private static bool Covers(double[] v, double[] regionLow, int objectives) {
     310      for (int i = 0; i < objectives - 1; i++) {
     311        if (v[i] > regionLow[i]) return false;
     312      }
     313      return true;
     314    }
     315
     316
     317  }
    377318}
    378319
Note: See TracChangeset for help on using the changeset viewer.