MultiDomainFunctionTest.h 17.7 KB
Newer Older
1
2
3
4
5
6
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//     NScD Oak Ridge National Laboratory, European Spallation Source
//     & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
7
#pragma once
8
9

#include "MantidAPI/FunctionDomain1D.h"
LamarMoore's avatar
LamarMoore committed
10
#include "MantidAPI/FunctionFactory.h"
11
12
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/IFunction1D.h"
LamarMoore's avatar
LamarMoore committed
13
14
#include "MantidAPI/JointDomain.h"
#include "MantidAPI/MultiDomainFunction.h"
15
16
17
#include "MantidAPI/ParamFunction.h"

#include <algorithm>
LamarMoore's avatar
LamarMoore committed
18
19
#include <boost/make_shared.hpp>
#include <cxxtest/TestSuite.h>
20
21
22
23

using namespace Mantid;
using namespace Mantid::API;

24
25
class MultiDomainFunctionTest_Function : public IFunction1D,
                                         public ParamFunction {
26
public:
27
28
29
  MultiDomainFunctionTest_Function() : IFunction1D(), ParamFunction() {
    this->declareParameter("A", 0);
    this->declareParameter("B", 0);
30
  }
31
  std::string name() const override {
32
33
34
    return "MultiDomainFunctionTest_Function";
  }

35
protected:
36
37
  void function1D(double *out, const double *xValues,
                  const size_t nData) const override {
38
39
40
    const double A = getParameter(0);
    const double B = getParameter(1);

41
    for (size_t i = 0; i < nData; ++i) {
42
43
44
45
      double x = xValues[i];
      out[i] = A + B * x;
    }
  }
46
47
  void functionDeriv1D(Jacobian *out, const double *xValues,
                       const size_t nData) override {
48
    for (size_t i = 0; i < nData; ++i) {
49
      double x = xValues[i];
50
51
      out->set(i, 0, x);
      out->set(i, 1, 1.0);
52
53
54
55
    }
  }
};

56
DECLARE_FUNCTION(MultiDomainFunctionTest_Function)
57

58
59
namespace {

60
class JacobianToTestNumDeriv : public Jacobian {
61
62
  size_t n[3];
  size_t np;
63

64
65
public:
  double off_diag;
66
  JacobianToTestNumDeriv() : np(2), off_diag(0.0) {
67
    // sizes of the three domains
68
    n[0] = 9;
69
70
71
    n[1] = 10;
    n[2] = 11;
  }
72
  void set(size_t iY, size_t iP, double value) override {
73
74
75
    // domain index the data point #iY comes from
    size_t jY = 2;
    size_t size = 0;
76
    for (size_t k = 0; k < 3; ++k) {
77
      size += n[k];
78
      if (iY < size) {
79
80
81
82
83
84
        jY = k;
        break;
      }
    }
    // domain index of function that has parameter #iP
    auto jP = iP / np;
85
86
    if (jY != jP)
      off_diag += value;
87
  }
88
  double get(size_t, size_t) override { return 0.0; }
89
  void zero() override {}
90
};
LamarMoore's avatar
LamarMoore committed
91
} // namespace
92

93
class MultiDomainFunctionTest : public CxxTest::TestSuite {
94
public:
95
96
  // This pair of boilerplate methods prevent the suite being created statically
  // This means the constructor isn't called when running other tests
97
98
99
100
  static MultiDomainFunctionTest *createSuite() {
    return new MultiDomainFunctionTest();
  }
  static void destroySuite(MultiDomainFunctionTest *suite) { delete suite; }
101

102
  MultiDomainFunctionTest() {
103
104
105
106
    multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());
    multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());
    multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());

107
108
    multi.getFunction(0)->setParameter("A", 0);
    multi.getFunction(0)->setParameter("B", 1);
109

110
111
    multi.getFunction(1)->setParameter("A", 1);
    multi.getFunction(1)->setParameter("B", 2);
112

113
114
    multi.getFunction(2)->setParameter("A", 2);
    multi.getFunction(2)->setParameter("B", 3);
