<commit_before>/*
 constraint.cpp

 Copyright (c) 2014, 2015, 2016 <NAME>

 This file is distributed under the terms of the MIT license.
 Please see the file 'LICENCE.txt' in the root directory 
 or http://opensource.org/licenses/mit-license.php for information.
*/

#include "constraint.h"
#include "combination.h"
#include "constants.h"
#include "error.h"
#include "fcs.h"
#include "cluster.h"
#include "mathfunctions.h"
#include "memory.h"
#include "rref.h"
#include "symmetry.h"
#include "system.h"
#include "timer.h"
#include "xml_parser.h"
#include <iostream>
#include <iomanip>
#include <boost/bimap.hpp>
#include <algorithm>
#include <unordered_set>
#include <map>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/xml_parser.hpp>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>

using namespace ALM_NS;

Constraint::Constraint()
{
    set_default_variables();
}

Constraint::~Constraint()
{
    deallocate_variables();
}

void Constraint::set_default_variables()
{
    constraint_mode = 11;
    rotation_axis = "";
    fix_harmonic = false;
    fix_cubic = false;
    constraint_algebraic = 1;
    fc2_file = "";
    fc3_file = "";
    exist_constraint = false;
    extra_constraint_from_symmetry = false;
    const_mat = nullptr;
    const_rhs = nullptr;
    const_symmetry = nullptr;
    const_fix = nullptr;
    const_relate = nullptr;
    const_relate_rotation = nullptr;
    index_bimap = nullptr;
    number_of_constraints = 0;
    tolerance_constraint = eps8;
}

void Constraint::deallocate_variables()
{
    if (const_symmetry) {
        deallocate(const_symmetry);
    }
    if (const_fix) {
        deallocate(const_fix);
    }
    if (const_relate) {
        deallocate(const_relate);
    }
    if (const_relate_rotation) {
        deallocate(const_relate_rotation);
    }
    if (index_bimap) {
        deallocate(index_bimap);
    }
    if (const_mat) {
        deallocate(const_mat);
    }
    if (const_rhs) {
        deallocate(const_rhs);
    }
}

