text:
/**
 * GeoDa TM, Copyright (C) 2011-2015 by Luc Anselin - all rights reserved
 *
 * This file is part of GeoDa.
 *
 * GeoDa is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * GeoDa 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <math.h>
#include <iostream>
#include <iomanip>
#include <utility> // std::pair
#include <boost/foreach.hpp>

#include <wx/wx.h>
#include <wx/xrc/xmlres.h>
#include <wx/dcclient.h>
#include <wx/wfstream.h>
#include <wx/txtstrm.h>

#include "../HighlightState.h"
#include "../GeneralWxUtils.h"
#include "../GeoDa.h"
#include "../logger.h"
#include "../Project.h"
#include "SimpleAxisCanvas.h"
#include "SimpleHistCanvas.h"
#include "CorrelogramView.h"

#ifdef __WIN32__
#ifndef isnan
#define isnan(x) _isnan(x)
#endif
#endif

BEGIN_EVENT_TABLE(CorrelogramFrame, TemplateFrame)
	EVT_MOUSE_EVENTS(CorrelogramFrame::OnMouseEvent)
	EVT_ACTIVATE(CorrelogramFrame::OnActivate)
END_EVENT_TABLE()

CorrelogramFrame::CorrelogramFrame(wxFrame *parent, Project* project,
                                   const wxString& title,
                                   const wxPoint& pos,
                                   const wxSize& size)
: TemplateFrame(parent, project, title, pos, size, wxDEFAULT_FRAME_STYLE),
correl_params_frame(0), panel(0),
panel_v_szr(0), bag_szr(0), top_h_sizer(0),
hist_plot(0), local_hl_state(0), message_win(0), project(project), shs_plot(0)
{
    wxLogMessage("Open CorrelogramFrame.");
	local_hl_state = new HighlightState();
	supports_timeline_changes = true;
	{
		std::vector<wxString> tm_strs;
		project->GetTableInt()->GetTimeStrings(tm_strs);
		var_man.ClearAndInit(tm_strs);
	}

	int width, height;
	GetClientSize(&width, &height);

	panel = new wxPanel(this);
	panel->SetBackgroundColour(*wxWHITE);
	SetBackgroundColour(*wxWHITE);
	panel->Bind(wxEVT_MOTION, &CorrelogramFrame::OnMouseEvent, this);

	bag_szr = new wxGridBagSizer(0, 0); // 0 vgap, 0 hgap

	panel_v_szr = new wxBoxSizer(wxVERTICAL);
	panel_v_szr->Add(bag_szr, 1, wxEXPAND);

	wxBoxSizer* panel_h_szr = new wxBoxSizer(wxHORIZONTAL);
	panel_h_szr->Add(panel_v_szr, 1, wxEXPAND);

	panel->SetSizer(panel_h_szr);

	// Top Sizer for Frame
	top_h_sizer = new wxBoxSizer(wxHORIZONTAL);
	top_h_sizer->Add(panel, 1, wxEXPAND|wxALL, 8);

	SetSizer(top_h_sizer);
	DisplayStatusBar(true);

	local_hl_state->SetSize(par.bins);

	UpdateDataMapFromVarMan();
	UpdateCorrelogramData();

	SetupPanelForNumVariables(var_man.GetVarsCount());

	Show(true);

	wxCommandEvent ev;
	OnShowCorrelParams(ev);
}

CorrelogramFrame::~CorrelogramFrame()
{
	if (correl_params_frame) {
		correl_params_frame->removeObserver(this);
		correl_params_frame->closeAndDeleteWhenEmpty();
	}
	if (local_hl_state) local_hl_state->closeAndDeleteWhenEmpty();
	if (HasCapture()) ReleaseMouse();
	DeregisterAsActive();
}

void CorrelogramFrame::OnMouseEvent(wxMouseEvent& event)
{
    if (event.RightDown()) {
        const wxPoint& pos = event.GetPosition();
    	// Workaround for right-click not changing window focus in OSX / wxW 3.0
    	wxActivateEvent ae(wxEVT_NULL, true, 0, wxActivateEvent::Reason_Mouse);
    	OnActivate(ae);

        OnRightClick(pos);
    }
}

void CorrelogramFrame::OnRightClick(const wxPoint& pos)
{
    wxMenu* optMenu;
    optMenu = wxXmlResource::Get()->LoadMenu("ID_CORRELOGRAM_MENU_OPTIONS");
    if (!optMenu) return;

    UpdateContextMenuItems(optMenu);
    PopupMenu(optMenu, pos);
    UpdateOptionMenuItems();

    wxMenuItem* save_menu = optMenu->FindItem(XRCID("ID_SAVE_CORRELOGRAM_STATS"));
    Connect(save_menu->GetId(), wxEVT_MENU,
            wxCommandEventHandler(CorrelogramFrame::OnSaveResult));
}

void CorrelogramFrame::OnSaveResult(wxCommandEvent& event)
{
    wxLogMessage("In CorrelogramFrame::OnSaveResult()");
    wxFileDialog
    saveFileDialog(this, _("Save Statistics file"), "", "",
                   "csv files (*.csv)|*.csv", wxFD_SAVE|wxFD_OVERWRITE_PROMPT);
    if (saveFileDialog.ShowModal() == wxID_CANCEL)
        return;

    wxFileOutputStream output_stream(saveFileDialog.GetPath());
    if (!output_stream.IsOk())
    {
        wxLogError("Cannot save current contents in file '%s'.", saveFileDialog.GetPath());
        return;
    }
    // write logReport to a text file
    wxTextOutputStream txt_out( output_stream );
    txt_out << "";

    vector<wxString> lbls;
    lbls.push_back("Autocorr.");
    lbls.push_back("Min");
    lbls.push_back("Max");
    lbls.push_back("# Pairs");

    wxString header = "";
    int total_pairs = 0;
    for (size_t i=0; i<cbins.size(); ++i) {
        header << "," << "bin-" << i;
        lbls[0] << "," << cbins[i].corr_avg;
        lbls[1] << "," << cbins[i].dist_min;
        lbls[2] << "," << cbins[i].dist_max;
        lbls[3] << "," << cbins[i].num_pairs;
        total_pairs += cbins[i].num_pairs;
    }

    txt_out << header << "\n";
    for (size_t i=0; i<lbls.size(); i++) {
        txt_out << lbls[i] << "\n";
    }

    double min = cbins[0].dist_min;
    double max = cbins[cbins.size()-1].dist_max;
    double range_left = 0, range_right=0;
    double est_dist = GetEstDistWithZeroAutocorr(range_left, range_right);

    txt_out << "\nSummary:\n";
    txt_out << "min dist, max dist, total # pairs, Autocorr.=0 range, est. distance\n";
    txt_out << min << "," << max << "," << total_pairs << ",";
    txt_out << "[" << range_left << " : " << range_right << "],";
    txt_out << est_dist << "\n";

    txt_out.Flush();
    output_stream.Close();
}

void CorrelogramFrame::OnActivate(wxActivateEvent& event)
{
	if (event.GetActive()) {
        wxLogMessage("In CorrelogramFrame::OnActivate()");
		RegisterAsActive("CorrelogramFrame", GetTitle());
	}
}

void CorrelogramFrame::MapMenus()
{
	wxMenuBar* mb = GdaFrame::GetGdaFrame()->GetMenuBar();
	// Map Options Menus
	wxMenu* optMenu;
	optMenu = wxXmlResource::Get()->LoadMenu("ID_CORRELOGRAM_MENU_OPTIONS");
	CorrelogramFrame::UpdateContextMenuItems(optMenu);

	GeneralWxUtils::ReplaceMenu(mb, "Options", optMenu);
	UpdateOptionMenuItems();
}

void CorrelogramFrame::UpdateOptionMenuItems()
{
	//TemplateFrame::UpdateOptionMenuItems(); // set common items first
	wxMenuBar* mb = GdaFrame::GetGdaFrame()->GetMenuBar();
	int menu = mb->FindMenu("Options");
	if (menu != wxNOT_FOUND) {
		CorrelogramFrame::UpdateContextMenuItems(mb->GetMenu(menu));
	}
}

void CorrelogramFrame::UpdateContextMenuItems(wxMenu* menu)
{
	// Update the checkmarks and enable/disable state for the
	// following menu items if they were specified for this particular
	// view in the xrc file.  Items that cannot be enable/disabled,
	// or are not checkable do not appear.
	TemplateFrame::UpdateContextMenuItems(menu); // set common items
}

void CorrelogramFrame::OnShowCorrelParams(wxCommandEvent& event)
{
    wxLogMessage("In CorrelogramFrame::OnShowCorrelParams()");
	if (correl_params_frame) {
		correl_params_frame->Iconize(false);
		correl_params_frame->Raise();
		correl_params_frame->SetFocus();
	} else {
		CorrelParams cp;
		cp.dist_metric = project->GetDefaultDistMetric();
		cp.dist_units = project->GetDefaultDistUnits();
		correl_params_frame = new CorrelParamsFrame(cp, var_man, project);
		correl_params_frame->registerObserver(this);
	}
}

void CorrelogramFrame::OnDisplayStatistics(wxCommandEvent& event)
{
    wxLogMessage("In CorrelogramFrame::OnDisplayStatistics()");
	UpdateOptionMenuItems();
}

/** Implementation of TableStateObserver interface */
void CorrelogramFrame::update(TableState* o)
{
	if (correl_params_frame)
        correl_params_frame->UpdateFromTable();
}

/** Implementation of TimeStateObserver interface */
void CorrelogramFrame::update(TimeState* o)
{
    bool has_time_var = false;
    for (data_map_type::iterator i=data_map.begin(); i!=data_map.end(); ++i) {
        if (i->second.size() > 1) {
            has_time_var = true;
            break;
        }
    }

    if (has_time_var) {
        var_man.UpdateGlobalTime(o->GetCurrTime());
        UpdateDataMapFromVarMan();
        UpdateCorrelogramData();
        SetupPanelForNumVariables(var_man.GetVarsCount());
        Refresh();
    }
}

void CorrelogramFrame::ReDraw()
{
	UpdateDataMapFromVarMan();
	UpdateCorrelogramData();
	SetupPanelForNumVariables(var_man.GetVarsCount());
	Refresh();
}

/** Implementation of CorrelParams interface */
void CorrelogramFrame::update(CorrelParamsObservable* o)
{
    par = o->GetCorrelParams();
    UpdateDataMapFromVarMan();
    UpdateCorrelogramData();
    SetupPanelForNumVariables(var_man.GetVarsCount());
    Refresh();
}

void CorrelogramFrame::notifyOfClosing(CorrelParamsObservable* o)
{
	correl_params_frame = 0;
    UpdateMessageWin();
}

void CorrelogramFrame::notifyNewHover(const std::vector<int>& hover_obs,
                                      int total_hover_obs)
{
	wxStatusBar* sb = GetStatusBar();
	if (!sb) return;

	wxString s;
	if (total_hover_obs > 0) {
		int id = hover_obs[0];
		s << "autocorrelation is " << cbins[id].corr_avg;
		s << " for obs within distance band ";
		s << cbins[id].dist_min << " to " << cbins[id].dist_max;
	}
	sb->SetStatusText(s);
}

/** Implementation of SimpleScatterPlotCanvasCbInt interface */
void CorrelogramFrame::notifyNewHistHover(const std::vector<int>& hover_obs,
                                          int total_hover_obs)
{
	wxStatusBar* sb = GetStatusBar();
	if (!sb) return;
	if (total_hover_obs == 0) {
		sb->SetStatusText("");
		return;
	}
	wxString s;
	int id = hover_obs[0];
	s << cbins[id].num_pairs << " pairs in distance band ";
	s << cbins[id].dist_min << " to " << cbins[id].dist_max;
	sb->SetStatusText(s);
}

/** At this time only one variable is supported, but this could easily
 be changed in the future, so will leave the num_vars parameter. */
void CorrelogramFrame::SetupPanelForNumVariables(int num_vars)
{
	if (!panel || !bag_szr)
        return;
	int num_top_rows = GenUtils::max<int>(1, num_vars);
	int num_rows_total = num_top_rows + 3;
	if (message_win) {
		message_win->Unbind(wxEVT_MOTION, &CorrelogramFrame::OnMouseEvent, this);
		message_win->Destroy();
		message_win = 0;
	}
	bag_szr->Clear();
	panel_v_szr->Remove(bag_szr); // bag_szr is deleted automatically
	bag_szr = new wxGridBagSizer(0, 0); // 0 vgap, 0 hgap
	for (size_t i=0, sz=scatt_plots.size(); i<sz; ++i) {
		if (scatt_plots[i]) {
			scatt_plots[i]->Unbind(wxEVT_MOTION,
                                   &CorrelogramFrame::OnMouseEvent,
                                   this);
			scatt_plots[i]->Destroy();
		}
	}
	scatt_plots.clear();
	if (hist_plot) {
		hist_plot->Unbind(wxEVT_MOTION, &CorrelogramFrame::OnMouseEvent, this);
		hist_plot->Destroy();
	}
	hist_plot = 0;
    if (shs_plot) {
        shs_plot->Unbind(wxEVT_MOTION, &CorrelogramFrame::OnMouseEvent, this);
        shs_plot->Destroy();
    }
    shs_plot = 0;
	for (size_t i=0, sz=vert_labels.size(); i<sz; ++i) {
		if (vert_labels[i])
            vert_labels[i]->Destroy();
	}
	vert_labels.clear();
	for (size_t i=0, sz=horiz_labels.size(); i<sz; ++i) {
		if (horiz_labels[i])
            horiz_labels[i]->Destroy();
	}
	horiz_labels.clear();

	wxString type_str;
	if (par.method == CorrelParams::ALL_PAIRS ||
			par.method == CorrelParams::ALL_PAIRS_THRESH) {
		type_str = " - all pairs";
	} else if (par.method == CorrelParams::RAND_SAMP ||
						 par.method == CorrelParams::RAND_SAMP_THRESH) {
		type_str = " - random";
	}
	if (par.method == CorrelParams::ALL_PAIRS_THRESH ||
			par.method == CorrelParams::RAND_SAMP_THRESH) {
		type_str << ", cutoff: " << GenUtils::DblToStr(par.threshold, 4);
	}

    bool valid_sampling = true;

	double freq_min = 0;
	double freq_max = 0;
	double dist_min = 0;
	double dist_max = 1;
	if (cbins.size() > 0) {
		dist_min = cbins[0].dist_min;
		dist_max = cbins[cbins.size()-1].dist_max;
	}
	std::vector<double> dist_vals(cbins.size());
	std::vector<double> freq_vals(cbins.size());
	for (size_t i=0; i<cbins.size(); ++i) {
		freq_vals[i] = (double) cbins[i].num_pairs;
		if (i == 0) {
			dist_vals[i] = cbins[i].dist_min;
		} else {
			dist_vals[i] = (double) cbins[i].dist_max;
		}
		if (freq_vals[i] > freq_max) freq_max = freq_vals[i];
	}
	if (num_vars < 1) {
		// add blank cell upper left-hand corner
		bag_szr->Add(50, 50, wxGBPosition(0, 0), wxGBSpan(1,1));
		message_win = new wxHtmlWindow(panel, wxID_ANY, wxDefaultPosition, wxSize(200,-1));
		message_win->Bind(wxEVT_RIGHT_UP, &CorrelogramFrame::OnMouseEvent, this);
        UpdateMessageWin();
		bag_szr->Add(message_win, wxGBPosition(0,1), wxGBSpan(1,1), wxEXPAND);
		SetTitle("Correlogram" + type_str);

	} else {
		for (int row=0; row<num_vars; ++row) {
			wxString nm=var_man.GetName(row);
			int tm=var_man.GetTime(row);
			double y_min = -1;
			double y_max = 1;
			if (row == 0) {
				SetTitle("Correlogram - " + var_man.GetNameWithTime(row) + type_str);
			}

			wxString title("Autocorr. of " + var_man.GetNameWithTime(row));

			std::vector<double> Y(cbins.size());
			std::vector<bool> Y_undef(cbins.size());

			for (size_t i=0; i<cbins.size(); ++i) {
				Y[i] = cbins[i].corr_avg;
                Y_undef[i] = false;
                if (isnan(Y[i])) {
                    valid_sampling = false;
                }
			}

			AxisScale v_axs;
			v_axs.ticks = 5;
			v_axs.data_min = -1;
			v_axs.data_max = 1;
			v_axs.scale_min = v_axs.data_min;
			v_axs.scale_max = v_axs.data_max;
			v_axs.scale_range = v_axs.data_max-v_axs.data_min;
			v_axs.tic_inc = v_axs.scale_range / (double) (v_axs.ticks-1);
			v_axs.p = 1;
			v_axs.tics.resize(v_axs.ticks);
			v_axs.tics_str.resize(v_axs.ticks);
			v_axs.tics_str_show.resize(v_axs.ticks);
			for (int i=0; i<v_axs.ticks; ++i) {
				double d = i;
				v_axs.tics[i] = v_axs.data_min + d*v_axs.tic_inc;
				std::stringstream s;
				s << v_axs.tics[i];
				v_axs.tics_str[i] = s.str();
				v_axs.tics_str_show[i] = true;
			}

			SimpleAxisCanvas* sa_can = 0;
			{
				sa_can = new SimpleAxisCanvas(panel, this, project,
											  local_hl_state,
											  Y, Y_undef, title, y_min, y_max,
											  false, // is horizontal ?
											  true, // show axes
											  false, // hide negative labels
											  false, // add auto padding min
											  false, // add auto padding max
											  -1, // number ticks
											  true, // force tick at min
											  true, // force tick at max
											  &v_axs,
                                              false, // non-standardized
											  wxDefaultPosition,
											  wxSize(50, -1));
				bag_szr->Add(sa_can, wxGBPosition(row, 0), wxGBSpan(1,1), wxEXPAND);
				vert_labels.push_back(sa_can);
			}

			std::vector<double> X(cbins.size());
			std::vector<bool> X_undef(cbins.size());
			for (size_t i=0; i<cbins.size(); ++i) {
				X[i] = cbins[i].dist_min + (cbins[i].dist_max - cbins[i].dist_min)/2.0;
                X_undef[i] = false;
			}

			SimpleScatterPlotCanvas* sp_can = 0;
			sp_can = new SimpleScatterPlotCanvas(panel, this, project,
												 local_hl_state, this,
												 X, Y, X_undef, Y_undef,
                                                 title, title,
												 dist_min, dist_max,
												 y_min, y_max,
												 false, false, true,
												 "ID_CORRELOGRAM_MENU_OPTIONS",
												 false, // show axes
												 true, // show horiz axis thru orig
												 false, // show vert axis thru orig
												 false, false,
												 valid_sampling, // show LOWESS fit
												 false);
            sp_can->SetFixedAspectRatioMode(false);
			sp_can->ChangeLoessParams(0.2,5,0.02);
			bag_szr->Add(sp_can, wxGBPosition(row, 1), wxGBSpan(1,1), wxEXPAND);
			scatt_plots.push_back(sp_can);
		}
	}

	{
		{
			AxisScale v_axs;
			v_axs.ticks = 4;
			v_axs.data_min = freq_min;
			v_axs.data_max = freq_max;
			v_axs.scale_min = v_axs.data_min;
			v_axs.scale_max = v_axs.data_max;
			v_axs.scale_range = v_axs.data_max-v_axs.data_min;
			v_axs.tic_inc = v_axs.scale_range / (double) (v_axs.ticks-1);
			v_axs.p = 1;
			v_axs.tics.resize(v_axs.ticks);
			v_axs.tics_str.resize(v_axs.ticks);
			v_axs.tics_str_show.resize(v_axs.ticks);
			for (int i=0; i<v_axs.ticks; ++i) {
				double d = i;
				v_axs.tics[i] = v_axs.data_min + d*v_axs.tic_inc;
				std::stringstream s;
				// round to nearest whole number
				s << setprecision(0) << fixed << v_axs.tics[i];
				v_axs.tics_str[i] = s.str();
				v_axs.tics_str_show[i] = true;
			}
			SimpleAxisCanvas* sa_can = 0;
			{
				sa_can = new SimpleAxisCanvas(panel, this, project,
											  local_hl_state,
											  freq_vals, "Frequency",
											  freq_min, freq_max,
											  false, // is horizontal ?
											  true, // show axes
											  false, // hide negative labels
											  false, // add auto padding min
											  false, // add auto padding max
											  -1, // number ticks
											  true, // force tick at min
											  true, // force tick at max
											  &v_axs,
                                              false, // non-standardized
											  wxDefaultPosition,
											  wxSize(50, -1));
				bag_szr->Add(sa_can, wxGBPosition(num_top_rows, 0),
										 wxGBSpan(1,1), wxEXPAND);
				vert_labels.push_back(sa_can);
			}

			SimpleBinsHistCanvas* sh_can = 0;
			std::vector<SimpleBin> hist_bins(cbins.size());
			for (size_t i=0; i<cbins.size(); ++i) {
				hist_bins[i].min = cbins[i].dist_min;
				hist_bins[i].max = cbins[i].dist_max;
				hist_bins[i].count = cbins[i].num_pairs;
			}
			sh_can = new SimpleBinsHistCanvas(panel, this, project,
											  local_hl_state,
											  this,
											  hist_bins, "", dist_min, dist_max,
											  "ID_CORRELOGRAM_MENU_OPTIONS",
											  false);
			bag_szr->Add(sh_can, wxGBPosition(num_top_rows, 1),
									 wxGBSpan(1,1), wxEXPAND);
			hist_plot = sh_can;
		}

		// add blank cell lower left-hand corner
		bag_szr->Add(50, 50, wxGBPosition(num_top_rows+1, 0), wxGBSpan(1,1));

		wxString title;
		if (par.dist_metric == WeightsMetaInfo::DM_arc) {
			title << "Arc Distance, ";
			if (par.dist_units == WeightsMetaInfo::DU_mile) {
				title << " mi";
			} else {
				title << " km";
			}
		} else {
			title << "Euclidean Distance";
		}

		AxisScale h_axs;
		h_axs.ticks = cbins.size()+1;
		h_axs.data_min = cbins[0].dist_min;
		h_axs.data_max = cbins[cbins.size()-1].dist_max;
		h_axs.scale_min = h_axs.data_min;
		h_axs.scale_max = h_axs.data_max;
		h_axs.scale_range = h_axs.data_max-h_axs.data_min;
		h_axs.tic_inc = h_axs.scale_range / (double) (h_axs.ticks-1);
		h_axs.p = 1;
		h_axs.tics.resize(h_axs.ticks);
		h_axs.tics_str.resize(h_axs.ticks);
		h_axs.tics_str_show.resize(h_axs.ticks);
		for (int i=0; i<h_axs.ticks; ++i) {
			double d = i;
			h_axs.tics[i] = h_axs.data_min + d*h_axs.tic_inc;
			stringstream ss;
			if (h_axs.tics[i] < 10000000) {
				ss << std::fixed << std::setprecision(1) << h_axs.tics[i];
				h_axs.tics_str[i] = ss.str();
			} else {
				ss << std::setprecision(1) << h_axs.tics[i];
				h_axs.tics_str[i] = ss.str();
			}
			h_axs.tics_str_show[i] = true;
		}

		SimpleAxisCanvas* sa_can = 0;
		sa_can = new SimpleAxisCanvas(panel, this, project, local_hl_state,
									dist_vals, title, dist_min, dist_max,
									true, // is horizontal ?
									true, // show axes
									false, // hide negative labels
									false, // add auto padding min
									false, // add auto padding max
									cbins.size(), // number ticks,
									true, // force tick at min
									true, // force tick at max
									&h_axs, // custom axis scale pointer
                                    false, // non-standardized
									wxDefaultPosition, wxSize(-1, 50));

		bag_szr->Add(sa_can, wxGBPosition(num_top_rows+1, 1),
								 wxGBSpan(1,1), wxEXPAND);
		horiz_labels.push_back(sa_can);

        // add blank cell lower left-hand corner
        bag_szr->Add(0, 0, wxGBPosition(num_top_rows+2, 0), wxGBSpan(1,1));
        bag_szr->Add(0, 0, wxGBPosition(num_top_rows+2, 1), wxGBSpan(1,1));

	}

	bag_szr->SetFlexibleDirection(wxBOTH);

	// first column
    if (bag_szr->IsColGrowable(0)) {
        bag_szr->RemoveGrowableCol(0);
    }
	// final row  - 1(axis)
	if (bag_szr->IsRowGrowable(num_rows_total-2)) {
		bag_szr->RemoveGrowableRow(num_rows_total-2);
	}

	// second column
    if (bag_szr->IsColGrowable(1)) {
        bag_szr->RemoveGrowableCol(1);
    }
	bag_szr->AddGrowableCol(1, 2);

	// all rows exluding last two
	for (int i=0; i<num_top_rows; ++i) {
        if (bag_szr->IsRowGrowable(i)) {
            bag_szr->RemoveGrowableRow(i);
        }
		bag_szr->AddGrowableRow(i, 1);
	}

	// second-to-last row
	if (bag_szr->IsRowGrowable(num_top_rows)) {
		bag_szr->RemoveGrowableRow(num_top_rows);
	}
	bag_szr->AddGrowableRow(num_top_rows, 1);

	// last row
   /*
	if (bag_szr->IsRowGrowable(num_top_rows+1)) {
		bag_szr->RemoveGrowableRow(num_top_rows+1);
	}
	bag_szr->AddGrowableRow(num_top_rows+1, 1);
	if (bag_szr->IsRowGrowable(num_top_rows+2)) {
		bag_szr->RemoveGrowableRow(num_top_rows+2);
	}
	bag_szr->AddGrowableRow(num_top_rows+2, 1);
	*/
	panel_v_szr->Add(bag_szr, 1, wxALL | wxEXPAND);

    vector<wxString> lbls;
    lbls.push_back("Autocorr.");
    lbls.push_back("Min");
    lbls.push_back("Max");
    lbls.push_back("# Pairs");
    vector<vector<double> > vals;
    vector<double> stats;

    int sum_pairs = 0;
    for (size_t i=0; i<cbins.size(); ++i) {
        vector<double> sub_vals;
        sub_vals.push_back(cbins[i].corr_avg);
        sub_vals.push_back(cbins[i].dist_min);
        sub_vals.push_back(cbins[i].dist_max);
        sub_vals.push_back(cbins[i].num_pairs);
        vals.push_back(sub_vals);
        sum_pairs += cbins[i].num_pairs;
    }

    double min = cbins[0].dist_min;
    double max = cbins[cbins.size()-1].dist_max;
    double range_left = 0, range_right=0;
    double est_dist = GetEstDistWithZeroAutocorr(range_left, range_right);
    stats.push_back(min);
    stats.push_back(max);
    stats.push_back(sum_pairs);
    stats.push_back(range_left);
    stats.push_back(range_right);
    stats.push_back(est_dist);

    SimpleHistStatsCanvas* shs_can = 0;
    shs_can = new SimpleHistStatsCanvas(panel, this, project, local_hl_state,
                                        lbls, vals, stats,
                                        "ID_CORRELOGRAM_MENU_OPTIONS",
                                        wxDefaultPosition, wxSize(-1, 90));
    shs_can->SetFixedAspectRatioMode(false);

    //bag_szr->Add(shs_can, wxGBPosition(num_top_rows+2, 1), wxGBSpan(1,1), wxEXPAND);

    shs_plot = shs_can;

    panel_v_szr->Add(shs_can, 0, wxLEFT | wxRIGHT | wxEXPAND);

	top_h_sizer->RecalcSizes();

    if (valid_sampling == false ) {
        wxString msg = _("The sample size for random sampling is too small.\nPlease increase the number of iterations.");
        wxString title = _("Insufficient Random Sampling");
        wxMessageDialog dlg (this, msg, title, wxOK | wxICON_WARNING);
        dlg.ShowModal();
    }
}

double CorrelogramFrame::GetEstDistWithZeroAutocorr(double& rng_left,
                                                    double& rng_right)
{
    double rst = -1;
    for (size_t i=0; i<cbins.size()-1; ++i) {
        double a1 = cbins[i].corr_avg;
        double a2 = cbins[i+1].corr_avg;
        bool cross_axis_down = a1 == abs(a1) && a2 != abs(a2);
        bool cross_axis_up = a1 != abs(a1) && a2 == abs(a2);

        if ( cross_axis_down || cross_axis_up ) {
            double d1 = (cbins[i].dist_max + cbins[i].dist_min) / 2.0;
            double d2 = (cbins[i+1].dist_max + cbins[i+1].dist_min) / 2.0;

            rng_left = d1;
            rng_right = d2;

            //(d2 - d) / (a2 -0) = (d2 - d1) / (a2 - a1) ;
            double d = d2 - (d2 - d1) / (a2 - a1) * a2;
            rst = d;
            break;
        }
    }
    return rst;
}

void CorrelogramFrame::UpdateMessageWin()
{
	if (!message_win) return;
	wxString s;
	s << "<!DOCTYPE html>";
	s << "<html>";
	s << "<head>";
	s << "  <style type=\"text/css\">";
	s << "  body {";
	s << "    font-family: \"Trebuchet MS\", Arial, Helvetica, sans-serif;";
	s << "    font-size: small;";
	s << "    color: blue;";
	s << "  }";

	s << "  h1 {";
	s << "    font-family: \"Trebuchet MS\", Arial, Helvetica, sans-serif;";
	s << "    color: blue;";
	s << "  }";
	s <<   "</style>";

	s << "</head>";
	s << "<body>";
	s << "<br /><br /><br />";
	//s << "  <h1 align=\"center\">Message</h1>";
	s << "<center><p>";
	s << "<font face=\"verdana,arial,sans-serif\" color=\"black\" size=\"4\">";

	int count = var_man.GetVarsCount();
	if (count == 0) {
		s << "Please right-click or use<br />";
		s << "<font color=\"blue\">Options > Change Parameters<br /></font>";
		s << "to specify variable and distance parameters.";
	} else if (count > 0) {
		s << "Variables specified: <br />";
		for (int i=0; i<count; ++i) {
			s << var_man.GetName(i);
			if (i+1 < count) s << "<br />";
		}
	}

	s << "  </font></p></center>";
	s << "</body>";
	s << "</html>";
	message_win->SetPage(s );
}

/** Adds/removes variables from data_map according to variables present
 in var_man. */
void CorrelogramFrame::UpdateDataMapFromVarMan()
{
	using namespace std;
	// get set of var_man names
	set<wxString> vm_nms;
	for (int i=0; i<var_man.GetVarsCount(); ++i) {
		vm_nms.insert(var_man.GetName(i));
	}

	// remove items from data_map not in vm_nms
	set<wxString> to_remove;
	for (data_map_type::iterator i=data_map.begin(); i!=data_map.end(); ++i) {
		wxString nm(i->first);
		if (vm_nms.find(nm) != vm_nms.end())
            continue;
		to_remove.insert(nm);
	}

	for (set<wxString>::iterator i=to_remove.begin(); i!=to_remove.end(); ++i) {
		data_map.erase(*i);
        data_undef_map.erase(*i);
	}

	// add items to data_map that are in vm_nms, but not currently in data_map
	set<wxString> to_add;
	for (set<wxString>::iterator i=vm_nms.begin(); i!=vm_nms.end(); ++i) {
		wxString nm(*i);
		if (data_map.find(nm) != data_map.end())
            continue;
		to_add.insert(nm);
	}

	TableInterface* table_int = project->GetTableInt();
	for (set<wxString>::iterator i=to_add.begin(); i!=to_add.end(); ++i) {
		wxString nm = (*i);
		int c_id = table_int->FindColId(nm);
		if (c_id < 0) {
			//LOG_MSG("Error, variable not found in table: " + nm);
			continue;
		}
		int tms = table_int->GetColTimeSteps(c_id);
        vec_vec_dbl_type vec_vec_data(tms);
        vec_vec_bool_type vec_vec_undef(tms);
		for (int t=0; t<tms; ++t) {
			table_int->GetColData(c_id, t, vec_vec_data[t]);
			table_int->GetColUndefined(c_id, t, vec_vec_undef[t]);
		}
		pair<wxString, vec_vec_dbl_type> p(nm, vec_vec_data);
		data_map.insert(p);

		pair<wxString, vec_vec_bool_type> p_undef(nm, vec_vec_undef);
		data_undef_map.insert(p_undef);
	}
}

/** Update histogram and update correlogram if data_map has at least
 one variable */
bool CorrelogramFrame::UpdateCorrelogramData()
{
	using namespace GenGeomAlgs;
	using namespace CorrelogramAlgs;
	bool success = false;
	std::vector<double> Z;
	std::vector<bool> Z_undef;
	if (var_man.GetVarsCount() > 0) {
		wxString nm = var_man.GetName(0);
		int tm = var_man.GetTime(0);
		wxString title(var_man.GetNameWithTime(0));
		const std::vector<double>& data(data_map[nm][tm]);
		const std::vector<bool>& data_undef(data_undef_map[nm][tm]);
		Z.resize(data.size());
		Z_undef.resize(data.size());
        for (size_t i=0, sz=data.size(); i<sz; ++i) {
            Z[i] = data[i];
            Z_undef[i] = data_undef[i];
        }
	}
	bool is_arc = par.dist_metric == WeightsMetaInfo::DM_arc;
	bool is_mi = par.dist_units == WeightsMetaInfo::DU_mile;
	double th_rad = par.threshold;
	if (is_arc) {
		th_rad = is_mi ? EarthMiToRad(par.threshold) : EarthKmToRad(par.threshold);
	}

	std::vector<wxRealPoint> pts;
	project->GetCentroids(pts);

	if (par.method == CorrelParams::ALL_PAIRS) {
		success = MakeCorrAllPairs(pts, Z, Z_undef, is_arc, par.bins, cbins);

	} else if (par.method == CorrelParams::ALL_PAIRS_THRESH) {
		if (is_arc) {
			success = MakeCorrThresh(project->GetUnitSphereRtree(), Z, Z_undef,
                                     th_rad, par.bins, cbins);
		} else {
			success = MakeCorrThresh(project->GetEucPlaneRtree(), Z, Z_undef,
                                     par.threshold, par.bins, cbins);
		}
	} else if (par.method == CorrelParams::RAND_SAMP) {
		success = MakeCorrRandSamp(pts, Z, Z_undef, is_arc, -1,
                                   par.bins,  par.max_iterations, cbins);
	}	else if (par.method == CorrelParams::RAND_SAMP_THRESH) {
		success = MakeCorrRandSamp(pts, Z, Z_undef, is_arc,
                                   (is_arc ? th_rad : par.threshold),
                                   par.bins,  par.max_iterations, cbins);
	}

    if (success == false) {
        wxString msg = "Please select another variable with values more suitable for computing a correlogram.";
        wxString title = "Variable Value Error";
        wxMessageDialog dlg (this, msg, title, wxOK | wxICON_ERROR);
        dlg.ShowModal();
        return success;
    }

	if (success && par.bins != local_hl_state->GetHighlightSize()) {
		local_hl_state->SetSize(par.bins);
	}

	if (success && is_arc) {
		for (size_t i=0; i<cbins.size(); ++i) {
			if (is_mi) {
				cbins[i].dist_min = EarthRadToMi(cbins[i].dist_min);
				cbins[i].dist_max = EarthRadToMi(cbins[i].dist_max);
			} else {
				cbins[i].dist_min = EarthRadToKm(cbins[i].dist_min);
				cbins[i].dist_max = EarthRadToKm(cbins[i].dist_max);
			}
		}
	}

	return success;
}

max_stars_repo_name: chenyoujie/GeoDa

max_issues_repo_name: chenyoujie/GeoDa

max_forks_repo_name: chenyoujie/GeoDa


text:
#include <stdio.h>
#include <omp.h>

#include <Eigen/Dense>
#include <Parareal/core.h>
#include <Parareal/parareal.h>
#include <Parareal/forward_euler.h>

#include "rhs_linear1d.h"

typedef Eigen::VectorXd Evec;
typedef Eigen::MatrixXd Emat;

int main(int argc, char **argv)
{
  int P = omp_get_max_threads();

  printf("Using %d processors, testing weak scaling. \n", P);
  printf("Threads\tTime\n");
  for (int k = 2; k < P; k++)
  {
    omp_set_num_threads(k);
    ode_system ode;
    ode.dimension = 1; ode.t_init = 0; ode.t_final = 0.5*k;
    ode.y0 = Evec(1); ode.y0(0) = 1;
    ode.f = std::function<int(double, Evec&, Evec&)>(&linear1d);

    time_stepper course; time_stepper fine;

    // Forward Euler
    course.dt = 0.5;
    course.F = std::function<int(ode_system&, double, Evec &)>(&forward_euler);
    course.F_allt = std::function<int(ode_system&, double, Emat &)>(&forward_euler_allt);
    fine.dt = .0000001;
    fine.F = std::function<int(ode_system&, double, Evec &)>(&forward_euler);
    fine.F_allt = std::function<int(ode_system&, double, Emat &)>(&forward_euler_allt);

    // Test Parareal's Strong Scaling

    int csteps = ode.num_steps(course.dt);
    double tt = 0.0;

    // Parareal Solver
    Emat yf(csteps, ode.dimension);
    tt = omp_get_wtime();
    parareal(ode, course, fine, 6, yf);
    tt = omp_get_wtime() - tt;

    printf("%d\t%f\n", k, tt);
  }
  return 0;
}

ext: cpp

max_stars_repo_name: abhijit-c/Parareal

max_issues_repo_name: abhijit-c/Parareal

max_forks_repo_name: abhijit-c/Parareal


text:
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
////using namespace boost::numeric::ublas;

extern "C"{
  #include <cblas.h>
}

#include <gsl/gsl_integration.h>

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cstring>
#include <complex>      // std::complex
using namespace std;

#ifndef _CPLXCONST_
#define _CPLXCONST_

const  complex<double> ZeroC (0.0,0.0);
const  complex<double> OneC (1.0,0.0);
const  complex<double> OneImC (0.0,1.0);

#endif

//#ifndef _MYEPSCONST_
// #define _MYEPSCONST_

// const double MyEPS(1.0e-13);

// #endif
// Declares the classes
////////////////////////////////////////////////
////////////// Class CNRGarray \\\\\\\\\\\\\\\\\
////////////////////////////////////////////////

#ifndef _CNRGARRAY_
#define _CNRGARRAY_

class CNRGarray {
/* Declares the class */

public:

  // Default constructor

  //CNRGarray():totalS(false){}
  CNRGarray():MyEPS(1.0e-18)
  {}
  // Constructor with NQNumbers defined
  CNRGarray(int Nqns){
    NQNumbers=Nqns;
  }
  // Default  Destructor
  ~CNRGarray(){}

  // Member variables (regular)
  int Nshell;
  int NQNumbers;

  // Member variables (STL vectors)
  vector<double> QNumbers;

  // Vectors that spam the full basis
  vector<double> dEigVec;
  vector<double> dEn;
  vector<complex<double> > cEigVec;

  vector<int> iDegen;

  vector<int> BlockBegEnd;

  // April 09: list of N+1 states generated
  // (to be used in DM-NRG)

  vector< vector<int> > ChildStates;

 // August 2011: mark kept states. (spans the full basis)
  vector<bool> Kept;

  // April 10: totalS qnumbers
  bool totalS; // tru -> there are totalS QNs
  vector<int> Sqnumbers; // position(s) of total S in QNumbers
  // Idea: This can be used in an updated version of BuildBasis...
  //bool ComplexEVec; // tru -> complex eigenvectors

  double MyEPS;

  // Member functions
  int NumBlocks();
  int Nstates();
  double EniBlockj(int i, int j);
  void SetDegen();
  void SetE0zero();

  void FilterQNumbers();

  void PrintQNumbers();

  void PrintEn();

  void PrintAll();

  void PrintBlockQNumbers(int iblock);

  void PrintBlockEn(int iblock);

  void PrintBlock(int iblock);

  void ClearAll();

  double GetQNumber(int iblock, int whichqn);

  int GetBlockLimit(int iblock, int whichlimit);
  int GetBlockSize(int iblock);

  int GetBlockSize(int iblock,bool kp);

  int GetBlockFromSt(int ist,int &pos_in_block);

  int GetBlockLimitEigv(int iblock, int whichlimit);

  boost::numeric::ublas::matrix<double> EigVec2BLAS(int iblock);

  // 2013
  boost::numeric::ublas::matrix<complex<double> > cEigVec2BLAS(int iblock);

  double GetEigVecComponent(int ist,int istbasis); // April 09

  double Ecut(int Ncut);

  double GetBlockEmin(int iblock);

  double GetBlockEmax(int iblock);

  int GetBlockEcutPos(int iblock,double Ecut);

  void SetKept(int Ncutoff);
  bool CheckKept(int ist, bool kp);
  bool CheckKept(int iblock, int istbl, bool kp);

  int GetBlockFromQNumbers(double *qnums);

  double GetQNumberFromSt(int ist, int whichqn);

  // May 08
  int GetiGS();
  void SaveQSParameters();

  // Sep 08
  void AddQNumber();

  // Feb 09
  void SaveBin(char arqname[]);
  void ReadBin(char arqname[]);

  void ReadBinInStream(ifstream &InFile);

  // April 09

  double PartitionFuncTeq0();

  double PartitionFunc(double betabar);

  double PartitionFuncDisc(double betabar);

  // September 09
  void SetEigVecToOne();

  // Nov 2010
  boost::numeric::ublas::matrix<double> ExpEi2BLAS(int iblock, double betabar);

  // March 2014
  bool CheckComplex();

protected :

};

#endif

////////////////////////////////////////////////
/////////// Class CNRGbasisarray \\\\\\\\\\\\\\\
////////////////////////////////////////////////
//
//  Subclass of CNRGarray
//

#ifndef _CNRGBASISARRAY_
#define _CNRGBASISARRAY_

class CNRGbasisarray: public CNRGarray{

public :

  // Default constructor
  CNRGbasisarray(){}

  CNRGbasisarray(int Nqns){
    NQNumbers=Nqns;
  }
  // Default  Destructor

   // Constructor that copies members from a CNRGarray object
   CNRGbasisarray(CNRGarray A1){
     CNRGarray::Nshell=A1.Nshell;
     CNRGarray::NQNumbers=A1.NQNumbers;
     CNRGarray::QNumbers=A1.QNumbers;
     CNRGarray::dEn=A1.dEn;
     CNRGarray::dEigVec=A1.dEigVec;
     CNRGarray::cEigVec=A1.cEigVec;
     CNRGarray::Kept=A1.Kept;
     CNRGarray::BlockBegEnd=A1.BlockBegEnd;
     CNRGarray::totalS=A1.totalS; // tru -> there are totalS QNs
     CNRGarray::Sqnumbers=A1.Sqnumbers; // position(s) of total S in QNumbers
   }

  ~CNRGbasisarray(){}

  vector<int> iType;
  vector<int> StCameFrom;

  // moved here

  // NEW
  vector<int> BlockBegEndBC;
  vector<int> BlockBegEndEigVec;

  // Newer: if the symmetry is reduced at a certain point, need to retain
  // old Qnumbers

  int NQNumbers_stcf;
  vector<double>  StCameFromQNumbers;

  // Member functions

  void ClearAll();
  void PrintAll();
  void PrintBlockBasis(int iblock);

  void PrintBasisAll(){
    for (int ibl=0;ibl<NumBlocks();ibl++)
      PrintBlockBasis(ibl);
  }

  void SetSCFfromdEn();

  int  GetBlockSizeBC(int iblock);

  int GetBlockFromStBC(int ist,int &pos_in_block);

  void SyncNRGarray(CNRGarray &A1){
     CNRGarray::Nshell=A1.Nshell;
     CNRGarray::NQNumbers=A1.NQNumbers;
     CNRGarray::QNumbers=A1.QNumbers;
     CNRGarray::dEn=A1.dEn;
     CNRGarray::dEigVec=A1.dEigVec;
     CNRGarray::cEigVec=A1.cEigVec;
     CNRGarray::Kept=A1.Kept;
     CNRGarray::BlockBegEnd=A1.BlockBegEnd;
     CNRGarray::totalS=A1.totalS; // tru -> there are totalS QNs
     CNRGarray::Sqnumbers=A1.Sqnumbers; // position(s) of total S in QNumbers
  }

  // New
  void FalseCut(CNRGarray* pA1){
    CNRGbasisarray::ClearAll();
    CNRGbasisarray::SyncNRGarray(*pA1);
    CNRGbasisarray::SetSCFfromdEn();
    CNRGbasisarray::BlockBegEndBC=pA1->BlockBegEnd;
    CNRGbasisarray::SetBlockBegEndEigVec();
  }

  void SetBlockBegEndEigVec();

  void RemoveBlock(int iblock); // StCameFrom keeps the pre-cut state indexes

  void RemoveState(int ist); // StCameFrom keeps the pre-cut state indexes

  void RemoveStatesFromBlock(int ibl, int ist1,int ist2);

  int GetBlockLimitEigv(int iblock, int whichlimit); // April 09, overloads

  double GetEigVecComponent(int ist,int istbasis); // April 09 overloads

  // Ubuntu's g++ does not like this...
  boost::numeric::ublas::matrix<double> EigVecCut2BLAS(int iblock);
  boost::numeric::ublas::matrix<complex<double> > cEigVecCut2BLAS(int iblock);
  // converts the non-square eigenvector matrix to BLAS format

  boost::numeric::ublas::matrix<double> EigVecCut2BLAS(int iblock, bool kp);
  // converts the non-square eigenvector matrix to BLAS format (kp states)

  // April 08
  double GetStCameFromQNumberFromSt(int ist, int whichqn);

  // Sets dEn as the energies of the previous states in StCameFrom
  // SetdEnPrevious();

  void CopyBlock(int iblock, int NoCopies);
  // Create copies of block states and set iDegen

  // Sep 08
  void SetLastQNumber(vector<double> InputVec);

  // Feb 09

  void SaveBin(char arqname[]);
  void ReadBin(char arqname[]);

  void ReadBinInStream(ifstream &InFile);

};

#endif

////////////////////////////////////////////////
/////////// Class CNRGmatrix     \\\\\\\\\\\\\\\
////////////////////////////////////////////////
//
//  Subclass of CNRGarray
//

#ifndef _CNRGMATRIX_
#define _CNRGMATRIX_

class CNRGmatrix: public CNRGarray{

public :

  // Default constructor (with initialization list)
  CNRGmatrix():UpperTriangular(false),NeedOld(false),SaveMatYN(false),CalcAvg(false),IsComplex(false),WignerEckartL(0.5)
  {}
  //CNRGmatrix(){}

   // Constructor that copies members from a CNRGarray object
   CNRGmatrix(CNRGbasisarray A1){
     UpperTriangular=false;
     SaveMatYN=false;
     CalcAvg=false;
     IsComplex=false;
     CNRGarray::Nshell=A1.Nshell;
     CNRGarray::NQNumbers=A1.NQNumbers;
     CNRGarray::QNumbers=A1.QNumbers;
     CNRGarray::BlockBegEnd=A1.BlockBegEnd;
     CNRGarray::Kept=A1.Kept;
     //     CNRGarray::sizeBlock=A1.sizeBlock;
     iType=A1.iType;
     CNRGarray::totalS=A1.totalS; // tru -> there are totalS QNs
     CNRGarray::Sqnumbers=A1.Sqnumbers; // position(s) of total S in QNumbers
   }
   // Default  Destructor
   ~CNRGmatrix(){}

  //////////////////////
  //   Static members //
  //////////////////////

  vector<double> MatEl;
  // new
  //template <typename AnyType>
  vector< complex<double> > MatElCplx;
  // Best solution: look for MatEl everywhere and
  // add a if (IsComplex)-> use MatElCmplx instead

  vector<int> MatBlockMap;
  vector<int> MatBlockBegEnd;

  // New: same as in CNRGbasisarray
  vector<int> iType;

  // NEW: indicate kept-kept,disc-disc, disc-kept,...
  vector<bool> MatKept;

  // Does this matrix needs the "old" one on updating?
  bool NeedOld;

  // Stores upper triangular part of the blocks only?
  bool UpperTriangular;

  // Save matrix?
  bool SaveMatYN;

  // CalcAvg?
  bool CalcAvg;
  char MatName[20];

  // Is complex?
  bool IsComplex;

  // Wigner-Eckart rank (usually 1/2 but it can be 1)
  // Needed when using SU(2) symmetry.
  double WignerEckartL;

  //////////////////////
  // Function members //
  //////////////////////

  void SyncNRGarray(CNRGarray A1){
     CNRGarray::Nshell=A1.Nshell;
     CNRGarray::NQNumbers=A1.NQNumbers;
     CNRGarray::QNumbers=A1.QNumbers;
     CNRGarray::BlockBegEnd=A1.BlockBegEnd;
     CNRGarray::Kept=A1.Kept;
     CNRGarray::totalS=A1.totalS; // tru -> there are totalS QNs
     CNRGarray::Sqnumbers=A1.Sqnumbers; // position(s) of total S in QNumbers
  }
  /// Same for basis array
  void SyncNRGarray(CNRGbasisarray A1){
     CNRGarray::Nshell=A1.Nshell;
     CNRGarray::NQNumbers=A1.NQNumbers;
     CNRGarray::QNumbers=A1.QNumbers;
     CNRGarray::BlockBegEnd=A1.BlockBegEnd;
     CNRGarray::Kept=A1.Kept;
     CNRGarray::totalS=A1.totalS; // tru -> there are totalS QNs
     CNRGarray::Sqnumbers=A1.Sqnumbers; // position(s) of total S in QNumbers
     iType=A1.iType;
  }
//   ////
//   void ReverseSyncNRGarray(CNRGarray &A1){
//      A1.Nshell=CNRGarray::Nshell;
//      A1.NQNumbers=CNRGarray::NQNumbers;
//      A1.QNumbers=CNRGarray::QNumbers;
//      A1.BlockBegEnd=CNRGarray::BlockBegEnd;
//   }
  /// Same for basis array
  void ReverseSyncNRGarray(CNRGbasisarray &A1){
     A1.Nshell=CNRGarray::Nshell;
     A1.NQNumbers=CNRGarray::NQNumbers;
     A1.QNumbers=CNRGarray::QNumbers;
     A1.BlockBegEnd=CNRGarray::BlockBegEnd;
     A1.Kept=CNRGarray::Kept;
     A1.totalS=CNRGarray::totalS; // tru -> there are totalS QNs
     A1.Sqnumbers=CNRGarray::Sqnumbers; // position(s) of total S in QNumbers
     A1.iType=iType;
  }

  bool ChkSync(CNRGarray* pA1){
    if ( (CNRGarray::Nshell==pA1->Nshell)&&
	 (CNRGarray::NQNumbers==pA1->NQNumbers)&&
	 (CNRGarray::QNumbers==pA1->QNumbers)&&
	 (CNRGarray::BlockBegEnd==pA1->BlockBegEnd)&&
	 (CNRGarray::Kept==pA1->Kept)&&
	 (CNRGarray::totalS==pA1->totalS)&&
	 (CNRGarray::Sqnumbers==pA1->Sqnumbers)
	 ){
      return(true);
    }
    else{return(false);}
  }
  /// Same for basis array
  bool ChkSync(CNRGbasisarray* pA1){
    if ( (CNRGarray::Nshell==pA1->Nshell)&&
	 (CNRGarray::NQNumbers==pA1->NQNumbers)&&
	 (CNRGarray::QNumbers==pA1->QNumbers)&&
	 (CNRGarray::BlockBegEnd==pA1->BlockBegEnd)&&
	 (CNRGarray::Kept==pA1->Kept)&&
	 (CNRGarray::totalS==pA1->totalS)&&
	 (CNRGarray::Sqnumbers==pA1->Sqnumbers)&&
	 (iType==pA1->iType) // new (Sep 09, checking)
	 ){
      return(true);
    }
    else{return(false);}
  }

  ////

  void CopyData(CNRGmatrix* pMat1){
    // Array data
    Nshell=pMat1->Nshell;
    NQNumbers=pMat1->NQNumbers;
    QNumbers=pMat1->QNumbers;
    BlockBegEnd=pMat1->BlockBegEnd;
    Kept=pMat1->Kept;
    iType=pMat1->iType;
    totalS=pMat1->totalS;
    Sqnumbers=pMat1->Sqnumbers;
    // Matrix data
    MatEl=pMat1->MatEl;
    MatElCplx=pMat1->MatElCplx;
    MatBlockMap=pMat1->MatBlockMap;
    MatBlockBegEnd=pMat1->MatBlockBegEnd;
    MatKept=pMat1->MatKept;
    UpperTriangular=pMat1->UpperTriangular;
    IsComplex=pMat1->IsComplex;
    WignerEckartL=pMat1->WignerEckartL;
  }

  int NumMatBlocks();

  int FindMatBlock(int iblock1, int iblock2);

  int GetMatBlockLimit(int iblock1, int iblock2, int whichlimit);
  int GetMatBlockLimit(int imatblock, int whichlimit);

  int GetMatBlockSize(int iblock1, int iblock2);
  int GetMatBlockSize(int imatblock );

  double GetBlockMatEl(int iblock1, int iblock2, int iel, int jel);
  double GetMatEl(int ist, int jst);
  complex<double> cGetBlockMatEl(int iblock1, int iblock2, int iel, int jel);
  complex<double> cGetMatEl(int ist, int jst);

  // April 2019
  double CalcTrace();

  int GetMatElPosition(int iblock1, int iblock2, int iel, int jel);

  // deprecated
  void PushBlockMatEl(double El, int iblock1, int iblock2, int iel, int jel);
  ///
  void PushMatEl(double El, int ist, int jst);
  void PushMatEl(complex<double> El, int ist, int jst);

  ////// TEMPLATES! Testing. Not easy.
  template<class T> void TemplateTest(T El){

    cout << " This is the element: " << El << endl;
      return;
  }

// Replace EXISTING matrix element ist,jst by El
  template <class T> void TPushMatEl(T El, int ist, int jst){
    int ibl,jbl;

    int iblock1=CNRGarray::GetBlockFromSt(ist,ibl);
    int iblock2=CNRGarray::GetBlockFromSt(jst,jbl);

    int iPos=CNRGmatrix::GetMatElPosition(iblock1,iblock2,ibl,jbl);

    // This does not compile!
//     if (iPos>-1){
//       if (IsComplex) MatElCplx[iPos]=El;
//       else MatEl[iPos]=El;
//     }

  }
  //// END TEMPLATE TESTING

  void FilterMap_SetBegEnd();

  //  double* Block2Array(int imatblock);

  void ClearAll();

  void PrintMatBlockQNumbers(int imatblock);

  void PrintMatBlock(int imatblock);
  void PrintMatBlock(int iblock1, int iblock2);

  void PrintAllBlocks();

  void DiagBlock(int iblock,
		 vector<double> &eigvalues,
		 vector<double> &eigvectors);
  void DiagBlock(int iblock,
		 vector<double> &eigvalues,
		 vector<complex<double> > &eigvectors);

  // Pointers to functions to be defined

  // Checks whether there are non-zero matrix elements
  bool (*CheckForMatEl)(CNRGbasisarray *pAeigCut, int iblock1, int iblock2);

  // Calculates matrix elements between iblock1 and iblock2
  double (*CalcMatEl)(CNRGbasisarray *pAbasis, CNRGbasisarray *pSingleSite,
		      int isti, int istj);

  complex<double> (*CalcMatElCplx)(CNRGbasisarray *pAbasis,
			       CNRGbasisarray *pSingleSite,
			       int isti, int istj);

  double (*CalcHNMatEl)(vector<double> Params,
		      CNRGbasisarray* pAbasis,
		      CNRGbasisarray* pSingleSite,
		      CNRGmatrix* MatArray,
		      int ist, int jst);

  complex<double> (*CalcHNMatElCplx)(vector<double> Params,
		      CNRGbasisarray* pAbasis,
		      CNRGbasisarray* pSingleSite,
		      CNRGmatrix* MatArray,
		      int ist, int jst);

  //
  // THIS one is looking good!!!
  //

  void DiagHN(vector<double> ParamsHN,
	      CNRGbasisarray* pAbasis,
	      CNRGbasisarray* pSingleSite,
	      CNRGmatrix* MatArray,
	      CNRGarray* pAeig,bool display=false);

  //
  // This one too!
  //

  void PutInRegularForm(int iblock);

  //
  // Save in Binary
  //

  void SaveBinary(char arqname[]);

  void SaveInOldFormat();

  // Mar 09

  void SaveBin(char arqname[]);
  void ReadBin(char arqname[]);

  void ReadBinInStream(ifstream &InFile);

  // April 09

  // Given a block ibl_base, and stores the connecting blocks
  // and respective BlockMat in the two-column array
  // Connect_BlockMatBlock
  void GetConnectingBlocks(int ibl_base,
			   vector < vector<int> > &Connect_BlockMatBlock,
			   bool ToBase=false);
  // May 09

  void GetBlocksFromMatBlock(int iMatBlock, int &ibl1, int &ibl2);

  boost::numeric::ublas::matrix<double> MatBlock2BLAS(int iMatBlock);

  boost::numeric::ublas::matrix<double> MatBlock2BLAS(int ibl1, int ibl2);
  boost::numeric::ublas::matrix<complex<double> > cMatBlock2BLAS(int ibl1, int ibl2);

  boost::numeric::ublas::matrix<double> RotateBlock2BLAS(CNRGbasisarray* pAcut, int iMatBlock);

  boost::numeric::ublas::matrix<double> RotateBlock2BLAS(CNRGbasisarray* pAcut, int ibl1, int ibl2);
  boost::numeric::ublas::matrix<complex<double> > cRotateBlock2BLAS(CNRGbasisarray* pAcut, int ibl1, int ibl2);

  boost::numeric::ublas::matrix<double> RotateBlock2BLAS_NoCut(CNRGarray* pAeig, int ibl1, int ibl2, bool forward);
  boost::numeric::ublas::matrix<complex<double> > cRotateBlock2BLAS_NoCut(CNRGarray* pAeig, int ibl1, int ibl2, bool forward);

  // August 2011

  boost::numeric::ublas::matrix<double> MatBlock2BLAS(int iMatBlock,
						      bool kp1, bool kp2);

  boost::numeric::ublas::matrix<double> MatBlock2BLAS(int ibl1, int ibl2, bool kp1, bool kp2);
  boost::numeric::ublas::matrix<complex<double> > cMatBlock2BLAS(int ibl1, int ibl2, bool kp1, bool kp2);

  boost::numeric::ublas::matrix<double> RotateBlock2BLAS(CNRGbasisarray* pAcut, int iMatBlock, bool kp1, bool kp2);

  boost::numeric::ublas::matrix<double> RotateBlock2BLAS(CNRGbasisarray* pAcut, int ibl1, int ibl2, bool kp1, bool kp2);

  // April 2019

  void MatBlock2STL(vector<double>& Mout, int iMatBlock);

  void MatBlock2STL(vector<double>& Mout, int ibl1, int ibl2);

  void MatBlock2STL(vector<double>& Mout, int iMatBlock, bool kp1, bool kp2);

  void MatBlock2STL(vector<double>& Mout, int ibl1, int ibl2, bool kp1, bool kp2);

  void cMatBlock2STL(vector<complex<double> >& cMout, int iMatBlock);

  void cMatBlock2STL(vector<complex<double> >& cMout, int ibl1, int ibl2);

  void cMatBlock2STL(vector<complex<double> >& cMout, int iMatBlock, bool kp1, bool kp2);

  void cMatBlock2STL(vector<complex<double> >& cMout, int ibl1, int ibl2, bool kp1, bool kp2);

  // Sep 09
  // Set as Diagonal matrix (diagonal, UpperTriang=false)
  void SetDiagonalMatrix(double El);
  void SetIdentityMatrix();

  // Uses it's own QNumbers and iType to check for mat els.
  bool SelfCheckForMatEl(int iblock1, int iblock2);
  // Sets Matrix as zero, based on CheckForMatEl (SelfCheckForMatEl)
  void SetZeroMatrix();

  // May 2011
  // IF (*CalcMatEl) and (*CheckForMatEl) are defined, then it sets the matrix
  void SetMatrix(CNRGbasisarray* pAbasis,
		 CNRGbasisarray* pSingleSite);

};

#endif

#ifndef _CNRGCHAIN_
#define _CNRGCHAIN_

class CNRGchain{

public :

  // Default constructor (with initialization list)

  CNRGchain():Lambda(2.5),z_twist(1.0),BandOffset(0.0),HybRef(1.0),HybRefIsSet(false),HybFuncIsSumDeltas(false),Fsq(2.0),DiscScheme(0)
  {}

  CNRGchain(double ThisLambda):z_twist(1.0),BandOffset(0.0),HybRef(1.0),HybRefIsSet(false),HybFuncIsSumDeltas(false),Fsq(2.0),DiscScheme(0){
    Lambda=ThisLambda;
  }

  CNRGchain(double ThisLambda, int Nsitesmax):z_twist(1.0),BandOffset(0.0),HybRef(1.0),HybRefIsSet(false),HybFuncIsSumDeltas(false),Fsq(2.0),DiscScheme(0)
  {
    Lambda=ThisLambda;
    //z_twist=1.0;
    SetChainWilson(Nsitesmax);
  }

  // Default  Destructor
  ~CNRGchain(){}

  // Static members
  //int Nshell;
  //vector<double> eps_chi;

  double Lambda;
  double z_twist;
  double BandOffset;

  vector<double> en;
  vector<double> chin;

  vector<double> *HybFuncParams; // pointer to an STL vector
  // It SHOULD NOT be a POINTER!! Should allocate this WITHIN THE OBJECT!
  // Leave it like that for now...

  //vector< vector<double> > HybFuncVecParams; // STL array

  double HybRef;
  bool HybRefIsSet;
  bool HybFuncIsSumDeltas;
  double Fsq;
  int DiscScheme; // 0- usual; 1- CampoOliveira

  // Function members

  //void SetChain(double Lambda);

  void PrintChain(int Nsites);

  void PrintAll();

  void SetChainWilson(int Nsitesmax);

  void ReadParams(char arqname[],int Nparams);

  void ReadParams(char arqname[]); // Reads the whole thing into the STL vector
                                   // LAST NUMBER: Nlines.

  //void ReadVecParams(char arqname[], int Ncols);

  // Hybridization function in a GSL format (need to be defined externally!)
  double (*HybFunction)(double omega,
			void *params);

  double (*HybFuncWithEn)(double omega,
  			  void *params);

  // Ok, this is a bummer but...
  // hybridizations (Delta(en)) AND en*Delta(en) (or Delta(en)/en)
  // have to be defined OUTSIDE as static members
  // in order to be passed in here.
  // So, for every static HybFunction (=Delta(e))
  // there HAS TO BE a HybFuncWithEn defined as well

  // Tried hard to pass non-static function to GSL...
  // ...doesn't work. Ditto (may 2010)
  //   double HFtimesEn(double omega,
  // 		   void *params);

  //   static double HFtimesEnCallBack(double omega,
  // 				  void* params);

  ////////////

  double GetEn(int Nsites);
  double GetChin(int Nsites);

  // GSL integration
  double GSL_Integration(gsl_function *pIntegrandGSL,
			 gsl_integration_workspace *w_gsl,
			 double x0, double x1);

  // Lanczos procedure
  void SetChainLanczos(int Nsitesmax);

  double GetHyb_w(double omega);

  double ScaleFactor(int n);

private:

  // GSL stuff
//   gsl_function HybFunGSL;
//   gsl_function EnHybF;
  gsl_integration_workspace *wsp_gsl;

  void SetHybRef();
  double CalcF2(gsl_integration_workspace *w_gsl);

  // Adds elements of *pFm and *pEm
  void CalcFabmEabm(int m,
		    gsl_integration_workspace *w_gsl,
		    char which,
		    vector<double> *pFm,
		    vector<double> *pEm);

  double RHS(int n, double eam, double um, double um1,
	     double en, double chinm1, bool scale=true);

};

// struct NRGchainGSLCallbackHolder
// {
//   //CNRGchain* cls;
//   //void* data;
//   CNRGchain cls;
//   void* data;

// };

#endif

/////////////////////////////////////////////////
////////////// Class CNRGthermo \\\\\\\\\\\\\\\\\
/////////////////////////////////////////////////

#ifndef _CNRGTHERMO_
#define _CNRGTHERMO_

class CNRGthermo{

public :

  //   // Static members

  // Input/Output
  char ThermoName[64];
  char ChainArqName[64];
  char ArqName[64];

  double dImpValue;
  bool CalcChain;
  int Nsite0Chain;

  double betabar;

   //
  vector<double> dChainValues;
  vector<double> dValues;
  vector<double> dTempValues;

   // pointer to function:
  double (*Calc)(vector<double> Params,
		 CNRGarray* pAeig,
		 int sqnumber,bool totalS);

//   // Function members

  void AddValue(vector<double> Params,
		CNRGarray* pAeig,
		int sqnumber,bool totalS,
		double Temp);

  int ReadChain();

  void ReadNChainValue(int Nsites,int Nsites0);

  void SaveNValue(int Nsites, int Nsites0);

  void SaveToFile();

  bool CheckForFile();

};
// end CNRGthermo
#endif

//////////////////////////////////////////////////////
////////////// Class CNRGCodeHandler \\\\\\\\\\\\\\\\\
//////////////////////////////////////////////////////

#ifndef _CNRGCODEHANDLER_
#define _CNRGCODEHANDLER_

class CNRGCodeHandler{

public :

  // Default Constructor
  CNRGCodeHandler():NumChannels(1),Lambda(2.5){}

  // Default  Destructor
  ~CNRGCodeHandler(){}

  // Static members
  char ModelOption[32];
  char Symmetry[32];
  char BandType[32];

  int ModelNo;
  int SymNo;
  int BandNo;

  // Files

  char ParamFileName[64];
  char BeginFileName[64];
  char EndFileName[64];
  char SaveArraysFileName[64];

  char LancInFileName[64];

  int NoInputParamsDouble;
  vector <double> dInitParams;
  //  vector <string> cParamsName;

  int calcdens;

  int Ncutoff;
  int Nsitesmax;
  int Nsites0;
  int UpdateBefCut;
  int Nsites;

  int NumNRGmats;
  int NumThermoMats;
  int NumChannels;
  int NopsSaved;

  double Lambda;
//  double Dband;
  double ChemPot;
  double DN;
  double betabar;

  double HalfLambdaFactor;

  // totalS QNumbers
  bool totalS;
  int Sqnumber;

  // instream (perhaps not such a good idea)
  std::ifstream CodeInFile;

  // CNRGarray objects (pointers)

  CNRGbasisarray* pAcut;
  CNRGbasisarray* pAbasis;
  CNRGmatrix* MatArray;

  vector<CNRGthermo> ThermoSTLArray; // New Sep 09

  bool SaveData;

  // Need an array of chains
  //CNRGchain chain;
  vector<CNRGchain> Chains;
  // I will have to add a general ztwist and DiscScheme
  //ThisCode.chain.z_twist
  double code_z_twist;
  int code_DiscScheme; // 0- usual; 1- CampoOliveira
  // ThisCode.chain.DiscScheme

  // Function members

  //void InitialSetUp();
  // Trying this
  void InitialSetUp(bool ReadParamsOnly=false);

  void ReadParams();

  void ReadParams(char FileName[], int NoAdditionalParams);

  void SetChain();

  void SaveBegFile();

  void WrapUp();

  void SetSingleSite(CNRGbasisarray* pSingleSite);

  void CalcThermo(CNRGthermo* ThermoArray,
		  CNRGarray* pAeig);

  void CalcStuff(CNRGarray* pAeig, CNRGbasisarray* pAeigCut);

  void DoesNothing();

  void ZeroParams();

  void SetTotS();

  //void SaveGenPars(char idname[]);
  void SaveGenPars();

  void SaveArrays();
  void SaveArrays(char idname[]);

  //void ReadGenPars(char idname[]);
  //void ReadGenPars();
  // Keep this until all ThisCodePars.dat have been updated :)
  void ReadGenPars(bool ReadModelBand=false);

  void ReadArrays();
  void ReadArrays(char idname[]);

  void SetCurrentDN();

  void SetNopsSaved();

  bool CheckFileExists(char arqname[]);

  // New: May 2010

  void ModelSwitch(  vector<int> &CommonQNs,
		     vector<int> &totSpos,
		     CNRGarray* pAeig,
		     CNRGbasisarray* pSingleSite,
		     CNRGmatrix* pHN,
		     vector<CNRGmatrix> &STLMatArray,
		     CNRGthermo* ThermoArray,
		     vector<double> &chi_m1);

  double chiN(int Nsites, double Lambda);

  void PrintSettings();

};

#endif

ext: hpp


text:
// Copyright (c) 2014, Im Kyeong-Hyeon
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice, this
//   list of conditions and the following disclaimer.
//
// * Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef IWL_DEFINES_HPP_
#define IWL_DEFINES_HPP_

#include <functional>
#include <algorithm>
#include <vector>
#include <string>
#include <array>
#include <memory>
#include <stdexcept>
#include <utility>
#include <type_traits>

#include <cstdint>
#include <cstring>
#include <cassert>

#include <boost/noncopyable.hpp>

// iwl namespace macro
#define BEGIN_IWL() \
    namespace iwl {

#define END_IWL() \
    }

#endif // IWL_DEFINES_HPP_


text:
/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, development version     *
*                (c) 2006-2017 INRIA, USTL, UJF, CNRS, MGH                    *
*                                                                             *
* This program is free software; you can redistribute it and/or modify it     *
* under the terms of the GNU Lesser General Public License as published by    *
* the Free Software Foundation; either version 2.1 of the License, or (at     *
* your option) any later version.                                             *
*                                                                             *
* This program 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 Lesser General Public License *
* for more details.                                                           *
*                                                                             *
* You should have received a copy of the GNU Lesser General Public License    *
* along with this program. If not, see <http://www.gnu.org/licenses/>.        *
*******************************************************************************
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: contact@sofa-framework.org                             *
******************************************************************************/
#include <SofaBaseTopology/EdgeSetTopologyModifier.h>
#include <sofa/core/visual/VisualParams.h>
#include <sofa/core/topology/TopologyChange.h>
#include <SofaBaseTopology/EdgeSetTopologyContainer.h>
#include <algorithm>
#include <functional>
#include <iostream>
#include <sofa/core/ObjectFactory.h>

// Use BOOST GRAPH LIBRARY :

#include <boost/config.hpp>
#include <iostream>
#include <vector>
#include <utility>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/connected_components.hpp>

#include <boost/graph/cuthill_mckee_ordering.hpp>
#include <boost/graph/bandwidth.hpp>

namespace sofa
{

namespace component
{

namespace topology
{
using namespace sofa::defaulttype;
SOFA_DECL_CLASS(EdgeSetTopologyModifier)
int EdgeSetTopologyModifierClass = core::RegisterObject("Edge set topology modifier")
        .add< EdgeSetTopologyModifier >();

using namespace std;
using namespace sofa::defaulttype;
using namespace sofa::core::topology;

void EdgeSetTopologyModifier::init()
{
    PointSetTopologyModifier::init();
    getContext()->get(m_container);
}

void EdgeSetTopologyModifier::addEdgeProcess(Edge e)
{
#ifndef NDEBUG
    // check if the 2 vertices are different
    if(e[0] == e[1])
    {
        serr << "Error: [EdgeSetTopologyModifier::addEdge] : invalid edge: "
                << e[0] << ", " << e[1] << sendl;

        return;
    }

    // check if there already exists an edge.
    // Important: getEdgeIndex creates the edge vertex shell array
    if(m_container->hasEdgesAroundVertex())
    {
        if(m_container->getEdgeIndex(e[0],e[1]) != -1)
        {
            serr << "Error: [EdgeSetTopologyModifier::addEdgesProcess] : Edge "
                    << e[0] << ", " << e[1] << " already exists." << sendl;
            return;
        }
    }
#endif
    if (m_container->hasEdgesAroundVertex())
    {
        const unsigned int edgeId = m_container->getNumberOfEdges();

        sofa::helper::vector< unsigned int > &shell0 = m_container->getEdgesAroundVertexForModification( e[0] );
        shell0.push_back(edgeId);

        sofa::helper::vector< unsigned int > &shell1 = m_container->getEdgesAroundVertexForModification( e[1] );
        shell1.push_back(edgeId);
    }

    helper::WriteAccessor< Data< sofa::helper::vector<Edge> > > m_edge = m_container->d_edge;
    m_edge.push_back(e);
}

void EdgeSetTopologyModifier::addEdgesProcess(const sofa::helper::vector< Edge > &edges)
{
    for (unsigned int i=0; i<edges.size(); ++i)
    {
        addEdgeProcess(edges[i]);
    }
}

void EdgeSetTopologyModifier::addEdgesWarning(const unsigned int nEdges)
{
    m_container->setEdgeTopologyToDirty();

    // Warning that edges just got created
    EdgesAdded *e = new EdgesAdded(nEdges);
    addTopologyChange(e);
}

void EdgeSetTopologyModifier::addEdgesWarning(const unsigned int nEdges,
        const sofa::helper::vector< Edge >& edgesList,
        const sofa::helper::vector< unsigned int >& edgesIndexList)
{
    m_container->setEdgeTopologyToDirty();

    // Warning that edges just got created
    EdgesAdded *e = new EdgesAdded(nEdges, edgesList, edgesIndexList);
    addTopologyChange(e);
}

void EdgeSetTopologyModifier::addEdgesWarning(const unsigned int nEdges,
        const sofa::helper::vector< Edge >& edgesList,
        const sofa::helper::vector< unsigned int >& edgesIndexList,
        const sofa::helper::vector< sofa::helper::vector< unsigned int > > & ancestors)
{
    m_container->setEdgeTopologyToDirty();

    // Warning that edges just got created
    EdgesAdded *e = new EdgesAdded(nEdges, edgesList, edgesIndexList, ancestors);
    addTopologyChange(e);
}

void EdgeSetTopologyModifier::addEdgesWarning(const unsigned int nEdges,
        const sofa::helper::vector< Edge >& edgesList,
        const sofa::helper::vector< unsigned int >& edgesIndexList,
        const sofa::helper::vector< sofa::helper::vector< unsigned int > > & ancestors,
        const sofa::helper::vector< sofa::helper::vector< double > >& baryCoefs)
{
    m_container->setEdgeTopologyToDirty();

    // Warning that edges just got created
    EdgesAdded *e = new EdgesAdded(nEdges, edgesList, edgesIndexList, ancestors, baryCoefs);
    addTopologyChange(e);
}

void EdgeSetTopologyModifier::addEdgesWarning(const unsigned int nEdges,
        const sofa::helper::vector< Edge >& edgesList,
        const sofa::helper::vector< unsigned int >& edgesIndexList,
        const sofa::helper::vector< core::topology::EdgeAncestorElem >& ancestorElems)
{
    m_container->setEdgeTopologyToDirty();

    sofa::helper::vector< sofa::helper::vector< unsigned int > > ancestors;
    sofa::helper::vector< sofa::helper::vector< double > > baryCoefs;
    ancestors.resize(nEdges);
    baryCoefs.resize(nEdges);

    for (unsigned int i=0; i<nEdges; ++i)
    {
        for (unsigned int j=0; j<ancestorElems[i].srcElems.size(); ++j)
        {
            sofa::core::topology::TopologyElemID src = ancestorElems[i].srcElems[j];
            if (src.type == sofa::core::topology::EDGE && src.index != sofa::core::topology::Topology::InvalidID)
                ancestors[i].push_back(src.index);
        }
        if (!ancestors[i].empty())
        {
            baryCoefs[i].insert(baryCoefs[i].begin(), ancestors[i].size(), 1.0/ancestors[i].size());
        }
    }

    // Warning that edges just got created
    EdgesAdded *e = new EdgesAdded(nEdges, edgesList, edgesIndexList, ancestorElems, ancestors, baryCoefs);
    addTopologyChange(e);
}

void EdgeSetTopologyModifier::removeEdgesWarning(sofa::helper::vector<unsigned int> &edges )
{
    m_container->setEdgeTopologyToDirty();

    // sort edges to remove in a descendent order
    std::sort( edges.begin(), edges.end(), std::greater<unsigned int>() );

    // Warning that these edges will be deleted
    EdgesRemoved *e = new EdgesRemoved(edges);
    addTopologyChange(e);

    //std::cout<<"addTopologyChange EdgesRemoved"<<std::endl;
}

void EdgeSetTopologyModifier::removeEdgesProcess(const sofa::helper::vector<unsigned int> &indices,
        const bool removeIsolatedItems)
{
    if(!m_container->hasEdges())	// this method should only be called when edges exist
    {
        sout << "Warning. [EdgeSetTopologyModifier::removeEdgesProcess] edge array is empty." << sendl;
        return;
    }

    if(removeIsolatedItems && !m_container->hasEdgesAroundVertex())
    {
        m_container->createEdgesAroundVertexArray();
    }

    sofa::helper::vector<unsigned int> vertexToBeRemoved;
    helper::WriteAccessor< Data< sofa::helper::vector<Edge> > > m_edge = m_container->d_edge;

    unsigned int lastEdgeIndex = m_container->getNumberOfEdges() - 1;
    for (unsigned int i=0; i<indices.size(); ++i, --lastEdgeIndex)
    {
        // now updates the shell information of the edge formely at the end of the array
        if(m_container->hasEdgesAroundVertex())
        {

            const Edge &e = m_edge[ indices[i] ];
            const Edge &q = m_edge[ lastEdgeIndex ];
            const unsigned int point0 = e[0], point1 = e[1];
            const unsigned int point2 = q[0], point3 = q[1];

            sofa::helper::vector< unsigned int > &shell0 = m_container->m_edgesAroundVertex[ point0 ];
            shell0.erase( std::remove( shell0.begin(), shell0.end(), indices[i] ), shell0.end() );
            if(removeIsolatedItems && shell0.empty())
            {
                vertexToBeRemoved.push_back(point0);
            }

            sofa::helper::vector< unsigned int > &shell1 = m_container->m_edgesAroundVertex[ point1 ];
            shell1.erase( std::remove( shell1.begin(), shell1.end(), indices[i] ), shell1.end() );
            if(removeIsolatedItems && shell1.empty())
            {
                vertexToBeRemoved.push_back(point1);
            }

            if(indices[i] < lastEdgeIndex)
            {
                //replaces the edge index oldEdgeIndex with indices[i] for the first vertex
                sofa::helper::vector< unsigned int > &shell2 = m_container->m_edgesAroundVertex[ point2 ];
                replace(shell2.begin(), shell2.end(), lastEdgeIndex, indices[i]);

                //replaces the edge index oldEdgeIndex with indices[i] for the second vertex
                sofa::helper::vector< unsigned int > &shell3 = m_container->m_edgesAroundVertex[ point3 ];
                replace(shell3.begin(), shell3.end(), lastEdgeIndex, indices[i]);
            }
        }

        // removes the edge from the edgelist
        m_edge[ indices[i] ] = m_edge[ lastEdgeIndex ]; // overwriting with last valid value.
        m_edge.resize( lastEdgeIndex ); // resizing to erase multiple occurence of the edge.
    }

    if (! vertexToBeRemoved.empty())
    {
        removePointsWarning(vertexToBeRemoved);
        // inform other objects that the points are going to be removed
        propagateTopologicalChanges();
        removePointsProcess(vertexToBeRemoved);
    }
}

void EdgeSetTopologyModifier::addPointsProcess(const unsigned int nPoints)
{
    // start by calling the parent's method.
    PointSetTopologyModifier::addPointsProcess( nPoints );

    if(m_container->hasEdgesAroundVertex())
        m_container->m_edgesAroundVertex.resize( m_container->getNbPoints() );
}

void EdgeSetTopologyModifier::removePointsProcess(const sofa::helper::vector<unsigned int> &indices,
        const bool removeDOF)
{
    // Note: edges connected to the points being removed are not removed here (this situation should not occur)

    if(m_container->hasEdges())
    {
        helper::WriteAccessor< Data< sofa::helper::vector<Edge> > > m_edge = m_container->d_edge;
        // forces the construction of the edge shell array if it does not exists
        if(!m_container->hasEdgesAroundVertex())
            m_container->createEdgesAroundVertexArray();

        unsigned int lastPoint = m_container->getNbPoints() - 1;
        for (unsigned int i=0; i<indices.size(); ++i, --lastPoint)
        {
            // updating the edges connected to the point replacing the removed one:
            // for all edges connected to the last point
            for (unsigned int j=0; j<m_container->m_edgesAroundVertex[lastPoint].size(); ++j)
            {
                const int edgeId = m_container->m_edgesAroundVertex[lastPoint][j];
                // change the old index for the new one
                if ( m_edge[ edgeId ][0] == lastPoint )
                    m_edge[ edgeId ][0] = indices[i];
                else
                    m_edge[ edgeId ][1] = indices[i];
            }

            // updating the edge shell itself (change the old index for the new one)
            m_container->m_edgesAroundVertex[ indices[i] ] = m_container->m_edgesAroundVertex[ lastPoint ];
        }

        m_container->m_edgesAroundVertex.resize( m_container->m_edgesAroundVertex.size() - indices.size() );
    }

    // Important : the points are actually deleted from the mechanical object's state vectors iff (removeDOF == true)
    // call the parent method.
    PointSetTopologyModifier::removePointsProcess( indices, removeDOF );
}

void EdgeSetTopologyModifier::renumberPointsProcess( const sofa::helper::vector<unsigned int> &index,
        const sofa::helper::vector<unsigned int> &inv_index,
        const bool renumberDOF)
{
    if(m_container->hasEdges())
    {
        helper::WriteAccessor< Data< sofa::helper::vector<Edge> > > m_edge = m_container->d_edge;
        if(m_container->hasEdgesAroundVertex())
        {
            // copy of the the edge vertex shell array
            sofa::helper::vector< sofa::helper::vector< unsigned int > > edgesAroundVertex_cp = m_container->getEdgesAroundVertexArray();

            for (unsigned int i=0; i<index.size(); ++i)
            {
                m_container->m_edgesAroundVertex[i] = edgesAroundVertex_cp[ index[i] ];
            }
        }

        for (unsigned int i=0; i<m_edge.size(); ++i)
        {
            const unsigned int p0 = inv_index[ m_edge[i][0]  ];
            const unsigned int p1 = inv_index[ m_edge[i][1]  ];

            // FIXME : Edges should not be flipped during simulations as it will break code such as FEM storing a rest shape.
            // Commented by pierre-jean.bensoussan@digital-trainers.com
            /*
            if(p0 < p1)
            {
            	m_container->m_edge[i][0] = p0;
            	m_container->m_edge[i][1] = p1;
            }
            else
            {
            	m_container->m_edge[i][0] = p1;
            	m_container->m_edge[i][1] = p0;
            }
            */

            m_edge[i][0] = p0;
            m_edge[i][1] = p1;
        }
    }

    // call the parent method
    PointSetTopologyModifier::renumberPointsProcess( index, inv_index, renumberDOF );
}

void EdgeSetTopologyModifier::swapEdgesProcess(const sofa::helper::vector< sofa::helper::vector< unsigned int > >& edgesPairs)
{
    if(!m_container->hasEdges())
        return;

    // first create the edges
    sofa::helper::vector< Edge > v;
    v.reserve(2*edgesPairs.size());

    sofa::helper::vector< unsigned int > edgeIndexList;
    edgeIndexList.reserve(2*edgesPairs.size());

    sofa::helper::vector<sofa::helper::vector<unsigned int> > ancestorsArray;
    ancestorsArray.reserve(edgesPairs.size());

    unsigned int nbEdges = m_container->getNumberOfEdges();

    for (unsigned int i=0; i<edgesPairs.size(); ++i)
    {
        const unsigned int i1 = edgesPairs[i][0];
        const unsigned int i2 = edgesPairs[i][1];

        const unsigned int p11 = m_container->getEdge(i1)[0];
        const unsigned int p12 = m_container->getEdge(i1)[1];
        const unsigned int p21 = m_container->getEdge(i2)[0];
        const unsigned int p22 = m_container->getEdge(i2)[1];

        const Edge e1(p11, p21), e2(p12, p22);

        v.push_back(e1);
        v.push_back(e2);
        edgeIndexList.push_back(nbEdges);
        edgeIndexList.push_back(nbEdges+1);
        nbEdges += 2;

        sofa::helper::vector<unsigned int> ancestors(2);
        ancestors[0] = i1;
        ancestors[1] = i2;
        ancestorsArray.push_back(ancestors);
    }

    addEdgesProcess( v );

    // now warn about the creation
    addEdgesWarning((unsigned int)v.size(), v, edgeIndexList, ancestorsArray);

    // now warn about the destruction of the old edges
    sofa::helper::vector< unsigned int > indices;
    indices.reserve(2*edgesPairs.size());
    for (unsigned int i=0; i<edgesPairs.size(); ++i)
    {
        indices.push_back( edgesPairs[i][0]  );
        indices.push_back( edgesPairs[i][1] );
    }
    removeEdgesWarning(indices );

    // propagate the warnings
    propagateTopologicalChanges();

    // now destroy the old edges.
    removeEdgesProcess( indices );
}

void EdgeSetTopologyModifier::fuseEdgesProcess(const sofa::helper::vector< sofa::helper::vector< unsigned int > >& edgesPairs,
        const bool removeIsolatedPoints)
{
    if(!m_container->hasEdges())
        return;

    // first create the edges
    sofa::helper::vector< Edge > v;
    v.reserve(edgesPairs.size());

    sofa::helper::vector< unsigned int > edgeIndexList;
    edgeIndexList.reserve(edgesPairs.size());

    sofa::helper::vector<sofa::helper::vector<unsigned int> > ancestorsArray;
    ancestorsArray.reserve(edgesPairs.size());

    unsigned int nbEdges=m_container->getNumberOfEdges();

    for (unsigned int i=0; i<edgesPairs.size(); ++i)
    {
        const unsigned int i1 = edgesPairs[i][0];
        const unsigned int i2 = edgesPairs[i][1];

        unsigned int p11 = m_container->getEdge(i1)[0];
        unsigned int p22 = m_container->getEdge(i2)[1];

        if(p11 == p22)
        {
            p11 = m_container->getEdge(i2)[0];
            p22 = m_container->getEdge(i1)[1];
        }

        const Edge e (p11, p22);
        v.push_back(e);

        edgeIndexList.push_back(nbEdges);
        nbEdges += 1;

        sofa::helper::vector<unsigned int> ancestors(2);
        ancestors[0] = i1;
        ancestors[1] = i2;
        ancestorsArray.push_back(ancestors);
    }

    addEdgesProcess( v );

    // now warn about the creation
    addEdgesWarning((unsigned int)v.size(), v, edgeIndexList, ancestorsArray);

    // now warn about the destruction of the old edges
    sofa::helper::vector< unsigned int > indices;
    indices.reserve(2*edgesPairs.size());
    for (unsigned int i=0; i<edgesPairs.size(); ++i)
    {
        indices.push_back( edgesPairs[i][0] );
        indices.push_back( edgesPairs[i][1] );
    }

    removeEdgesWarning( indices );

    // propagate the warnings
    propagateTopologicalChanges();

    // now destroy the old edges.
    removeEdgesProcess( indices, removeIsolatedPoints );
}

void EdgeSetTopologyModifier::splitEdgesProcess(sofa::helper::vector<unsigned int> &indices,
        const bool removeIsolatedPoints)
{
    if(!m_container->hasEdges())
        return;

    sofa::helper::vector< sofa::helper::vector< double > > defaultBaryCoefs(indices.size());

    sofa::helper::vector< sofa::helper::vector< unsigned int > > v(indices.size());

    sofa::helper::vector< Edge >  edges;
    edges.reserve(2*indices.size());

    sofa::helper::vector< unsigned int >  edgesIndex;
    edgesIndex.reserve(2*indices.size());

    unsigned int nbEdges = m_container->getNumberOfEdges();

    for (unsigned int i=0; i<indices.size(); ++i)
    {
        const unsigned int p1 = m_container->getEdge( indices[i] )[0];
        const unsigned int p2 = m_container->getEdge( indices[i] )[1];

        // Adding the new point
        v[i].resize(2);
        v[i][0] = p1;
        v[i][1] = p2;

        // Adding the new Edges
        const Edge e1( p1, m_container->getNbPoints() + i );
        const Edge e2( m_container->getNbPoints() + i, p2 );
        edges.push_back( e1 );
        edges.push_back( e2 );
        edgesIndex.push_back(nbEdges++);
        edgesIndex.push_back(nbEdges++);

        defaultBaryCoefs[i].resize(2, 0.5f);
    }

    addPointsProcess((unsigned int)indices.size());

    addEdgesProcess( edges );

    // warn about added points and edges
    addPointsWarning((unsigned int)indices.size(), v, defaultBaryCoefs);

    addEdgesWarning((unsigned int)edges.size(), edges, edgesIndex);

    // warn about old edges about to be removed
    removeEdgesWarning( indices );

    propagateTopologicalChanges();

    // Removing the old edges
    removeEdgesProcess( indices, removeIsolatedPoints );
}

void EdgeSetTopologyModifier::splitEdgesProcess(sofa::helper::vector<unsigned int> &indices,
        const sofa::helper::vector< sofa::helper::vector< double > >& baryCoefs,
        const bool removeIsolatedPoints)
{
    if(!m_container->hasEdges())
        return;

    sofa::helper::vector< sofa::helper::vector< unsigned int > > v(indices.size());

    sofa::helper::vector< Edge >  edges;
    edges.reserve(2*indices.size());

    sofa::helper::vector< unsigned int >  edgesIndex;
    edgesIndex.reserve(2*indices.size());

    unsigned int nbEdges = m_container->getNumberOfEdges();

    for (unsigned int i=0; i<indices.size(); ++i)
    {
        const unsigned int p1 = m_container->getEdge( indices[i] )[0];
        const unsigned int p2 = m_container->getEdge( indices[i] )[1];

        // Adding the new point
        v[i].resize(2);
        v[i][0] = p1;
        v[i][1] = p2;

        // Adding the new Edges
        const Edge e1( p1, m_container->getNbPoints() + i );
        const Edge e2( m_container->getNbPoints() + i, p2 );
        edges.push_back( e1 );
        edges.push_back( e2 );
        edgesIndex.push_back(nbEdges++);
        edgesIndex.push_back(nbEdges++);
    }

    addPointsProcess((unsigned int)indices.size());

    addEdgesProcess( edges );

    // warn about added points and edges
    addPointsWarning((unsigned int)indices.size(), v, baryCoefs);

    addEdgesWarning((unsigned int)edges.size(), edges, edgesIndex);

    // warn about old edges about to be removed
    removeEdgesWarning( indices );

    propagateTopologicalChanges();

    // Removing the old edges
    removeEdgesProcess( indices, removeIsolatedPoints );
}

void EdgeSetTopologyModifier::removeEdges(const sofa::helper::vector< unsigned int >& edgeIds,
        const bool removeIsolatedPoints, const bool resetTopoChange)
{
    sofa::helper::vector<unsigned int> edgeIds_filtered;
    for (unsigned int i = 0; i < edgeIds.size(); i++)
    {
        if( edgeIds[i] >= m_container->getNumberOfEdges())
            std::cout << "Error: EdgeSetTopologyModifier::removeEdges: edges: "<< edgeIds[i] <<" is out of bound and won't be removed." << std::endl;
        else
            edgeIds_filtered.push_back(edgeIds[i]);
    }

    /// add the topological changes in the queue
    removeEdgesWarning(edgeIds_filtered);
    // inform other objects that the edges are going to be removed
    if (resetTopoChange)
        propagateTopologicalChanges();
    else
        propagateTopologicalChangesWithoutReset();
    // now destroy the old edges.
    removeEdgesProcess( edgeIds_filtered, removeIsolatedPoints );

    m_container->checkTopology();
}

void EdgeSetTopologyModifier::removeItems(const sofa::helper::vector< unsigned int >& items)
{
    removeEdges(items);
}

void EdgeSetTopologyModifier::renumberPoints( const sofa::helper::vector<unsigned int> &index,
        const sofa::helper::vector<unsigned int> &inv_index)
{
    /// add the topological changes in the queue
    renumberPointsWarning(index, inv_index);
    // inform other objects that the triangles are going to be removed
    propagateTopologicalChanges();
    // now renumber the points
    renumberPointsProcess(index, inv_index);

    m_container->checkTopology();
}

void EdgeSetTopologyModifier::addEdges(const sofa::helper::vector< Edge >& edges)
{
    unsigned int nEdges = m_container->getNumberOfEdges();

    /// actually add edges in the topology container
    addEdgesProcess(edges);

    sofa::helper::vector<unsigned int> edgesIndex;
    edgesIndex.reserve(edges.size());

    for (unsigned int i=0; i<edges.size(); ++i)
    {
        edgesIndex.push_back(nEdges+i);
    }

    // add topology event in the stack of topological events
    addEdgesWarning((unsigned int)edges.size(), edges, edgesIndex);

    // inform other objects that the edges are already added
    propagateTopologicalChanges();
}

void EdgeSetTopologyModifier::addEdges(const sofa::helper::vector< Edge >& edges,
        const sofa::helper::vector< sofa::helper::vector< unsigned int > > & ancestors,
        const sofa::helper::vector< sofa::helper::vector< double > >& baryCoefs)
{
    unsigned int nEdges = m_container->getNumberOfEdges();

    /// actually add edges in the topology container
    addEdgesProcess(edges);

    sofa::helper::vector<unsigned int> edgesIndex;
    edgesIndex.reserve(edges.size());

    for (unsigned int i=0; i<edges.size(); ++i)
    {
        edgesIndex.push_back(nEdges+i);
    }

    // add topology event in the stack of topological events
    addEdgesWarning((unsigned int)edges.size(), edges, edgesIndex, ancestors, baryCoefs);

    // inform other objects that the edges are already added
    propagateTopologicalChanges();
}

void EdgeSetTopologyModifier::addEdges(const sofa::helper::vector< Edge >& edges,
        const sofa::helper::vector< core::topology::EdgeAncestorElem >& ancestorElems)
{
    unsigned int nEdges = m_container->getNumberOfEdges();

    assert(ancestorElems.size() == edges.size());

    /// actually add edges in the topology container
    addEdgesProcess(edges);

    sofa::helper::vector<unsigned int> edgesIndex;
    edgesIndex.resize(edges.size());
    for (unsigned int i=0; i<edges.size(); ++i)
    {
        edgesIndex[i] = nEdges+i;
    }

    // add topology event in the stack of topological events
    addEdgesWarning((unsigned int)edges.size(), edges, edgesIndex, ancestorElems);

    // inform other objects that the edges are already added
    propagateTopologicalChanges();
}

void EdgeSetTopologyModifier::swapEdges(const sofa::helper::vector< sofa::helper::vector< unsigned int > >& edgesPairs)
{
    swapEdgesProcess(edgesPairs);
    m_container->checkTopology();
}

void EdgeSetTopologyModifier::fuseEdges(const sofa::helper::vector< sofa::helper::vector< unsigned int > >& edgesPairs, const bool removeIsolatedPoints)
{
    fuseEdgesProcess(edgesPairs, removeIsolatedPoints);
    m_container->checkTopology();
}

void EdgeSetTopologyModifier::splitEdges( sofa::helper::vector<unsigned int> &indices,
        const bool removeIsolatedPoints)
{
    splitEdgesProcess(indices, removeIsolatedPoints);
    m_container->checkTopology();
}

void EdgeSetTopologyModifier::splitEdges( sofa::helper::vector<unsigned int> &indices,
        const sofa::helper::vector< sofa::helper::vector< double > >& baryCoefs,
        const bool removeIsolatedPoints)
{
    splitEdgesProcess(indices, baryCoefs, removeIsolatedPoints);
    m_container->checkTopology();
}

// Give the optimal vertex permutation according to the Reverse CuthillMckee algorithm (use BOOST GRAPH LIBRAIRY)
void EdgeSetTopologyModifier::resortCuthillMckee(sofa::helper::vector<int>& inverse_permutation)
{
    using namespace boost;
    using namespace std;
    typedef adjacency_list<vecS, vecS, undirectedS,
            property<vertex_color_t, default_color_type,
            property<vertex_degree_t,int> > > Graph;
    typedef graph_traits<Graph>::vertex_descriptor Vertex;
    typedef graph_traits<Graph>::vertices_size_type size_type;

    Graph G;

    const sofa::helper::vector<Edge> &ea=m_container->getEdgeArray();

    for (unsigned int k=0; k<ea.size(); ++k)
    {
        add_edge(ea[k][0], ea[k][1], G);
    }

    inverse_permutation.resize(num_vertices(G));

    property_map<Graph, vertex_index_t>::type index_map = get(vertex_index, G);

    sout << "original bandwidth: " << bandwidth(G) << sendl;

    std::vector<Vertex> inv_perm(num_vertices(G));
    std::vector<size_type> perm(num_vertices(G));

    //reverse cuthill_mckee_ordering
    cuthill_mckee_ordering(G, inv_perm.rbegin());

    //sout << "Reverse Cuthill-McKee ordering:" << sendl;
    //sout << "  ";
    unsigned int ind_i = 0;
    for (std::vector<Vertex>::const_iterator it = inv_perm.begin();
            it != inv_perm.end(); ++it)
    {
        //sout << index_map[*it] << " ";
        inverse_permutation[ind_i++] = (int)index_map[*it];
    }
    //sout << sendl;

    for (size_type c=0; c!=inv_perm.size(); ++c)
        perm[index_map[inv_perm[c]]] = c;

    sout << "  bandwidth: "
            << bandwidth(G, make_iterator_property_map(&perm[0], index_map, perm[0]))
            << sendl;
}

void EdgeSetTopologyModifier::movePointsProcess (const sofa::helper::vector <unsigned int>& id,
        const sofa::helper::vector< sofa::helper::vector< unsigned int > >& ancestors,
        const sofa::helper::vector< sofa::helper::vector< double > >& coefs,
        const bool moveDOF)
{
    SOFA_UNUSED(moveDOF);
    unsigned int nbrVertex = (unsigned int)id.size();
    bool doublet;
    sofa::helper::vector<unsigned int> edgesAroundVertex2Move;
    sofa::helper::vector< Edge > edgeArray;

    // Step 1/4 - Creating trianglesAroundVertex to moved due to moved points:
    for (unsigned int i = 0; i<nbrVertex; ++i)
    {
        const sofa::helper::vector <unsigned int>& edgesAroundVertex = m_container->getEdgesAroundVertexArray()[ id[i] ];

        for (unsigned int j = 0; j<edgesAroundVertex.size(); ++j)
        {
            doublet = false;

            for (unsigned int k =0; k<edgesAroundVertex2Move.size(); ++k) //Avoid double
            {
                if (edgesAroundVertex2Move[k] == edgesAroundVertex[j])
                {
                    doublet = true;
                    break;
                }
            }

            if(!doublet)
                edgesAroundVertex2Move.push_back (edgesAroundVertex[j]);

        }
    }

    std::sort( edgesAroundVertex2Move.begin(), edgesAroundVertex2Move.end(), std::greater<unsigned int>() );

    // Step 2/4 - Create event to delete all elements before moving and propagate it:
    EdgesMoved_Removing *ev1 = new EdgesMoved_Removing (edgesAroundVertex2Move);
    this->addTopologyChange(ev1);
    propagateTopologicalChanges();

    // Step 3/4 - Physically move all dof:
    PointSetTopologyModifier::movePointsProcess (id, ancestors, coefs);

    // Step 4/4 - Create event to recompute all elements concerned by moving and propagate it:

    // Creating the corresponding array of Triangles for ancestors
    for (unsigned int i = 0; i<edgesAroundVertex2Move.size(); i++)
        edgeArray.push_back (m_container->getEdgeArray()[ edgesAroundVertex2Move[i] ]);

    EdgesMoved_Adding *ev2 = new EdgesMoved_Adding (edgesAroundVertex2Move, edgeArray);
    this->addTopologyChange(ev2); // This event should be propagated with global workflow
}

bool EdgeSetTopologyModifier::removeConnectedComponents(unsigned int elemID)
{
    if(!m_container)
    {
        serr << "TopologyContainer pointer is empty." << sendl;
        return false;
    }

    sofa::helper::vector <unsigned int> elems = m_container->getConnectedElement(elemID);
    this->removeItems(elems);
    return true;
}

bool EdgeSetTopologyModifier::removeConnectedElements(unsigned int elemID)
{
    if(!m_container)
    {
        serr << "TopologyContainer pointer is empty." << sendl;
        return false;
    }

    sofa::helper::vector <unsigned int> elems = m_container->getElementAroundElement(elemID);
    this->removeItems(elems);
    return true;
}

bool EdgeSetTopologyModifier::removeIsolatedElements()
{
    return this->removeIsolatedElements(0);
}

bool EdgeSetTopologyModifier::removeIsolatedElements(unsigned int scaleElem)
{
    if(!m_container)
    {
        serr << "TopologyContainer pointer is empty." << sendl;
        return false;
    }

    unsigned int nbr = this->m_container->getNumberOfElements();
    sofa::helper::vector <unsigned int> elemAll = m_container->getConnectedElement(0);
    sofa::helper::vector <unsigned int> elem, elemMax, elemToRemove;

    if (nbr == elemAll.size()) // nothing to do
        return true;

    elemMax = elemAll;

    if (scaleElem == 0) //remove all isolated elements
        scaleElem = nbr;

    while (elemAll.size() < nbr)
    {
        std::sort(elemAll.begin(), elemAll.end());
        EdgeID other_edgeID = (EdgeID)elemAll.size();

        for (EdgeID i = 0; i<elemAll.size(); ++i)
            if (elemAll[i] != i)
            {
                other_edgeID = i;
                break;
            }

        elem = this->m_container->getConnectedElement(other_edgeID);
        elemAll.insert(elemAll.begin(), elem.begin(), elem.end());

        if (elemMax.size() < elem.size())
        {
            if (elemMax.size() <= scaleElem)
                elemToRemove.insert(elemToRemove.begin(), elemMax.begin(), elemMax.end());

            elemMax = elem;
        }
        else
        {
            if (elem.size() <= scaleElem)
                elemToRemove.insert(elemToRemove.begin(), elem.begin(), elem.end());
        }
    }

    this->removeItems(elemToRemove);

    return true;
}

void EdgeSetTopologyModifier::propagateTopologicalEngineChanges()
{
    if (m_container->beginChange() == m_container->endChange()) // nothing to do if no event is stored
        return;

    if (!m_container->isEdgeTopologyDirty()) // edge Data has not been touched
        return PointSetTopologyModifier::propagateTopologicalEngineChanges();

    sofa::helper::list <sofa::core::topology::TopologyEngine *>::iterator it;
    for ( it = m_container->m_enginesList.begin(); it!=m_container->m_enginesList.end(); ++it)
    {
        sofa::core::topology::TopologyEngine* topoEngine = (*it);
        if (topoEngine->isDirty())
        {
#ifndef NDEBUG
            std::cout << "EdgeSetTopologyModifier::performing: " << topoEngine->getName() << std::endl;
#endif
            topoEngine->requestUpdate();
        }
    }

    m_container->cleanEdgeTopologyFromDirty();
    PointSetTopologyModifier::propagateTopologicalEngineChanges();
}

} // namespace topology

} // namespace component

} // namespace sofa

