dense_train_tools.hpp

Go to the documentation of this file.

#ifndef FSTDENSETRAINTOOLS_H
#define FSTDENSETRAINTOOLS_H
 /*
=========================================================================
Copyright and license:
  * FST4 software (with the exception of any bundled or externally linked
    third-party components) is Copyright (C) 2000-present Institute of
    Information Theory and Automation of the Czech Academy of Sciences (UTIA).
  * FST4 is free software: you can redistribute it and/or modify it under
    the terms of the GNU Affero General Public License as published by the
    Free Software Foundation, version 3. The full text of the license is in
    the file LICENSE in the distribution root and at
    https://www.gnu.org/licenses/agpl-3.0.html.
  * FST4 is distributed in the hope that it will be useful, but WITHOUT ANY
    WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public
    License for more details.
  * Closed-source use without the AGPL's copyleft obligations requires a
    commercial license from UTIA; in case of interest write to
    utia@utia.cas.cz. LICENSING.md in the distribution root explains both
    options, the contributor license agreement that every contribution
    requires, and the contacts (technical questions and pull requests:
    somol@utia.cas.cz, with "FST" in the subject line).
========================================================================== */
 
#include <cmath>
#include <cstddef>
#include <memory>
#include <vector>
#include "data_accessor.hpp"
 
namespace FST {

template<typename REALTYPE, typename IDXTYPE, typename DIMTYPE, class DATAACCESSOR, class AFTERPATTERN>
std::size_t standardize_design(const std::shared_ptr<DATAACCESSOR> &da, const DIMTYPE classes,
        const std::size_t d, const std::size_t zstride, const std::vector<DIMTYPE> &index,
        std::vector<REALTYPE> &Z, std::vector<REALTYPE> &mean, std::vector<REALTYPE> &scale,
        AFTERPATTERN after_pattern)
{
        const DIMTYPE stride=da->getNoOfFeatures();
        mean.assign(d,0.0); scale.assign(d,0.0);
        Z.clear();
        typename DATAACCESSOR::PPattern p;
        IDXTYPE s,i;
        std::size_t N=0;
        for(DIMTYPE c=0;c<classes;c++)
        {
                da->setClass(c);
                for(const auto &blk : da->blocks(TRAIN)) for(i=0,p=blk.data,s=blk.size;i<s;i++)
                {
                        const auto *pat=&p[(std::size_t)i*stride];
                        for(std::size_t f=0;f<d;f++) {
                                const REALTYPE v=(REALTYPE)pat[index[f]];
                                Z.push_back(v);
                                mean[f]+=v;
                        }
                        after_pattern(pat,c);
                        ++N;
                }
        }
        if(N==0) return 0;
        for(std::size_t f=0;f<d;f++) mean[f]/=(REALTYPE)N;
        // second pass over the stored raw values: centered sums of squares. The one-pass form
        // E[x^2]-mean^2 cancels catastrophically once |mean| >> sd (features offset by 1e8 changed
        // logreg/linsvm/rbfsvm/ridge results although standardization is translation-invariant);
        // the design matrix is in memory anyway, so the extra pass costs one sweep
        for(std::size_t n=0;n<N;n++)
                for(std::size_t f=0;f<d;f++) {const REALTYPE dev=Z[n*zstride+f]-mean[f]; scale[f]+=dev*dev;}
        for(std::size_t f=0;f<d;f++)
        {
                const REALTYPE var=scale[f]/(REALTYPE)N;
                scale[f]=(var>0.0)?std::sqrt(var):1.0; // constant feature -> scale 1, stays inert
        }
        for(std::size_t n=0;n<N;n++)
                for(std::size_t f=0;f<d;f++) Z[n*zstride+f]=(Z[n*zstride+f]-mean[f])/scale[f];
        return N;
}

template<typename REALTYPE>
bool cholesky_factor(std::vector<REALTYPE> &A, const std::size_t m)
{
        for(std::size_t j=0;j<m;j++)
        {
                REALTYPE diag=A[j*m+j];
                for(std::size_t k=0;k<j;k++) diag-=A[j*m+k]*A[j*m+k];
                if(!(diag>0.0) || !std::isfinite(diag)) return false;
                diag=std::sqrt(diag);
                A[j*m+j]=diag;
                for(std::size_t i=j+1;i<m;i++)
                {
                        REALTYPE v=A[i*m+j];
                        for(std::size_t k=0;k<j;k++) v-=A[i*m+k]*A[j*m+k];
                        A[i*m+j]=v/diag;
                }
        }
        return true;
}

template<typename REALTYPE>
void cholesky_solve(const std::vector<REALTYPE> &A, const std::size_t m,
        const REALTYPE *rhs, const bool negate_rhs, REALTYPE *out)
{
        for(std::size_t i=0;i<m;i++)
        {
                REALTYPE v=negate_rhs?-rhs[i]:rhs[i];
                for(std::size_t k=0;k<i;k++) v-=A[i*m+k]*out[k];
                out[i]=v/A[i*m+i];
        }
        for(std::size_t ii=m;ii>0;ii--)
        {
                const std::size_t i=ii-1;
                REALTYPE v=out[i];
                for(std::size_t k=i+1;k<m;k++) v-=A[k*m+i]*out[k];
                out[i]=v/A[i*m+i];
        }
}
 
} // namespace FST
#endif // FSTDENSETRAINTOOLS_H