<commit_before>/*
Copyright (c) 2015-2017 by the parties listed in the AUTHORS file.
All rights reserved.  Use of this source code is governed by
a BSD-style license that can be found in the LICENSE file.
*/

#include <toast_tod_internal.hpp>

#include <sstream>
#include <iostream>
#include <iomanip>


void toast::filter::polyfilter(
    const long order, double **signals, uint8_t *flags,
    const size_t n, const size_t nsignal,
    const long *starts, const long *stops, const size_t nscan ) {

    // Process the signals, one subscan at a time.  There is only one
    // flag vector, because the flags must be identical to apply the
    // same template matrix.

    char upper = 'U';
    char lower = 'L';
    char notrans = 'N';
    char trans = 'T';
    char diag = 'U';
    char nodiag = 'N';
    double fzero = 0.0;
    double fone = 1.0;
    int one = 1;
    int zero = 0;

    int norder = order + 1;

#pragma omp parallel for schedule(static)
    for ( int iscan=0; iscan<nscan; ++iscan ) {
        int start = starts[iscan];
        int stop = stops[iscan];
        if ( start < 0 ) start = 0;
        if ( stop > n-1 ) stop = n-1;
        if ( stop < start ) continue;
        int scanlen = stop - start + 1;

        // Build the template matrix

        double *templates = static_cast< double* >(
            toast::mem::aligned_alloc (
                scanlen*norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double dx = 2. / scanlen;
        double xstart = 0.5*dx - 1;
        for ( int i=0; i<scanlen; ++i ) {
            if ( flags[start+i] ) continue;
            double x = xstart + i*dx;
            int offset = i*norder;
            if ( norder > 0 ) templates[offset++] = 1;
            if ( norder > 1 ) templates[offset++] = x;
            for ( int iorder=1; iorder<norder-1; ++iorder ) {
                templates[offset] =
                    ((2*iorder+1)*x*templates[offset-1]
                     - iorder*templates[offset-2]) / (iorder+1);
                ++offset;
            }
        }

        double *cov = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        // invcov = templates x templates.T

        toast::lapack::syrk( &upper, &notrans, &norder, &scanlen, &fone,
                             templates, &norder, &fzero, cov,
                             &norder );

        // invert cov = invcov^-1 using LU decomposition

        int info = 0;
        toast::lapack::potrf( &upper, &norder, cov, &norder, &info );
        if ( info == 0 ) {
            toast::lapack::potri( &upper, &norder, cov, &norder, &info );
        }

        if ( info ) {
            // The matrix is singular.  Raise the appropriate quality
            // flags.
            for ( int i=0; i<scanlen; ++i ) {
                flags[start+i] = 255;
            }
            continue;
        }

        // Symmetrize for dgemv later on (workaround for dtrmv issue)

        for ( int row=0; row<norder; ++row ) {
          for ( int col=row+1; col<norder; ++col ) {
            cov[row*norder + col] = cov[col*norder + row];
          }
        }

        // Filter every signal

        double *proj = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double *coeff = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double *noise = static_cast< double* >(
            toast::mem::aligned_alloc (
                scanlen*sizeof(double), toast::mem::SIMD_ALIGN ) );

        for ( int isignal=0; isignal<nsignal; ++isignal ) {
            double *signal = signals[isignal] + start;

            // proj = templates x signal

            toast::lapack::gemv( &notrans, &norder, &scanlen, &fone, templates,
                                 &norder, signal, &one, &fzero, proj, &one );

            // coeff = cov x proj

            // For whatever reason, trmv refused to yield the right answer...
            // Use general matrix vector multiply instead

            //toast::lapack::trmv( &upper, &trans, &nodiag, &norder, cov,
            //                     &norder, coeff, &one );

            toast::lapack::gemv( &trans, &norder, &norder, &fone,
                                 cov, &norder, proj, &one,
                                 &fzero, coeff, &one );

            // noise = templates.T x coeff

            toast::lapack::gemv( &trans, &norder, &scanlen, &fone, templates,
                                 &norder, coeff, &one, &fzero, noise, &one );

            // Subtract noise

            for ( int i=0; i<scanlen; ++i ) signal[i] -= noise[i];
        }

        // Free workspace

        toast::mem::aligned_free( templates );
        toast::mem::aligned_free( cov );
        toast::mem::aligned_free( proj );
        toast::mem::aligned_free( coeff );
        toast::mem::aligned_free( noise );

    }

    return;
}
<commit_msg>Modify the polynomial filter to handle flagging more intuitively.<commit_after>/*
Copyright (c) 2015-2017 by the parties listed in the AUTHORS file.
All rights reserved.  Use of this source code is governed by
a BSD-style license that can be found in the LICENSE file.
*/

#include <toast_tod_internal.hpp>

#include <sstream>
#include <iostream>
#include <iomanip>


void toast::filter::polyfilter(
    const long order, double **signals, uint8_t *flags,
    const size_t n, const size_t nsignal,
    const long *starts, const long *stops, const size_t nscan ) {

    // Process the signals, one subscan at a time.  There is only one
    // flag vector, because the flags must be identical to apply the
    // same template matrix.

    char upper = 'U';
    char lower = 'L';
    char notrans = 'N';
    char trans = 'T';
    char diag = 'U';
    char nodiag = 'N';
    double fzero = 0.0;
    double fone = 1.0;
    int one = 1;
    int zero = 0;

    int norder = order + 1;

#pragma omp parallel for schedule(static)
    for ( int iscan=0; iscan<nscan; ++iscan ) {
        int start = starts[iscan];
        int stop = stops[iscan];
        if ( start < 0 ) start = 0;
        if ( stop > n-1 ) stop = n-1;
        if ( stop < start ) continue;
        int scanlen = stop - start + 1;

        // Build the full template matrix used to clean the signal.
        // We subtract the template value even from flagged samples to
        // support point source masking etc.

        double *full_templates = static_cast< double* >(
            toast::mem::aligned_alloc (
                scanlen*norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double dx = 2. / scanlen;
        double xstart = 0.5*dx - 1;
        for ( int i=0; i<scanlen; ++i ) {
            double x = xstart + i*dx;
            int offset = i*norder;
            if ( norder > 0 ) full_templates[offset++] = 1;
            if ( norder > 1 ) full_templates[offset++] = x;
            for ( int iorder=1; iorder<norder-1; ++iorder ) {
                full_templates[offset] =
                    ((2*iorder+1)*x*full_templates[offset-1]
                     - iorder*full_templates[offset-2]) / (iorder+1);
                ++offset;
            }
        }

        // Assemble the flagged template matrix used in the linear regression

        double *templates = static_cast< double* >(
            toast::mem::aligned_alloc (
                scanlen*norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        for ( int i=0; i<scanlen; ++i ) {
            if ( flags[start+i] ) continue;
            for ( int offset=i*norder; offset<(i+1)*norder; ++offset )
                templates[offset] = full_templates[offset];
        }

        double *cov = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        // invcov = templates x templates.T

        toast::lapack::syrk( &upper, &notrans, &norder, &scanlen, &fone,
                             templates, &norder, &fzero, cov,
                             &norder );

        // invert cov = invcov^-1 using LU decomposition

        int info = 0;
        toast::lapack::potrf( &upper, &norder, cov, &norder, &info );
        if ( info == 0 ) {
            toast::lapack::potri( &upper, &norder, cov, &norder, &info );
        }

        if ( info ) {
            // The matrix is singular.  Raise the appropriate quality
            // flags.
            for ( int i=0; i<scanlen; ++i ) {
                flags[start+i] = 255;
            }
            continue;
        }

        // Symmetrize for dgemv later on (workaround for dtrmv issue)

        for ( int row=0; row<norder; ++row ) {
          for ( int col=row+1; col<norder; ++col ) {
            cov[row*norder + col] = cov[col*norder + row];
          }
        }

        // Filter every signal

        double *proj = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double *coeff = static_cast< double* >(
            toast::mem::aligned_alloc (
                norder*sizeof(double), toast::mem::SIMD_ALIGN ) );

        double *noise = static_cast< double* >(
            toast::mem::aligned_alloc (
                scanlen*sizeof(double), toast::mem::SIMD_ALIGN ) );

        for ( int isignal=0; isignal<nsignal; ++isignal ) {
            double *signal = signals[isignal] + start;

            // proj = templates x signal

            toast::lapack::gemv( &notrans, &norder, &scanlen, &fone, templates,
                                 &norder, signal, &one, &fzero, proj, &one );

            // coeff = cov x proj

            // For whatever reason, trmv refused to yield the right answer...
            // Use general matrix vector multiply instead

            //toast::lapack::trmv( &upper, &trans, &nodiag, &norder, cov,
            //                     &norder, coeff, &one );

            toast::lapack::gemv( &trans, &norder, &norder, &fone,
                                 cov, &norder, proj, &one,
                                 &fzero, coeff, &one );

            // noise = templates.T x coeff

            toast::lapack::gemv( &trans, &norder, &scanlen, &fone,
                                 full_templates, &norder, coeff, &one, &fzero,
                                 noise, &one );

            // Subtract noise

            for ( int i=0; i<scanlen; ++i ) signal[i] -= noise[i];
        }

        // Free workspace

        toast::mem::aligned_free( full_templates );
        toast::mem::aligned_free( templates );
        toast::mem::aligned_free( cov );
        toast::mem::aligned_free( proj );
        toast::mem::aligned_free( coeff );
        toast::mem::aligned_free( noise );

    }

    return;
}

{'fraction_non_alphanumeric': 0.1271326362135388, 'fraction_numerical': 0.008989176297926987, 'mean_word_length': 1.756066734074823, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7466482', 'n_tokens_mistral': 3398, 'n_tokens_neox': 3287, 'n_words': 922}