tsgLoadUnstructuredPoints.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 __TASMANIAN_ADDONS_LOADUNSTRUCTURED_HPP
#define __TASMANIAN_ADDONS_LOADUNSTRUCTURED_HPP

 
#include "tsgLoadNeededValues.hpp"

 
namespace TasGrid{

inline bool hasGPUBasis(TasmanianSparseGrid const &grid){
    return AccelerationMeta::isAvailable(accel_gpu_cuda)
        and not (grid.isLocalPolynomial() and ((grid.getOrder() < 0) or (grid.getOrder() > 2)))
        and not (grid.isWavelet() and grid.getOrder() == 3);
}

template<typename scalar_type>
void generateCoefficientsGPU(double const data_points[], int num_data, scalar_type model_values[],
                             double tolerance, TasmanianSparseGrid &grid){
    AccelerationContext const *acceleration = grid.getAccelerationContext();
    int num_outputs = grid.getNumOutputs();
    int num_points = grid.getNumPoints();
    int num_equations = (tolerance > 0.0) ? num_data + num_points : num_data;
 
    GpuVector<scalar_type> basis_matrix(acceleration, num_equations, num_points);
    grid.evaluateHierarchicalFunctionsGPU(data_points, num_data, reinterpret_cast<double*>(basis_matrix.data()));
 
    if (tolerance > 0.0){
        double correction = std::sqrt(tolerance);
        constexpr long long esize = (std::is_same<scalar_type, double>::value) ? 1 : 2;
        long long num_total = static_cast<long long>(num_points) * static_cast<long long>(num_points) * esize;
        TasGpu::fillDataGPU(acceleration, 0.0, num_total, 1,
                            reinterpret_cast<double*>(basis_matrix.data() + Utils::size_mult(num_data, num_points)));
 
        long long stride = static_cast<long long>(num_points + 1) * esize;
        TasGpu::fillDataGPU(acceleration, correction, num_points, stride,
                            reinterpret_cast<double*>(basis_matrix.data() + Utils::size_mult(num_data, num_points)));
 
        TasGpu::fillDataGPU(acceleration, 0.0, esize * num_outputs * num_points, 1,
                            reinterpret_cast<double*>(model_values + Utils::size_mult(num_data, num_outputs)));
    }
 
    TasmanianDenseSolver::solvesLeastSquaresGPU(acceleration, num_equations, num_points,
                                                basis_matrix.data(), num_outputs, model_values);
}

template<typename scalar_type>
Data2D<scalar_type> generateCoefficients(double const data_points[], int num_data, double const model_values[], double tolerance, TasmanianSparseGrid &grid){
    int num_dimensions = grid.getNumDimensions();
    int num_outputs = grid.getNumOutputs();
    int num_points = grid.getNumPoints();
    int num_equations = (tolerance > 0.0) ? num_data + num_points : num_data;
 
    AccelerationContext const *acceleration = grid.getAccelerationContext();
 
    Data2D<scalar_type> basis_matrix(num_points, num_equations, 0.0);
    Data2D<scalar_type> coefficients(num_outputs, num_equations, 0.0);
 
    if (acceleration->mode == accel_gpu_magma and hasGPUBasis(grid)){
        int mem_usage = 268435456 / ((std::is_same<double, scalar_type>::value) ? 1 : 2); // 2GB of double numbers
        if (grid.getGPUMemory(acceleration->device) < 2048) mem_usage /= 2; // use minimum of 1GB
        int num_batch = std::min(num_data, mem_usage / num_points);
        if (num_batch < 128){
            // cannot use the GPU with less than 128 entries, might as well use the CPU
            grid.evaluateHierarchicalFunctions(data_points, num_data, reinterpret_cast<double*>(basis_matrix.data()));
        }else{
            acceleration->setDevice();
            GpuVector<double> gpu_points(acceleration, num_dimensions, num_batch);
            GpuVector<scalar_type> gpu_matrix(acceleration, num_points, num_batch);
            for(int i = 0; i < num_data; i += num_batch){
                int num_this_batch = std::min(num_batch, num_data - i);
                gpu_points.load(acceleration, Utils::size_mult(num_this_batch, num_dimensions),
                                data_points + Utils::size_mult(i, num_dimensions));
                grid.evaluateHierarchicalFunctionsGPU(gpu_points.data(), num_this_batch, reinterpret_cast<double*>(gpu_matrix.data()));
                gpu_matrix.unload(acceleration, Utils::size_mult(num_this_batch, num_points), basis_matrix.getStrip(i));
            }
        }
    }else{
        grid.evaluateHierarchicalFunctions(data_points, num_data, reinterpret_cast<double*>(basis_matrix.data()));
    }
 
    if (tolerance > 0.0){
        double correction = std::sqrt(tolerance);
        for(int i=0; i<grid.getNumPoints(); i++)
            basis_matrix.getStrip(i + num_data)[i] = correction;
    }
 
    auto icoeff = coefficients.begin();
    for(size_t i=0; i<Utils::size_mult(num_data, grid.getNumOutputs()); i++)
        *icoeff++ = model_values[i];
 
    TasmanianDenseSolver::solvesLeastSquares(acceleration, num_equations, grid.getNumPoints(),
                                             basis_matrix.data(), grid.getNumOutputs(), coefficients.data());
 
    return coefficients;
}

template<typename scalar_type>
inline void loadUnstructuredDataL2tmpl(double const data_points[], int num_data, double const model_values[],
                                       double tolerance, TasmanianSparseGrid &grid){
 
    if (grid.empty()) throw std::runtime_error("Cannot use loadUnstructuredDataL2() with an empty grid.");
    if (grid.getNumNeeded() != 0)
        grid.mergeRefinement();
 
    AccelerationContext const *acceleration = grid.getAccelerationContext();
    if (acceleration->mode == accel_none)
        throw std::runtime_error("The loadUnstructuredDataL2() method cannot be used with acceleration mode accel_none.");
 
    int num_dimensions = grid.getNumDimensions();
    int num_outputs = grid.getNumOutputs();
    int num_points = grid.getNumPoints();
    int num_equations = (tolerance > 0.0) ? num_data + num_points : num_data;
 
    Data2D<scalar_type> coefficients =
        [&]()->Data2D<scalar_type>{
            if (acceleration->mode == accel_gpu_cuda and hasGPUBasis(grid)){
                acceleration->setDevice();
                GpuVector<double> gpu_points(acceleration, num_dimensions, num_data, data_points);
                GpuVector<scalar_type> gpu_values(acceleration, num_outputs, num_equations);
                TasGpu::load_n(acceleration, model_values, Utils::size_mult(num_outputs, num_data), gpu_values.data());
                generateCoefficientsGPU<scalar_type>(gpu_points.data(), num_data, gpu_values.data(), tolerance, grid);
                return Data2D<scalar_type>(num_outputs, num_equations, gpu_values.unload(acceleration));
            }else{
                return generateCoefficients<scalar_type>(data_points, num_data, model_values, tolerance, grid);
            }
        }();
 
    // do the real-complex split (used in Fourier grids)
    if (std::is_same<scalar_type, std::complex<double>>::value){
        std::vector<double> real_coeffs(Utils::size_mult(2 * grid.getNumOutputs(), grid.getNumPoints()));
        auto icoeff = coefficients.begin();
        for(size_t i=0; i<Utils::size_mult(grid.getNumOutputs(), grid.getNumPoints()); i++)
            real_coeffs[i] = std::real(*icoeff++);
        icoeff = coefficients.begin();
        for(size_t i=Utils::size_mult(grid.getNumOutputs(), grid.getNumPoints()); i<Utils::size_mult(2 * grid.getNumOutputs(), grid.getNumPoints()); i++)
            real_coeffs[i] = std::imag(*icoeff++);
        grid.setHierarchicalCoefficients(real_coeffs.data());
    }else{
        grid.setHierarchicalCoefficients(reinterpret_cast<double*>(coefficients.data()));
    }
}

inline void loadUnstructuredDataL2(double const data_points[], int num_data, double const model_values[],
                                   double tolerance, TasmanianSparseGrid &grid){
    if (grid.isFourier()){
        loadUnstructuredDataL2tmpl<std::complex<double>>(data_points, num_data, model_values, tolerance, grid);
    }else{
        loadUnstructuredDataL2tmpl<double>(data_points, num_data, model_values, tolerance, grid);
    }
}

inline void loadUnstructuredDataL2(std::vector<double> const &data_points, std::vector<double> const &model_values,
                                   double tolerance, TasmanianSparseGrid &grid){
    if (grid.empty()) throw std::runtime_error("Cannot use loadUnstructuredDataL2() with an empty grid.");
    int num_data = static_cast<int>(data_points.size() / grid.getNumDimensions());
    if (model_values.size() < Utils::size_mult(num_data, grid.getNumOutputs()))
        throw std::runtime_error("In loadUnstructuredDataL2(), provided more points than data.");
    loadUnstructuredDataL2(data_points.data(), num_data, model_values.data(), tolerance, grid);
}
 
}
 
#endif