27 #ifndef CDPL_UTIL_DGCOORDINATESGENERATOR_HPP
28 #define CDPL_UTIL_DGCOORDINATESGENERATOR_HPP
35 #include <boost/random/mersenne_twister.hpp>
36 #include <boost/random/uniform_real.hpp>
37 #include <boost/random/uniform_int.hpp>
52 template <std::
size_t Dim,
typename T,
typename Derived>
60 typedef std::vector<DistanceConstraint> DistanceConstraintList;
90 std::size_t point1Idx;
91 std::size_t point2Idx;
136 template <
typename CoordsArray>
137 void generate(std::size_t num_points, CoordsArray& coords);
139 template <
typename CoordsArray>
154 std::size_t getNumVolumeConstraints()
const;
156 template <
typename CoordsArray>
157 void embedCoords(std::size_t num_points, CoordsArray& coords);
159 template <
typename CoordsArray>
160 void adjCoordsForDistanceConstraint(CoordsArray& coords,
const ValueType& lambda, std::size_t constr_idx)
const;
162 template <
typename Vec>
165 template <
typename CoordsArray>
166 void adjCoordsForVolumeConstraint(CoordsArray& coords,
const ValueType& lambda, std::size_t constr_idx)
const;
168 template <
typename Vec>
169 static ValueType calcDiffVectorAndSquaredDist(
const Vec& pt1_pos,
const Vec& pt2_pos,
ValueType diff[]);
171 typedef boost::random::mt11213b RandNumEngine;
173 std::size_t numCycles;
174 double cycleStepCountFactor;
177 DistanceConstraintList distConstraints;
178 RandNumEngine randomEngine;
181 template <std::
size_t Dim,
typename T,
typename Derived>
184 template <std::
size_t Dim,
typename T,
typename Derived>
187 template <std::
size_t Dim,
typename T,
typename Derived>
190 template <std::
size_t Dim,
typename T,
typename Derived>
194 template <std::
size_t Dim,
typename T,
typename Derived>
202 template <std::
size_t Dim,
typename T>
210 template <
typename T>
218 class VolumeConstraint;
221 typedef std::vector<VolumeConstraint> VolumeConstraintList;
228 class VolumeConstraint
248 std::size_t point1Idx;
249 std::size_t point2Idx;
250 std::size_t point3Idx;
251 std::size_t point4Idx;
279 template <
typename CoordsArray>
283 template <
typename CoordsArray>
284 void adjCoordsForVolumeConstraint(CoordsArray& coords,
const ValueType& lambda, std::size_t constr_idx)
const;
286 template <
typename Vec>
287 void adjCoordsForConstraint(
const VolumeConstraint& constr, Vec& pt1_pos, Vec& pt2_pos, Vec& pt3_pos,
288 Vec& pt4_pos,
const ValueType& lambda)
const;
290 template <
typename Vec>
291 static void calcDiffVector(
const Vec& pt1_pos,
const Vec& pt2_pos,
ValueType diff[]);
293 VolumeConstraintList volConstraints;
304 template <std::
size_t Dim,
typename T,
typename Derived>
306 const ValueType& lb,
const ValueType& ub):
308 point2Idx(pt2_idx), lowerBound(lb), upperBound(ub)
311 template <std::
size_t Dim,
typename T,
typename Derived>
317 template <std::
size_t Dim,
typename T,
typename Derived>
323 template <std::
size_t Dim,
typename T,
typename Derived>
329 template <std::
size_t Dim,
typename T,
typename Derived>
335 template <std::
size_t Dim,
typename T,
typename Derived>
338 if (point1Idx < constr.point1Idx)
341 if (point1Idx == constr.point1Idx)
342 return (point2Idx < constr.point2Idx);
350 template <std::
size_t Dim,
typename T,
typename Derived>
352 numCycles(DEF_NUM_CYCLES), cycleStepCountFactor(DEF_CYCLE_STEP_COUNT_FACTOR), startLearningRate(DEF_START_LEARNING_RATE),
353 learningRateDecr(DEF_LEARNING_RATE_DECREMENT), randomEngine(170375)
356 template <std::
size_t Dim,
typename T,
typename Derived>
358 numCycles(gen.numCycles), cycleStepCountFactor(gen.cycleStepCountFactor),
359 startLearningRate(gen.startLearningRate), learningRateDecr(gen.learningRateDecr),
360 distConstraints(gen.distConstraints),
361 randomEngine(gen.randomEngine)
364 template <std::
size_t Dim,
typename T,
typename Derived>
371 numCycles = gen.numCycles;
372 cycleStepCountFactor = gen.cycleStepCountFactor;
373 startLearningRate = gen.startLearningRate;
374 learningRateDecr = gen.learningRateDecr;
375 distConstraints = gen.distConstraints;
376 randomEngine = gen.randomEngine;
381 template <std::
size_t Dim,
typename T,
typename Derived>
384 distConstraints.clear();
387 template <std::
size_t Dim,
typename T,
typename Derived>
389 const ValueType& lb,
const ValueType& ub)
394 template <std::
size_t Dim,
typename T,
typename Derived>
397 return distConstraints.size();
400 template <std::
size_t Dim,
typename T,
typename Derived>
404 if (idx >= distConstraints.size())
405 throw Base::IndexError(
"DGCoordinatesGeneratorBase: constraint index out of bounds");
407 return distConstraints[idx];
410 template <std::
size_t Dim,
typename T,
typename Derived>
414 if (idx >= distConstraints.size())
415 throw Base::IndexError(
"DGCoordinatesGeneratorBase: constraint index out of bounds");
417 return distConstraints[idx];
420 template <std::
size_t Dim,
typename T,
typename Derived>
423 if (idx >= distConstraints.size())
424 throw Base::IndexError(
"DGCoordinatesGeneratorBase: constraint index out of bounds");
426 distConstraints.erase(distConstraints.begin() + idx);
429 template <std::
size_t Dim,
typename T,
typename Derived>
432 if ((it - distConstraints.begin()) >= distConstraints.size())
433 throw Base::IndexError(
"DGCoordinatesGeneratorBase: constraint iterator out of bounds");
435 distConstraints.erase(it);
438 template <std::
size_t Dim,
typename T,
typename Derived>
442 return distConstraints.begin();
445 template <std::
size_t Dim,
typename T,
typename Derived>
449 return distConstraints.end();
452 template <std::
size_t Dim,
typename T,
typename Derived>
456 return distConstraints.begin();
459 template <std::
size_t Dim,
typename T,
typename Derived>
463 return distConstraints.end();
466 template <std::
size_t Dim,
typename T,
typename Derived>
469 numCycles = num_cycles;
472 template <std::
size_t Dim,
typename T,
typename Derived>
475 cycleStepCountFactor = fact;
478 template <std::
size_t Dim,
typename T,
typename Derived>
481 startLearningRate = rate;
484 template <std::
size_t Dim,
typename T,
typename Derived>
487 learningRateDecr = decr;
490 template <std::
size_t Dim,
typename T,
typename Derived>
496 template <std::
size_t Dim,
typename T,
typename Derived>
499 return cycleStepCountFactor;
502 template <std::
size_t Dim,
typename T,
typename Derived>
506 return startLearningRate;
509 template <std::
size_t Dim,
typename T,
typename Derived>
513 return learningRateDecr;
516 template <std::
size_t Dim,
typename T,
typename Derived>
519 randomEngine.seed(seed);
522 template <std::
size_t Dim,
typename T,
typename Derived>
523 template <
typename CoordsArray>
527 ValueType error = ValueType();
528 ValueType pos_diff[Dim];
530 for (
typename DistanceConstraintList::const_iterator it = distConstraints.begin(), end = distConstraints.end(); it != end; ++it) {
533 ValueType dist_2 = calcDiffVectorAndSquaredDist(coords[constr.getPoint1Index()], coords[constr.getPoint2Index()], pos_diff);
534 ValueType dist = std::sqrt(dist_2);
535 ValueType lb = constr.getLowerBound();
536 ValueType ub = constr.getUpperBound();
538 if (dist >= lb && dist <= ub)
542 ValueType tmp = (dist_2 - lb * lb) / (0.000001 + lb * lb);
547 ValueType tmp = (dist_2 - ub * ub) / (0.000001 + ub * ub);
556 template <std::
size_t Dim,
typename T,
typename Derived>
562 template <std::
size_t Dim,
typename T,
typename Derived>
565 std::sort(distConstraints.begin(), distConstraints.end());
568 template <std::
size_t Dim,
typename T,
typename Derived>
569 template <
typename CoordsArray>
572 embedCoords(num_points, coords);
575 template <std::
size_t Dim,
typename T,
typename Derived>
576 template <
typename CoordsArray>
579 std::size_t num_dist_constrs = distConstraints.size();
580 std::size_t num_vol_constrs =
static_cast<Derived&
>(*this).getNumVolumeConstraints();
582 if ((num_dist_constrs + num_vol_constrs) == 0)
585 std::size_t num_steps = std::size_t((num_dist_constrs + num_vol_constrs) * cycleStepCountFactor);
586 ValueType lambda = startLearningRate;
588 if (num_dist_constrs > 0 && num_vol_constrs > 0) {
589 boost::random::uniform_int_distribution<std::size_t> constr_sd(0, num_dist_constrs + num_vol_constrs - 1);
591 for (std::size_t i = 0; i < numCycles; i++, lambda -= learningRateDecr) {
592 for (std::size_t j = 0; j < num_steps; j++) {
593 std::size_t constr_idx = constr_sd(randomEngine);
595 if (constr_idx < num_dist_constrs)
596 adjCoordsForDistanceConstraint(coords, lambda, constr_idx);
598 static_cast<Derived&
>(*this).template adjCoordsForVolumeConstraint<CoordsArray>(coords, lambda, constr_idx - num_dist_constrs);
605 if (num_dist_constrs > 0) {
606 boost::random::uniform_int_distribution<std::size_t> constr_sd(0, num_dist_constrs - 1);
608 for (std::size_t i = 0; i < numCycles; i++, lambda -= learningRateDecr)
609 for (std::size_t j = 0; j < num_steps; j++)
610 adjCoordsForDistanceConstraint(coords, lambda, constr_sd(randomEngine));
615 boost::random::uniform_int_distribution<std::size_t> constr_sd(0, num_vol_constrs - 1);
617 for (std::size_t i = 0; i < numCycles; i++, lambda -= learningRateDecr)
618 for (std::size_t j = 0; j < num_steps; j++)
619 static_cast<Derived&
>(*this).template adjCoordsForVolumeConstraint<CoordsArray>(coords, lambda, constr_sd(randomEngine));
622 template <std::
size_t Dim,
typename T,
typename Derived>
623 template <
typename CoordsArray>
625 std::size_t constr_idx)
const
629 adjCoordsForConstraint(constr, coords[constr.getPoint1Index()], coords[constr.getPoint2Index()], lambda);
632 template <std::
size_t Dim,
typename T,
typename Derived>
633 template <
typename Vec>
635 Vec& pt2_pos,
const ValueType& lambda)
const
637 ValueType pos_diff[Dim];
638 ValueType dist = std::sqrt(calcDiffVectorAndSquaredDist(pt1_pos, pt2_pos, pos_diff));
640 ValueType ub = constr.getUpperBound();
641 ValueType lb = constr.getLowerBound();
643 if (dist >= lb && dist <= ub)
646 ValueType bound = (dist > ub ? ub : lb);
647 ValueType factor = lambda / 2 * (bound - dist) / (0.000001 + dist);
649 for (std::size_t i = 0; i < Dim; i++) {
650 ValueType pos_delta = pos_diff[i] * factor;
652 pt1_pos[i] -= pos_delta;
653 pt2_pos[i] += pos_delta;
657 template <std::
size_t Dim,
typename T,
typename Derived>
658 template <
typename CoordsArray>
660 std::size_t constr_idx)
const
663 template <std::
size_t Dim,
typename T,
typename Derived>
664 template <
typename Vec>
668 ValueType dist_2 = ValueType();
670 for (std::size_t i = 0; i < Dim; i++) {
671 diff[i] = pt2_pos[i] - pt1_pos[i];
672 dist_2 += diff[i] * diff[i];
681 template <
typename T>
683 std::size_t pt3_idx, std::size_t pt4_idx,
684 const ValueType& lb,
const ValueType& ub):
686 point2Idx(pt2_idx), point3Idx(pt3_idx), point4Idx(pt4_idx), lowerBound(lb), upperBound(ub)
689 template <
typename T>
695 template <
typename T>
701 template <
typename T>
707 template <
typename T>
713 template <
typename T>
720 template <
typename T>
730 template <
typename T>
733 volConstraints.clear();
736 template <
typename T>
738 std::size_t pt4_idx,
const ValueType& lb,
const ValueType& ub)
740 volConstraints.push_back(VolumeConstraint(pt1_idx, pt2_idx, pt3_idx, pt4_idx, lb, ub));
743 template <
typename T>
746 return volConstraints.size();
749 template <
typename T>
753 if (idx >= volConstraints.size())
754 throw Base::IndexError(
"DGCoordinatesGenerator: constraint index out of bounds");
756 return volConstraints[idx];
759 template <
typename T>
762 if (idx >= volConstraints.size())
763 throw Base::IndexError(
"DGCoordinatesGenerator: constraint index out of bounds");
765 volConstraints.erase(volConstraints.begin() + idx);
768 template <
typename T>
771 if ((it - volConstraints.begin()) >= volConstraints.size())
772 throw Base::IndexError(
"DGCoordinatesGenerator: constraint iterator out of bounds");
774 volConstraints.erase(it);
777 template <
typename T>
781 return volConstraints.begin();
784 template <
typename T>
788 return volConstraints.end();
791 template <
typename T>
795 return volConstraints.begin();
798 template <
typename T>
802 return volConstraints.end();
805 template <
typename T>
806 template <
typename CoordsArray>
810 ValueType error = ValueType();
815 for (
typename VolumeConstraintList::const_iterator it = volConstraints.begin(), end = volConstraints.end(); it != end; ++it) {
816 const VolumeConstraint& constr = *it;
818 calcDiffVector(coords[constr.getPoint4Index()], coords[constr.getPoint1Index()], v_41);
819 calcDiffVector(coords[constr.getPoint4Index()], coords[constr.getPoint2Index()], v_42);
820 calcDiffVector(coords[constr.getPoint4Index()], coords[constr.getPoint3Index()], v_43);
822 ValueType vol = (v_41[0] * (v_42[1] * v_43[2] - v_42[2] * v_43[1]) - v_41[1] * (v_42[0] * v_43[2] - v_42[2] * v_43[0]) + v_41[2] * (v_42[0] * v_43[1] - v_42[1] * v_43[0])) / 6;
824 ValueType lb = constr.getLowerBound();
825 ValueType ub = constr.getUpperBound();
827 if (vol >= lb && vol <= ub)
831 ValueType tmp = (vol - lb);
836 ValueType tmp = (vol - ub);
845 template <
typename T>
846 template <
typename CoordsArray>
849 const VolumeConstraint& constr = volConstraints[constr_idx];
851 adjCoordsForConstraint(constr, coords[constr.getPoint1Index()], coords[constr.getPoint2Index()],
852 coords[constr.getPoint3Index()], coords[constr.getPoint4Index()], lambda);
855 template <
typename T>
856 template <
typename Vec>
858 Vec& pt4_pos,
const ValueType& lambda)
const
860 Vec* pt_pos[4] = {&pt1_pos, &pt2_pos, &pt3_pos, &pt4_pos};
865 calcDiffVector(*pt_pos[3], *pt_pos[0], v_41);
866 calcDiffVector(*pt_pos[3], *pt_pos[1], v_42);
867 calcDiffVector(*pt_pos[3], *pt_pos[2], v_43);
871 g[0][0] = (v_42[1] * v_43[2] - v_42[2] * v_43[1]) / 6;
872 g[0][1] = -(v_42[0] * v_43[2] - v_42[2] * v_43[0]) / 6;
873 g[0][2] = (v_42[0] * v_43[1] - v_42[1] * v_43[0]) / 6;
875 ValueType vol = (v_41[0] * g[0][0] + v_41[1] * g[0][1] + v_41[2] * g[0][2]);
876 ValueType ub = constr.getUpperBound();
877 ValueType lb = constr.getLowerBound();
879 if (vol >= lb && vol <= ub)
882 g[1][0] = (v_41[2] * v_43[1] - v_41[1] * v_43[2]) / 6;
883 g[1][1] = (v_41[0] * v_43[2] - v_41[2] * v_43[0]) / 6;
884 g[1][2] = (v_41[1] * v_43[0] - v_41[0] * v_43[1]) / 6;
886 g[2][0] = (v_41[1] * v_42[2] - v_41[2] * v_42[1]) / 6;
887 g[2][1] = (v_41[2] * v_42[0] - v_41[0] * v_42[2]) / 6;
888 g[2][2] = (v_41[0] * v_42[1] - v_41[1] * v_42[0]) / 6;
890 g[3][0] = -g[0][0] - g[1][0] - g[2][0];
891 g[3][1] = -g[0][1] - g[1][1] - g[2][1];
892 g[3][2] = -g[0][2] - g[1][2] - g[2][2];
894 ValueType g_len2_sum = ValueType();
896 for (std::size_t i = 0; i < 4; i++)
897 g_len2_sum += g[i][0] * g[i][0] + g[i][1] * g[i][1] + g[i][2] * g[i][2];
899 ValueType bound = (vol < lb ? lb : ub);
900 ValueType fact = lambda * (bound - vol) / g_len2_sum;
902 for (std::size_t i = 0; i < 4; i++)
903 for (std::size_t j = 0; j < 3; j++)
904 (*pt_pos[i])[j] += fact * g[i][j];
907 template <
typename T>
908 template <
typename Vec>
911 for (std::size_t i = 0; i < 3; i++)
912 diff[i] = pt2_pos[i] - pt1_pos[i];
917 #endif // CDPL_UTIL_DGCOORDINATESGENERATOR_HPP