Status.hpp

Go to the documentation of this file.

/*
 *
 * Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
 * All rights reserved.
 *
 * 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. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *    This product includes software developed by the Delft University of Technology.
 * 4. Neither the name of the Delft University of Technology 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 LAURENS VAN DER MAATEN ''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 LAURENS VAN DER MAATEN 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.
 *
 */
 
#ifndef QDTSNE_STATUS_HPP
#define QDTSNE_STATUS_HPP
 
#include <vector>
#include <cmath>
#include <algorithm>
#include <type_traits>
#include <limits>
 
#include "SPTree.hpp"
#include "Options.hpp"
#include "utils.hpp"
 
namespace qdtsne {
 
template<int num_dim_, typename Index_, typename Float_> 
class Status {
public:
    Status(NeighborList<Index_, Float_> neighbors, Options options) :
        my_neighbors(std::move(neighbors)),
        my_dY(my_neighbors.size() * num_dim_), 
        my_uY(my_neighbors.size() * num_dim_), 
        my_gains(my_neighbors.size() * num_dim_, 1.0), 
        my_pos_f(my_neighbors.size() * num_dim_), 
        my_neg_f(my_neighbors.size() * num_dim_), 
        my_tree(my_neighbors.size(), options.max_depth),
        my_options(std::move(options))
    {
        if (options.num_threads > 1) {
            my_parallel_buffer.resize(my_neighbors.size());
        }
    }
private:
    NeighborList<Index_, Float_> my_neighbors; 
    std::vector<Float_> my_dY, my_uY, my_gains, my_pos_f, my_neg_f;
 
    internal::SPTree<num_dim_, Float_> my_tree;
    std::vector<Float_> my_parallel_buffer; // Buffer to hold parallel-computed results prior to reduction.
 
    Options my_options;
    int my_iter = 0;
 
    typename decltype(my_tree)::LeafApproxWorkspace my_leaf_workspace;
 
public:
    int iteration() const {
        return my_iter;
    }
 
    int max_iterations() const {
        return my_options.max_iterations;
    }
 
    size_t num_observations() const {
        return my_neighbors.size();
    }
 
#ifndef NDEBUG
    const auto& get_neighbors() const {
        return my_neighbors;
    }
#endif
 
public:
    void run(Float_* Y, int limit) {
        Float_ multiplier = (my_iter < my_options.stop_lying_iter ? my_options.exaggeration_factor : 1);
        Float_ momentum = (my_iter < my_options.mom_switch_iter ? my_options.start_momentum : my_options.final_momentum);
 
        for(; my_iter < limit; ++my_iter) {
            // Stop lying about the P-values after a while, and switch momentum
            if (my_iter == my_options.stop_lying_iter) {
                multiplier = 1;
            }
            if (my_iter == my_options.mom_switch_iter) {
                momentum = my_options.final_momentum;
            }
 
            iterate(Y, multiplier, momentum);
        }
    }
 
    void run(Float_* Y) {
        run(Y, my_options.max_iterations);
    }
 
private:
    static Float_ sign(Float_ x) { 
        constexpr Float_ zero = 0;
        constexpr Float_ one = 1;
        return (x == zero ? zero : (x < zero ? -one : one));
    }
 
