Skip to content

GitLab

Commit be7bd3ec authored by Sergey Kireev's avatar Sergey Kireev
Browse files

abs/des implemented in cpp

parent 4bbd2143
Hide whitespace changes
Inline Side-by-side
nle_cpp/src/main.cpp
View file @ be7bd3ec
......@@ -30,6 +30,9 @@ void count(const Parameters &p) {
std::shared_ptr<CA> ca(CA::create(nx,ny,p.get_nsize()));
ca->set_parameters(p.get_en(),p.get_el(),p.get_enl(),p.get_kT(),p.get_mu(),p.get_evap(),p.get_cond(),p.get_pgroup());
ca->set_parameter_abs(p.get_abs());
ca->set_parameter_des(p.get_des());
ca->set_parameter_minndes(p.get_minndes());
ca->set_test_move_NL_dH(p.is_check());
if (p.is_load_from_file()) {
......
nle_cpp/src/nle_ca.h
View file @ be7bd3ec
......@@ -23,7 +23,8 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
const int nsize; // size of N side
int nnsize; // number of cells in N
int ngroups; // number of groups, set by find_groups()
double en,el,enl,kT,mu,evap,cond,pgroup;
int minndes; // minimum number of neighbors to prevent desorb N
double en,el,enl,kT,mu,evap,cond,absn,desn,pgroup;
std::vector<int> sN,xN,yN;
std::vector<int> oppdir;
......@@ -48,7 +49,11 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
kT = 0.836;
mu = 0.0;
evap = 0.0;
cond = 0.0;
absn = 0.0;
desn = 0.0;
pgroup = 0.0;
minndes = 0;
test_move_NL_dH = false;
//printf("ndirs: %d nneighbors: %d\n",ndirs,nneighbors);
//printf("dirs.count %d hood.count: %d\n",ca::dirs.count,ca::hood.count);
......@@ -87,6 +92,10 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
void set_parameter_cond(double val) { cond = val; }
void set_parameter_pgroup(double val) { pgroup = val; }
void set_parameter_abs(double val) { absn = val; }
void set_parameter_des(double val) { desn = val; }
void set_parameter_minndes(int val) { minndes = val; }
int xy2s(int i,int j) const { return (i+1) + (j+1)*1024; }
int s2x(int s) const { return s%1024-1; }
int s2y(int s) const { return s/1024-1; }
......@@ -152,12 +161,15 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
for (int d=0;d<ca::hood.count;d++) {
const int state = ca::getn(x,y,d);
if (isN(state)) {
const int bx = s2x(state);
const int by = s2y(state);
if (!wasn || bx != pbx || by != pby) count_n++;
wasn = true;
pbx = bx;
pby = by;
if (nsize == 1) count_n++;
else {
const int bx = s2x(state);
const int by = s2y(state);
if (!wasn || bx != pbx || by != pby) count_n++;
wasn = true;
pbx = bx;
pby = by;
}
}
else {
wasn = false;
......@@ -172,6 +184,8 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
return { count_n, count_l, count_e };
}
virtual int count_N_neighbors(int bx,int by) { exit(-1); return 0; }
bool try_swap_LE(const int lx,const int ly,const int ex,const int ey) {
//printf("try_swap_LE(%d,%d, %d,%d) : %d\n",lx,ly,ex,ey,ca::hood.count); fflush(stdout);
const auto cl = neighborhood(lx,ly);
......@@ -199,6 +213,9 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
return ca::rnd01() < exp(-dH/kT);
}
bool may_absorb(int bx,int by) { return is_filled_L(bx,by); }
bool may_desorb(int bx,int by) { return count_N_neighbors(bx,by) < minndes; }
double count_dH_move_NL(int bx,int by,int dir) {
double dH = 0;
for (int i=0;i<st[dir].size();i++) {
......@@ -321,7 +338,7 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
virtual void rule(int i,int j) {
ca::dirs.fix(i,j,ca::nx,ca::ny);
int state = ca::get(i,j);
if (isE(state)) {
if (isE(state)) { // E
if (cond > 0.0 && ca::rnd01() < cond) { // condensation: E -> L
if (try_cond_L(i,j)) {
ca::get(i,j) = L;
......@@ -340,14 +357,21 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
// }
//}
}
else if (state == L) {
if (evap > 0.0 && ca::rnd01() < evap) { // evaporation: L -> E
else if (isL(state)) { // L
const double rndL = absn + evap > 0.0 ? ca::rnd01() : 0.0;
if (rndL < absn) { // absorption: L..L -> N
const int bx = i - xN[0];
const int by = j - yN[0];
if (may_absorb(bx,by))
put_N(bx,by);
}
else if (rndL < absn + evap) { // evaporation: L -> E
if (try_evap_L(i,j)) {
ca::get(i,j) = E;
}
}
else { // movement: L <-> E
int dir = ca::dirs.rnd();
const int dir = ca::dirs.rnd();
int newi = i;
int newj = j;
ca::dirs.move(newi,newj,dir); ca::dirs.fix(newi,newj,ca::nx,ca::ny);
......@@ -359,24 +383,32 @@ class NLE_CA_base : public CA2D<int,dirs_type,hood_type> {
}
}
}
else if (pgroup > 0 && ca::rnd01() < pgroup) { // group movement N <-> L
find_groups();
int dir = ca::dirs.rnd();
int label = ca::group[ca::xy1d(i,j)];
if (may_move_group_NL(label,dir)) {
int size = ca::groupsize[label];
double p = 1.0/size;
if (ca::rnd01() < p)
do_move_group_NL(label,dir);
else { // N
const int bx = i - s2x(state);
const int by = j - s2y(state);
const double rndN = desn + pgroup > 0.0 ? ca::rnd01() : 0.0;
if (rndN < desn) { // desorption: N -> L..L
if (may_desorb(bx,by))
for (int k=0;k<sN.size();k++)
ca::getfixd(bx+xN[k],by+yN[k]) = L;
}
else if (rndN < desn + pgroup) { // group movement N <-> L
find_groups();
const int dir = ca::dirs.rnd();
const int label = ca::group[ca::xy1d(i,j)];
if (may_move_group_NL(label,dir)) {
const int size = ca::groupsize[label];
const double p = 1.0/size;
if (ca::rnd01() < p)
do_move_group_NL(label,dir);
}
}
else { // movement: N <-> L
const int dir = ca::dirs.rnd();
if (may_move_NL(bx,by,dir))
if (try_move_NL(bx,by,dir))
do_move_NL(bx,by,dir);
}
}
else { // movement: N <-> L
int bx = i - s2x(state);
int by = j - s2y(state);
int dir = ca::dirs.rnd();
if (may_move_NL(bx,by,dir))
if (try_move_NL(bx,by,dir))
do_move_NL(bx,by,dir);
}
}
......@@ -805,6 +837,20 @@ class NLE_sqr4_CA : public NLE_sqr_CA<hood_vonNeumann> {
return { count_n, count_l, count_e };
}
virtual int count_N_neighbors(int bx,int by) {
if (nsize == 1) return neighborhood(bx,by)[0];
const int sLT = xy2s(0,0);
const int sLB = xy2s(0,nsize-1);
const int sRT = xy2s(nsize-1,0);
const int sRB = xy2s(nsize-1,nsize-1);
//int[] test_count_nle_on_line(int x,int y,int dx,int dy,int n,int initial_state,int final_state)
const auto c0 = test_count_nle_on_line(bx,by-1,1,0,nsize,sLB,sRB);
const auto c1 = test_count_nle_on_line(bx+nsize,by,0,1,nsize,sLT,sLB);
const auto c2 = test_count_nle_on_line(bx,by+nsize,1,0,nsize,sLT,sRT);
const auto c3 = test_count_nle_on_line(bx-1,by,0,1,nsize,sRT,sRB);
return c0[0] + c1[0] + c2[0] + c3[0];
}
std::array<int,3> test_count_nle_on_line_l(int x,int y,int dx,int dy,int n,int initial_state,int final_state) {
int count_n = 0;
int count_l = 0;
......@@ -1147,6 +1193,25 @@ class NLE_sqr8_CA : public NLE_sqr_CA<hood_Moore> {
return { snn, snl, sll };
}
virtual int count_N_neighbors(int bx,int by) {
if (nsize == 1) return neighborhood(bx,by)[0];
const int sLT = xy2s(0,0);
const int sLB = xy2s(0,nsize-1);
const int sRT = xy2s(nsize-1,0);
const int sRB = xy2s(nsize-1,nsize-1);
//int[] test_count_nle_on_line(int x,int y,int dx,int dy,int n,int initial_state,int final_state)
const auto c0 = test_count_nle_on_line(bx-1,by-1,1,0,nsize+2,sLB,sRB);
const auto c1 = test_count_nle_on_line(bx+nsize,by-1,0,1,nsize+2,sLT,sLB);
const auto c2 = test_count_nle_on_line(bx-1,by+nsize,1,0,nsize+2,sLT,sRT);
const auto c3 = test_count_nle_on_line(bx-1,by-1,0,1,nsize+2,sRT,sRB);
int snn = c0[0] + c1[0] + c2[0] + c3[0];
{ const int s = ca::getfixd(bx-1 ,by-1 ); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx+nsize,by-1 ); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx-1 ,by+nsize); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx+nsize,by+nsize); if (isN(s)) snn--; }
return snn;
}
public:
static NLE_sqr8_CA* create(int nx,int ny,int ns = 1) {
......@@ -1550,6 +1615,32 @@ class NLE_hex_CA : public NLE_CA_base<hood_hex> {
return { snn, snl, sll };
}
virtual int count_N_neighbors(int bx,int by) {
if (nsize == 1) return neighborhood(bx,by)[0];
const int ns = nsize;
const int sTL = xy2s(ns-1,0);
const int sTR = xy2s(2*ns-2,0);
const int sRR = xy2s(2*ns-2,ns-1);
const int sBR = xy2s(ns-1,2*ns-2);
const int sBL = xy2s(0,2*ns-2);
const int sLL = xy2s(0,ns-1);
//int[] test_count_nle_on_line(int x,int y,int dx,int dy,int n,int initial_state,int final_state)
const auto c0 = test_count_nle_on_line(bx+ ns-1,by -1, 1, 0,nsize+1,sBL,sBR);
const auto c1 = test_count_nle_on_line(bx+2*ns-1,by -1, 0, 1,nsize+1,sLL,sBL);
const auto c2 = test_count_nle_on_line(bx+2*ns-1,by+ ns-1,-1, 1,nsize+1,sTL,sLL);
const auto c3 = test_count_nle_on_line(bx+ ns-1,by+2*ns-1,-1, 0,nsize+1,sTR,sTL);
const auto c4 = test_count_nle_on_line(bx -1,by+2*ns-1, 0,-1,nsize+1,sRR,sTR);
const auto c5 = test_count_nle_on_line(bx+ -1,by+ ns-1, 1,-1,nsize+1,sBR,sRR);
int snn = c0[0] + c1[0] + c2[0] + c3[0] + c4[0] + c5[0];
{ const int s = ca::getfixd(bx+ ns-1,by -1); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx+2*ns-1,by -1); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx+2*ns-1,by+ ns-1); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx+ ns-1,by+2*ns-1); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx -1,by+2*ns-1); if (isN(s)) snn--; }
{ const int s = ca::getfixd(bx -1,by+ ns-1); if (isN(s)) snn--; }
return snn;
}
virtual void do_move_NL(int bx,int by,int dir) {
XY nb(bx,by);
ca::dirs.move(nb,dir);
......@@ -1685,7 +1776,6 @@ class NLE_hex_CA : public NLE_CA_base<hood_hex> {
}
}
void print_neighbourhood(int bx,int by,int dir) {
int x0 = -1, x1 = 2*nsize;
int y0 = -1, y1 = 2*nsize;
......
nle_cpp/src/param.cpp
View file @ be7bd3ec
......@@ -73,6 +73,18 @@
#define COND 0.0
#endif
#ifndef ABSORB
#define ABSORB 0.0
#endif
#ifndef DESORB
#define DESORB 0.0
#endif
#ifndef MINNDES
#define MINNDES 0
#endif
#ifndef PGROUP
#define PGROUP 0.0
#endif
......@@ -142,6 +154,9 @@ Parameters::Parameters():
default_mu (MU),
default_evaporation(EVAP),
default_condensation(COND),
default_absorption(ABSORB),
default_desorption(DESORB),
default_minndes(MINNDES),
default_pgroup(PGROUP),
default_long_output(LONG),
default_datafile_prefix(PREFIX),
......@@ -167,6 +182,7 @@ void Parameters::reset() {
nsize = default_nsize;
long_output = default_long_output;
initial_step = default_initial_step;
minndes = default_minndes;
fraction_n = default_fraction_n;
fraction_l = default_fraction_l;
......@@ -178,6 +194,8 @@ void Parameters::reset() {
mu = default_mu;
evaporation = default_evaporation;
condensation = default_condensation;
absorption = default_absorption;
desorption = default_desorption;
pgroup = default_pgroup;
datafile_prefix = default_datafile_prefix;
......@@ -236,9 +254,12 @@ bool Parameters::read_from_cmdline(int argc,char *argv[]) {
std::cout << " -enl <value> - set value of parameter 'enl' (default: " << default_enl << ")" << std::endl;
std::cout << " -kt <value> - set value of parameter 'kT' (default: " << default_kT << ")" << std::endl;
std::cout << " -mu <value> - set value of parameter 'mu' (default: " << default_mu << ")" << std::endl;
std::cout << " -evap <value> - set evaporation probability: 'L' -> 'E' (default: " << default_evaporation << ")" << std::endl;
std::cout << " -cond <value> - set condensation probability: 'E' -> 'L' (default: " << default_evaporation << ")" << std::endl;
std::cout << " -evap <value> - set L evaporation probability: 'L' -> 'E' (default: " << default_evaporation << ")" << std::endl;
std::cout << " -cond <value> - set L condensation probability: 'E' -> 'L' (default: " << default_evaporation << ")" << std::endl;
std::cout << " -ec <value> - set evaporation and condensation probability" << std::endl;
std::cout << " -absorb <value> - set N absorption probability: 'L' -> 'N' (default: " << default_absorption << ")" << std::endl;
std::cout << " -desorb <value> - set N desorption probability: 'N' -> 'L' (default: " << default_desorption << ")" << std::endl;
std::cout << " -minndes <count> - set minimal number of neighboring N to enable desorption (default: " << default_minndes << ")" << std::endl;
std::cout << " -g <value> - set probability of group movement (default: " << default_pgroup << ")" << std::endl;
std::cout << " -ofp <string> - set output file prefix (default: '" << default_datafile_prefix << "')" << std::endl;
std::cout << " -ofs <string> - set output file suffix (default: '" << default_datafile_suffix << "')" << std::endl;
......@@ -377,6 +398,24 @@ bool Parameters::read_from_cmdline(int argc,char *argv[]) {
evaporation = condensation = std::stod(argv[cnt]);
cnt++;
}
else if (strcmp("-abs",argv[cnt])==0 || strcmp("-absorb",argv[cnt])==0) {
cnt++;
if (cnt == argc) { std::cout << "Error: Option '-absorb' requires floating point parameter\n"; return false; }
absorption = std::stod(argv[cnt]);
cnt++;
}
else if (strcmp("-des",argv[cnt])==0 || strcmp("-desorb",argv[cnt])==0) {
cnt++;
if (cnt == argc) { std::cout << "Error: Option '-desorb' requires floating point parameter\n"; return false; }
desorption = std::stod(argv[cnt]);
cnt++;
}
else if (strcmp("-minndes",argv[cnt])==0) {
cnt++;
if (cnt == argc) { std::cout << "Error: Option '-minndes' requires integer parameter\n"; return false; }
minndes = std::stoi(argv[cnt]);
cnt++;
}
else if (strcmp("-ell",argv[cnt])==0 || strcmp("-el",argv[cnt])==0) {
cnt++;
if (cnt == argc) { std::cout << "Error: Option '-ell' requires floating point parameter\n"; return false; }
......@@ -537,6 +576,8 @@ bool Parameters::verify() const {
}
if (evaporation<0.0 || evaporation>1.0) { std::cout << "Error: wrong 'evaporation' parameter value: " << evaporation << std::endl; result = false; }
if (condensation<0.0 || condensation>1.0) { std::cout << "Error: wrong 'condensation' parameter value: " << condensation << std::endl; result = false; }
if (absorption<0.0 || absorption>1.0) { std::cout << "Error: wrong 'absorption' parameter value: " << absorption << std::endl; result = false; }
if (desorption<0.0 || desorption>1.0) { std::cout << "Error: wrong 'desorption' parameter value: " << desorption << std::endl; result = false; }
if (pgroup<0.0 || pgroup>1.0) { std::cout << "Error: wrong group movement probability value: " << pgroup << std::endl; result = false; }
......@@ -558,8 +599,14 @@ void Parameters::print() const {
else std::cout << " ";
std::cout << "En: " << enn << " Enl: " << enl << " El: " << ell
<< " kT: " << kT << " mu: " << mu;
if (evaporation == condensation) std::cout << " evap/cond: " << evaporation;
else std::cout << " evap: " << evaporation << " cond: " << condensation;
if (evaporation != 0.0 || condensation != 0.0)
if (evaporation == condensation) std::cout << " evap/cond: " << evaporation;
else std::cout << " evap: " << evaporation << " cond: " << condensation;
if (absorption > 0.0 || desorption > 0.0) {
if (absorption == desorption) std::cout << " abs/des: " << absorption;
else std::cout << " abs: " << absorption << " des: " << desorption;
if (desorption > 0) std::cout << " (min:" << minndes << ")";
}
if (pgroup > 0) std::cout << " G: " << pgroup;
if (long_output) { std::cout << std::endl;
if (load_file) std::cout << "Loading file: " << input_filename << std::endl;
......
nle_cpp/src/param.h
View file @ be7bd3ec
......@@ -15,6 +15,7 @@ class Parameters {
int default_check;
int default_fill_debug;
int default_gstmin;
int default_minndes;
double default_fraction_n;
double default_fraction_l;
......@@ -27,6 +28,8 @@ class Parameters {
double default_evaporation;
double default_condensation;
double default_pgroup;
double default_absorption;
double default_desorption;
std::string default_datafile_prefix;
std::string default_datafile_suffix;
......@@ -43,6 +46,8 @@ class Parameters {
double evaporation;
double condensation;
double pgroup;
double absorption;
double desorption;
int nx,ny;
int nt,ns,np;
......@@ -52,6 +57,7 @@ class Parameters {
int squhex;
int save_groups;
int gstmin;
int minndes;
std::string datafile_prefix;
std::string datafile_suffix;
......@@ -89,6 +95,8 @@ public:
double get_mu() const { return mu; }
double get_evap() const { return evaporation; }
double get_cond() const { return condensation; }
double get_abs() const { return absorption; }
double get_des() const { return desorption; }
double get_pgroup() const { return pgroup; }
int get_nx() const { return nx; }
......@@ -98,6 +106,7 @@ public:
int get_np() const { return np; }
int get_nsize() const { return nsize; }
int get_initial_step() const { return initial_step; }
int get_minndes() const { return minndes; }
std::string get_prefix() const { return datafile_prefix; }
std::string get_suffix() const { return datafile_suffix; }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment