tsgConstructSurrogate.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_ADDONS_CONSTRUCT_SURROGATE_HPP
#define __TASMANIAN_ADDONS_CONSTRUCT_SURROGATE_HPP
 
#include "tsgCandidateManager.hpp"
 
namespace TasGrid{
 
using ModelSignature = std::function<void(std::vector<double> const &x, std::vector<double> &y, size_t thread_id)>;
 
constexpr bool mode_parallel = true;
 
constexpr bool mode_sequential = false;
 
constexpr bool with_initial_guess = true;
 
constexpr bool no_initial_guess = false;
 
template<bool parallel_construction, bool use_initial_guess>
void constructCommon(ModelSignature model,
                     size_t max_num_points, size_t num_parallel_jobs, size_t max_samples_per_job,
                     TasmanianSparseGrid &grid,
                     std::function<std::vector<double>(TasmanianSparseGrid &)> candidates,
                     std::string const &checkpoint_filename){
 
    num_parallel_jobs   = std::max(size_t(1), num_parallel_jobs);
    max_samples_per_job = std::max(size_t(1), max_samples_per_job);
    size_t num_dimensions = (size_t) grid.getNumDimensions();
    size_t num_outputs    = (size_t) grid.getNumOutputs();
    CandidateManager manager(num_dimensions, max_samples_per_job); // keeps track of started and ordered samples
    CompleteStorage complete(num_dimensions); // temporarily stores complete samples (batch loading is faster)
 
    std::string filename = checkpoint_filename;
    std::string filename_old = checkpoint_filename + "_old";
 
    if (!filename.empty()){ // recover from an existing checkpoint
        std::ifstream infile(filename, std::ios::binary);
        try{ // attempt to recover from filename
            if (!infile.good()) throw std::runtime_error("missing main checkpoint");
            grid.read(infile, mode_binary);
            complete.read(infile);
        }catch(std::runtime_error &){
            // main file is missing or is corrupt, try the older version
            std::ifstream oldfile(filename_old, std::ios::binary);
            try{
                if (!oldfile.good()) throw std::runtime_error("missing main checkpoint");
                grid.read(oldfile, mode_binary);
                complete.read(oldfile);
            }catch(std::runtime_error &){
                // nothing could be recovered, start over from the current grid
            }
        }
    }
 
    if (!filename.empty()){ // initial checkpoint
        std::ofstream ofs(filename, std::ios::binary);
        grid.write(ofs, mode_binary); // write grid to current
        complete.write(ofs);
    }
 
    // prepare several commonly used steps
    auto checkpoint = [&]()->void{ // keeps two saved states for the constructed grid
        if (!filename.empty()){
            { // copy current into old and write to current
                std::ifstream current_state(filename, std::ios::binary);
                std::ofstream previous_state(filename, std::ios::binary);
                previous_state << current_state.rdbuf();
            }
            std::ofstream ofs(filename, std::ios::binary);
            grid.write(ofs, mode_binary); // write grid to current
            complete.write(ofs);
        }
    };
 
    auto load_complete = [&]()->void{ // loads any complete points, does nothing if getNumStored() is zero
        if (complete.getNumStored() > 0)
            complete.load(grid);
    };
 
    auto refresh_candidates = [&]()->void{ // loads complete and asks for new candidates
        load_complete(); // always load before computing new candidates
        manager = candidates(grid); // get new candidates
    };
 
    size_t total_num_launched = complete.getNumStored() + grid.getNumLoaded(); // count all launched jobs, including the ones already complete
    auto checkout_sample = [&]()->std::vector<double>{ // get the "most-important" point that has not started yet
        auto x = manager.next(max_num_points - total_num_launched);
        if (x.empty()){ // did not find a job, maybe we need to refresh the candidates
            refresh_candidates();
            x = manager.next(max_num_points - total_num_launched); // if this is empty, then we have exhausted all possible candidates
        }
        return x;
    };
 
    // load the initial guess into y (is using initial guess), otherwise set y to the correct size
    auto set_initial_guess = [&](std::vector<double> const &x, std::vector<double> &y)->void{
        if (use_initial_guess){
            if (grid.getNumLoaded()) grid.evaluateBatch(x, y);
            else y.clear();
        }else{
            y.resize(num_outputs * (x.size() / num_dimensions));
        }
    };
 
    if (!grid.isUsingConstruction()) // the procedure assumes dynamic construction
        grid.beginConstruction();
 
    refresh_candidates();
 
    if (parallel_construction == mode_parallel){ // parallel version
        // allocate space for all x and y pairs, will be filled by workers and processed by main
        std::vector<std::vector<double>> x(num_parallel_jobs),
                                         y(num_parallel_jobs, std::vector<double>(max_samples_per_job * num_outputs));
 
        std::vector<int> work_flag(num_parallel_jobs);
        constexpr int flag_done = 0;
        constexpr int flag_computing = 1;
        constexpr int flag_shutdown = 2;
 
        std::condition_variable until_someone_done;
        std::condition_variable until_new_job;
        std::mutex access_count_done;
        int count_done = 0;
 
        // lambda that will handle the work
        auto do_work = [&](size_t thread_id)->void{
 
            int my_flag = flag_computing;
            while(my_flag == flag_computing){
                model(x[thread_id], y[thread_id], thread_id); // does the model evaluations
 
                { // must guarantee sync between work_flag and count_done, use a lock
                    std::lock_guard<std::mutex> lock(access_count_done);
                    work_flag[thread_id] = flag_done;
                    count_done++;
                }
                until_someone_done.notify_one(); // just finished some work, notify the main thread
 
                { // wait till the main thread gives us an new piece of work
                    std::unique_lock<std::mutex> lock(access_count_done);
                    until_new_job.wait(lock, [&]()->bool{ return (work_flag[thread_id] != flag_done); });
                    my_flag = work_flag[thread_id];
                }
            }
        };
 
        // launch initial set of jobs
        std::vector<std::thread> workers(num_parallel_jobs);
        for(size_t id=0; id<num_parallel_jobs; id++){
            x[id] = manager.next(max_num_points - total_num_launched);
            if (!x[id].empty()){
                total_num_launched += x[id].size() / num_dimensions;
                set_initial_guess(x[id], y[id]);
                work_flag[id] = flag_computing;
                workers[id] = std::thread(do_work, id);
            }else{
                work_flag[id] = flag_shutdown; // not enough samples, cancel the thread
            }
        }
 
        auto collect_finished = [&]()->bool{
            bool any_done = false;
            for(size_t id=0; id<num_parallel_jobs; id++){
                if (work_flag[id] == flag_done){
                    if (!x.empty()){ // shouldn't be empty
                        complete.add(x[id], y[id]);
                        manager.complete(x[id]);
                        any_done = true;
                    }
                    if ((grid.getNumLoaded() < 1000) || (double(complete.getNumStored()) / double(grid.getNumLoaded()) > 0.2))
                        load_complete(); // move from complete into the grid
 
                    if (total_num_launched < max_num_points){
                        // refresh the candidates if enough of the current candidates have completed
                        if (double(manager.getNumDone()) / double(manager.getNumCandidates()) > 0.2)
                            refresh_candidates();
 
                        x[id] = checkout_sample(); // if necessary this will call refresh_candidates()
 
                        if (!x[id].empty()){ // if empty, then we have finished all possible candidates (reached tolerance)
                            total_num_launched += x[id].size() / num_dimensions;
                            set_initial_guess(x[id], y[id]);
                            work_flag[id] = flag_computing;
                        }else{
                            work_flag[id] = flag_shutdown; // not enough samples, cancel the thread
                        }
                    }else{
                        work_flag[id] = flag_shutdown; // reached the budget, shutdown the thread
                    }
                }
            }
            return any_done;
        };
 
        while(manager.getNumRunning() > 0){ // main loop
            {   // lock access to the count_done variable
                std::unique_lock<std::mutex> lock(access_count_done);
                // unlock and wait until some else increments the "done" count
                until_someone_done.wait(lock, [&]()->bool{ return (count_done > 0); });
                // the lock is back on at this point, process the completed samples, reset the count and go back to waiting
                count_done = 0;
                if (collect_finished()) checkpoint(); // if new samples were computed, save the state
            } // unlock the access_count_done and notify that we have loaded new jobs
            // without the unlock, the threads will wake up but will not be able to read the worker flags
            until_new_job.notify_all();
        }
 
        load_complete(); // flush completed jobs
 
        for(auto &w : workers) if (w.joinable()) w.join(); // join all threads
 
    }else{
        std::vector<double> x(grid.getNumDimensions()), y( grid.getNumOutputs());
 
        while((total_num_launched < max_num_points) && (manager.getNumCandidates() > 0)){
            x = manager.next(max_num_points - total_num_launched);
            if (x.empty()){ // need more candidates
                refresh_candidates();
                x = manager.next(max_num_points - total_num_launched); // if this is empty, then we have exhausted the candidates
            }
            if (!x.empty()){ // could be empty if there are no more candidates
                total_num_launched += x.size() / num_dimensions;
                set_initial_guess(x, y);
                model(x, y, 0); // compute a sample
                complete.add(x, y);
                manager.complete(x);
 
                // the fist thousand points can be loaded one at a time, then add when % increase of the grid is achieved
                if ((grid.getNumLoaded() < 1000) || (double(complete.getNumStored()) / double(grid.getNumLoaded()) > 0.2))
                    load_complete(); // also does checkpoint save
                // if done with the top % of the grid points, recompute the candidates
                if (double(manager.getNumDone()) / double(manager.getNumCandidates()) > 0.2)
                    refresh_candidates();
            }
            checkpoint();
        }
 
        load_complete(); // flush completed jobs
    }
}
 
template<bool parallel_construction = TasGrid::mode_parallel, bool initial_guess = no_initial_guess>
void constructSurrogate(ModelSignature model,
                        size_t max_num_points, size_t num_parallel_jobs, size_t max_samples_per_job,
                        TasmanianSparseGrid &grid,
                        double tolerance, TypeRefinement criteria, int output = -1,
                        std::vector<int> const &level_limits = std::vector<int>(),
                        std::string const &checkpoint_filename = std::string()){
    if (!grid.isLocalPolynomial() && !grid.isWavelet()) throw std::runtime_error("ERROR: construction (with tolerance and criteria) called for a grid that is not local polynomial or wavelet.");
    constructCommon<parallel_construction, initial_guess>
                                          (model, max_num_points, num_parallel_jobs, max_samples_per_job, grid,
                                           [&](TasmanianSparseGrid &g)->std::vector<double>{
                                               return g.getCandidateConstructionPoints(tolerance, criteria, output, level_limits);
                                           }, checkpoint_filename);
}
 
template<bool parallel_construction = TasGrid::mode_parallel, bool initial_guess = no_initial_guess>
void constructSurrogate(ModelSignature model,
                        size_t max_num_points, size_t num_parallel_jobs, size_t max_samples_per_job,
                        TasmanianSparseGrid &grid,
                        TypeDepth type, std::vector<int> const &anisotropic_weights = std::vector<int>(),
                        std::vector<int> const &level_limits = std::vector<int>(),
                        std::string const &checkpoint_filename = std::string()){
    constructCommon<parallel_construction, initial_guess>
                                           (model, max_num_points, num_parallel_jobs, max_samples_per_job, grid,
                                           [&](TasmanianSparseGrid &g)->std::vector<double>{
                                               return g.getCandidateConstructionPoints(type, anisotropic_weights, level_limits);
                                           }, checkpoint_filename);
}
 
template<bool parallel_construction = TasGrid::mode_parallel, bool initial_guess = no_initial_guess>
void constructSurrogate(ModelSignature model,
                        size_t max_num_points, size_t num_parallel_jobs, size_t max_samples_per_job,
                        TasmanianSparseGrid &grid,
                        TypeDepth type, int output, std::vector<int> const &level_limits = std::vector<int>(),
                        std::string const &checkpoint_filename = std::string()){
    constructCommon<parallel_construction, initial_guess>
                                           (model, max_num_points, num_parallel_jobs, max_samples_per_job, grid,
                                           [&](TasmanianSparseGrid &g)->std::vector<double>{
                                               return g.getCandidateConstructionPoints(type, output, level_limits);
                                           }, checkpoint_filename);
}
 
}
 
#endif