<commit_before>namespace sirius
{

class SHT
{
    private:

        int lmax_;

        int lmmax_;

        int num_points_;

        mdarray<double, 2> coord_;

        mdarray<double, 2> tp_;

        std::vector<double> w_;

        /// backward transformation from Ylm to spherical coordinates
        mdarray<complex16, 2> ylm_backward_;
        
        /// forward transformation from spherical coordinates to Ylm
        mdarray<complex16, 2> ylm_forward_;
        
        /// backward transformation from Rlm to spherical coordinates
        mdarray<double, 2> rlm_backward_;

        /// forward transformation from spherical coordinates to Rlm
        mdarray<double, 2> rlm_forward_;

        mdarray<complex16, 2> ylm_dot_rlm_;

        mdarray<int, 2> l_m_by_lm_;

    public:
        
        /// Default constructor
        SHT()
        {

        }

        void set_lmax(int lmax__)
        {
            lmax_ = lmax__;
            lmmax_ = (lmax_ + 1) * (lmax_ + 1);

            num_points_ = Lebedev_Laikov_npoint(2 * lmax_);
            //num_points_ = lmmax_;
            
            std::vector<double> x(num_points_);
            std::vector<double> y(num_points_);
            std::vector<double> z(num_points_);

            coord_.set_dimensions(3, num_points_);
            coord_.allocate();

            tp_.set_dimensions(2, num_points_);
            tp_.allocate();

            w_.resize(num_points_);

            Lebedev_Laikov_sphere(num_points_, &x[0], &y[0], &z[0], &w_[0]);
            //uniform_coverage();

            ylm_backward_.set_dimensions(lmmax_, num_points_);
            ylm_backward_.allocate();

            ylm_forward_.set_dimensions(num_points_, lmmax_);
            ylm_forward_.allocate();

            rlm_backward_.set_dimensions(lmmax_, num_points_);
            rlm_backward_.allocate();

            rlm_forward_.set_dimensions(num_points_, lmmax_);
            rlm_forward_.allocate();

            for (int itp = 0; itp < num_points_; itp++)
            {
                coord_(0, itp) = x[itp];
                coord_(1, itp) = y[itp];
                coord_(2, itp) = z[itp];
                
                double vs[3];

                spherical_coordinates(&coord_(0, itp), vs);

                //double t = tp_(0, itp);
                //double p = tp_(1, itp);

                //coord_(0, itp) = sin(t) * cos(p);
                //coord_(1, itp) = sin(t) * sin(p);
                //coord_(2, itp) = cos(t);

                spherical_harmonics(lmax_, vs[1], vs[2], &ylm_backward_(0, itp));
                spherical_harmonics(lmax_, vs[1], vs[2], &rlm_backward_(0, itp));

                for (int lm = 0; lm < lmmax_; lm++)
                {
                    ylm_forward_(itp, lm) = conj(ylm_backward_(lm, itp)) * w_[itp] * fourpi;
                    rlm_forward_(itp, lm) = rlm_backward_(lm, itp) * w_[itp] * fourpi;
                    //ylm_forward_(lm, itp) = ylm_backward_(lm, itp);
                    //rlm_forward_(lm, itp) = rlm_backward_(lm, itp);
                }
            }

            //linalg<lapack>::invert_ge(&ylm_forward_(0, 0), lmmax_);
            //linalg<lapack>::invert_ge(&rlm_forward_(0, 0), lmmax_);

            double dr = 0;
            double dy = 0;

            for (int lm = 0; lm < lmmax_; lm++)
            {
                for (int lm1 = 0; lm1 < lmmax_; lm1++)
                {
                    double t = 0;
                    complex16 zt(0, 0);
                    for (int itp = 0; itp < num_points_; itp++)
                    {
                        zt += ylm_forward_(itp, lm) * ylm_backward_(lm1, itp);
                        t += rlm_forward_(itp, lm) * rlm_backward_(lm1, itp);
                    }
                    
                    if (lm == lm1) 
                    {
                        zt -= 1.0;
                        t -= 1.0;
                    }<commit_msg>small update of SHT<commit_after>namespace sirius
{

class SHT
{
    private:

        int lmax_;

        int lmmax_;

        int num_points_;

        mdarray<double, 2> coord_;

        mdarray<double, 2> tp_;

        std::vector<double> w_;

        /// backward transformation from Ylm to spherical coordinates
        mdarray<complex16, 2> ylm_backward_;
        
        /// forward transformation from spherical coordinates to Ylm
        mdarray<complex16, 2> ylm_forward_;
        
        /// backward transformation from Rlm to spherical coordinates
        mdarray<double, 2> rlm_backward_;

        /// forward transformation from spherical coordinates to Rlm
        mdarray<double, 2> rlm_forward_;

        mdarray<complex16, 2> ylm_dot_rlm_;

        mdarray<int, 2> l_m_by_lm_;
        
        int mesh_type_;

    public:
        
        /// Default constructor
        SHT() : mesh_type_(0)
        {

        }

        void set_lmax(int lmax__)
        {
            lmax_ = lmax__;
            lmmax_ = (lmax_ + 1) * (lmax_ + 1);
            
            if (mesh_type_ == 0) num_points_ = Lebedev_Laikov_npoint(2 * lmax_);
            if (mesh_type_ == 1) num_points_ = lmmax_;
            
            std::vector<double> x(num_points_);
            std::vector<double> y(num_points_);
            std::vector<double> z(num_points_);

            coord_.set_dimensions(3, num_points_);
            coord_.allocate();

            tp_.set_dimensions(2, num_points_);
            tp_.allocate();

            w_.resize(num_points_);

            if (mesh_type_ == 0) Lebedev_Laikov_sphere(num_points_, &x[0], &y[0], &z[0], &w_[0]);
            if (mesh_type_ == 1) uniform_coverage();

            ylm_backward_.set_dimensions(lmmax_, num_points_);
            ylm_backward_.allocate();

            ylm_forward_.set_dimensions(num_points_, lmmax_);
            ylm_forward_.allocate();

            rlm_backward_.set_dimensions(lmmax_, num_points_);
            rlm_backward_.allocate();

            rlm_forward_.set_dimensions(num_points_, lmmax_);
            rlm_forward_.allocate();

            for (int itp = 0; itp < num_points_; itp++)
            {
                if (mesh_type_ == 0)
                {
                    coord_(0, itp) = x[itp];
                    coord_(1, itp) = y[itp];
                    coord_(2, itp) = z[itp];
                    
                    double vs[3];

                    spherical_coordinates(&coord_(0, itp), vs);
                    spherical_harmonics(lmax_, vs[1], vs[2], &ylm_backward_(0, itp));
                    spherical_harmonics(lmax_, vs[1], vs[2], &rlm_backward_(0, itp));
                    for (int lm = 0; lm < lmmax_; lm++)
                    {
                        ylm_forward_(itp, lm) = conj(ylm_backward_(lm, itp)) * w_[itp] * fourpi;
                        rlm_forward_(itp, lm) = rlm_backward_(lm, itp) * w_[itp] * fourpi;
                    }
                }
                if (mesh_type_ == 1)
                {
                    double t = tp_(0, itp);
                    double p = tp_(1, itp);

                    coord_(0, itp) = sin(t) * cos(p);
                    coord_(1, itp) = sin(t) * sin(p);
                    coord_(2, itp) = cos(t);

                    spherical_harmonics(lmax_, t, p, &ylm_backward_(0, itp));
                    spherical_harmonics(lmax_, t, p, &rlm_backward_(0, itp));

                    for (int lm = 0; lm < lmmax_; lm++)
                    {
                        ylm_forward_(lm, itp) = ylm_backward_(lm, itp);
                        rlm_forward_(lm, itp) = rlm_backward_(lm, itp);
                    }
                }
            }

            if (mesh_type_ == 1)
            {
                linalg<lapack>::invert_ge(&ylm_forward_(0, 0), lmmax_);
                linalg<lapack>::invert_ge(&rlm_forward_(0, 0), lmmax_);
            }

            double dr = 0;
            double dy = 0;

            for (int lm = 0; lm < lmmax_; lm++)
            {
                for (int lm1 = 0; lm1 < lmmax_; lm1++)
                {
                    double t = 0;
                    complex16 zt(0, 0);
                    for (int itp = 0; itp < num_points_; itp++)
                    {
                        zt += ylm_forward_(itp, lm) * ylm_backward_(lm1, itp);
                        t += rlm_forward_(itp, lm) * rlm_backward_(lm1, itp);
                    }
                    
                    if (lm == lm1) 
                    {
                        zt -= 1.0;
                        t -= 1.0;
                    }
{'fraction_non_alphanumeric': 0.11421884142773552, 'fraction_numerical': 0.01650087770626097, 'mean_word_length': 1.3267084127416282, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7403850', 'n_tokens_mistral': 2861, 'n_tokens_neox': 2628, 'n_words': 540}