    void iterate(Float_* Y, Float_ multiplier, Float_ momentum) {
        compute_gradient(Y, multiplier);
 
        // Update gains
        size_t ngains = my_gains.size(); 
        for (size_t i = 0; i < ngains; ++i) {
            Float_& g = my_gains[i];
            constexpr Float_ lower_bound = 0.01;
            constexpr Float_ to_add = 0.2;
            constexpr Float_ to_mult = 0.8;
            g = std::max(lower_bound, sign(my_dY[i]) != sign(my_uY[i]) ? (g + to_add) : (g * to_mult));
        }
 
        // Perform gradient update (with momentum and gains)
        for (size_t i = 0; i < ngains; ++i) {
            my_uY[i] = momentum * my_uY[i] - my_options.eta * my_gains[i] * my_dY[i];
            Y[i] += my_uY[i];
        }
 
        // Make solution zero-mean
        size_t N = num_observations();
        for (int d = 0; d < num_dim_; ++d) {
            auto start = Y + d;
 
            // Compute means from column-major coordinates.
            Float_ sum = 0;
            for (size_t i = 0; i < N; ++i, start += num_dim_) {
                sum += *start;
            }
            sum /= N;
 
            start = Y + d;
            for (size_t i = 0; i < N; ++i, start += num_dim_) {
                *start -= sum;
            }
        }
 
        return;
    }
 
private:
    void compute_gradient(const Float_* Y, Float_ multiplier) {
        my_tree.set(Y);
        compute_edge_forces(Y, multiplier);
 
        size_t N = num_observations();
        std::fill(my_neg_f.begin(), my_neg_f.end(), 0);
 
        Float_ sum_Q = compute_non_edge_forces();
 
        // Compute final t-SNE gradient
        size_t ntotal = N * static_cast<size_t>(num_dim_);
        for (size_t i = 0; i < ntotal; ++i) {
            my_dY[i] = my_pos_f[i] - (my_neg_f[i] / sum_Q);
        }
    }
 
    void compute_edge_forces(const Float_* Y, Float_ multiplier) {
        std::fill(my_pos_f.begin(), my_pos_f.end(), 0);
        size_t N = num_observations();
 
        parallelize(my_options.num_threads, N, [&](int, size_t start, size_t length) -> void {
            for (size_t n = start, end = start + length; n < end; ++n) {
                const auto& current = my_neighbors[n];
                size_t offset = n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
                const Float_* self = Y + offset;
                Float_* pos_out = my_pos_f.data() + offset;
 
                for (const auto& x : current) {
                    Float_ sqdist = 0; 
                    const Float_* neighbor = Y + static_cast<size_t>(x.first) * num_dim_; // cast to avoid overflow.
                    for (int d = 0; d < num_dim_; ++d) {
                        Float_ delta = self[d] - neighbor[d];
                        sqdist += delta * delta;
                    }
 
                    const Float_ mult = multiplier * x.second / (static_cast<Float_>(1) + sqdist);
                    for (int d = 0; d < num_dim_; ++d) {
                        pos_out[d] += mult * (self[d] - neighbor[d]);
                    }
                }
            }
        });
 
        return;
    }
 
    Float_ compute_non_edge_forces() {
        size_t N = num_observations();
 
        if (my_options.leaf_approximation) {
            my_tree.compute_non_edge_forces_for_leaves(my_options.theta, my_leaf_workspace, my_options.num_threads);
        }
 
        if (my_options.num_threads > 1) {
            // Don't use reduction methods, otherwise we get numeric imprecision
            // issues (and stochastic results) based on the order of summation.
            parallelize(my_options.num_threads, N, [&](int, size_t start, size_t length) -> void {
                for (size_t n = start, end = start + length; n < end; ++n) {
                    auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
                    if (my_options.leaf_approximation) {
                        my_parallel_buffer[n] = my_tree.compute_non_edge_forces_from_leaves(n, neg_ptr, my_leaf_workspace);
                    } else {
                        my_parallel_buffer[n] = my_tree.compute_non_edge_forces(n, my_options.theta, neg_ptr);
                    }
                }
            });
            return std::accumulate(my_parallel_buffer.begin(), my_parallel_buffer.end(), static_cast<Float_>(0));
        }
 
        Float_ sum_Q = 0;
        for (size_t n = 0; n < N; ++n) {
            auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
            if (my_options.leaf_approximation) {
                sum_Q += my_tree.compute_non_edge_forces_from_leaves(n, neg_ptr, my_leaf_workspace);
            } else {
                sum_Q += my_tree.compute_non_edge_forces(n, my_options.theta, neg_ptr);
            }
        }
        return sum_Q;
    }
};
 
}
 
#endif