tsgHierarchyManipulator.hpp

/*
 * 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 __TSG_HIERARCHY_MANIPULATOR_HPP
#define __TSG_HIERARCHY_MANIPULATOR_HPP
 
#include "tsgRuleLocalPolynomial.hpp"
#include "tsgIndexManipulator.hpp"


 
namespace TasGrid{

namespace HierarchyManipulations{

template<RuleLocal::erule effrule>
Data2D<int> computeDAGup(MultiIndexSet const &mset);
Data2D<int> computeDAGup(MultiIndexSet const &mset, RuleLocal::erule effrule);
template<RuleLocal::erule effrule>
Data2D<int> computeDAGup(MultiIndexSet const &mset, bool &is_complete);

template<RuleLocal::erule effrule>
Data2D<int> computeDAGDown(MultiIndexSet const &mset);

template<RuleLocal::erule effrule>
std::vector<int> computeLevels(MultiIndexSet const &mset);
std::vector<int> computeLevels(MultiIndexSet const &mset, RuleLocal::erule effrule);

template<RuleLocal::erule effrule>
void completeToLower(MultiIndexSet const &mset, MultiIndexSet &refined);

template<RuleLocal::erule effrule, typename callable_method>
void touchAllImmediateRelatives(std::vector<int> &point, MultiIndexSet const &mset, callable_method apply){
    int max_kids = RuleLocal::getMaxNumKids<effrule>();
    for(auto &v : point){
        int save = v; // replace one by one each index of p with either parent or kid
 
        // check the parents
        v = RuleLocal::getParent<effrule>(save);
        if (v > -1){
            int parent_index = mset.getSlot(point);
            if (parent_index > -1)
                apply(parent_index);
        }
 
        v = RuleLocal::getStepParent<effrule>(save);
        if (v > -1){
            int parent_index = mset.getSlot(point);
            if (parent_index > -1)
                apply(parent_index);
        }
 
        for(int k=0; k<max_kids; k++){
            v = RuleLocal::getKid<effrule>(save, k);
            if (v > -1){
                int kid_index = mset.getSlot(point);
                if (kid_index > -1)
                    apply(kid_index);
            }
        }
 
        v = save; // restore the original index for the next iteration
    }
}

template<RuleLocal::erule effrule>
MultiIndexSet getLevelZeroPoints(size_t num_dimensions){
    int num_parents = 0;
    while(RuleLocal::getParent<effrule>(num_parents) == -1) num_parents++;
    return MultiIndexManipulations::generateFullTensorSet(std::vector<int>(num_dimensions, num_parents));
}

template<RuleLocal::erule effrule>
MultiIndexSet getLargestConnected(MultiIndexSet const &current, MultiIndexSet const &candidates){
    if (candidates.empty()) return MultiIndexSet();
    auto num_dimensions = candidates.getNumDimensions();
 
    // always consider the points without parents
    MultiIndexSet level_zero = getLevelZeroPoints<effrule>(num_dimensions);
 
    MultiIndexSet result; // keep track of the cumulative result
    MultiIndexSet total = current; // forms a working copy of the entire merged graph
 
    // do not consider the points already included in total, complexity is level_zero.getNumIndexes()
    if (!total.empty()) level_zero = level_zero - total;
 
    if (level_zero.getNumIndexes() > 0){ // level zero nodes are missing from current
        Data2D<int> roots(num_dimensions, 0);
        for(int i=0; i<level_zero.getNumIndexes(); i++){
            std::vector<int> p(level_zero.getIndex(i), level_zero.getIndex(i) + num_dimensions);
            if (!candidates.missing(p))
                roots.appendStrip(p);
        }
 
        result = MultiIndexSet(roots);
        total += result;
    }
 
    if (total.empty()) return MultiIndexSet(); // current was empty and no roots could be added
 
    int max_kids      = RuleLocal::getMaxNumKids<effrule>();
    int max_relatives = RuleLocal::getMaxNumParents<effrule>() + max_kids;
    Data2D<int> update;
    do{
        update = Data2D<int>(num_dimensions, 0);
 
        for(int i=0; i<total.getNumIndexes(); i++){
            std::vector<int> relative(total.getIndex(i), total.getIndex(i) + num_dimensions);
            for(auto &r : relative){
                int k = r; // save the value
                for(int j=0; j<max_relatives; j++){
                    r = (j < max_kids) ? RuleLocal::getKid<effrule>(k, j)
                                       : ((j - max_kids == 0) ? RuleLocal::getParent<effrule>(k) : RuleLocal::getStepParent<effrule>(k));
                    if ((r != -1) && !candidates.missing(relative) && total.missing(relative))
                        update.appendStrip(relative);
                }
                r = k;
            }
        }
 
        if (update.getNumStrips() > 0){
            MultiIndexSet update_set(update);
            result += update_set;
            total  += update_set;
        }
    }while(update.getNumStrips() > 0);
 
    return result;
}
 

template<typename T>
std::vector<Data2D<T>> splitByLevels(size_t stride, typename std::vector<T>::const_iterator ibegin, typename std::vector<T>::const_iterator iend, std::vector<int>::const_iterator ilevels){
    size_t top_level = (size_t) *std::max_element(ilevels, ilevels + std::distance(ibegin, iend) / stride);
 
    std::vector<Data2D<T>> split(top_level + 1, Data2D<T>(stride, 0));
 
    for( struct { std::vector<int>::const_iterator il; typename std::vector<T>::const_iterator idata;} v = {ilevels, ibegin};
         v.idata != iend;
         v.il++, std::advance(v.idata, stride))
        split[*v.il].appendStrip(v.idata);
 
    return split;
}

inline std::vector<Data2D<int>> splitByLevels(MultiIndexSet const &mset, std::vector<int> const &levels){
    return splitByLevels<int>(mset.getNumDimensions(), mset.begin(), mset.end(), levels.begin());
}

template<typename T>
inline std::vector<Data2D<T>> splitByLevels(Data2D<T> const &data, std::vector<int> const &levels){
    return splitByLevels<T>(data.getStride(), data.begin(), data.end(), levels.begin());
}

inline std::vector<Data2D<double>> splitByLevels(StorageSet const &stortage, std::vector<int> const &levels){
    return splitByLevels<double>(stortage.getNumOutputs(), stortage.begin(), stortage.end(), levels.begin());
}

class SplitDirections{
public:
    SplitDirections(const MultiIndexSet &points);
    ~SplitDirections(){}

    int getNumJobs() const{ return (int) job_pnts.size(); }
    int getJobDirection(int job) const{ return job_directions[job]; }
    int getJobNumPoints(int job) const{ return (int) job_pnts[job].size(); }
    const int* getJobPoints(int job) const{ return job_pnts[job].data(); }
    int getMaxNumPoints() const { return (job_pnts.size() > 0) ? (int) std::max_element(job_pnts.begin(), job_pnts.end(),
        [&](std::vector<int> const &a, std::vector<int> const& b)->bool{ return (a.size() < b.size()); })->size() : 0; }
 
private:
    std::vector<int> job_directions;
    std::vector<std::vector<int>> job_pnts;
};
 
}
 
}
 
#endif