115

116
117
118
    domain.addDomain(boost::make_shared<FunctionDomain1DVector>(0, 1, 9));
    domain.addDomain(boost::make_shared<FunctionDomain1DVector>(1, 2, 10));
    domain.addDomain(boost::make_shared<FunctionDomain1DVector>(2, 3, 11));
119
120
  }

121
122
123
124
125
126
127
128
  void test_calc_domain0_only() {
    multi.setDomainIndex(0, 0);
    multi.setDomainIndices(1, std::vector<size_t>());
    multi.setDomainIndices(2, std::vector<size_t>());
    // multi.setDomainIndex(1,1);
    // multi.setDomainIndex(2,2);

    // std::vector<size_t> ii;
129
130
    // ii.emplace_back(0);
    // ii.emplace_back(1);
131
132
    // multi.setDomainIndices(1,ii);
    // ii.clear();
133
134
    // ii.emplace_back(0);
    // ii.emplace_back(2);
135
    // multi.setDomainIndices(2,ii);
136
137

    FunctionValues values(domain);
138
    multi.function(domain, values);
139
140
141

    const double A = multi.getFunction(0)->getParameter("A");
    const double B = multi.getFunction(0)->getParameter("B");
142
143
144
    const FunctionDomain1D &d =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
    for (size_t i = 0; i < 9; ++i) {
145
146
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i]);
    }
147
    for (size_t i = 9; i < 19; ++i) {
148
149
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
150
    for (size_t i = 19; i < 30; ++i) {
151
152
153
154
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
  }

155
156
157
158
  void test_calc_domain1_only() {
    multi.setDomainIndex(0, 1);
    multi.setDomainIndices(1, std::vector<size_t>());
    multi.setDomainIndices(2, std::vector<size_t>());
159
160

    FunctionValues values(domain);
161
    multi.function(domain, values);
162
163
164

    const double A = multi.getFunction(0)->getParameter("A");
    const double B = multi.getFunction(0)->getParameter("B");
165
166
167
    const FunctionDomain1D &d =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 0; i < 9; ++i) {
168
169
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
170
171
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i - 9]);
172
    }
173
    for (size_t i = 19; i < 30; ++i) {
174
175
176
177
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
  }

178
179
180
181
  void test_calc_domain2_only() {
    multi.setDomainIndex(0, 2);
    multi.setDomainIndices(1, std::vector<size_t>());
    multi.setDomainIndices(2, std::vector<size_t>());
182
183

    FunctionValues values(domain);
184
    multi.function(domain, values);
185
186
187

    const double A = multi.getFunction(0)->getParameter("A");
    const double B = multi.getFunction(0)->getParameter("B");
188
189
190
    const FunctionDomain1D &d =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 0; i < 9; ++i) {
191
192
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
193
    for (size_t i = 9; i < 19; ++i) {
194
195
      TS_ASSERT_EQUALS(values.getCalculated(i), 0);
    }
196
197
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i - 19]);
198
199
200
    }
  }

201
  void test_calc_all_domains() {
202
    multi.clearDomainIndices();
203
204
    multi.setDomainIndices(1, std::vector<size_t>());
    multi.setDomainIndices(2, std::vector<size_t>());
205
206

    FunctionValues values(domain);
207
    multi.function(domain, values);
208
209
210

    const double A = multi.getFunction(0)->getParameter("A");
    const double B = multi.getFunction(0)->getParameter("B");
211
212
213
    const FunctionDomain1D &d0 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
    for (size_t i = 0; i < 9; ++i) {
214
215
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
    }
216
217
218
219
    const FunctionDomain1D &d1 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i - 9]);
220
    }
221
222
223
224
    const FunctionDomain1D &d2 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i - 19]);
225
226
227
    }
  }

