ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4171 from akohlmey/sort-vector

Browse Source 
Add special variable functions sort() and rsort() for sorting vectors by value
This commit is contained in:
Stan Moore
2024-05-28 16:22:08 -06:00
committed by GitHub
parent dedcfa157c 3bc367e0b0
commit e8ee0d9e71
6 changed files with 137 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
doc/src/variable.rst
View File

@ -67,7 +67,7 @@ Syntax
bound(group,dir,region), gyration(group,region), ke(group,reigon),
angmom(group,dim,region), torque(group,dim,region),
inertia(group,dimdim,region), omega(group,dim,region)
special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)
special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), rsort(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)
feature functions = is_available(category,feature), is_active(category,feature), is_defined(category,id)
atom value = id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i]
atom vector = id, mass, type, mol, radius, q, x, y, z, vx, vy, vz, fx, fy, fz
@ -547,7 +547,7 @@ variables.
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Region functions | count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR) |
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Special functions | sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label) |
| Special functions | sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), rsort(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label) |
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
| Feature functions | is_available(category,feature), is_active(category,feature), is_defined(category,id) |
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
@ -913,23 +913,27 @@ Special Functions
Special functions take specific kinds of arguments, meaning their
arguments cannot be formulas themselves.
The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
each take 1 argument which is of the form "c_ID" or "c_ID[N]" or
"f_ID" or "f_ID[N]" or "v_name". The first two are computes and the
second two are fixes; the ID in the reference should be replaced by
the ID of a compute or fix defined elsewhere in the input script. The
compute or fix must produce either a global vector or array. If it
produces a global vector, then the notation without "[N]" should be
used. If it produces a global array, then the notation with "[N]"
should be used, when N is an integer, to specify which column of the
global array is being referenced. The last form of argument "v_name"
is for a vector-style variable where "name" is replaced by the name of
the variable.
The sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), and
rsort(x) functions each take 1 argument which is of the form "c_ID" or
"c_ID[N]" or "f_ID" or "f_ID[N]" or "v_name". The first two are
computes and the second two are fixes; the ID in the reference should be
replaced by the ID of a compute or fix defined elsewhere in the input
script. The compute or fix must produce either a global vector or
array. If it produces a global vector, then the notation without "[N]"
should be used. If it produces a global array, then the notation with
"[N]" should be used, where N is an integer, to specify which column of
the global array is being referenced. The last form of argument
"v_name" is for a vector-style variable where "name" is replaced by the
name of the variable.
These functions operate on a global vector of inputs and reduce it to
a single scalar value. This is analogous to the operation of the
:doc:`compute reduce <compute_reduce>` command, which performs similar
operations on per-atom and local vectors.
The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
operate on a global vector of inputs and reduce it to a single scalar
value. This is analogous to the operation of the :doc:`compute reduce
<compute_reduce>` command, which performs similar operations on per-atom
and local vectors.
The sort(x) and rsort(x) functions operate on a global vector of inputs
and return a global vector of the same length.
The sum() function calculates the sum of all the vector elements. The
min() and max() functions find the minimum and maximum element
@ -953,6 +957,12 @@ of points, equally spaced by 1 in their x coordinate: (1,V1), (2,V2),
length N. The returned value is the slope of the line. If the line
has a single point or is vertical, it returns 1.0e20.
.. versionadded:: TBD
The sort(x) and rsort(x) functions sort the data of the input vector by
their numeric value: sort(x) sorts in ascending order, rsort(x) sorts
in descending order.
The gmask(x) function takes 1 argument which is a group ID. It
can only be used in atom-style variables. It returns a 1 for
atoms that are in the group, and a 0 for atoms that are not.

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -3252,6 +3252,7 @@ rRESPA
Rsi
Rso
Rspace
rsort
rsq
rst
rstyle

8
src/info.cpp
View File

