<commit_before>/* 
 * File:   imc.cc
 * Author: ruehle
 *
 * Created on March 10, 2009, 3:43 PM
 */

#include <iostream>
#include <fstream>
#include <iomanip>
#include <boost/lexical_cast.hpp>
#include <tools/rangeparser.h>
#include "beadlist.h"
#include "nblistgrid.h"
#include "imc.h"
#include "imcio.h"

Imc::Imc()
   : _write_every(0), _do_blocks(false), _do_imc(false)
{
}

Imc::~Imc()
{
}

// begin the coarse graining process
// here the data structures are prepared to handle all the data
void Imc::BeginCG(Topology *top, Topology *top_atom) {
    // do some output
    cout << "begin to calculate inverse monte carlo parameters\n";
    cout << "# of bonded interactions: " << _bonded.size() << endl;
    cout << "# of non-bonded interactions: " << _nonbonded.size() << endl;
    
   // we didn't process any frames so far
    _nframes = 0;
    _nblock = 0;
    
   // initialize non-bonded structures
   for (list<Property*>::iterator iter = _nonbonded.begin();
        iter != _nonbonded.end(); ++iter) {
            interaction_t *i = AddInteraction(*iter);
            // calculate normalization factor for rdf
            BeadList beads1, beads2;
        
            beads1.Generate(*top, (*iter)->get("type1").value());
            beads2.Generate(*top, (*iter)->get("type2").value());

            if((*iter)->get("type1").value() ==  (*iter)->get("type2").value())
                i->_norm = top->BoxVolume()/(4.*M_PI* i->_step * beads1.size()*(beads2.size()-1.)/2.);
            else
                i->_norm = top->BoxVolume()/(4.*M_PI* i->_step * beads1.size()*beads2.size());

            i->_is_bonded = false;
   }
   
    // initialize non-bonded structures
   for (list<Property*>::iterator iter = _bonded.begin();
        iter != _bonded.end(); ++iter) {
            interaction_t *i = AddInteraction(*iter);
            i->_is_bonded = true;
   }
   
   // initialize the group structures
    if(_do_imc)
        InitializeGroups();
};

// create an entry for interactions
Imc::interaction_t *Imc::AddInteraction(Property *p)
{
    string name = p->get("name").value();
    string group;
    if(_do_imc)
        group = p->get("imc.group").value();
    else
        group = "none";
    
    interaction_t *i = new interaction_t;    
    _interactions[name] = i;
    getGroup(group)->_interactions.push_back(i);
    
    i->_step = p->get("step").as<double>();
    i->_min = p->get("min").as<double>();
    i->_max = p->get("max").as<double>();
    i->_norm = 1.0;
    i->_p=p;
    
    // initialize the current and average histogram
    int n = (int)((i->_max - i->_min) / i->_step + 1.000000001);

    i->_current.Initialize(i->_min, i->_max+i->_step, n);
    i->_average.Initialize(i->_min, i->_max+i->_step, n);
    
    return i;
}

// end of trajectory, post processing data
void Imc::EndCG()
{
    if(!_do_blocks) {
        WriteDist();
        if(_do_imc)
            WriteIMCData();
    }
    // clear interactions and groups
    _interactions.clear();
    _groups.clear();
}

// load options from xml file
void Imc::LoadOptions(const string &file) {
    load_property_from_xml(_options, file);
    _bonded = _options.Select("cg.bonded");
    _nonbonded = _options.Select("cg.non-bonded");
        
}

// evaluate current conformation
void Imc::EvalConfiguration(Topology *top, Topology *top_atom) {
    _nframes++;
    // process non-bonded interactions
    DoNonbonded(top);
    // process bonded interactions
    DoBonded(top);
    // update correlation matrices
    if(_do_imc)
        DoCorrelations();
    
    if(_write_every != 0) {
        if((_nframes % _write_every)==0) {
            _nblock++;
            string suffix = string("_") + boost::lexical_cast<string>(_nblock);
            WriteDist(suffix);
            WriteIMCData(suffix);
            WriteIMCBlock(suffix);
            if(_do_blocks)
                ClearAverages();
        }
    }
            
}

void Imc::ClearAverages()
{
    map<string, interaction_t *>::iterator ic_iter;
    map<string, group_t *>::iterator group_iter;
    
    _nframes = 0;
    for (ic_iter = _interactions.begin(); ic_iter != _interactions.end(); ++ic_iter)
        ic_iter->second->_average.Clear();    
    
    for (group_iter = _groups.begin(); group_iter != _groups.end(); ++group_iter)
        group_iter->second->_corr.clear();      
}

// process non-bonded interactions for current frame
void Imc::DoNonbonded(Topology *top)
{
    for (list<Property*>::iterator iter = _nonbonded.begin();
            iter != _nonbonded.end(); ++iter) {
        string name = (*iter)->get("name").value();
        
        interaction_t &i = *_interactions[name];
        
        // generate the bead lists
        BeadList beads1, beads2;
        
        beads1.Generate(*top, (*iter)->get("type1").value());
        beads2.Generate(*top, (*iter)->get("type2").value());
        
        // generate the neighbour list
        NBListGrid nb;
        nb.setCutoff(i._max + i._step);
        
        // is it same types or different types?
        if((*iter)->get("type1").value() == (*iter)->get("type2").value())
            nb.Generate(beads1);
        else
            nb.Generate(beads1, beads2);
        
        // clear the current histogram
        i._current.Clear();
        
        // process all pairs
        NBListGrid::iterator pair_iter;
        for(pair_iter = nb.begin(); pair_iter!=nb.end();++pair_iter) {
                i._current.Process((*pair_iter)->dist());            
        }
        
        // update the average
        i._average.data().y() = (((double)_nframes-1.0)*i._average.data().y() 
                + i._current.data().y())/(double)_nframes;
    }    
}

// process non-bonded interactions for current frame
void Imc::DoBonded(Topology *top)<commit_msg>switched to old nblist till bugs in grid search are fixed<commit_after>/* 
 * File:   imc.cc
 * Author: ruehle
 *
 * Created on March 10, 2009, 3:43 PM
 */

#include <iostream>
#include <fstream>
#include <iomanip>
#include <boost/lexical_cast.hpp>
#include <tools/rangeparser.h>
#include "beadlist.h"
#include "nblistgrid.h"
#include "imc.h"
#include "imcio.h"

Imc::Imc()
   : _write_every(0), _do_blocks(false), _do_imc(false)
{
}

Imc::~Imc()
{
}

// begin the coarse graining process
// here the data structures are prepared to handle all the data
void Imc::BeginCG(Topology *top, Topology *top_atom) {
    // do some output
    cout << "begin to calculate inverse monte carlo parameters\n";
    cout << "# of bonded interactions: " << _bonded.size() << endl;
    cout << "# of non-bonded interactions: " << _nonbonded.size() << endl;
    
   // we didn't process any frames so far
    _nframes = 0;
    _nblock = 0;
    
   // initialize non-bonded structures
   for (list<Property*>::iterator iter = _nonbonded.begin();
        iter != _nonbonded.end(); ++iter) {
            interaction_t *i = AddInteraction(*iter);
            // calculate normalization factor for rdf
            BeadList beads1, beads2;
        
            beads1.Generate(*top, (*iter)->get("type1").value());
            beads2.Generate(*top, (*iter)->get("type2").value());

            if((*iter)->get("type1").value() ==  (*iter)->get("type2").value())
                i->_norm = top->BoxVolume()/(4.*M_PI* i->_step * beads1.size()*(beads2.size()-1.)/2.);
            else
                i->_norm = top->BoxVolume()/(4.*M_PI* i->_step * beads1.size()*beads2.size());

            i->_is_bonded = false;
   }
   
    // initialize non-bonded structures
   for (list<Property*>::iterator iter = _bonded.begin();
        iter != _bonded.end(); ++iter) {
            interaction_t *i = AddInteraction(*iter);
            i->_is_bonded = true;
   }
   
   // initialize the group structures
    if(_do_imc)
        InitializeGroups();
};

// create an entry for interactions
Imc::interaction_t *Imc::AddInteraction(Property *p)
{
    string name = p->get("name").value();
    string group;
    if(_do_imc)
        group = p->get("imc.group").value();
    else
        group = "none";
    
    interaction_t *i = new interaction_t;    
    _interactions[name] = i;
    getGroup(group)->_interactions.push_back(i);
    
    i->_step = p->get("step").as<double>();
    i->_min = p->get("min").as<double>();
    i->_max = p->get("max").as<double>();
    i->_norm = 1.0;
    i->_p=p;
    
    // initialize the current and average histogram
    int n = (int)((i->_max - i->_min) / i->_step + 1.000000001);

    i->_current.Initialize(i->_min, i->_max+i->_step, n);
    i->_average.Initialize(i->_min, i->_max+i->_step, n);
    
    return i;
}

// end of trajectory, post processing data
void Imc::EndCG()
{
    if(!_do_blocks) {
        WriteDist();
        if(_do_imc)
            WriteIMCData();
    }
    // clear interactions and groups
    _interactions.clear();
    _groups.clear();
}

// load options from xml file
void Imc::LoadOptions(const string &file) {
    load_property_from_xml(_options, file);
    _bonded = _options.Select("cg.bonded");
    _nonbonded = _options.Select("cg.non-bonded");
        
}

// evaluate current conformation
void Imc::EvalConfiguration(Topology *top, Topology *top_atom) {
    _nframes++;
    // process non-bonded interactions
    DoNonbonded(top);
    // process bonded interactions
    DoBonded(top);
    // update correlation matrices
    if(_do_imc)
        DoCorrelations();
    
    if(_write_every != 0) {
        if((_nframes % _write_every)==0) {
            _nblock++;
            string suffix = string("_") + boost::lexical_cast<string>(_nblock);
            WriteDist(suffix);
            WriteIMCData(suffix);
            WriteIMCBlock(suffix);
            if(_do_blocks)
                ClearAverages();
        }
    }
            
}

void Imc::ClearAverages()
{
    map<string, interaction_t *>::iterator ic_iter;
    map<string, group_t *>::iterator group_iter;
    
    _nframes = 0;
    for (ic_iter = _interactions.begin(); ic_iter != _interactions.end(); ++ic_iter)
        ic_iter->second->_average.Clear();    
    
    for (group_iter = _groups.begin(); group_iter != _groups.end(); ++group_iter)
        group_iter->second->_corr.clear();      
}

// process non-bonded interactions for current frame
void Imc::DoNonbonded(Topology *top)
{
    for (list<Property*>::iterator iter = _nonbonded.begin();
            iter != _nonbonded.end(); ++iter) {
        string name = (*iter)->get("name").value();
        
        interaction_t &i = *_interactions[name];
        
        // generate the bead lists
        BeadList beads1, beads2;
        
        beads1.Generate(*top, (*iter)->get("type1").value());
        beads2.Generate(*top, (*iter)->get("type2").value());
        
        // generate the neighbour list
        NBList nb;
        nb.setCutoff(i._max + i._step);
        
        // is it same types or different types?
        if((*iter)->get("type1").value() == (*iter)->get("type2").value())
            nb.Generate(beads1);
        else
            nb.Generate(beads1, beads2);
        
        // clear the current histogram
        i._current.Clear();
        
        // process all pairs
        NBList::iterator pair_iter;
        for(pair_iter = nb.begin(); pair_iter!=nb.end();++pair_iter) {
                i._current.Process((*pair_iter)->dist());            
        }
        
        // update the average
        i._average.data().y() = (((double)_nframes-1.0)*i._average.data().y() 
                + i._current.data().y())/(double)_nframes;
    }    
}

// process non-bonded interactions for current frame
void Imc::DoBonded(Topology *top)
{'fraction_non_alphanumeric': 0.1601448026202379, 'fraction_numerical': 0.009653508015859335, 'mean_word_length': 2.395668715247293, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 119, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 0, 'substrings_counts': 2, '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': '2592598', 'n_tokens_mistral': 3874, 'n_tokens_neox': 3716, 'n_words': 773}