228
229
230
  void test_set_wrong_index() {
    multi.setDomainIndices(1, std::vector<size_t>());
    multi.setDomainIndices(2, std::vector<size_t>());
231
232
233

    FunctionValues values(domain);

234
    multi.setDomainIndex(0, 3);
Samuel Jones's avatar
Samuel Jones committed
235
236
    TS_ASSERT_THROWS(multi.function(domain, values),
                     const std::invalid_argument &);
237

238
    multi.setDomainIndex(0, 4);
Samuel Jones's avatar
Samuel Jones committed
239
240
    TS_ASSERT_THROWS(multi.function(domain, values),
                     const std::invalid_argument &);
241
242
  }

243
244
  void test_calc() {
    multi.setDomainIndex(0, 0);
245
    std::vector<size_t> ii;
246
247
    ii.emplace_back(0);
    ii.emplace_back(1);
248
    multi.setDomainIndices(1, ii);
249
    ii.clear();
250
251
    ii.emplace_back(0);
    ii.emplace_back(2);
252
    multi.setDomainIndices(2, ii);
253
254

    FunctionValues values(domain);
255
256
257
258
259
260
261
262
263
264
265
    multi.function(domain, values);

    double A = multi.getFunction(0)->getParameter("A") +
               multi.getFunction(1)->getParameter("A") +
               multi.getFunction(2)->getParameter("A");
    double B = multi.getFunction(0)->getParameter("B") +
               multi.getFunction(1)->getParameter("B") +
               multi.getFunction(2)->getParameter("B");
    const FunctionDomain1D &d0 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
    for (size_t i = 0; i < 9; ++i) {
266
267
268
269
270
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
    }

    A = multi.getFunction(1)->getParameter("A");
    B = multi.getFunction(1)->getParameter("B");
271
272
273
274
    const FunctionDomain1D &d1 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i - 9]);
275
276
277
278
    }

    A = multi.getFunction(2)->getParameter("A");
    B = multi.getFunction(2)->getParameter("B");
279
280
281
282
    const FunctionDomain1D &d2 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i - 19]);
283
284
285
    }
  }

286
  void test_attribute() {
287
    multi.clearDomainIndices();
288
289
290
    multi.setLocalAttributeValue(0, "domains", "i");
    multi.setLocalAttributeValue(1, "domains", "0,1");
    multi.setLocalAttributeValue(2, "domains", "0,2");
291
292

    FunctionValues values(domain);
293
294
295
296
297
298
299
300
301
302
303
    multi.function(domain, values);

    double A = multi.getFunction(0)->getParameter("A") +
               multi.getFunction(1)->getParameter("A") +
               multi.getFunction(2)->getParameter("A");
    double B = multi.getFunction(0)->getParameter("B") +
               multi.getFunction(1)->getParameter("B") +
               multi.getFunction(2)->getParameter("B");
    const FunctionDomain1D &d0 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
    for (size_t i = 0; i < 9; ++i) {
304
305
306
307
308
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
    }

    A = multi.getFunction(1)->getParameter("A");
    B = multi.getFunction(1)->getParameter("B");
309
310
311
312
    const FunctionDomain1D &d1 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i - 9]);
313
314
315
316
    }

    A = multi.getFunction(2)->getParameter("A");
    B = multi.getFunction(2)->getParameter("B");
317
318
319
320
    const FunctionDomain1D &d2 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i - 19]);
321
322
323
    }
  }

324
  void test_attribute_domain_range() {
325
    multi.clearDomainIndices();
326
    multi.setLocalAttributeValue(0, "domains", "0-2");
327
    return;
328
329
    multi.setLocalAttributeValue(1, "domains", "i");
    multi.setLocalAttributeValue(2, "domains", "i");
330
331

    FunctionValues values(domain);
332
    multi.function(domain, values);
333
334
335

    double A = multi.getFunction(0)->getParameter("A");
    double B = multi.getFunction(0)->getParameter("B");
336
337
338
    const FunctionDomain1D &d0 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
    for (size_t i = 0; i < 9; ++i) {
339
340
341
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
    }

342
343
344
345
346
347
348
349
    A = multi.getFunction(0)->getParameter("A") +
        multi.getFunction(1)->getParameter("A");
    B = multi.getFunction(0)->getParameter("B") +
        multi.getFunction(1)->getParameter("B");
    const FunctionDomain1D &d1 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i - 9]);