void Constraint::setup(const System *system,
                       const Fcs *fcs,
                       const Cluster *cluster,
                       const Symmetry *symmetry,
                       const int linear_model,
                       const int mirror_image_conv,
                       const int verbosity,
                       Timer *timer)
{
    timer->start_clock("constraint");

    if (verbosity > 0) {
        std::cout << " CONSTRAINT" << std::endl;
        std::cout << " ==========" << std::endl << std::endl;
    }

    constraint_algebraic = constraint_mode / 10;
    constraint_mode = constraint_mode % 10;
    const auto maxorder = cluster->get_maxorder();

    if (linear_model == 2) {
        if (constraint_mode > 1) {
            warn("Constraint::setup", "Sorry, only ICONST = 11 is supported \n"
                                      "                      when LMODEL = enet. We set ICONST = 11 in this run.\n");
            constraint_mode = 1;
        }
        constraint_algebraic = 1;
    }

    switch (constraint_mode) {
        case 0: // do nothing
            impose_inv_T = false;
            impose_inv_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 0: Constraint for translational/rotational invariance" << std::endl;
                std::cout << "              will NOT be considered." << std::endl;
            }
            break;
        case 1:
            impose_inv_T = true;
            impose_inv_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 1: Constraints for translational invariance" << std::endl;
                std::cout << "              will be considered." << std::endl;
            }
            break;
        case 2:
            impose_inv_T = true;
            impose_inv_R = true;
            exclude_last_R = true;
            if (verbosity > 0) {
                std::cout << "  ICONST = 2: Constraints for translational and rotational invariance" << std::endl;
                std::cout << "              will be considered. Axis of rotation is " << rotation_axis << std::endl;
                std::cout << "              Rotational invariance of the maximum order will be neglected" << std::endl;
            }
            break;
        case 3:
            impose_inv_T = true;
            impose_inv_R = true;
            exclude_last_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 3: Constraints for translational and rotational invariance" << std::endl;
                std::cout << "              will be considered. Axis of rotation is " << rotation_axis << std::endl;
            }
            break;
        default:
            exit("Constraint::setup", "invalid constraint_mode", constraint_mode);
            break;
    }

    if (fcs->get_forceconstant_basis() == "Lattice" && impose_inv_R) {
        exit("Constraint::setup()", "Sorry, rotational invariance with FCSYM_BASIS = Lattice is "
                                    "not supported.\n Use FCSYM_BASIS = Cartesian instead.");
    }

    if (verbosity > 0) std::cout << std::endl;

    if (const_symmetry) {
        deallocate(const_symmetry);
    }
    if (const_fix) {
        deallocate(const_fix);
    }
    if (const_relate) {
        deallocate(const_relate);
    }
    if (index_bimap) {
        deallocate(index_bimap);
    }

    allocate(const_fix, maxorder);
    allocate(const_relate, maxorder);
    allocate(index_bimap, maxorder);
    allocate(const_symmetry, maxorder);
    allocate(const_translation, maxorder);
    allocate(const_rotation_self, maxorder);
    allocate(const_rotation_cross, maxorder);
    allocate(const_self, maxorder);

    for (auto order = 0; order < maxorder; ++order) {
        const_translation[order].clear();
        const_rotation_self[order].clear();
        const_rotation_cross[order].clear();
        const_self[order].clear();
        const_fix[order].clear();
    }

    if (fix_harmonic) {

        if (verbosity > 0) {
            std::cout << "  FC2XML is given : Harmonic force constants will be " << std::endl;
            std::cout << "                    fixed to the values given in " << fc2_file << std::endl;
            std::cout << std::endl;
        }

        fix_forceconstants_to_file(0,
                                   symmetry,
                                   fcs,
                                   fc2_file,
                                   const_fix[0]);
    }

    fix_cubic = fix_cubic & (cluster->get_maxorder() > 1);
    if (fix_cubic) {

        if (verbosity > 0) {
            std::cout << "  FC3XML is given : Cubic force constants will be " << std::endl;
            std::cout << "                    fixed to the values given in " << fc3_file << std::endl;
            std::cout << std::endl;
        }

        fix_forceconstants_to_file(1,
                                   symmetry,
                                   fcs,
                                   fc3_file,
                                   const_fix[1]);
    }

    generate_symmetry_constraint_in_cartesian(system->get_supercell().number_of_atoms,
                                              symmetry,
                                              cluster,
                                              fcs,
                                              verbosity);


    extra_constraint_from_symmetry = false;

    for (auto order = 0; order < cluster->get_maxorder(); ++order) {
        if (!const_symmetry[order].empty()) extra_constraint_from_symmetry = true;
    }

    if (impose_inv_T) {
        // const_translation is updated.
        generate_translational_constraint(system->get_supercell(),
                                          symmetry,
                                          cluster,
                                          fcs,
                                          mirror_image_conv,
                                          verbosity);
    }

    if (impose_inv_R) {
        generate_rotational_constraint(system,
                                       symmetry,
                                       cluster,
                                       fcs,
                                       verbosity,
                                       tolerance_constraint);
    }

    // Merge intra-order constrants and do reduction

    for (auto order = 0; order < maxorder; ++order) {

        const auto nparam = fcs->get_nequiv()[order].size();

        const_self[order].reserve(
                const_translation[order].size()
                + const_rotation_self[order].size()
                + const_symmetry[order].size());

        // The order of const_symmetry and const_translation
        // should not be changed since the rref_sparse is sensitive to
        // the numerical accuracy of the input matrix.
        const_self[order].insert(const_self[order].end(),
                                 const_symmetry[order].begin(),
                                 const_symmetry[order].end());

        const_self[order].insert(const_self[order].end(),
                                 const_translation[order].begin(),
                                 const_translation[order].end());

        const_self[order].insert(const_self[order].end(),
                                 const_rotation_self[order].begin(),
                                 const_rotation_self[order].end());

//        rref_sparse(nparam, const_self[order], tolerance_constraint);
        rref_sparse(nparam, const_self[order], tolerance_constraint);
    }

    get_mapping_constraint(maxorder,
                           fcs->get_nequiv(),
                           const_self,
                           const_fix,
                           const_relate,
                           index_bimap);

#ifdef _DEBUG
    for (auto order = 0; order < maxorder; ++order) {
        std::cout << "const_relate:\n";

        for (const auto &it: const_relate[order]) {
            std::cout << std::setw(5) << it.p_index_target;
            std::cout << " : ";
            for (auto m = 0; m < it.alpha.size(); ++m) {
                std::cout << "(" << std::setw(15) << it.alpha[m] << ", " << std::setw(5) << it.p_index_orig[m] << ") ";
            }
            std::cout << std::endl;
        }

        std::cout << "\nindex_bimap:\n";

        for (const auto &it: index_bimap[order]) {
            std::cout << std::setw(5) << it.left << " <--> " << it.right << std::endl;
        }

    }

#endif

    if (!constraint_algebraic) {

        size_t Pmax = 0;
        size_t nparams = 0;
        for (auto order = 0; order < maxorder; ++order) {
            Pmax += const_self[order].size()
                    + const_rotation_cross[order].size();
        }
        if (fix_harmonic) {
            Pmax -= const_self[0].size();
            Pmax += fcs->get_nequiv()[0].size();
        }
        if (fix_cubic) {
            Pmax -= const_self[1].size();
            Pmax += fcs->get_nequiv()[1].size();
        }
        for (auto order = 0; order < maxorder; ++order) {
            nparams += fcs->get_nequiv()[order].size();
        }

        if (const_mat) {
            deallocate(const_mat);
        }
        allocate(const_mat, Pmax, nparams);

        if (const_rhs) {
            deallocate(const_rhs);
        }
        allocate(const_rhs, Pmax);

        // const_mat and const_rhs are updated.
        number_of_constraints = calc_constraint_matrix(maxorder,
                                                       fcs->get_nequiv(),
                                                       nparams);
    }

    exist_constraint
            = impose_inv_T
              || fix_harmonic
              || fix_cubic
              || extra_constraint_from_symmetry;

    if (verbosity > 0) {
        if (exist_constraint) {

            int order;

            if (impose_inv_T || impose_inv_R) {
                std::cout << "  Number of constraints [T-inv, R-inv (self), R-inv (cross)]:" << std::endl;
                for (order = 0; order < maxorder; ++order) {
                    std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                    std::cout << " " << std::setw(6) << const_translation[order].size();
                    std::cout << std::setw(5) << const_rotation_self[order].size();
                    std::cout << std::setw(5) << const_rotation_cross[order].size();
                    std::cout << std::endl;
                }
                std::cout << std::endl;
            }

            if (extra_constraint_from_symmetry) {
                std::cout << "  There are constraints from crystal symmetry." << std::endl;
                std::cout << "  The number of such constraints for each order:" << std::endl;
                for (order = 0; order < maxorder; ++order) {
                    std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                    std::cout << " " << std::setw(6) << const_symmetry[order].size();
                    std::cout << std::endl;
                }
                std::cout << std::endl;
            }

            if (extra_constraint_from_symmetry) {
                std::cout << "  Constraints of T-inv, R-inv (self), and those from crystal symmetry are merged."
                          << std::endl;
            } else {
                std::cout << "  Constraints of T-inv and R-inv (self) are merged." << std::endl;
            }
            std::cout << "  If there are redundant constraints, they are removed in this process." << std::endl;
            std::cout << std::endl;
            std::cout << "  Number of inequivalent constraints (self, cross) : " << std::endl;

            for (order = 0; order < maxorder; ++order) {
                std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                std::cout << " " << std::setw(6) << const_self[order].size();
                std::cout << std::setw(5) << const_rotation_cross[order].size();
                std::cout << std::endl;
            }
            std::cout << std::endl;

            if (constraint_algebraic) {

                std::cout << "  ICONST >= 10 : Constraints will be considered algebraically."
                          << std::endl << std::endl;

                if (impose_inv_R) {
                    std::cout << "  WARNING : Inter-order constraints for rotational invariance will be neglected."
                              << std::endl;
                }

                for (order = 0; order < maxorder; ++order) {
                    std::cout << "  Number of free" << std::setw(9) << cluster->get_ordername(order)
                              << " FCs : " << index_bimap[order].size() << std::endl;
                }
                std::cout << std::endl;

            } else {

                std::cout << "  Total number of constraints = " << number_of_constraints << std::endl << std::endl;

            }
        }
        timer->print_elapsed();
        std::cout << " -------------------------------------------------------------------" << std::endl;
        std::cout << std::endl;
    }


    deallocate(const_translation);
    const_translation = nullptr;
    deallocate(const_rotation_self);
    const_rotation_self = nullptr;
    deallocate(const_rotation_cross);
    const_rotation_cross = nullptr;
    deallocate(const_self);
    const_self = nullptr;

    timer->stop_clock("constraint");
}

size_t Constraint::calc_constraint_matrix(const int maxorder,
                                          const std::vector<size_t> *nequiv,
                                          const size_t nparams) const
{
    size_t i, j;
    int order;
    double *arr_tmp;
    std::vector<ConstraintClass> const_total;

    const_total.clear();
    allocate(arr_tmp, nparams);

    size_t nshift = 0;

    for (order = 0; order < maxorder; ++order) {
        const auto nelems = nequiv[order].size();

        if (const_fix[order].empty()) {
            for (auto &p: const_self[order]) {
                for (i = 0; i < nparams; ++i) arr_tmp[i] = 0.0;
                for (const auto &it: p) {
                    arr_tmp[nshift + it.first] = it.second;
                }
                const_total.emplace_back(nparams, arr_tmp);
            }
        }
        nshift += nelems;
    }


    const auto nconst1 = const_total.size();

    // Inter-order constraints
    size_t nshift2 = 0;
    for (order = 0; order < maxorder; ++order) {
        if (order > 0) {
            if (const_fix[order - 1].empty() && const_fix[order].empty()) {
                for (auto &p: const_rotation_cross[order]) {
                    for (i = 0; i < nparams; ++i) arr_tmp[i] = 0.0;<commit_msg>get_mapping_constraint should be called after generating<commit_after>/*
 constraint.cpp

 Copyright (c) 2014, 2015, 2016 <NAME>

 This file is distributed under the terms of the MIT license.
 Please see the file 'LICENCE.txt' in the root directory 
 or http://opensource.org/licenses/mit-license.php for information.
*/

#include "constraint.h"
#include "combination.h"
#include "constants.h"
#include "error.h"
#include "fcs.h"
#include "cluster.h"
#include "mathfunctions.h"
#include "memory.h"
#include "rref.h"
#include "symmetry.h"
#include "system.h"
#include "timer.h"
#include "xml_parser.h"
#include <iostream>
#include <iomanip>
#include <boost/bimap.hpp>
#include <algorithm>
#include <unordered_set>
#include <map>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/xml_parser.hpp>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>

using namespace ALM_NS;

Constraint::Constraint()
{
    set_default_variables();
}

Constraint::~Constraint()
{
    deallocate_variables();
}

void Constraint::set_default_variables()
{
    constraint_mode = 11;
    rotation_axis = "";
    fix_harmonic = false;
    fix_cubic = false;
    constraint_algebraic = 1;
    fc2_file = "";
    fc3_file = "";
    exist_constraint = false;
    extra_constraint_from_symmetry = false;
    const_mat = nullptr;
    const_rhs = nullptr;
    const_symmetry = nullptr;
    const_fix = nullptr;
    const_relate = nullptr;
    const_relate_rotation = nullptr;
    index_bimap = nullptr;
    number_of_constraints = 0;
    tolerance_constraint = eps8;
}

void Constraint::deallocate_variables()
{
    if (const_symmetry) {
        deallocate(const_symmetry);
    }
    if (const_fix) {
        deallocate(const_fix);
    }
    if (const_relate) {
        deallocate(const_relate);
    }
    if (const_relate_rotation) {
        deallocate(const_relate_rotation);
    }
    if (index_bimap) {
        deallocate(index_bimap);
    }
    if (const_mat) {
        deallocate(const_mat);
    }
    if (const_rhs) {
        deallocate(const_rhs);
    }
}

void Constraint::setup(const System *system,
                       const Fcs *fcs,
                       const Cluster *cluster,
                       const Symmetry *symmetry,
                       const int linear_model,
                       const int mirror_image_conv,
                       const int verbosity,
                       Timer *timer)
{
    timer->start_clock("constraint");

    if (verbosity > 0) {
        std::cout << " CONSTRAINT" << std::endl;
        std::cout << " ==========" << std::endl << std::endl;
    }

    constraint_algebraic = constraint_mode / 10;
    constraint_mode = constraint_mode % 10;
    const auto maxorder = cluster->get_maxorder();

    if (linear_model == 2) {
        if (constraint_mode > 1) {
            warn("Constraint::setup", "Sorry, only ICONST = 11 is supported \n"
                                      "                      when LMODEL = enet. We set ICONST = 11 in this run.\n");
            constraint_mode = 1;
        }
        constraint_algebraic = 1;
    }

    switch (constraint_mode) {
        case 0: // do nothing
            impose_inv_T = false;
            impose_inv_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 0: Constraint for translational/rotational invariance" << std::endl;
                std::cout << "              will NOT be considered." << std::endl;
            }
            break;
        case 1:
            impose_inv_T = true;
            impose_inv_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 1: Constraints for translational invariance" << std::endl;
                std::cout << "              will be considered." << std::endl;
            }
            break;
        case 2:
            impose_inv_T = true;
            impose_inv_R = true;
            exclude_last_R = true;
            if (verbosity > 0) {
                std::cout << "  ICONST = 2: Constraints for translational and rotational invariance" << std::endl;
                std::cout << "              will be considered. Axis of rotation is " << rotation_axis << std::endl;
                std::cout << "              Rotational invariance of the maximum order will be neglected" << std::endl;
            }
            break;
        case 3:
            impose_inv_T = true;
            impose_inv_R = true;
            exclude_last_R = false;
            if (verbosity > 0) {
                std::cout << "  ICONST = 3: Constraints for translational and rotational invariance" << std::endl;
                std::cout << "              will be considered. Axis of rotation is " << rotation_axis << std::endl;
            }
            break;
        default:
            exit("Constraint::setup", "invalid constraint_mode", constraint_mode);
            break;
    }

    if (fcs->get_forceconstant_basis() == "Lattice" && impose_inv_R) {
        exit("Constraint::setup()", "Sorry, rotational invariance with FCSYM_BASIS = Lattice is "
                                    "not supported.\n Use FCSYM_BASIS = Cartesian instead.");
    }

    if (verbosity > 0) std::cout << std::endl;

    if (const_symmetry) {
        deallocate(const_symmetry);
    }
    if (const_fix) {
        deallocate(const_fix);
    }
    if (const_relate) {
        deallocate(const_relate);
    }
    if (index_bimap) {
        deallocate(index_bimap);
    }

    allocate(const_fix, maxorder);
    allocate(const_relate, maxorder);
    allocate(index_bimap, maxorder);
    allocate(const_symmetry, maxorder);
    allocate(const_translation, maxorder);
    allocate(const_rotation_self, maxorder);
    allocate(const_rotation_cross, maxorder);
    allocate(const_self, maxorder);

    for (auto order = 0; order < maxorder; ++order) {
        const_translation[order].clear();
        const_rotation_self[order].clear();
        const_rotation_cross[order].clear();
        const_self[order].clear();
        const_fix[order].clear();
    }

    if (fix_harmonic) {

        if (verbosity > 0) {
            std::cout << "  FC2XML is given : Harmonic force constants will be " << std::endl;
            std::cout << "                    fixed to the values given in " << fc2_file << std::endl;
            std::cout << std::endl;
        }

        fix_forceconstants_to_file(0,
                                   symmetry,
                                   fcs,
                                   fc2_file,
                                   const_fix[0]);
    }

    fix_cubic = fix_cubic & (cluster->get_maxorder() > 1);
    if (fix_cubic) {

        if (verbosity > 0) {
            std::cout << "  FC3XML is given : Cubic force constants will be " << std::endl;
            std::cout << "                    fixed to the values given in " << fc3_file << std::endl;
            std::cout << std::endl;
        }

        fix_forceconstants_to_file(1,
                                   symmetry,
                                   fcs,
                                   fc3_file,
                                   const_fix[1]);
    }

    generate_symmetry_constraint_in_cartesian(system->get_supercell().number_of_atoms,
                                              symmetry,
                                              cluster,
                                              fcs,
                                              verbosity);


    extra_constraint_from_symmetry = false;

    for (auto order = 0; order < cluster->get_maxorder(); ++order) {
        if (!const_symmetry[order].empty()) extra_constraint_from_symmetry = true;
    }

    if (impose_inv_T) {
        // const_translation is updated.
        generate_translational_constraint(system->get_supercell(),
                                          symmetry,
                                          cluster,
                                          fcs,
                                          mirror_image_conv,
                                          verbosity);
    }

    if (impose_inv_R) {
        generate_rotational_constraint(system,
                                       symmetry,
                                       cluster,
                                       fcs,
                                       verbosity,
                                       tolerance_constraint);
    }

    // Merge intra-order constrants and do reduction

    for (auto order = 0; order < maxorder; ++order) {

        const auto nparam = fcs->get_nequiv()[order].size();

        const_self[order].reserve(
                const_translation[order].size()
                + const_rotation_self[order].size()
                + const_symmetry[order].size());

        // The order of const_symmetry and const_translation
        // should not be changed since the rref_sparse is sensitive to
        // the numerical accuracy of the input matrix.
        const_self[order].insert(const_self[order].end(),
                                 const_symmetry[order].begin(),
                                 const_symmetry[order].end());

        const_self[order].insert(const_self[order].end(),
                                 const_translation[order].begin(),
                                 const_translation[order].end());

        const_self[order].insert(const_self[order].end(),
                                 const_rotation_self[order].begin(),
                                 const_rotation_self[order].end());

//        rref_sparse(nparam, const_self[order], tolerance_constraint);
        rref_sparse(nparam, const_self[order], tolerance_constraint);
    }


#ifdef _DEBUG
    for (auto order = 0; order < maxorder; ++order) {
        std::cout << "const_relate:\n";

        for (const auto &it: const_relate[order]) {
            std::cout << std::setw(5) << it.p_index_target;
            std::cout << " : ";
            for (auto m = 0; m < it.alpha.size(); ++m) {
                std::cout << "(" << std::setw(15) << it.alpha[m] << ", " << std::setw(5) << it.p_index_orig[m] << ") ";
            }
            std::cout << std::endl;
        }

        std::cout << "\nindex_bimap:\n";

        for (const auto &it: index_bimap[order]) {
            std::cout << std::setw(5) << it.left << " <--> " << it.right << std::endl;
        }

    }

#endif

    if (!constraint_algebraic) {

        size_t Pmax = 0;
        size_t nparams = 0;
        for (auto order = 0; order < maxorder; ++order) {
            Pmax += const_self[order].size()
                    + const_rotation_cross[order].size();
        }
        if (fix_harmonic) {
            Pmax -= const_self[0].size();
            Pmax += fcs->get_nequiv()[0].size();
        }
        if (fix_cubic) {
            Pmax -= const_self[1].size();
            Pmax += fcs->get_nequiv()[1].size();
        }
        for (auto order = 0; order < maxorder; ++order) {
            nparams += fcs->get_nequiv()[order].size();
        }

        if (const_mat) {
            deallocate(const_mat);
        }
        allocate(const_mat, Pmax, nparams);

        if (const_rhs) {
            deallocate(const_rhs);
        }
        allocate(const_rhs, Pmax);

        // const_mat and const_rhs are updated.
        number_of_constraints = calc_constraint_matrix(maxorder,
                                                       fcs->get_nequiv(),
                                                       nparams);
    }

    get_mapping_constraint(maxorder,
                           fcs->get_nequiv(),
                           const_self,
                           const_fix,
                           const_relate,
                           index_bimap);

    exist_constraint
            = impose_inv_T
              || fix_harmonic
              || fix_cubic
              || extra_constraint_from_symmetry;

    if (verbosity > 0) {
        if (exist_constraint) {

            int order;

            if (impose_inv_T || impose_inv_R) {
                std::cout << "  Number of constraints [T-inv, R-inv (self), R-inv (cross)]:" << std::endl;
                for (order = 0; order < maxorder; ++order) {
                    std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                    std::cout << " " << std::setw(6) << const_translation[order].size();
                    std::cout << std::setw(5) << const_rotation_self[order].size();
                    std::cout << std::setw(5) << const_rotation_cross[order].size();
                    std::cout << std::endl;
                }
                std::cout << std::endl;
            }

            if (extra_constraint_from_symmetry) {
                std::cout << "  There are constraints from crystal symmetry." << std::endl;
                std::cout << "  The number of such constraints for each order:" << std::endl;
                for (order = 0; order < maxorder; ++order) {
                    std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                    std::cout << " " << std::setw(6) << const_symmetry[order].size();
                    std::cout << std::endl;
                }
                std::cout << std::endl;
            }

            if (extra_constraint_from_symmetry) {
                std::cout << "  Constraints of T-inv, R-inv (self), and those from crystal symmetry are merged."
                          << std::endl;
            } else {
                std::cout << "  Constraints of T-inv and R-inv (self) are merged." << std::endl;
            }
            std::cout << "  If there are redundant constraints, they are removed in this process." << std::endl;
            std::cout << std::endl;
            std::cout << "  Number of inequivalent constraints (self, cross) : " << std::endl;

            for (order = 0; order < maxorder; ++order) {
                std::cout << "   " << std::setw(8) << cluster->get_ordername(order);
                std::cout << " " << std::setw(6) << const_self[order].size();
                std::cout << std::setw(5) << const_rotation_cross[order].size();
                std::cout << std::endl;
            }
            std::cout << std::endl;

            if (constraint_algebraic) {

                std::cout << "  ICONST >= 10 : Constraints will be considered algebraically."
                          << std::endl << std::endl;

                if (impose_inv_R) {
                    std::cout << "  WARNING : Inter-order constraints for rotational invariance will be neglected."
                              << std::endl;
                }

                for (order = 0; order < maxorder; ++order) {
                    std::cout << "  Number of free" << std::setw(9) << cluster->get_ordername(order)
                              << " FCs : " << index_bimap[order].size() << std::endl;
                }
                std::cout << std::endl;

            } else {

                std::cout << "  Total number of constraints = " << number_of_constraints << std::endl << std::endl;

            }
        }
        timer->print_elapsed();
        std::cout << " -------------------------------------------------------------------" << std::endl;
        std::cout << std::endl;
    }


    deallocate(const_translation);
    const_translation = nullptr;
    deallocate(const_rotation_self);
    const_rotation_self = nullptr;
    deallocate(const_rotation_cross);
    const_rotation_cross = nullptr;
    deallocate(const_self);
    const_self = nullptr;

    timer->stop_clock("constraint");
}

size_t Constraint::calc_constraint_matrix(const int maxorder,
                                          const std::vector<size_t> *nequiv,
                                          const size_t nparams) const
{
    size_t i, j;
    int order;
    double *arr_tmp;
    std::vector<ConstraintClass> const_total;

    const_total.clear();
    allocate(arr_tmp, nparams);

    size_t nshift = 0;

    for (order = 0; order < maxorder; ++order) {
        const auto nelems = nequiv[order].size();
        if (const_fix[order].empty()) {
            for (auto &p: const_self[order]) {
                for (i = 0; i < nparams; ++i) arr_tmp[i] = 0.0;
                for (const auto &it: p) {
                    arr_tmp[nshift + it.first] = it.second;
                }
                const_total.emplace_back(nparams, arr_tmp);
            }
        }
        nshift += nelems;
    }

    const auto nconst1 = const_total.size();

    // Inter-order constraints
    size_t nshift2 = 0;
    for (order = 0; order < maxorder; ++order) {
        if (order > 0) {
            if (const_fix[order - 1].empty() && const_fix[order].empty()) {
                for (auto &p: const_rotation_cross[order]) {
                    for (i = 0; i < nparams; ++i) arr_tmp[i] = 0.0;
{'fraction_non_alphanumeric': 0.12406977436446985, 'fraction_numerical': 0.006370185152110496, 'mean_word_length': 1.5808558039486824, 'pattern_counts': {'":': 0, '<': 541, '<?xml version=': 0, '>': 95, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 102, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2925276', 'n_tokens_mistral': 9769, 'n_tokens_neox': 9232, 'n_words': 2006}