SumErrSqRadialPeak.h 3.89 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
//
// Struct SumErrSqRadialPeak
//   A simple expression computing the sum of all peak errors (given as
//   first and second argument) for a range of peaks (third argument and fourth argument)
//   according to
//
//   \f[
//     result = \frac{1}{n} * \sqrt{\sum_{i=start}^end (measurement_i - value_i)^2}
//   \f]
//
// Copyright (c) 2017 - 2018, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
// and the paper
// "Matching of turn pattern measurements for cyclotrons using multiobjective optimization"
// (https://doi.org/10.1103/PhysRevAccelBeams.22.064602)
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
adelmann's avatar
adelmann committed
30 31 32 33 34
#ifndef __SUMERRSQRADIALPEAK_H__
#define __SUMERRSQRADIALPEAK_H__

#include <string>

35
#include "boost/type_traits/remove_cv.hpp"
adelmann's avatar
adelmann committed
36 37 38 39 40 41 42 43 44 45 46 47 48
#include "boost/variant/get.hpp"
#include "boost/variant/variant.hpp"
#include "boost/smart_ptr.hpp"

#include "Util/Types.h"
#include "Util/PeakReader.h"
#include "Expression/Parser/function.hpp"

struct SumErrSqRadialPeak {

    static const std::string name;

    Expressions::Result_t operator()(client::function::arguments_t args) {
49 50 51 52
        if (args.size() != 4) {
            throw OptPilotException("SumErrSqRadialPeak::operator()",
                                    "sumErrSqRadialPeak expects 4 arguments, " + std::to_string(args.size()) + " given");
        }
adelmann's avatar
adelmann committed
53 54 55 56 57

        meas_filename_ = boost::get<std::string>(args[0]);
        sim_filename_  = boost::get<std::string>(args[1]);
        begin_         = boost::get<double>(args[2]);
        end_           = boost::get<double>(args[3]);
58

adelmann's avatar
adelmann committed
59
        bool is_valid = true;
60

adelmann's avatar
adelmann committed
61 62 63 64 65
        boost::scoped_ptr<PeakReader> meas_peaks(new PeakReader(meas_filename_));
        boost::scoped_ptr<PeakReader> sim_peaks(new PeakReader(sim_filename_));
        try {
            sim_peaks->parseFile();
            meas_peaks->parseFile();
66

adelmann's avatar
adelmann committed
67 68 69
            if ( end_ < begin_ || end_ < 0 || begin_ < 0 )
                throw OptPilotException("SumErrSqRadialPeak::operator()",
                                        "Error check turn number range");
70

adelmann's avatar
adelmann committed
71 72 73 74 75 76 77
        } catch (OptPilotException &ex) {
            std::cout << "Caught exception: " << ex.what() << std::endl;
            is_valid = false;
        }

        double sum = 0;
        int nPeaks = end_ - begin_ + 1;
78

adelmann's avatar
adelmann committed
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
        for (int turn = begin_; turn < end_ + 1; ++turn) {
            double sim_value = 0.0, meas_value = 0.0;
            try {
                sim_peaks->getPeak(turn, sim_value);
                meas_peaks->getPeak(turn, meas_value);
            } catch(OptPilotException &e) {
                std::cout << "Exception while getting value "
                          << "from peak file: " << e.what()
                          << std::endl;
                is_valid = false;
            }
            double val = meas_value - sim_value;
            sum += val * val;
        }

        return boost::make_tuple(std::sqrt(sum) / (double)nPeaks, is_valid);
    }

private:
    std::string meas_filename_;
    std::string sim_filename_;
    int begin_;
    int end_;

    // define a mapping to arguments in argument vector
    boost::tuple<std::string, std::string, int, int> argument_types;
gsell's avatar
gsell committed
105 106
    // :FIXME: remove unused enum
#if 0
adelmann's avatar
adelmann committed
107 108 109 110 111 112
    enum {
          meas_filename
        , sim_filename
        , begin
        , end
    } argument_type_id;
gsell's avatar
gsell committed
113
#endif
adelmann's avatar
adelmann committed
114 115 116
};

#endif
gsell's avatar
gsell committed
117 118 119 120
// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c
// c-basic-offset: 4
gsell's avatar
gsell committed
121 122 123
// indent-tabs-mode: nil
// require-final-newline: nil
// End: