egs_scoring.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ scoring headers
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:
#
###############################################################################
*/


#ifndef EGS_SCORING_
#define EGS_SCORING_

#include "egs_libconfig.h"
#include "egs_functions.h"
#include "egs_math.h"

#include <iostream>
using namespace std;

class EGS_EXPORT EGS_ScoringSingle {

public:

    EGS_ScoringSingle() : sum(0), sum2(0), tmp(0), current_ncase(0) {};

    inline void score(unsigned short ncase, EGS_Float f) {
        if (ncase == current_ncase) {
            tmp += f;
        }
        else {
            finishCase(ncase,f);
        }
    };

    inline void finishCase(unsigned short new_case, EGS_Float new_result) {
        current_ncase = new_case;
        sum += tmp;
        sum2 += tmp*tmp;
        tmp = new_result;
    };

    EGS_Float currentScore() const {
        return tmp;
    };

    void      currentScore(EGS_Float &s, unsigned short &ncase) const {
        s = tmp;
        ncase = current_ncase;
    };

    void      currentScore(double &s, double &s2) {
        s=sum;
        s2=sum2;
    };

    void      currentResult(EGS_I64 ncase, double &r, double &dr) {
        r = sum + tmp;
        dr = sum2 + tmp*tmp;
        r /= ncase;
        dr /= ncase;
        dr -= r*r;
        if (dr > 0) {
            dr = sqrt(dr/(ncase-1));
        }
    };

    bool storeState(ostream &data) {
        //sum += tmp; sum2 += tmp*tmp; tmp = 0;
        //data << current_ncase << "  " << sum << "  " << sum2 << endl;
        data << current_ncase << "  " << sum+tmp << "  " << sum2+tmp *tmp
             << "\n";
        return data.good();
    };

    bool setState(istream &data) {
        data >> current_ncase >> sum >> sum2;
        tmp = 0;
        return data.good();
    };

    void reset() {
        current_ncase = 0;
        tmp = 0;
        sum = 0;
        sum2 = 0;
    };

    EGS_ScoringSingle &operator+=(const EGS_ScoringSingle &x) {
        sum += tmp + x.sum + x.tmp;
        sum2 += tmp*tmp + x.sum2 + x.tmp*x.tmp;
        current_ncase = 0;
        tmp = 0;
        return *this;
    };


protected:

    double         sum;
    double         sum2;
    EGS_Float      tmp;
    unsigned short current_ncase;

};

class EGS_EXPORT EGS_ScoringArray {

public:

    EGS_ScoringArray(int N);

    ~EGS_ScoringArray();

    void setHistory(EGS_I64 ncase);

    inline void score(int ireg, EGS_Float f) {
        result[ireg].score(current_ncase_short,f);
    };

    EGS_Float currentScore(int ireg) const {
        return result[ireg].currentScore();
    };

    EGS_Float thisHistoryScore(int ireg) const {
        EGS_Float res;
        unsigned short nc;
        result[ireg].currentScore(res,nc);
        return nc == current_ncase_short ? res : 0;
    };

    void currentScore(int ireg, double &s, double &s2) {
        result[ireg].currentScore(s,s2);
    };

    void currentResult(int ireg, double &r, double &dr) {
        result[ireg].currentResult(current_ncase,r,dr);
    };

    void reportResults(double norm, const char *title, bool relative_error,
                       const char *format = 0);

    bool storeState(ostream &data) {
        data << nreg << "  " << current_ncase_short << "\n";
        if (!egsStoreI64(data,current_ncase)) {
            return false;
        }
        if (!egsStoreI64(data,current_ncase_65536)) {
            return false;
        }
        data << "\n";
        for (int j=0; j<nreg; j++) {
            if (!result[j].storeState(data)) {
                return false;
            }
        }
        return true;
    };

    bool setState(istream &data) {
        int nreg1;
        data >> nreg1 >> current_ncase_short;
        if (!data.good() || nreg1 < 1) {
            return false;
        }
        if (!egsGetI64(data,current_ncase)) {
            return false;
        }
        if (!egsGetI64(data,current_ncase_65536)) {
            return false;
        }
        if (nreg1 != nreg) {
            if (nreg > 0) {
                delete [] result;
            }
            nreg = nreg1;
            result = new EGS_ScoringSingle [nreg];
        }
        for (int j=0; j<nreg; j++) {
            if (!result[j].setState(data)) {
                return false;
            }
        }
        return true;
    };

    void reset() {
        current_ncase = 0;
        current_ncase_65536 = 0;
        current_ncase_short = 0;
        for (int j=0; j<nreg; j++) {
            result[j].reset();
        }
    };

    EGS_ScoringArray &operator+=(const EGS_ScoringArray &x) {
        current_ncase += x.current_ncase;
        current_ncase_65536 = current_ncase >> 16;
        EGS_I64 aux = current_ncase - (current_ncase_65536 << 16);
        current_ncase_short = (unsigned short) aux;
        for (int j=0; j<nreg; j++) {
            result[j] += x.result[j];
        }
        return *this;
    };

    int bins() const {
        return nreg;
    };

    int regions() const {
        return nreg;
    };

protected:

    EGS_I64           current_ncase;
    EGS_I64           current_ncase_65536;
    int               nreg;
    EGS_ScoringSingle *result;
    unsigned short    current_ncase_short;

};

#endif