max_stars_repo_name: sofa-framework/issofa

max_issues_repo_name: sofa-framework/issofa

max_forks_repo_name: sofa-framework/issofa


text:
/* Copyright 2017 Battelle Energy Alliance, LLC

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/
#include "randomClass.h"
#include <boost/random/mersenne_twister.hpp>

class RandomClassImpl {
  /**
   * This class create the instance of a mersenne-twister random number generator (BOOST)
   */
public:
  boost::random::mt19937 _backend;
};

RandomClass::RandomClass() : _rng(new RandomClassImpl()), _range(_rng->_backend.max() - _rng->_backend.min()) {
}

void RandomClass::seed(unsigned int seed) {
    _rng->_backend.seed(seed);
  }

double RandomClass::random() {
    return (_rng->_backend()-_rng->_backend.min())/_range;
  }

RandomClass::~RandomClass(){
  delete _rng;
}

ext: cxx

max_stars_repo_name: milljm/raven

max_issues_repo_name: milljm/raven

max_forks_repo_name: milljm/raven


text:
#include "interpreter/exception.h"
#include "interpreter/handler.h"
#include "interpreter/ffi_conversion.h"

#include <boost/python.hpp>
#include <boost/python/numpy.hpp>

namespace Asampl::Interpreter::Handler {

namespace py = boost::python;
namespace np = py::numpy;

class PythonHandler;

class PythonHandlerDownload : public IHandlerContextDownload {
public:
    PythonHandlerDownload(PythonHandler& handler);

    void push(const std::vector<uint8_t>& data) override;
    DownloadResponse download() override;

private:
    py::object context_;
};

class PythonHandlerUpload : public IHandlerContextUpload {
public:
    PythonHandlerUpload(PythonHandler& handler);

    UploadResponse pull() override;
    void upload(const ValuePtr& value) override;

private:
    py::object context_;
};

class PythonHandler : public IHandler {
public:
    PythonHandler(std::filesystem::path path)
    {
        path.replace_extension("py");
        py::exec_file(path.c_str(), ns);
        download_class = ns["Download"];
        upload_class = ns["Upload"];
    }

    std::unique_ptr<IHandlerContextDownload> open_download() override {
        return std::unique_ptr<PythonHandlerDownload>(new PythonHandlerDownload(*this));
    }

    std::unique_ptr<IHandlerContextUpload> open_upload() override {
        throw InterpreterException("FFI upload is not supported yet");
    }

private:
    py::dict ns;
    py::object download_class;
    py::object upload_class;

private:
    friend class PythonHandlerDownload;
    friend class PythonHandlerUpload;
};

PythonHandlerDownload::PythonHandlerDownload(PythonHandler& handler)
    : context_{handler.download_class()}
{
}

void PythonHandlerDownload::push(const std::vector<uint8_t>& data) {
    auto ptr = reinterpret_cast<uint8_t*>(malloc(data.size()));
    std::copy(data.begin(), data.end(), ptr);

    auto array = np::from_data(
        ptr,
        np::dtype::get_builtin<uint8_t>(),
        py::make_tuple(data.size()), py::make_tuple(1),
        py::object());

    context_.attr("push")(array);
}

DownloadResponse PythonHandlerDownload::download() {
    py::object response = context_.attr("download")();
    if (response.contains("status")) {
        if (response["status"] == "fatal") {
            throw InterpreterException(py::extract<std::string>(response["error"]));
        }
        if (response["status"] == "not_ready") {
            return DownloadResponse::new_not_ready();
        }
        if (response["status"] == "out_of_data") {
            return DownloadResponse::new_out_of_data();
        }
        throw InterpreterException("Unknown response status");
    }

    const double timestamp = py::extract<double>(response["timestamp"]);
    auto value = convert_from_python(response["value"]);
    return DownloadResponse::new_valid(std::move(value), timestamp);
}

PythonHandlerUpload::PythonHandlerUpload(PythonHandler& handler)
    : context_(handler.upload_class())
{
}

UploadResponse PythonHandlerUpload::pull() {
    py::object response = context_.attr("pull")();
    if (auto extract = py::extract<np::ndarray>(response); extract.check()) {
        np::ndarray array = extract();
        if (array.get_nd() != 1) {
            throw InterpreterException("Wrong number of dimensions. Expected 1");
        }

        std::vector<uint8_t> data;
        data.resize(array.shape(0));
        std::copy(array.get_data(), array.get_data() + array.shape(0), data.data());

        return UploadResponse::new_valid(std::move(data));
    } else if (response.is_none()) {
        return UploadResponse::new_not_ready();
    } else {
        throw InterpreterException("Unexpected return value from handler");
    }
}

void PythonHandlerUpload::upload(const ValuePtr& value) {
    context_.attr("upload")(convert_to_python(value));
}

}

max_stars_repo_name: Asampl-development-team/Asampl

max_issues_repo_name: Asampl-development-team/Asampl

max_forks_repo_name: Asampl-development-team/Asampl