Skip to content

seed(n) variable function #28

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
g0mb4 wants to merge 9 commits into
base: master
Choose a base branch
from
Open

seed(n) variable function #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
005fb3f
variable: parse randomseed() function.
g0mb4 Mar 24, 2023
639ad73
variable: implement randomseed().
g0mb4 Mar 24, 2023
4fd218f
variable: clean up randomseed().
g0mb4 Mar 24, 2023
c5eabe0
Variable: move srand() to LAMMPS initialization.
g0mb4 Mar 27, 2023
f9e9252
doc: Add randomseed().
g0mb4 Mar 27, 2023
fa19e6e
variable: fix range in random_seed().
g0mb4 Mar 27, 2023
86e3f95
Variable: fix typo.
g0mb4 Mar 27, 2023
c920d5c
variable: rework seed(n).
g0mb4 May 3, 2023
32b6580
Merge branch 'CFDEMproject:master' into random-seed
g0mb4 Jul 11, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion doc/variable.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,8 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
math functions = sqrt(x), exp(x), ln(x), log(x), abs(x),
sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x),
random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x)
ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z),
seed(n)
group functions = count(group), mass(group), charge(group),
xcm(group,dim), vcm(group,dim), fcm(group,dim),
bound(group,xmin), gyration(group), ke(group),
Expand Down Expand Up @@ -553,6 +554,12 @@ the {start} keyword of the "run"_run.html command. See the
"thermo_style"_thermo_style.html keyword elaplong =
timestep-startstep.

The seed(n) function generates a prime number between 10000 and
INT_MAX (2147483647). The output of this function can be used to be the seed of
particletemplate, particledistribution, etc.
If the argument is 0, a random prime number is generated, else the n-th greater,
than 10000 prime (n = 1 -> 10007, n = 2 -> 10009, etc) is the output.

:line

Group and Region Functions :h4
Expand Down
11 changes: 11 additions & 0 deletions src/lammps.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,8 @@
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.

Tóth János (MATE, Gödöllő)
------------------------------------------------------------------------- */

#include <mpi.h>
Expand Down Expand Up @@ -86,6 +88,9 @@
#include "error.h"
#include "granular_styles.h"

#include <stdlib.h>
#include <time.h>

using namespace LAMMPS_NS;

/* ----------------------------------------------------------------------
Expand Down Expand Up @@ -559,6 +564,12 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
input->one(cmd);
error->done();
}

// used in Variable::random_seed()
if(me == 0){
// NOTE: better random generator for random_seed()?
srand(time(NULL));
Copy link
Contributor

@Azrael3000 Azrael3000 Apr 24, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will cause random numbers to be truly random, which is something you don't want as it destroys repeatability.

Copy link
Contributor Author

@g0mb4 g0mb4 Apr 24, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is my point. If I want to repeat my simulation, I can use a predefined seed, but if I want to repeat an "experiment" with "random" particle positions, I can use randomseed().

Copy link
Contributor

@schrummy14 schrummy14 Apr 24, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could create a list of, say 500, prime numbers given an initial seed. If someone provides 0, you get the current behavior, and for any other number you get a list of semi-random numbers that is repeatable. Then every call to randomseed simply increments the location in the pre-compiled list.
You would also be able to add additional checks that make sure each prime number is unique.

}
}

/* ----------------------------------------------------------------------
Expand Down
155 changes: 151 additions & 4 deletions src/variable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.

Tóth János (MATE, Gödöllő)
------------------------------------------------------------------------- */

#include <cmath>
Expand Down Expand Up @@ -69,6 +71,7 @@
#include "memory.h"
#include "error.h"
#include "force.h"
#include "math_extra_liggghts.h"

#if defined(_WIN32) || defined(_WIN64)
#include <windows.h>
Expand Down Expand Up @@ -99,7 +102,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE,
VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,SEED};

// customize by adding a special function

Expand Down Expand Up @@ -1673,7 +1676,8 @@ double Variable::evaluate(char *str, Tree **tree)
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z),
vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z),
gmask(x),rmask(x),grmask(x,y)
gmask(x),rmask(x),grmask(x,y),
seed(n)
---------------------------------------------------------------------- */

double Variable::collapse_tree(Tree *tree)
Expand Down Expand Up @@ -2107,6 +2111,17 @@ double Variable::collapse_tree(Tree *tree)
return tree->value;
}

if (tree->type == SEED) {
int ivalue1 = static_cast<int> (collapse_tree(tree->left));

if (ivalue1 < 0)
error->one(FLERR,"seed() must be called with a non-negative argument (>= 0)");

tree->type = VALUE;
tree->value = generate_seed(ivalue1);
return tree->value;
}

// mask functions do not become a single collapsed value

if (tree->type == GMASK) return 0.0;
Expand All @@ -2124,7 +2139,8 @@ double Variable::collapse_tree(Tree *tree)
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z),
vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z),
gmask(x),rmask(x),grmask(x,y)
gmask(x),rmask(x),grmask(x,y),
seed(n)
---------------------------------------------------------------------- */

double Variable::eval_tree(Tree *tree, int i)
Expand Down Expand Up @@ -2368,6 +2384,15 @@ double Variable::eval_tree(Tree *tree, int i)
return arg;
}

