tsgIndexManipulator.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_INDEX_MANIPULATOR_HPP
#define __TSG_INDEX_MANIPULATOR_HPP
 
#include "tsgIndexSets.hpp"
#include "tsgOneDimensionalWrapper.hpp"


 
namespace TasGrid{

namespace MultiIndexManipulations{

inline MultiIndexSet generateFullTensorSet(std::vector<int> const &num_entries){
    size_t num_dimensions = num_entries.size();
    int num_total = 1;
    for(auto &l : num_entries) num_total *= l;
    std::vector<int> indexes(Utils::size_mult(num_dimensions, num_total));
    auto iter = indexes.rbegin();
    for(int i=num_total-1; i>=0; i--){
        int t = i;
        auto l = num_entries.rbegin();
        // in order to generate indexes in the correct order, the for loop must go backwards
        for(size_t j = 0; j<num_dimensions; j++){
            *iter++ = (t % *l);
            t /= *l++;
        }
    }
    return MultiIndexSet(num_dimensions, std::move(indexes));
}

inline MultiIndexSet generateLowerMultiIndexSet(size_t num_dimensions, std::function<bool(const std::vector<int> &index)> inside){
    size_t c = num_dimensions -1;
    bool is_in = true;
    std::vector<int> root(num_dimensions, 0);
    std::vector<int> indexes;
    while(is_in || (c > 0)){
        if (is_in){
            indexes.insert(indexes.end(), root.begin(), root.end());
            c = num_dimensions-1;
            root[c]++;
        }else{
            std::fill(root.begin() + c, root.end(), 0);
            root[--c]++;
        }
        is_in = inside(root);
    }
    return MultiIndexSet(num_dimensions, std::move(indexes));
}

void completeSetToLower(MultiIndexSet &set);

struct ProperWeights{
    ProperWeights(size_t num_dimensions, TypeDepth type, std::vector<int> const &weights){
        contour = OneDimensionalMeta::getControurType(type);
        if (weights.empty()) linear = std::vector<int>(num_dimensions, 1);
        if (contour == type_level){ // linear weights only
            if (!weights.empty()){
                linear = weights; // simplest type
            }
        }else if (contour == type_curved){
            if (weights.empty()){
                curved = std::vector<double>(num_dimensions, 0.0); // cannot define default log-correction
            }else{
                linear = std::vector<int>(weights.begin(), weights.begin() + num_dimensions);
                curved = std::vector<double>(num_dimensions);
                std::transform(weights.begin() + num_dimensions, weights.end(), curved.begin(), [&](int i)->double{ return (double) i; });
            }
        }else{ // hyperbolic
            if (weights.empty()){
                curved = std::vector<double>(num_dimensions, 1.0);
            }else{
                linear = std::vector<int>(num_dimensions, 1); // used only for offset normalization
                double exponent_normalization = (double) *std::min_element(weights.begin(), weights.end());
                curved = std::vector<double>(num_dimensions);
                std::transform(weights.begin(), weights.end(), curved.begin(),
                               [&](int i)->double{ return ((double) i) / exponent_normalization; });
            }
        }
    }
    bool provenLower() const{
        if (contour == type_curved)
            for(size_t i=0; i<linear.size(); i++)
                if ((double)linear[i] + curved[i] < 0) return false;
        return true;
    }
    int minLinear() const{ return *std::min_element(linear.begin(), linear.end()); }
    size_t getNumDimensions() const{ return linear.size(); }

    TypeDepth contour;
    std::vector<int> linear;
    std::vector<double> curved;
};

template<typename CacheType, TypeDepth contour, bool isotropic>
std::vector<std::vector<CacheType>> generateLevelWeightsCache(ProperWeights const &weights, std::function<int(int i)> rule_exactness, int offset){
    size_t num_dimensions = weights.getNumDimensions();
    std::vector<std::vector<CacheType>> cache(num_dimensions);
    std::vector<int> exactness_cache;
 
    // if we know the sized here, then assign, otherwise dynamically build later
    if (isotropic){
        for(auto &vec : cache) vec.reserve((size_t) offset);
        exactness_cache.resize((size_t) offset);
        exactness_cache[0] = 0;
        for(int i=1; i<offset; i++)
            exactness_cache[i] = 1 + rule_exactness(i - 1);
    }else{
        exactness_cache.push_back(0); // exactness always starts from 0
    }
 
    for(size_t j=0; j<num_dimensions; j++){
        int wl = weights.linear[j]; // linear weights
        double wc = (contour == type_level) ? 0.0 : weights.curved[j]; // curved weights
 
        size_t i = 0; // keep track of the index
        CacheType w = (CacheType) ((contour == type_hyperbolic) ? 1 : 0);
        cache[j].push_back(w); // initial entry
        do{ // accumulate the cache
            i++;
            if (!isotropic && (i >= exactness_cache.size()))
                exactness_cache.push_back(1 + rule_exactness((int)(i - 1)));
 
            int e = exactness_cache[i];
 
            if (contour == type_level){
                w = (CacheType)(wl * e);
            }else if (contour == type_curved){
                w = (CacheType)(wl * e) + (CacheType)(wc * std::log1p((CacheType)e));
            }else{ // must be hyperbolic
                w = (CacheType)(pow((CacheType) (1 + e), wc));
            }
 
            cache[j].push_back(w);
        // isotropic mode works until offset, anisotropic mode works until weight reaches the offset
        }while( (!isotropic && (std::ceil(w) <= (CacheType) offset)) || (isotropic && (i + 1 < (size_t) offset)) );
    }
 
    return cache;
}

template<typename CacheType, TypeDepth contour>
inline CacheType getIndexWeight(int const index[], std::vector<std::vector<CacheType>> const &cache){
    CacheType w = (contour == type_hyperbolic) ? 1 : 0;
    for(size_t j=0; j<cache.size(); j++)
        if (contour == type_hyperbolic)
            w *= cache[j][index[j]];
        else
            w += cache[j][index[j]];
    return w;
}

MultiIndexSet selectTensors(size_t num_dimensions, int offset, TypeDepth type, std::function<int(int i)> rule_exactness,
                            std::vector<int> const &anisotropic_weights, std::vector<int> const &level_limits);

std::vector<int> computeLevels(MultiIndexSet const &mset);

std::vector<int> getMaxIndexes(const MultiIndexSet &mset);

Data2D<int> computeDAGup(MultiIndexSet const &mset);

MultiIndexSet selectFlaggedChildren(const MultiIndexSet &mset, const std::vector<bool> &flagged, const std::vector<int> &level_limits);

MultiIndexSet generateNestedPoints(const MultiIndexSet &tensors, std::function<int(int)> getNumPoints);

MultiIndexSet generateNonNestedPoints(const MultiIndexSet &tensors, const OneDimensionalWrapper &wrapper);

void resortIndexes(const MultiIndexSet &iset, std::vector<std::vector<int>> &map, std::vector<std::vector<int>> &lines1d);

template<bool nested>
std::vector<int> referencePoints(const int levels[], const OneDimensionalWrapper &wrapper, const MultiIndexSet &points){
    size_t num_dimensions = (size_t) points.getNumDimensions();
    std::vector<int> num_points(num_dimensions);
    int num_total = 1; // this will be a subset of all points, no danger of overflow
    for(size_t j=0; j<num_dimensions; j++) num_points[j] = wrapper.getNumPoints(levels[j]);
    for(auto n : num_points) num_total *= n;
 
    std::vector<int> refs(num_total);
    std::vector<int> p(num_dimensions);
 
    for(int i=0; i<num_total; i++){
        int t = i;
        auto n = num_points.rbegin();
        for(int j=(int) num_dimensions-1; j>=0; j--){
            p[j] = (nested) ? t % *n : wrapper.getPointIndex(levels[j], t % *n);
            t /= *n++;
        }
        refs[i] = points.getSlot(p);
    }
    return refs;
}

std::vector<int> computeTensorWeights(MultiIndexSet const &mset);

inline MultiIndexSet createActiveTensors(const MultiIndexSet &mset, const std::vector<int> &weights){
    size_t num_dimensions = mset.getNumDimensions();
    size_t nz_weights = 0;
    for(auto w: weights) if (w != 0) nz_weights++;
 
    std::vector<int> indexes(nz_weights * num_dimensions);
    nz_weights = 0;
    auto iter = indexes.begin();
    auto iset = mset.begin();
    for(auto w: weights){
        if (w != 0){
            std::copy_n(iset, num_dimensions, iter);
            std::advance(iter, num_dimensions);
        }
        std::advance(iset, num_dimensions);
    }
 
    return MultiIndexSet(num_dimensions, std::move(indexes));
}

inline void computeActiveTensorsWeights(MultiIndexSet const &tensors, MultiIndexSet &active_tensors, std::vector<int> &active_w){
    std::vector<int> tensors_w = MultiIndexManipulations::computeTensorWeights(tensors);
    active_tensors = MultiIndexManipulations::createActiveTensors(tensors, tensors_w);
 
    active_w = std::vector<int>();
    active_w.reserve(active_tensors.getNumIndexes());
    for(auto w : tensors_w) if (w != 0) active_w.push_back(w);
}

MultiIndexSet createPolynomialSpace(const MultiIndexSet &tensors, std::function<int(int)> exactness);

inline bool isLowerComplete(std::vector<int> const &point, MultiIndexSet const &mset, std::vector<int> &scratch){
    std::copy(point.begin(), point.end(), scratch.begin());
    for(int &d : scratch){
        if (d > 0){
            d--;
            if (mset.missing(scratch)) return false;
            d++;
        }
    }
    return true;
}

inline MultiIndexSet getLargestCompletion(MultiIndexSet const &current, MultiIndexSet const &candidates){
    if (candidates.empty()) return MultiIndexSet();
    auto num_dimensions = candidates.getNumDimensions();
    MultiIndexSet result; // start with an empty set
    if (current.empty()){
        if (candidates.missing(std::vector<int>(num_dimensions, 0))){
            return MultiIndexSet(); // current is empty, 0-th index is the only thing that can be added
        }else{
            result = MultiIndexSet(num_dimensions, std::vector<int>(num_dimensions, 0));
        }
    }
 
    std::vector<int> kid(num_dimensions);
    std::vector<int> scratch(num_dimensions);
    bool loopon = true;
    while(loopon){
        Data2D<int> update(num_dimensions, 0);
 
        MultiIndexSet total = current;
        if (!result.empty()) total += result;
        for(int i=0; i<total.getNumIndexes(); i++){
            std::copy_n(total.getIndex(i), num_dimensions, kid.begin());
            for(int &k : kid){
                k++; // construct the kid in the new direction
                if (!candidates.missing(kid) && result.missing(kid) && isLowerComplete(kid, total, scratch))
                    update.appendStrip(kid);
                k--;
            }
        }
 
        loopon = (update.getNumStrips() > 0);
        if (loopon) result += update;
    }
    return result;
}
 

template<bool limited>
MultiIndexSet addExclusiveChildren(const MultiIndexSet &tensors, const MultiIndexSet &exclude, const std::vector<int> level_limits){
    int num_dimensions = (int) tensors.getNumDimensions();
    Data2D<int> tens(num_dimensions, 0);
    std::vector<int> scratch(num_dimensions);
    for(int i=0; i<tensors.getNumIndexes(); i++){ // add the new tensors (so long as they are not included in the initial grid)
        const int *t = tensors.getIndex(i);
        std::vector<int> kid(t, t + num_dimensions);
        auto ilimit = level_limits.begin();
        for(auto &k : kid){
            k++;
            if (exclude.missing(kid) && tensors.missing(kid)){ // if the kid is not to be excluded and if not included in the current set
                if (isLowerComplete(kid, tensors, scratch)){
                    if (limited){
                        if ((*ilimit == -1) || (k <= *ilimit))
                            tens.appendStrip(kid);
                        ilimit++;
                    }else{
                        tens.appendStrip(kid);
                    }
                }
            }
            k--;
        }
    }
 
    return MultiIndexSet(tens);
}

template<class IndexList, class RuleLike, class OutputIteratorLike>
OutputIteratorLike indexesToNodes(IndexList const &list, RuleLike const &rule, OutputIteratorLike nodes){
    return std::transform(list.begin(), list.end(), nodes, [&](int i)->double{ return rule.getNode(i); });
}

template<class IteratorLike, class RuleLike, class OutputIteratorLike>
OutputIteratorLike indexesToNodes(IteratorLike ibegin, size_t num_entries, RuleLike const &rule, OutputIteratorLike nodes){
    return std::transform(ibegin, ibegin + num_entries, nodes, [&](int i)->double{ return rule.getNode(i); });
}

template<class IndexList, class RuleLike>
std::vector<double> getIndexesToNodes(IndexList const &list, RuleLike const &rule){
    std::vector<double> result(std::distance(list.begin(), list.end()));
    indexesToNodes(list, rule, result.begin());
    return result;
}

template<class IteratorLike, class RuleLike>
std::vector<double> getIndexesToNodes(IteratorLike ibegin, size_t num_entries, RuleLike const &rule){
    std::vector<double> result(num_entries);
    indexesToNodes(ibegin, num_entries, rule, result.begin());
    return result;
}

std::vector<int> inferAnisotropicWeights(AccelerationContext const *acceleration, TypeOneDRule rule, TypeDepth depth,
                                         MultiIndexSet const &points, std::vector<double> const &coefficients, double tol);
 
}
 
}
 
#endif