BoundaryConstraint.cpp 7.75 KB
Newer Older
1
2
3
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
4
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
5
6
#include "MantidAPI/Expression.h"
#include "MantidAPI/ConstraintFactory.h"
7
#include "MantidAPI/IFunction.h"
8
#include "MantidKernel/Logger.h"
9
10
#include <boost/lexical_cast.hpp>
#include <sstream>
11
#include <iostream>
12

13
14
namespace Mantid {
namespace CurveFitting {
15
namespace Constraints {
16
17
18
19
namespace {
/// static logger
Kernel::Logger g_log("BoundaryConstraint");
}
20

21
DECLARE_CONSTRAINT(BoundaryConstraint)
22

23
// using namespace Kernel;
24
25
using namespace API;

26
27
/// Default constructor
BoundaryConstraint::BoundaryConstraint()
28
29
    : API::IConstraint(), m_penaltyFactor(1000.0), m_hasLowerBound(false),
      m_hasUpperBound(false), m_lowerBound(DBL_MAX), m_upperBound(-DBL_MAX) {}
30
31
32
33

/// Constructor with no boundary arguments
/// @param paramName :: The parameter name
BoundaryConstraint::BoundaryConstraint(const std::string &paramName)
34
    : API::IConstraint(), m_penaltyFactor(1000.0), m_hasLowerBound(false),
35
      m_hasUpperBound(false) {UNUSED_ARG(paramName);}
36

37
/** Constructor with boundary arguments
38
39
40
41
 * @param fun :: The function
 * @param paramName :: The parameter name
 * @param lowerBound :: The lower bound
 * @param upperBound :: The upper bound
42
43
 * @param isDefault :: Flag to mark as default the value of an object associated
 * with this reference:
Jose Borreguero's avatar
Jose Borreguero committed
44
 *  a tie or a constraint.
45
 */
46
47
48
49
BoundaryConstraint::BoundaryConstraint(API::IFunction *fun,
                                       const std::string paramName,
                                       const double lowerBound,
                                       const double upperBound, bool isDefault)
50
51
    : m_penaltyFactor(1000.0), m_hasLowerBound(true), m_hasUpperBound(true),
      m_lowerBound(lowerBound), m_upperBound(upperBound) {
52
  reset(fun, fun->parameterIndex(paramName), isDefault);
53
54
}

55
56
57
BoundaryConstraint::BoundaryConstraint(API::IFunction *fun,
                                       const std::string paramName,
                                       const double lowerBound, bool isDefault)
58
59
    : m_penaltyFactor(1000.0), m_hasLowerBound(true), m_hasUpperBound(false),
      m_lowerBound(lowerBound), m_upperBound(-DBL_MAX) {
60
  reset(fun, fun->parameterIndex(paramName), isDefault);
61
62
}

63
/** Initialize the constraint from an expression.
64
65
 * @param fun :: The function
 * @param expr :: The initializing expression which must look like this:
66
67
 * " 10 < Sigma < 20 " or
 * " Sigma > 20 "
68
69
 * @param isDefault :: Flag to mark as default the value of an object associated
 * with this reference:
Jose Borreguero's avatar
Jose Borreguero committed
70
 *  a tie or a constraint.
71
 */
72
73
74
75
void BoundaryConstraint::initialize(API::IFunction *fun,
                                    const API::Expression &expr,
                                    bool isDefault) {
  if (expr.size() < 2 || expr.name() != "==") {
76
77
78
79
    g_log.error("Wrong initialization expression");
    throw std::invalid_argument("Wrong initialization expression");
  }
  clearBounds();
80
  const Expression &terms(expr);
81
82
83
84
85

  std::vector<double> values(3);
  int ilow = -1;
  int ihi = -1;
  std::string parName;
86
  for (size_t i = 0; i < terms.size(); i++) {
87
    std::string name = terms[i].str();
88
    try {
89
90
91
      double d = boost::lexical_cast<double>(name);
      values[i] = d;
      std::string op = terms[i].operator_name();
92
93
94
      if (op.empty()) {
        op = terms[i + 1].operator_name();
        if (op[0] == '<') {
95
          ilow = static_cast<int>(i);
96
        } else if (op[0] == '>') {
97
          ihi = static_cast<int>(i);
98
        } else {
99
          g_log.error("Unknown operator in initialization expression");
100
101
          throw std::invalid_argument(
              "Unknown operator in initialization expression");
102
        }
103
104
105
      } // if empty
      else {
        if (op[0] == '<') {
106
          ihi = static_cast<int>(i);
107
        } else if (op[0] == '>') {
108
          ilow = static_cast<int>(i);
109
        } else {
110
          g_log.error("Unknown operator in initialization expression");
111
112
          throw std::invalid_argument(
              "Unknown operator in initialization expression");
113
        }
114
115
116
      } // if not empty
    } catch (boost::bad_lexical_cast &) {
      if (!parName.empty()) {
117
118
119
120
121
        g_log.error("Non-numeric value for a bound");
        throw std::invalid_argument("Non-numeric value for a bound");
      }
      parName = name;
    }
122
  } // for i
123

124
  try {
125
    size_t i = fun->parameterIndex(parName);
126
127
128
129
    reset(fun, i, isDefault);
  } catch (...) {
    g_log.error() << "Parameter " << parName << " not found in function "
                  << fun->name() << '\n';
130
131
132
    throw;
  }

133
  if (ilow >= 0) {
134
135
    setLower(values[ilow]);
  }
136
  if (ihi >= 0) {
137
138
139
140
    setUpper(values[ihi]);
  }
}

141
142
/** Set penalty factor
 *
143
 *  @param c :: penalty factor
144
 */
145
146
147
148
149
void BoundaryConstraint::setPenaltyFactor(const double &c) {
  if (c <= 0.0) {
    g_log.warning()
        << "Penalty factor <= 0 selected for boundary constraint."
        << " Only positive penalty factor allowed. Penalty factor set to 1";
150
151
    m_penaltyFactor = 1;
  }
152
  { m_penaltyFactor = c; }
153
154
}

155
156
157
158
void BoundaryConstraint::setParamToSatisfyConstraint() {
  if (!(m_hasLowerBound || m_hasUpperBound)) {
    g_log.warning()
        << "No bounds have been set on BoundaryConstraint for parameter "
159
        << parameterName() << ". Therefore"
160
        << " this constraint serves no purpose!";
161
    return;
162
163
  }

164
  double paramValue = getParameter();
165
166

  if (m_hasLowerBound)
167
    if (paramValue < m_lowerBound)
168
      setParameter(m_lowerBound, false);
169
  if (m_hasUpperBound)
170
    if (paramValue > m_upperBound)
171
      setParameter(m_upperBound, false);
172
173
}

174
175
176
177
double BoundaryConstraint::check() {
  if (!(m_hasLowerBound || m_hasUpperBound)) {
    g_log.warning()
        << "No bounds have been set on BoundaryConstraint for parameter "
178
        << parameterName() << ". Therefore"
179
        << " this constraint serves no purpose!";
180
181
182
    return 0.0;
  }

183
  double paramValue = getParameter();
184
185
186

  double penalty = 0.0;

187
  if (m_hasLowerBound)
188
189
    if (paramValue < m_lowerBound) {
      double dp = m_lowerBound - paramValue;
190
191
      penalty = m_penaltyFactor * dp * dp;
    }
192
  if (m_hasUpperBound)
193
194
    if (paramValue > m_upperBound) {
      double dp = paramValue - m_upperBound;
195
196
      penalty = m_penaltyFactor * dp * dp;
    }
197
198
199
200

  return penalty;
}

201
double BoundaryConstraint::checkDeriv() {
202
203
  double penalty = 0.0;

204
205
206
207
  if (/*m_activeParameterIndex < 0 ||*/ !(m_hasLowerBound || m_hasUpperBound)) {
    // no point in logging any warning here since checkDeriv() will always be
    // called after
    // check() is called
208
    return penalty;
209
  }
210

211
  double paramValue = getParameter();
212
213

  if (m_hasLowerBound)
214
215
    if (paramValue < m_lowerBound) {
      double dp = m_lowerBound - paramValue;
216
217
218
      penalty = 2 * m_penaltyFactor * dp;
    }
  if (m_hasUpperBound)
219
220
    if (paramValue > m_upperBound) {
      double dp = paramValue - m_upperBound;
221
222
223
224
225
226
      penalty = 2 * m_penaltyFactor * dp;
    }

  return penalty;
}

227
double BoundaryConstraint::checkDeriv2() {
228
229
  double penalty = 0.0;

230
231
232
233
  if (/*m_activeParameterIndex < 0 ||*/ !(m_hasLowerBound || m_hasUpperBound)) {
    // no point in logging any warning here since checkDeriv() will always be
    // called after
    // check() is called
234
    return penalty;
235
  }
236

237
  double paramValue = getParameter();
238
239

  if (m_hasLowerBound)
240
    if (paramValue < m_lowerBound)
241
      penalty = 2 * m_penaltyFactor;
242
  if (m_hasUpperBound)
243
    if (paramValue > m_upperBound)
244
      penalty = 2 * m_penaltyFactor;
245
246
247
248

  return penalty;
}

249
std::string BoundaryConstraint::asString() const {
250
  std::ostringstream ostr;
251
  if (m_hasLowerBound) {
252
253
    ostr << m_lowerBound << '<';
  }
254
  ostr << parameterName();
255
256
  if (m_hasUpperBound) {
    ostr << '<' << m_upperBound;
257
258
259
260
  }
  return ostr.str();
}

261
} // namespace Constraints
262
263
} // namespace CurveFitting
} // namespace Mantid