@ -110,12 +110,6 @@ static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES
using namespace LAMMPS_NS;
// must match enumerator in variable.h
static const char *varstyles[] = {
"index", "loop", "world", "universe", "uloop", "string", "getenv",
"file", "atomfile", "format", "equal", "atom", "vector", "python",
"timer", "internal", "(unknown)"};
static const char *mapstyles[] = { "none", "array", "hash", "yes" };
static const char *commstyles[] = { "brick", "tiled" };
@ -1401,7 +1395,7 @@ std::string Info::get_variable_info(int num) {
std::string text;
int ndata = 1;
text = fmt::format("Variable[{:3d}]: {:16} style = {:16} def =", num,
std::string(names[num]) + ',', std::string(varstyles[style[num]]) + ',');
std::string(names[num]) + ',', Variable::varstyles[style[num]] + ',');
if (style[num] == Variable::INTERNAL) {
text += fmt::format("{:.8}\n",input->variable->dvalue[num]);
return text;

129
src/variable.cpp
View File

@ -40,9 +40,11 @@
#include "fmt/ranges.h"
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstring>
#include <functional>
#include <unordered_map>
using namespace LAMMPS_NS;
@ -54,6 +56,12 @@ static constexpr int MAXLINE = 256;
static constexpr int CHUNK = 1024;
static constexpr int MAXFUNCARG = 6;
// must match enumerator in variable.h
const std::vector<std::string> Variable::varstyles = {
"index", "loop", "world", "universe", "uloop", "string", "getenv",
"file", "atomfile", "format", "equal", "atom", "vector", "python",
"timer", "internal", "(unknown)"};
static inline double MYROUND(double a) { return ((a - floor(a)) >= 0.5) ? ceil(a) : floor(a); }
enum{ARG,OP};
@ -73,7 +81,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
// customize by adding a special function
enum{SUM,XMIN,XMAX,AVE,TRAP,SLOPE};
enum{SUM,XMIN,XMAX,AVE,TRAP,SLOPE,SORT,RSORT,NOVECTOR};
static constexpr double BIG = 1.0e20;
@ -2279,7 +2287,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
if (math_function(word,contents,tree,treestack,ntreestack,argstack,nargstack,ivar));
else if (group_function(word,contents,tree,treestack,ntreestack,argstack,nargstack,ivar));
else if (special_function(word,contents,tree,treestack,ntreestack,argstack,nargstack,ivar));
else if (special_function(std::string(word),contents,tree,treestack,ntreestack,argstack,nargstack,ivar,str,i,ptr));
else if (feature_function(word,contents,tree,treestack,ntreestack,argstack,nargstack,ivar));
else print_var_error(FLERR,fmt::format("Invalid math/group/special/feature function '{}()' "
"in variable formula", word),ivar);
@ -4261,32 +4269,34 @@ Region *Variable::region_function(char *id, int ivar)
extract_setting(x),label2type(x,y),is_typelabel(x,y)
------------------------------------------------------------------------- */
int Variable::special_function(char *word, char *contents, Tree **tree, Tree **treestack,
int &ntreestack, double *argstack, int &nargstack, int ivar)
// to simplify finding matches and assigning constants for functions operating on vectors
static const std::unordered_map<std::string,int> special_function_map = {
{"sum", SUM}, {"min", XMIN}, {"max", XMAX}, {"ave", AVE}, {"trap", TRAP}, {"slope", SLOPE},
{"sort", SORT}, {"rsort", RSORT}, {"gmask", NOVECTOR}, {"rmask", NOVECTOR}, {"grmask", NOVECTOR},
{"next", NOVECTOR}, {"is_file", NOVECTOR}, {"is_os", NOVECTOR}, {"extract_setting", NOVECTOR},
{"label2type", NOVECTOR}, {"is_typelabel", NOVECTOR} };
int Variable::special_function(const std::string &word, char *contents, Tree **tree,
Tree **treestack, int &ntreestack, double *argstack,
int &nargstack, int ivar, char *str, int &istr, char *&ptr)
{
double sx,sxx;
double value,sy,sxy;
// word is not a match to any special function
if (strcmp(word,"sum") != 0 && strcmp(word,"min") && strcmp(word,"max") != 0 &&
strcmp(word,"ave") != 0 && strcmp(word,"trap") != 0 && strcmp(word,"slope") != 0 &&
strcmp(word,"gmask") != 0 && strcmp(word,"rmask") != 0 && strcmp(word,"grmask") != 0 &&
strcmp(word,"next") != 0 && strcmp(word,"is_file") != 0 && strcmp(word,"is_os") != 0 &&
strcmp(word,"extract_setting") != 0 && strcmp(word,"label2type") != 0 &&
strcmp(word,"is_typelabel") != 0)
return 0;
// return if "word" is not a match to any special function
if (special_function_map.find(word) == special_function_map.end()) return 0;
// process label2type() separately b/c its label arg can have commas in it
if (strcmp(word,"label2type") == 0 || strcmp(word,"is_typelabel") == 0) {
if ((word == "label2type") || (word == "is_typelabel")) {
if (!atom->labelmapflag)
print_var_error(FLERR,fmt::format("Cannot use {}() function without a labelmap",word),ivar);
std::string contents_copy(contents);
auto pos = contents_copy.find_first_of(',');
if (pos == std::string::npos) {
if (strcmp(word,"label2type") == 0) {
if (word == "label2type") {
print_var_error(FLERR, fmt::format("Invalid label2type({}) function in variable formula",
contents_copy), ivar);
} else {
@ -4313,7 +4323,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
print_var_error(FLERR, fmt::format("Invalid kind {} in {}() in variable", kind, word),ivar);
}
if (strcmp(word,"label2type") == 0) {
if (word == "label2type") {
if (value == -1)
print_var_error(FLERR, fmt::format("Invalid {} type label {} in label2type() in variable",
kind, typestr), ivar);
@ -4340,22 +4350,13 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
char *args[MAXFUNCARG];
int narg = parse_args(contents,args);
// special functions that operate on global vectors
// special functions that operate on global vectors are NOT mapped to NOVECTOR
if (strcmp(word,"sum") == 0 || strcmp(word,"min") == 0 ||
strcmp(word,"max") == 0 || strcmp(word,"ave") == 0 ||
strcmp(word,"trap") == 0 || strcmp(word,"slope") == 0) {
int method = 0;
if (strcmp(word,"sum") == 0) method = SUM;
else if (strcmp(word,"min") == 0) method = XMIN;
else if (strcmp(word,"max") == 0) method = XMAX;
else if (strcmp(word,"ave") == 0) method = AVE;
else if (strcmp(word,"trap") == 0) method = TRAP;
else if (strcmp(word,"slope") == 0) method = SLOPE;
int method = special_function_map.find(word)->second;
if (method != NOVECTOR) {
if (narg != 1)
print_var_error(FLERR,"Invalid special function in variable formula",ivar);
print_var_error(FLERR,fmt::format("Invalid special function {}() in variable formula", word),ivar);
Compute *compute = nullptr;
Fix *fix = nullptr;
@ -4473,6 +4474,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
if (method == SLOPE) sx = sxx = sy = sxy = 0.0;
else if (method == XMIN) value = BIG;
else if (method == XMAX) value = -BIG;
std::vector<double> unsorted;
if (compute) {
double *vec;
@ -4481,6 +4483,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
else vec = nullptr;
} else vec = compute->vector;
if ((method == SORT) || (method == RSORT)) unsorted.reserve(nvec);
int j = 0;
for (int i = 0; i < nvec; i++) {
if (method == SUM) value += vec[j];
@ -4488,6 +4491,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
else if (method == XMAX) value = MAX(value,vec[j]);
else if (method == AVE) value += vec[j];
else if (method == TRAP) value += vec[j];
else if ((method == SORT) || (method == RSORT)) unsorted.push_back(vec[j]);
else if (method == SLOPE) {
sx += (double)i;
sy += vec[j];
@ -4509,6 +4513,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
else if (method == XMAX) value = MAX(value,one);
else if (method == AVE) value += one;
else if (method == TRAP) value += one;
else if ((method == SORT) || (method == RSORT)) unsorted.push_back(one);
else if (method == SLOPE) {
sx += (double)i;
sy += one;
@ -4534,6 +4539,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
else if (method == XMAX) value = MAX(value,one);
else if (method == AVE) value += one;
else if (method == TRAP) value += one;
else if ((method == SORT) || (method == RSORT)) unsorted.push_back(one);
else if (method == SLOPE) {
sx += (double) i;
sy += one;
@ -4553,18 +4559,55 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
else value = BIG;
}
// save value in tree or on argstack
if ((method == SORT) || (method == RSORT)) {
if (method == SORT) std::sort(unsorted.begin(), unsorted.end(), std::less<double>());
if (method == RSORT) std::sort(unsorted.begin(), unsorted.end(), std::greater<double>());
if (tree) {
auto newtree = new Tree();
newtree->type = VALUE;
newtree->value = value;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value;
if (tree) {
double *newvec;
memory->create(newvec,nvec,"variable:values");
for (int m = 0; m < nvec; m++)
newvec[m] = unsorted[m];
auto newtree = new Tree();
newtree->type = VECTORARRAY;
newtree->array = newvec;
newtree->nvector = nvec;
newtree->nstride = 1;
newtree->selfalloc = 1;
treestack[ntreestack++] = newtree;
} else {
// process one pair of trailing brackets
// point istr beyond last bracket
// nbracket = # of bracket pairs
// index = int inside each bracket pair, vector index
if (str[istr] == '[') {
ptr = &str[istr];
int index = int_between_brackets(ptr,1);
if (index > nvec)
print_var_error(FLERR, fmt::format("index {} exceeds vector size of {}\n", index, nvec),ivar);
istr = ptr-str+1;
argstack[nargstack++] = unsorted[index-1];
}
}
} else {
// save value in tree or on argstack
if (tree) {
auto newtree = new Tree();
newtree->type = VALUE;
newtree->value = value;
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value;
}
// mask special functions
} else if (strcmp(word,"gmask") == 0) {
} else if (word == "gmask") {
if (tree == nullptr)
print_var_error(FLERR,"Gmask function in equal-style variable formula",ivar);
if (narg != 1)
@ -4579,7 +4622,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
newtree->ivalue = group->bitmask[igroup];
treestack[ntreestack++] = newtree;
} else if (strcmp(word,"rmask") == 0) {
} else if (word == "rmask") {
if (tree == nullptr)
print_var_error(FLERR,"Rmask function in equal-style variable formula",ivar);
if (narg != 1)
@ -4593,7 +4636,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
newtree->region = region;
treestack[ntreestack++] = newtree;
} else if (strcmp(word,"grmask") == 0) {
} else if (word == "grmask") {
if (tree == nullptr)
print_var_error(FLERR,"Grmask function in equal-style variable formula",ivar);
if (narg != 2)
@ -4613,7 +4656,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
// special function for file-style or atomfile-style variables
} else if (strcmp(word,"next") == 0) {
} else if (word == "next") {
if (narg != 1)
print_var_error(FLERR,"Invalid special function in variable formula",ivar);
@ -4664,7 +4707,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
} else print_var_error(FLERR,"Invalid variable style in special function next",ivar);
} else if (strcmp(word,"is_file") == 0) {
} else if (word == "is_file") {
if (narg != 1)
print_var_error(FLERR,"Invalid is_file() function in variable formula",ivar);
@ -4681,7 +4724,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value;
} else if (strcmp(word,"is_os") == 0) {
} else if (word == "is_os") {
if (narg != 1) print_var_error(FLERR,"Invalid is_os() function in variable formula",ivar);
value = utils::strmatch(platform::os_info(), args[0]) ? 1.0 : 0.0;
@ -4694,7 +4737,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value;
} else if (strcmp(word,"extract_setting") == 0) {
} else if (word == "extract_setting") {
if (narg != 1) print_var_error(FLERR,"Invalid extract_setting() function in variable formula",ivar);
value = lammps_extract_setting(lmp, args[0]);

9
src/variable.h
View File

@ -56,7 +56,7 @@ class Variable : protected Pointers {
int nvar; // # of defined variables
char **names; // name of each variable
// must match "varstyles" array in info.cpp
// must match "varstyles" array in variables.cpp, UNKNOWN must be last.
enum {
INDEX,
LOOP,
@ -73,9 +73,11 @@ class Variable : protected Pointers {
VECTOR,
PYTHON,
TIMER,
INTERNAL
INTERNAL,
UNKNOWN
};
static constexpr int VALUELENGTH = 64;
static const std::vector<std::string> varstyles;
private:
int me;
@ -141,7 +143,8 @@ class Variable : protected Pointers {
int math_function(char *, char *, Tree **, Tree **, int &, double *, int &, int);
int group_function(char *, char *, Tree **, Tree **, int &, double *, int &, int);
Region *region_function(char *, int);
int special_function(char *, char *, Tree **, Tree **, int &, double *, int &, int);
int special_function(const std::string &, char *, Tree **, Tree **, int &, double *, int &,
int, char *, int &, char *&);
int feature_function(char *, char *, Tree **, Tree **, int &, double *, int &, int);
void peratom2global(int, char *, double *, int, tagint, Tree **, Tree **, int &, double *, int &);
void custom2global(int *, double *, int, tagint, Tree **, Tree **, int &, double *, int &);

15
unittest/commands/test_variables.cpp
View File

@ -337,6 +337,13 @@ TEST_F(VariableTest, Expressions)
command("variable vec1 vector \"[-2, 0, 1,2 ,3, 5 , 7\n]\"");
command("variable vec2 vector v_vec1*0.5");
command("variable vec3 equal v_vec2[3]");
command("variable vec4 vector '[1, 5, 2.5, -10, -5, 20, 120, 4, 3, 3]'");
command("variable sort vector sort(v_vec4)");
command("variable rsrt vector rsort(v_vec4)");
command("variable max2 equal sort(v_vec4)[2]");
command("variable rmax equal rsort(v_vec4)[1]");
command("variable xxxl equal rsort(v_vec4)[11]");
command("variable isrt vector sort(v_one)");
variable->set("dummy index 1 2");
END_HIDE_OUTPUT();
@ -366,6 +373,10 @@ TEST_F(VariableTest, Expressions)
EXPECT_THAT(variable->retrieve("vec1"), StrEq("[-2,0,1,2,3,5,7]"));
EXPECT_THAT(variable->retrieve("vec2"), StrEq("[-1,0,0.5,1,1.5,2.5,3.5]"));
ASSERT_DOUBLE_EQ(variable->compute_equal("v_vec3"), 0.5);
EXPECT_THAT(variable->retrieve("sort"), StrEq("[-10,-5,1,2.5,3,3,4,5,20,120]"));
EXPECT_THAT(variable->retrieve("rsrt"), StrEq("[120,20,5,4,3,3,2.5,1,-5,-10]"));
ASSERT_DOUBLE_EQ(variable->compute_equal("v_max2"), -5);
ASSERT_DOUBLE_EQ(variable->compute_equal("v_rmax"), 120);
TEST_FAILURE(".*ERROR: Variable six: Invalid thermo keyword 'XXX' in variable formula.*",
command("print \"${six}\""););
@ -377,6 +388,10 @@ TEST_F(VariableTest, Expressions)
command("print \"${err2}\""););
TEST_FAILURE(".*ERROR on proc 0: Variable err3: Invalid power expression in variable formula.*",
command("print \"${err3}\""););
TEST_FAILURE(".*ERROR: Variable one: Mis-matched special function variable in variable formula.*",
command("print \"${isrt}\""););
TEST_FAILURE(".*ERROR: Variable vec4: index 11 exceeds vector size of 10.*",
command("print \"${xxxl}\""););
}
TEST_F(VariableTest, Functions)