if (tree->type == SEED){
int ivalue1 = static_cast<int> (eval_tree(tree->left,i));

if (ivalue1 < 0)
error->one(FLERR,"seed() must be called with a non-negative argument (>= 0)");

return generate_seed(ivalue1);
}

if (tree->type == GMASK) {
if (atom->mask[i] & tree->ivalue1) return 1.0;
else return 0.0;
Expand Down Expand Up @@ -2496,7 +2521,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
strcmp(word,"ramp") && strcmp(word,"stagger") &&
strcmp(word,"logfreq") && strcmp(word,"stride") &&
strcmp(word,"vdisplace") &&
strcmp(word,"swiggle") && strcmp(word,"cwiggle"))
strcmp(word,"swiggle") && strcmp(word,"cwiggle") &&
strcmp(word,"seed"))
return 0;

// parse contents for arg1,arg2,arg3 separated by commas
Expand Down Expand Up @@ -2819,6 +2845,19 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
argstack[nargstack++] = value;
}

} else if (strcmp(word,"seed") == 0) {
if (narg != 1)
error->all(FLERR,"seed() must be called with 1 argument");
if (tree) newtree->type = SEED;
else {
int ivalue1 = static_cast<int> (value1);

if(ivalue1 < 0)
error->all(FLERR,"seed() must be called with a non-negative argument (>= 0)");

argstack[nargstack++] = generate_seed(ivalue1);
}
}

delete [] arg1;
Expand Down Expand Up @@ -3931,6 +3970,114 @@ unsigned int Variable::data_mask(char *str)
return datamask;
}

/* ----------------------------------------------------------------------
generates a random prime number between 10'000 and INT_MAX
in a form of 6k +- 1
------------------------------------------------------------------------- */
static int generate_random_seed(){
static const int MIN_K = 1678; // 6*1678-1 = 10067, prime
static const int MAX_K = 357913904; // 6*357913904-1 = 2147483423, prime

static const int RANGE_K = MAX_K - MIN_K + 1;

int p = 0;

// it will find a prime
while(true){
int k = (rand() % RANGE_K) + MIN_K;
p = 6 * k - 1;

if(MathExtraLiggghts::isPrime(p)){
break;
}

p += 2;

if(MathExtraLiggghts::isPrime(p)){
break;
}
}

return p;
}

/*
first 200 prime greater than 10'000:
*/
static const int PRIME_TABLE_COUNT = 200;
static const int PRIME_TABLE[PRIME_TABLE_COUNT] = {
10007, 10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099,
10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, 10181,
10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, 10273, 10289,
10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, 10357, 10369, 10391,
10399, 10427, 10429, 10433, 10453, 10457, 10459, 10463, 10477, 10487, 10499,
10501, 10513, 10529, 10531, 10559, 10567, 10589, 10597, 10601, 10607, 10613,
10627, 10631, 10639, 10651, 10657, 10663, 10667, 10687, 10691, 10709, 10711,
10723, 10729, 10733, 10739, 10753, 10771, 10781, 10789, 10799, 10831, 10837,
10847, 10853, 10859, 10861, 10867, 10883, 10889, 10891, 10903, 10909, 10937,
10939, 10949, 10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057,
11059, 11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149,
11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, 11257,
11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, 11351, 11353,
11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, 11447, 11467, 11471,
11483, 11489, 11491, 11497, 11503, 11519, 11527, 11549, 11551, 11579, 11587,
11593, 11597, 11617, 11621, 11633, 11657, 11677, 11681, 11689, 11699, 11701,
11717, 11719, 11731, 11743, 11777, 11779, 11783, 11789, 11801, 11807, 11813,
11821, 11827, 11831, 11833, 11839, 11863, 11867, 11887, 11897, 11903, 11909,
11923, 11927
};

static int generate_nth_prime_impl(int n){
int p = PRIME_TABLE[PRIME_TABLE_COUNT - 1];
n -= PRIME_TABLE_COUNT;

while(n > 0){
p += 2;

if(MathExtraLiggghts::isPrime(p)){
n -= 1;
}
}

return p;
}

/* ----------------------------------------------------------------------
gives the n-th prime number between 10'000 and INT_MAX

n p
1 10007
2 10009
3 10037
...
------------------------------------------------------------------------- */
static int generate_nth_prime(int n){
// fast
if(n <= PRIME_TABLE_COUNT){
return PRIME_TABLE[n-1];
// small
} else {
return generate_nth_prime_impl(n);
}
}

int Variable::generate_seed(int n){

int seed = 0;

if(me == 0){
if (n == 0){
seed = generate_random_seed();
} else {
seed = generate_nth_prime(n);
}
}

MPI_Bcast(&seed,1,MPI_INT,0,world);

return seed;
}

/* ----------------------------------------------------------------------
class to read variable values from a file
for flag = SCALARFILE, reads one value per line
Expand Down
4 changes: 4 additions & 0 deletions src/variable.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.

Tóth János (MATE, Gödöllő)
------------------------------------------------------------------------- */

#ifndef LMP_VARIABLE_H
Expand Down Expand Up @@ -128,6 +130,8 @@ class Variable : protected Pointers {
int inumeric(char *);
char *find_next_comma(char *);
void print_tree(Tree *, int);

int generate_seed(int);
};

class VarReader : protected Pointers {
Expand Down