tsgDreamSample.hpp

Go to the documentation of this file.

/*
 * Copyright (c) 2017, Miroslav Stoyanov
 *
 * This file is part of
 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
 *    and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
 *    or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
 * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
 * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
 * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
 */
 
#ifndef __TASMANIAN_DREAM_SAMPLE_HPP
#define __TASMANIAN_DREAM_SAMPLE_HPP
 
#include "tsgDreamState.hpp"
#include "tsgDreamCoreRandom.hpp"
#include "tsgDreamCorePDF.hpp"
 
namespace TasDREAM{
 
using DreamDomain = std::function<bool(std::vector<double> const &x)>;
 
inline DreamDomain hypercube(std::vector<double> const &lower, std::vector<double> const &upper){
    return [=](const std::vector<double> &x)->bool{
        auto il = lower.begin(), iu = upper.begin();
        for(auto v : x) if ((v < *il++) || (v > *iu++)) return false;
        return true;
    };
}
 
inline void no_update(std::vector<double> &){}
 
inline double const_one(){ return 1.0; }
 
template<int weight_percent>
double const_percent(){ return ((double) weight_percent) / 100.0; }
 
inline void uniform_prior(TypeSamplingForm, const std::vector<double> &, std::vector<double> &values){ values.clear(); }
 
using DreamPDF = std::function<void(const std::vector<double> &candidates, std::vector<double> &values)>;
 
using DreamModel = std::function<void(const std::vector<double> &candidates, std::vector<double> &outputs)>;
 
using DreamLikelihood = std::function<void(TypeSamplingForm form, const std::vector<double> &model_outputs, std::vector<double> &likely)>;
 
using DreamPrior = std::function<void(TypeSamplingForm form, const std::vector<double> &candidates, std::vector<double> &values)>;
 
using DreamMergedLikelyModel = std::function<void(const std::vector<double> &candidates, std::vector<double> &values)>;
 
template<TypeSamplingForm form = regform>
DreamPDF posterior(DreamModel model,
          DreamLikelihood likelihood,
          DreamPrior prior){
    return [=](const std::vector<double> &candidates, std::vector<double> &values)->void{
        std::vector<double> model_outs;
        model(candidates, model_outs);
        likelihood(form, model_outs, values);
 
        std::vector<double> prior_vals(values.size());
        prior(form, candidates, prior_vals);
 
        auto iv = values.begin();
        if (form == regform){
            for(auto p : prior_vals) *iv++ *= p;
        }else{
            for(auto p : prior_vals) *iv++ += p;
        }
    };
}
 
template<TypeSamplingForm form = regform>
DreamPDF posterior(DreamMergedLikelyModel likelihood_model,
          DreamPrior prior){
    return [=](const std::vector<double> &candidates, std::vector<double> &values)->void{
        likelihood_model(candidates, values);
 
        std::vector<double> prior_vals(values.size());
        prior(form, candidates, prior_vals);
 
        auto iv = values.begin();
        if (form == regform){
            for(auto p : prior_vals) *iv++ *= p;
        }else{
            for(auto p : prior_vals) *iv++ += p;
        }
    };
}
 
 
template<TypeSamplingForm form = regform>
void SampleDREAM(int num_burnup, int num_collect,
                 DreamPDF probability_distribution,
                 DreamDomain inside,
                 TasmanianDREAM &state,
                 std::function<void(std::vector<double> &x)> independent_update = no_update,
                 std::function<double(void)> differential_update = const_one,
                 std::function<double(void)> get_random01 = tsgCoreUniform01){
 
    size_t num_chains = (size_t) state.getNumChains(), num_dimensions = (size_t) state.getNumDimensions();
    double unitlength = (double) num_chains;
 
    if (num_chains == 0) return; // no sampling with a null state
 
    if (!state.isStateReady()) throw std::runtime_error("ERROR: DREAM sampling requires that the setState() has been called first on the TasmanianDREAM.");
 
    if (!state.isPDFReady()) // initialize probability density (if not initialized already)
        state.setPDFvalues(probability_distribution);
 
    if (num_collect > 0) // pre-allocate memory for the new history
        state.expandHistory(num_collect);
 
    int total_iterations = std::max(num_burnup, 0) + std::max(num_collect, 0);
    for(int t = 0; t < total_iterations; t++){
        std::vector<double> candidates, values;
        candidates.reserve(num_chains * num_dimensions);
        values.reserve(num_chains);
 
        std::vector<bool> valid(num_chains, true); // keep track whether the samples need to be evaluated
 
        for(size_t i=0; i<num_chains; i++){
            std::vector<double> propose(num_dimensions);
 
            size_t jindex = (size_t) (get_random01() * unitlength);
            size_t kindex = (size_t) (get_random01() * unitlength);
            if (jindex >= num_chains) jindex = num_chains - 1; // this is needed in case get_random01() returns 1
            if (kindex >= num_chains) jindex = num_chains - 1;
 
            state.getIJKdelta(i, jindex, kindex, differential_update(), propose); // propose = s_i + w ( s_k - s_j)
            independent_update(propose); // propose += correction
 
            if (inside(propose)){
                candidates.insert(candidates.end(), propose.begin(), propose.end());
                values.resize(values.size() + 1);
            }else{
                valid[i] = false;
            }
        }
 
        if (!candidates.empty()) // block the pathological case of all proposals leaving the domain
            probability_distribution(candidates, values);
 
        std::vector<double> new_state(num_chains * num_dimensions), new_values(num_chains);
 
        auto icand = candidates.begin(); // loop over all candidates and values, accept or reject
        auto ival = values.begin();
 
        size_t accepted = 0;
 
        for(size_t i=0; i<num_chains; i++){
            bool keep_new = valid[i]; // if not valid, automatically reject
            if (valid[i]){ // apply random test
                if (*ival > state.getPDFvalue(i)){ // if the new value has higher probability, automatically accept
                    keep_new = true;
                }else{
                    if (form == regform){
                        keep_new = (*ival / state.getPDFvalue(i) >= get_random01()); // keep if the new value has higher probability
                    }else{
                        keep_new = (*ival - state.getPDFvalue(i) >= log(get_random01()));
                    }
                    //std::cout << "Trsh = " << *ival / state.getPDFvalue(i) << "   " << ((keep_new) ? "Accept" : "Reject") << std:: endl;
                }
            }
 
            if (keep_new){
                std::copy_n(icand, num_dimensions, new_state.begin() + i * num_dimensions);
                new_values[i] = *ival;
                accepted++; // accepted one more proposal
            }else{ // reject and reuse the old state
                state.getChainState((int) i, &*(new_state.begin() + i * num_dimensions));
                new_values[i] = state.getPDFvalue(i);
            }
 
            if (valid[i]){ // kept or rejected, if this sample was valid then move to the next sample in the list
                std::advance(icand, num_dimensions);
                ival++;
            }
        }
 
        state.setState(new_state);
        state.setPDFvalues(new_values);
 
        if (t >= num_burnup)
            state.saveStateHistory(accepted);
    }
}
 
 
template<TypeSamplingForm form = regform>
void SampleDREAM(int num_burnup, int num_collect,
                 DreamPDF probability_distribution,
                 DreamDomain inside,
                 TasmanianDREAM &state,
                 TypeDistribution dist, double magnitude,
                 std::function<double(void)> differential_update = const_one,
                 std::function<double(void)> get_random01 = tsgCoreUniform01){
    if (dist == dist_uniform){
        SampleDREAM<form>(num_burnup, num_collect, probability_distribution, inside, state,
                         [&](std::vector<double> &x)->void{ applyUniformUpdate(x, magnitude, get_random01); }, differential_update, get_random01);
    }else if (dist == dist_gaussian){
        SampleDREAM<form>(num_burnup, num_collect, probability_distribution, inside, state,
                         [&](std::vector<double> &x)->void{ applyGaussianUpdate(x, magnitude, get_random01); }, differential_update, get_random01);
    }else{ // assuming none
        SampleDREAM<form>(num_burnup, num_collect, probability_distribution, inside, state, no_update, differential_update, get_random01);
    }
}
 
}
 
#endif