350
351
    }

352
353
354
355
356
357
358
359
    A = multi.getFunction(0)->getParameter("A") +
        multi.getFunction(2)->getParameter("A");
    B = multi.getFunction(0)->getParameter("B") +
        multi.getFunction(2)->getParameter("B");
    const FunctionDomain1D &d2 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i - 19]);
360
361
362
    }
  }

363
364
365
366
367
368
369
370
  void test_attribute_in_FunctionFactory() {
    std::string ini =
        "composite=MultiDomainFunction;"
        "name=MultiDomainFunctionTest_Function,A=0,B=1,$domains=i;"
        "name=MultiDomainFunctionTest_Function,A=1,B=2,$domains=(0,1);"
        "name=MultiDomainFunctionTest_Function,A=2,B=3,$domains=(0,2)";
    auto mfun = boost::dynamic_pointer_cast<CompositeFunction>(
        FunctionFactory::Instance().createInitialized(ini));
371
372

    FunctionValues values(domain);
373
374
375
376
377
378
379
380
381
382
    mfun->function(domain, values);

    double A = mfun->getFunction(0)->getParameter("A") +
               mfun->getFunction(1)->getParameter("A") +
               mfun->getFunction(2)->getParameter("A");
    double B = mfun->getFunction(0)->getParameter("B") +
               mfun->getFunction(1)->getParameter("B") +
               mfun->getFunction(2)->getParameter("B");
    const FunctionDomain1D &d0 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(0));
383
    double checksum = 0;
384
    for (size_t i = 0; i < 9; ++i) {
385
386
387
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
      checksum += values.getCalculated(i);
    }
388
    TS_ASSERT_DIFFERS(checksum, 0);
389
390
391
392

    checksum = 0;
    A = mfun->getFunction(1)->getParameter("A");
    B = mfun->getFunction(1)->getParameter("B");
393
394
395
396
    const FunctionDomain1D &d1 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(1));
    for (size_t i = 9; i < 19; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i - 9]);
397
398
      checksum += values.getCalculated(i);
    }
399
    TS_ASSERT_DIFFERS(checksum, 0);
400
401
402
403

    checksum = 0;
    A = mfun->getFunction(2)->getParameter("A");
    B = mfun->getFunction(2)->getParameter("B");
404
405
406
407
    const FunctionDomain1D &d2 =
        static_cast<const FunctionDomain1D &>(domain.getDomain(2));
    for (size_t i = 19; i < 30; ++i) {
      TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i - 19]);
408
409
      checksum += values.getCalculated(i);
    }
410
    TS_ASSERT_DIFFERS(checksum, 0);
411
412
  }

413
  void test_derivatives_for_tied_parameters() {
414
    multi.clearDomainIndices();
415
416
417
    multi.setDomainIndex(0, 0);
    multi.setDomainIndex(1, 1);
    multi.setDomainIndex(2, 2);
418
419
    {
      JacobianToTestNumDeriv jacobian;
420
421
      multi.functionDeriv(domain, jacobian);
      TS_ASSERT_EQUALS(jacobian.off_diag, 0.0);
422
    }
423
    multi.setAttributeValue("NumDeriv", true);
424
425
    {
      JacobianToTestNumDeriv jacobian;
426
427
      multi.functionDeriv(domain, jacobian);
      TS_ASSERT_EQUALS(jacobian.off_diag, 0.0);
428
    }
429
430
431
    multi.tie("f1.A", "f0.A");
    multi.tie("f2.A", "f0.A");
    multi.setAttributeValue("NumDeriv", false);
432
433
    {
      JacobianToTestNumDeriv jacobian;
434
435
      multi.functionDeriv(domain, jacobian);
      TS_ASSERT_EQUALS(jacobian.off_diag, 0.0);
436
    }
437
    multi.setAttributeValue("NumDeriv", true);
438
439
    {
      JacobianToTestNumDeriv jacobian;
440
441
      multi.functionDeriv(domain, jacobian);
      TS_ASSERT_DIFFERS(jacobian.off_diag, 0.0);
442
443
444
    }
  }

