366 lines
13 KiB
C++
366 lines
13 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's C++ library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: Root_Density.cc
|
|
|
|
Purpose: defines Root_Density & Root_Density_Set, and their member functions
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
|
|
Root_Density::Root_Density ()
|
|
: Npts(1), lambdamax(0.0), lambda(Vect_DP(0.0, Npts)), dlambda(Vect_DP(0.0, Npts)), value(Vect_DP(0.0, Npts)),
|
|
prev_value(Vect_DP(0.0, Npts)), diff(1.0e+6), value_infty_set(false), value_infty(0.0)
|
|
{
|
|
}
|
|
|
|
Root_Density::Root_Density (int Npts_ref, DP lambdamax_ref)
|
|
: Npts(Npts_ref), lambdamax(lambdamax_ref), lambda(Vect_DP(Npts)), dlambda(Vect_DP(0.0, Npts)), value(Vect_DP(0.0, Npts)),
|
|
prev_value(Vect_DP(0.0, Npts)), diff(1.0e+6), value_infty_set(false), value_infty(0.0)
|
|
{
|
|
for (int i = 0; i < value.size(); ++i) {
|
|
lambda[i] = lambdamax * (-(Npts - 1.0) + 2*i)/Npts;
|
|
dlambda[i] = 2.0 * lambdamax/Npts;
|
|
}
|
|
}
|
|
|
|
Root_Density& Root_Density::operator= (const Root_Density& RefDensity)
|
|
{
|
|
if (this != &RefDensity) {
|
|
Npts = RefDensity.Npts;
|
|
lambdamax = RefDensity.lambdamax;
|
|
lambda = RefDensity.lambda;
|
|
dlambda = RefDensity.dlambda;
|
|
value = RefDensity.value;
|
|
prev_value = RefDensity.prev_value;
|
|
diff = RefDensity.diff;
|
|
value_infty_set = RefDensity.value_infty_set;
|
|
value_infty = RefDensity.value_infty;
|
|
|
|
}
|
|
return(*this);
|
|
}
|
|
|
|
DP Root_Density::Return_Value (DP lambda_ref)
|
|
{
|
|
// This function returns a value for epsilon at any real lambda
|
|
// using simple linear interpolation.
|
|
// Degree 3 polynomical also programmed in, but commented out: no improvement.
|
|
|
|
DP answer = 0.0;
|
|
if (fabs(lambda_ref) >= fabs(lambda[0])) {
|
|
if (value_infty_set) answer = value_infty;
|
|
else ABACUSerror("Need to set asymptotics of Root_Density !");
|
|
}
|
|
|
|
else { // try to find the i such that lambda[i] <= lambda_ref < lambda[i+1]
|
|
|
|
int index = (Npts - 1)/2;
|
|
int indexstep = (Npts - 1)/4 + 1;
|
|
|
|
while (indexstep >= 1) {
|
|
|
|
if ( // if is "lower": we go up
|
|
lambda_ref >= lambda[index + 1]) {
|
|
index += indexstep;
|
|
}
|
|
|
|
else if ( // if is "higher" or equal: we go down
|
|
lambda[index] > lambda_ref) {
|
|
index -= indexstep;
|
|
}
|
|
|
|
index = ABACUS::max(0, index);
|
|
index = ABACUS::min(Npts - 2, index);
|
|
|
|
if (indexstep == 1) indexstep--;
|
|
else indexstep = (indexstep + 1)/2;
|
|
|
|
} // while ...
|
|
|
|
if (index < 0 || index >= Npts || lambda[index] > lambda_ref || lambda[index + 1] < lambda_ref) {
|
|
cout << "Seeking index: " << index << "\t" << lambda[index] << "\t <=? " << lambda_ref
|
|
<< "\t<? " << lambda[index + 1] << endl;
|
|
ABACUSerror("Calculating index wrong in Root_Density::Evaluate.");
|
|
}
|
|
|
|
answer = ((value[index] * (lambda[index+1] - lambda_ref)
|
|
+ value[index + 1] * (lambda_ref - lambda[index]))/(lambda[index+1] - lambda[index]));
|
|
}
|
|
|
|
return(answer);
|
|
}
|
|
|
|
void Root_Density::Set_Asymptotics (DP value_infty_ref)
|
|
{
|
|
value_infty = value_infty_ref;
|
|
value_infty_set = true;
|
|
}
|
|
|
|
Root_Density Root_Density::Compress_and_Match_Densities (DP comp_factor)
|
|
{
|
|
// Returns a 'compressed' version of the density, using 1/comp_factor as many points.
|
|
|
|
int Npts_used = int(2.0 * lambdamax/(dlambda[0] * comp_factor));
|
|
|
|
Root_Density compressed_density(Npts_used, lambdamax);
|
|
|
|
compressed_density.Set_Asymptotics (value_infty);
|
|
|
|
for (int i = 0; i < compressed_density.Npts; ++i)
|
|
compressed_density.value[i] = (*this).Return_Value (compressed_density.lambda[i]);
|
|
|
|
return(compressed_density);
|
|
}
|
|
|
|
void Root_Density::Save (const char* outfile_Cstr)
|
|
{
|
|
ofstream outfile;
|
|
outfile.open(outfile_Cstr);
|
|
outfile.precision(16);
|
|
|
|
for (int i = 0; i < Npts; ++i) {
|
|
if (i > 0) outfile << endl;
|
|
outfile << setw(20) << lambda[i] << "\t" << setw(20) << value[i];
|
|
}
|
|
|
|
outfile.close();
|
|
}
|
|
|
|
|
|
Root_Density_Set::Root_Density_Set () : ntypes(1), epsilon(Vect<Root_Density> (ntypes)), Npts_total(0), diff(1.0e+6)
|
|
{
|
|
}
|
|
|
|
Root_Density_Set::Root_Density_Set (int ntypes_ref, int Npts_ref, DP lambdamax_ref)
|
|
: ntypes(ntypes_ref), epsilon(Vect<Root_Density> (ntypes_ref)), Npts_total(ntypes_ref * Npts_ref), diff(1.0e+6)
|
|
{
|
|
for (int n = 0; n < ntypes; ++n) epsilon[n] = Root_Density(Npts_ref, lambdamax_ref);
|
|
}
|
|
|
|
Root_Density_Set::Root_Density_Set (int ntypes_ref, Vect_INT Npts_ref, Vect_DP lambdamax_ref)
|
|
: ntypes(ntypes_ref), epsilon(Vect<Root_Density> (ntypes_ref)), Npts_total(Npts_ref.sum()), diff(1.0e+6)
|
|
{
|
|
if (Npts_ref.size() != ntypes_ref || lambdamax_ref.size() != ntypes_ref)
|
|
ABACUSerror("Wrong vector sizes in Root_Density_Set.");
|
|
for (int n = 0; n < ntypes; ++n) epsilon[n] = Root_Density(Npts_ref[n], lambdamax_ref[n]);
|
|
}
|
|
|
|
Root_Density_Set& Root_Density_Set::operator= (const Root_Density_Set& RefSet)
|
|
{
|
|
if (this != &RefSet) {
|
|
ntypes = RefSet.ntypes;
|
|
epsilon = RefSet.epsilon;
|
|
Npts_total = RefSet.Npts_total;
|
|
diff = RefSet.diff;
|
|
}
|
|
return(*this);
|
|
}
|
|
|
|
void Root_Density_Set::Insert_new_function (DP asymptotic_value)
|
|
{
|
|
// This function extends a set by adding one epsilon_n function on top
|
|
|
|
Root_Density_Set Updated_Set (ntypes + 1, 10, 10.0); // last two parameters are meaningless
|
|
for (int n = 0; n < ntypes; ++n) Updated_Set.epsilon[n] = epsilon[n];
|
|
|
|
Updated_Set.epsilon[ntypes] = Root_Density (50, epsilon[ntypes - 1].lambdamax);
|
|
Updated_Set.epsilon[ntypes].Set_Asymptotics (asymptotic_value);
|
|
|
|
for (int i = 0; i < Updated_Set.epsilon[ntypes].Npts; ++i)
|
|
Updated_Set.epsilon[ntypes].value[i] = Updated_Set.epsilon[ntypes].value_infty;
|
|
|
|
ntypes = Updated_Set.ntypes;
|
|
epsilon = Updated_Set.epsilon;
|
|
Npts_total+= Updated_Set.epsilon[ntypes - 1].Npts;
|
|
}
|
|
|
|
void Root_Density_Set::Extend_limits (Vect<bool> need_to_extend_limit)
|
|
{
|
|
// Extend the limits of integration at each level, according to boolean
|
|
|
|
// The function extends the limits by 10% on both sides, putting the
|
|
// extra values to value_infty.
|
|
|
|
if (need_to_extend_limit.size() != epsilon.size())
|
|
ABACUSerror("Wrong size need_to_extend_limit boolean in Extend_limits.");
|
|
|
|
Vect_INT nr_new_points_needed(0, ntypes);
|
|
int total_nr_new_points_added = 0;
|
|
DP dlambda_used = 0.0;
|
|
for (int n = 0; n < ntypes; ++n) {
|
|
if (need_to_extend_limit[n]) {
|
|
|
|
Root_Density epsilon_n_before_update = epsilon[n];
|
|
|
|
// Determine the dlambda to be used:
|
|
dlambda_used = epsilon[n].dlambda[0];
|
|
|
|
// How many new points do we add ? Say 5\% on each side:
|
|
nr_new_points_needed[n] = ABACUS::max(1, epsilon[n].Npts/20);
|
|
|
|
epsilon[n] = Root_Density(epsilon_n_before_update.Npts + 2* nr_new_points_needed[n],
|
|
epsilon_n_before_update.lambdamax + nr_new_points_needed[n] * dlambda_used);
|
|
epsilon[n].Set_Asymptotics(epsilon_n_before_update.value_infty);
|
|
|
|
for (int i = 0; i < nr_new_points_needed[n]; ++i) {
|
|
epsilon[n].lambda[i] = epsilon_n_before_update.lambda[0] - (nr_new_points_needed[n] - i) * dlambda_used;
|
|
epsilon[n].dlambda[i] = dlambda_used;
|
|
epsilon[n].value[i] = epsilon_n_before_update.value_infty;
|
|
}
|
|
|
|
for (int i = 0; i < epsilon_n_before_update.Npts; ++i) {
|
|
epsilon[n].lambda[i + nr_new_points_needed[n] ] = epsilon_n_before_update.lambda[i];
|
|
epsilon[n].dlambda[i + nr_new_points_needed[n] ] = epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].value[i + nr_new_points_needed[n] ] = epsilon_n_before_update.value[i];
|
|
}
|
|
|
|
for (int i = 0; i < nr_new_points_needed[n]; ++i) {
|
|
epsilon[n].lambda[i + epsilon_n_before_update.Npts + nr_new_points_needed[n] ]
|
|
= epsilon_n_before_update.lambda[epsilon_n_before_update.Npts - 1] + (i+1.0) * dlambda_used;
|
|
epsilon[n].dlambda[i + epsilon_n_before_update.Npts + nr_new_points_needed[n] ] = dlambda_used;
|
|
epsilon[n].value[i + epsilon_n_before_update.Npts + nr_new_points_needed[n] ] = epsilon_n_before_update.value_infty;
|
|
}
|
|
|
|
total_nr_new_points_added += 2 * nr_new_points_needed[n];
|
|
|
|
} // if (need
|
|
} // for n
|
|
|
|
Npts_total += total_nr_new_points_added;
|
|
|
|
// Done !
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
void Root_Density_Set::Insert_new_points (Vect<Vect<bool> > need_new_point_around)
|
|
{
|
|
// need_new_point_around specifies whether a new point needs to be inserted around existing points.
|
|
|
|
// Count the number of new points needed per type:
|
|
Vect_INT nr_new_points_needed(0, ntypes);
|
|
int total_nr_new_points_needed = 0;
|
|
for (int n = 0; n < ntypes; ++n) {
|
|
if (need_new_point_around[n].size() != epsilon[n].Npts)
|
|
ABACUSerror("Wrong size need_new_point_around boolean in Insert_new_points.");
|
|
for (int i = 0; i < epsilon[n].Npts; ++i)
|
|
if (need_new_point_around[n][i]) nr_new_points_needed[n]++;
|
|
total_nr_new_points_needed += nr_new_points_needed[n];
|
|
}
|
|
|
|
// Working version using non-equispaced points
|
|
// Now update all data via interpolation:
|
|
for (int n = 0; n < ntypes; ++n) {
|
|
Root_Density epsilon_n_before_update = epsilon[n];
|
|
epsilon[n] = Root_Density(epsilon_n_before_update.Npts + nr_new_points_needed[n], epsilon_n_before_update.lambdamax);
|
|
epsilon[n].Set_Asymptotics(epsilon_n_before_update.value_infty);
|
|
int nr_pts_added_n = 0;
|
|
for (int i = 0; i < epsilon_n_before_update.Npts; ++i) {
|
|
if (!need_new_point_around[n][i]) {
|
|
epsilon[n].lambda[i + nr_pts_added_n] = epsilon_n_before_update.lambda[i];
|
|
epsilon[n].dlambda[i + nr_pts_added_n] = epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].value[i + nr_pts_added_n] = epsilon_n_before_update.value[i];
|
|
}
|
|
else if (need_new_point_around[n][i]) {
|
|
epsilon[n].lambda[i + nr_pts_added_n] = epsilon_n_before_update.lambda[i] - 0.25 * epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].dlambda[i + nr_pts_added_n] = 0.5 * epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].value[i + nr_pts_added_n] = epsilon_n_before_update.Return_Value(epsilon[n].lambda[i + nr_pts_added_n]);
|
|
nr_pts_added_n++;
|
|
epsilon[n].lambda[i + nr_pts_added_n] = epsilon_n_before_update.lambda[i] + 0.25 * epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].dlambda[i + nr_pts_added_n] = 0.5 * epsilon_n_before_update.dlambda[i];
|
|
epsilon[n].value[i + nr_pts_added_n] = epsilon_n_before_update.Return_Value(epsilon[n].lambda[i + nr_pts_added_n]);
|
|
}
|
|
}
|
|
if (nr_pts_added_n != nr_new_points_needed[n]) {
|
|
cout << nr_pts_added_n << "\t" << nr_new_points_needed[n] << endl;
|
|
ABACUSerror("Wrong counting of new points in Insert_new_points.");
|
|
}
|
|
|
|
} // for n
|
|
|
|
Npts_total += total_nr_new_points_needed;
|
|
|
|
// Done !
|
|
|
|
return;
|
|
}
|
|
|
|
DP Root_Density_Set::Return_Value (int n_ref, DP lambda_ref)
|
|
{
|
|
// Returns a value, no matter what !
|
|
|
|
if (n_ref < ntypes) return(epsilon[n_ref].Return_Value(lambda_ref));
|
|
|
|
else // assume asymptotic form of epsilon, proportional to n
|
|
return(epsilon[ntypes - 1].Return_Value(lambda_ref) * n_ref/(ntypes - 1.0));
|
|
|
|
}
|
|
|
|
Root_Density_Set Root_Density_Set::Return_Compressed_and_Matched_Set (DP comp_factor)
|
|
{
|
|
// Returns a set with 1/comp_factor as many points at each level
|
|
|
|
if (comp_factor >= 2.0)
|
|
ABACUSerror("Compression factor too large in Return_Compressed_and_Matched_Set, numerical instability will occur.");
|
|
|
|
Vect_INT nrpts_comp (ntypes);
|
|
Vect_DP lambdamax_comp (ntypes);
|
|
for (int n = 0; n < ntypes; ++n) {
|
|
nrpts_comp[n] = int(2.0 * epsilon[n].lambdamax/(epsilon[n].dlambda[0] * comp_factor));
|
|
lambdamax_comp[n] = epsilon[n].lambdamax;
|
|
}
|
|
|
|
Root_Density_Set Compressed_and_Matched_Set (ntypes, nrpts_comp, lambdamax_comp);
|
|
|
|
for (int n = 0; n < ntypes; ++n)
|
|
Compressed_and_Matched_Set.epsilon[n] = (*this).epsilon[n].Compress_and_Match_Densities (comp_factor);
|
|
|
|
return(Compressed_and_Matched_Set);
|
|
}
|
|
|
|
void Root_Density_Set::Match_Densities (Root_Density_Set& RefSet)
|
|
{
|
|
// matched densities to those in RefSet
|
|
|
|
for (int n = 0; n < ntypes; ++n)
|
|
for (int i = 0; i < epsilon[n].Npts; ++i)
|
|
epsilon[n].value[i] = RefSet.epsilon[n].Return_Value(epsilon[n].lambda[i]);
|
|
}
|
|
|
|
void Root_Density_Set::Save (const char* outfile_Cstr)
|
|
{
|
|
ofstream outfile;
|
|
outfile.open(outfile_Cstr);
|
|
outfile.precision(16);
|
|
|
|
// Determine what the maximal nr of pts is:
|
|
int Npts_n_max = 0;
|
|
for (int n = 0; n < ntypes; ++n) Npts_n_max = ABACUS::max(Npts_n_max, epsilon[n].Npts);
|
|
|
|
for (int i = 0; i < Npts_n_max; ++i) {
|
|
if (i > 0) outfile << endl;
|
|
for (int n = 0; n < ntypes; ++n) (i < epsilon[n].Npts) ?
|
|
(outfile << epsilon[n].lambda[i] << "\t" << epsilon[n].value[i] << "\t")
|
|
: (outfile << 0 << "\t" << 0 << "\t");
|
|
}
|
|
|
|
outfile.close();
|
|
}
|
|
|
|
|
|
} // namespace ABACUS
|