445
446
447
448
449
450
  void test_clone_preserves_domains() {
    const auto copy = multi.clone();
    TS_ASSERT_EQUALS(copy->getNumberDomains(), multi.getNumberDomains());
  }

  void test_string_representation() {
Roman Tolchenov's avatar
Roman Tolchenov committed
451
452
453
454
455
456
    const std::string expected = "composite=MultiDomainFunction,NumDeriv=true;"
                                 "name=MultiDomainFunctionTest_Function,A=0,B="
                                 "1,$domains=i;name=MultiDomainFunctionTest_"
                                 "Function,A=0,B=2,$domains=i;name="
                                 "MultiDomainFunctionTest_Function,A=0,B=3,$"
                                 "domains=i;ties=(f1.A=f0.A,f2.A=f0.A)";
457
458
459
460
    TS_ASSERT_EQUALS(multi.asString(), expected);
    TS_ASSERT_EQUALS(multi.asString(), multi.clone()->asString());
  }

Roman Tolchenov's avatar
Roman Tolchenov committed
461
462
463
464
465
466
467
468
469
470
471
  void test_equivalent_functions() {
    std::string ini =
        "composite=MultiDomainFunction;"
        "name=MultiDomainFunctionTest_Function,A=1,B=2,$domains=i;"
        "name=MultiDomainFunctionTest_Function,A=3,B=4,$domains=i;ties=(f1.A="
        "f0.B)";
    auto mfun = boost::dynamic_pointer_cast<CompositeFunction>(
        FunctionFactory::Instance().createInitialized(ini));

    auto eqFuns = mfun->createEquivalentFunctions();
    TS_ASSERT_EQUALS(eqFuns.size(), 2);
472
473
474
475
    TS_ASSERT_EQUALS(eqFuns[0]->asString(),
                     "name=MultiDomainFunctionTest_Function,A=1,B=2");
    TS_ASSERT_EQUALS(eqFuns[1]->asString(),
                     "name=MultiDomainFunctionTest_Function,A=2,B=4");
Roman Tolchenov's avatar
Roman Tolchenov committed
476
477
  }

478
479
480
481
482
483
484
485
486
487
  void test_string_nested_composite() {
    std::string ini = "composite=MultiDomainFunction;"
                      "(composite=CompositeFunction,$domains=i;"
                      "name=MultiDomainFunctionTest_Function; "
                      "(name=MultiDomainFunctionTest_Function;name="
                      "MultiDomainFunctionTest_Function));"
                      "(composite=CompositeFunction,$domains=i;"
                      "name=MultiDomainFunctionTest_Function; "
                      "(name=MultiDomainFunctionTest_Function;name="
                      "MultiDomainFunctionTest_Function))";
Roman Tolchenov's avatar
Roman Tolchenov committed
488
489
    auto f = FunctionFactory::Instance().createInitialized(ini);
    auto g = FunctionFactory::Instance().createInitialized(f->asString());
490
    TS_ASSERT_EQUALS(g->asString(),
Roman Tolchenov's avatar
Roman Tolchenov committed
491
                     "composite=MultiDomainFunction,NumDeriv=true;(composite="
492
493
494
495
496
497
498
499
500
501
                     "CompositeFunction,NumDeriv=false,$domains=i;name="
                     "MultiDomainFunctionTest_Function,A=0,B=0;(name="
                     "MultiDomainFunctionTest_Function,A=0,B=0;name="
                     "MultiDomainFunctionTest_Function,A=0,B=0));(composite="
                     "CompositeFunction,NumDeriv=false,$domains=i;name="
                     "MultiDomainFunctionTest_Function,A=0,B=0;(name="
                     "MultiDomainFunctionTest_Function,A=0,B=0;name="
                     "MultiDomainFunctionTest_Function,A=0,B=0))");
  }

502
503
504
505
private:
  MultiDomainFunction multi;
  JointDomain domain;
};