re-use existing potential file, fix screen/log output, plug memory leaks
This commit is contained in:
Axel Kohlmeyer
File diff suppressed because it is too large
Load Diff
1
examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs
Symbolic link
1
examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs
Symbolic link
@ -0,0 +1 @@
|
||||
../../../potentials/Zr_mm.eam.fs
|
||||
@ -1,3 +0,0 @@
|
||||
pair_style hybrid/overlay zero 9.0 eam/fs
|
||||
pair_coeff * * zero
|
||||
pair_coeff * * eam/fs dir.inputs/Zr_2.eam.fs Zr
|
||||
@ -27,7 +27,7 @@ change_box all triclinic
|
||||
|
||||
pair_style hybrid/overlay zero 9.0 eam/fs
|
||||
pair_coeff * * zero
|
||||
pair_coeff * * eam/fs Zr_2.eam.fs Zr
|
||||
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
|
||||
|
||||
timestep 0.002
|
||||
|
||||
|
||||
@ -43,7 +43,8 @@ Changing box ...
|
||||
|
||||
pair_style hybrid/overlay zero 9.0 eam/fs
|
||||
pair_coeff * * zero
|
||||
pair_coeff * * eam/fs Zr_2.eam.fs Zr
|
||||
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
|
||||
Reading eam/fs potential file Zr_mm.eam.fs with DATE: 2007-06-11
|
||||
|
||||
timestep 0.002
|
||||
|
||||
@ -60,6 +61,48 @@ compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1
|
||||
|
||||
# 2nd step : perform dimension reduction + logistic regression
|
||||
compute slcsa all slcsa/atom 8 4 dir.slcsa/mean_descriptor.dat dir.slcsa/lda_scalings.dat dir.slcsa/lr_decision.dat dir.slcsa/lr_bias.dat dir.slcsa/mahalanobis_file.dat c_bnnn[*]
|
||||
Files used:
|
||||
database mean descriptor: dir.slcsa/mean_descriptor.dat
|
||||
lda scalings : dir.slcsa/lda_scalings.dat
|
||||
lr decision : dir.slcsa/lr_decision.dat
|
||||
lr bias : dir.slcsa/lr_bias.dat
|
||||
maha stats : dir.slcsa/mahalanobis_file.dat
|
||||
For class 0 maha threshold = 5.054
|
||||
mean B:
|
||||
-23.8329
|
||||
4.6638
|
||||
3.9805
|
||||
icov:
|
||||
1.1377 0.1077 -0.0171
|
||||
0.1077 0.8846 -0.2577
|
||||
-0.0171 -0.2577 0.6783
|
||||
For class 1 maha threshold = 5.234
|
||||
mean B:
|
||||
-21.2853
|
||||
-6.1583
|
||||
1.7948
|
||||
icov:
|
||||
1.7124 0.0341 0.1966
|
||||
0.0341 0.6453 0.288
|
||||
0.1966 0.288 1.8991
|
||||
For class 2 maha threshold = 5.036
|
||||
mean B:
|
||||
-23.1593
|
||||
1.3059
|
||||
-5.7549
|
||||
icov:
|
||||
0.7496 -0.0806 -0.1101
|
||||
-0.0806 1.1178 0.1667
|
||||
-0.1101 0.1667 0.6711
|
||||
For class 3 maha threshold = 7.994
|
||||
mean B:
|
||||
68.1971
|
||||
0.1604
|
||||
-0.0067
|
||||
icov:
|
||||
0.9663 -0.1846 0.6622
|
||||
-0.1846 8.2371 0.9841
|
||||
0.6622 0.9841 5.9601
|
||||
|
||||
#dump d1 all custom ${freqdump} slcsa_demo.dump id x y z c_slcsa[*]
|
||||
|
||||
@ -101,37 +144,37 @@ Neighbor list info ...
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 31.9 | 31.9 | 31.9 Mbytes
|
||||
Step PotEng KinEng Temp c_max_slcsa[1] c_max_slcsa[2] c_max_slcsa[3] c_max_slcsa[4] c_max_slcsa[5] c_min_slcsa[1] c_min_slcsa[2] c_min_slcsa[3] c_min_slcsa[4] c_min_slcsa[5]
|
||||
0 -139689.93 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
|
||||
50 -138733.53 1165.0291 417.29048 8.2778907 18.221432 3.019507 113.65218 2 5.2735872 10.834132 0.042702874 108.2674 2
|
||||
100 -138821.33 1442.648 516.7281 8.3357676 17.470311 2.7073063 113.39999 2 6.145345 12.307441 0.034372881 108.61824 2
|
||||
150 -138607.9 1435.9634 514.33383 8.3817341 17.984657 2.892662 113.40287 2 5.3635301 10.064171 0.030283213 107.5829 2
|
||||
200 -138494.56 1560.9093 559.08699 8.5615388 18.429835 2.8130845 114.31517 2 5.4711815 11.113891 0.074490923 107.87418 2
|
||||
Loop time of 36.3263 on 1 procs for 200 steps with 21600 atoms
|
||||
0 -143297.23 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
|
||||
50 -142154.08 1007.7164 360.9442 8.8091564 19.23244 4.2093382 113.87959 2 5.0327148 9.6817454 0.02610585 106.71863 2
|
||||
100 -142365.33 1406.6559 503.83647 8.6272189 17.908949 2.9294666 113.75167 2 6.2058895 11.913521 0.033775944 108.66893 2
|
||||
150 -142188.18 1432.0075 512.91691 8.6441961 18.176321 2.9277374 114.27958 2 5.5899425 10.521867 0.014919473 108.14526 2
|
||||
200 -142000.4 1481.7247 530.72462 8.5895692 18.65646 3.1725758 114.55015 2 5.5955774 10.776385 0.061469343 108.35384 2
|
||||
Loop time of 36.3759 on 1 procs for 200 steps with 21600 atoms
|
||||
|
||||
Performance: 0.951 ns/day, 25.227 hours/ns, 5.506 timesteps/s, 118.922 katom-step/s
|
||||
Performance: 0.950 ns/day, 25.261 hours/ns, 5.498 timesteps/s, 118.760 katom-step/s
|
||||
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 9.012 | 9.012 | 9.012 | 0.0 | 24.81
|
||||
Neigh | 0.67881 | 0.67881 | 0.67881 | 0.0 | 1.87
|
||||
Comm | 0.045066 | 0.045066 | 0.045066 | 0.0 | 0.12
|
||||
Output | 26.422 | 26.422 | 26.422 | 0.0 | 72.74
|
||||
Modify | 0.14726 | 0.14726 | 0.14726 | 0.0 | 0.41
|
||||
Other | | 0.02121 | | | 0.06
|
||||
Pair | 9.0837 | 9.0837 | 9.0837 | 0.0 | 24.97
|
||||
Neigh | 0.52896 | 0.52896 | 0.52896 | 0.0 | 1.45
|
||||
Comm | 0.045416 | 0.045416 | 0.045416 | 0.0 | 0.12
|
||||
Output | 26.548 | 26.548 | 26.548 | 0.0 | 72.98
|
||||
Modify | 0.1493 | 0.1493 | 0.1493 | 0.0 | 0.41
|
||||
Other | | 0.02088 | | | 0.06
|
||||
|
||||
Nlocal: 21600 ave 21600 max 21600 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 36652 ave 36652 max 36652 min
|
||||
Nghost: 36674 ave 36674 max 36674 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 2.5961e+06 ave 2.5961e+06 max 2.5961e+06 min
|
||||
Neighs: 2.61729e+06 ave 2.61729e+06 max 2.61729e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 5.19109e+06 ave 5.19109e+06 max 5.19109e+06 min
|
||||
FullNghs: 5.24007e+06 ave 5.24007e+06 max 5.24007e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 5191091
|
||||
Ave neighs/atom = 240.32829
|
||||
Neighbor list builds = 5
|
||||
Total # of neighbors = 5240069
|
||||
Ave neighs/atom = 242.59579
|
||||
Neighbor list builds = 4
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:43
|
||||
|
||||
@ -43,7 +43,8 @@ Changing box ...
|
||||
|
||||
pair_style hybrid/overlay zero 9.0 eam/fs
|
||||
pair_coeff * * zero
|
||||
pair_coeff * * eam/fs Zr_2.eam.fs Zr
|
||||
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
|
||||
Reading eam/fs potential file Zr_mm.eam.fs with DATE: 2007-06-11
|
||||
|
||||
timestep 0.002
|
||||
|
||||
@ -60,6 +61,48 @@ compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1
|
||||
|
||||
# 2nd step : perform dimension reduction + logistic regression
|
||||
compute slcsa all slcsa/atom 8 4 dir.slcsa/mean_descriptor.dat dir.slcsa/lda_scalings.dat dir.slcsa/lr_decision.dat dir.slcsa/lr_bias.dat dir.slcsa/mahalanobis_file.dat c_bnnn[*]
|
||||
Files used:
|
||||
database mean descriptor: dir.slcsa/mean_descriptor.dat
|
||||
lda scalings : dir.slcsa/lda_scalings.dat
|
||||
lr decision : dir.slcsa/lr_decision.dat
|
||||
lr bias : dir.slcsa/lr_bias.dat
|
||||
maha stats : dir.slcsa/mahalanobis_file.dat
|
||||
For class 0 maha threshold = 5.054
|
||||
mean B:
|
||||
-23.8329
|
||||
4.6638
|
||||
3.9805
|
||||
icov:
|
||||
1.1377 0.1077 -0.0171
|
||||
0.1077 0.8846 -0.2577
|
||||
-0.0171 -0.2577 0.6783
|
||||
For class 1 maha threshold = 5.234
|
||||
mean B:
|
||||
-21.2853
|
||||
-6.1583
|
||||
1.7948
|
||||
icov:
|
||||
1.7124 0.0341 0.1966
|
||||
0.0341 0.6453 0.288
|
||||
0.1966 0.288 1.8991
|
||||
For class 2 maha threshold = 5.036
|
||||
mean B:
|
||||
-23.1593
|
||||
1.3059
|
||||
-5.7549
|
||||
icov:
|
||||
0.7496 -0.0806 -0.1101
|
||||
-0.0806 1.1178 0.1667
|
||||
-0.1101 0.1667 0.6711
|
||||
For class 3 maha threshold = 7.994
|
||||
mean B:
|
||||
68.1971
|
||||
0.1604
|
||||
-0.0067
|
||||
icov:
|
||||
0.9663 -0.1846 0.6622
|
||||
-0.1846 8.2371 0.9841
|
||||
0.6622 0.9841 5.9601
|
||||
|
||||
#dump d1 all custom ${freqdump} slcsa_demo.dump id x y z c_slcsa[*]
|
||||
|
||||
@ -101,37 +144,37 @@ Neighbor list info ...
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 11.26 | 11.45 | 11.64 Mbytes
|
||||
Step PotEng KinEng Temp c_max_slcsa[1] c_max_slcsa[2] c_max_slcsa[3] c_max_slcsa[4] c_max_slcsa[5] c_min_slcsa[1] c_min_slcsa[2] c_min_slcsa[3] c_min_slcsa[4] c_min_slcsa[5]
|
||||
0 -139689.93 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
|
||||
50 -138733.53 1165.0291 417.29048 8.2778907 18.221432 3.019507 113.65218 2 5.2735872 10.834132 0.042702874 108.2674 2
|
||||
100 -138821.33 1442.648 516.7281 8.3357676 17.470311 2.7073063 113.39999 2 6.145345 12.307441 0.034372881 108.61824 2
|
||||
150 -138607.9 1435.9634 514.33383 8.3817341 17.984657 2.892662 113.40287 2 5.3635301 10.064171 0.030283213 107.5829 2
|
||||
200 -138494.56 1560.9093 559.08699 8.5615388 18.429835 2.8130845 114.31517 2 5.4711815 11.113891 0.074490923 107.87418 2
|
||||
Loop time of 9.89929 on 4 procs for 200 steps with 21600 atoms
|
||||
0 -143297.23 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
|
||||
50 -142154.08 1007.7164 360.9442 8.8091564 19.23244 4.2093382 113.87959 2 5.0327148 9.6817454 0.02610585 106.71863 2
|
||||
100 -142365.33 1406.6559 503.83647 8.6272189 17.908949 2.9294666 113.75167 2 6.2058895 11.913521 0.033775944 108.66893 2
|
||||
150 -142188.18 1432.0075 512.91691 8.6441961 18.176321 2.9277374 114.27958 2 5.5899425 10.521867 0.014919473 108.14526 2
|
||||
200 -142000.4 1481.7247 530.72462 8.5895692 18.65646 3.1725758 114.55015 2 5.5955774 10.776385 0.061469343 108.35384 2
|
||||
Loop time of 9.81677 on 4 procs for 200 steps with 21600 atoms
|
||||
|
||||
Performance: 3.491 ns/day, 6.875 hours/ns, 20.203 timesteps/s, 436.395 katom-step/s
|
||||
99.8% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
Performance: 3.521 ns/day, 6.817 hours/ns, 20.373 timesteps/s, 440.063 katom-step/s
|
||||
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 2.61 | 2.618 | 2.622 | 0.3 | 26.45
|
||||
Neigh | 0.19338 | 0.19618 | 0.19967 | 0.5 | 1.98
|
||||
Comm | 0.051086 | 0.055818 | 0.065206 | 2.4 | 0.56
|
||||
Output | 6.9725 | 6.9725 | 6.9725 | 0.0 | 70.43
|
||||
Modify | 0.04635 | 0.046617 | 0.046825 | 0.1 | 0.47
|
||||
Other | | 0.01018 | | | 0.10
|
||||
Pair | 2.6508 | 2.6589 | 2.6698 | 0.5 | 27.09
|
||||
Neigh | 0.1516 | 0.15276 | 0.15406 | 0.3 | 1.56
|
||||
Comm | 0.047132 | 0.058969 | 0.066095 | 3.2 | 0.60
|
||||
Output | 6.8886 | 6.8886 | 6.8886 | 0.0 | 70.17
|
||||
Modify | 0.046437 | 0.04661 | 0.046825 | 0.1 | 0.47
|
||||
Other | | 0.01091 | | | 0.11
|
||||
|
||||
Nlocal: 5400 ave 5431 max 5376 min
|
||||
Histogram: 1 1 0 0 0 0 1 0 0 1
|
||||
Nghost: 12889.2 ave 12909 max 12847 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 1 2
|
||||
Neighs: 649026 ave 652559 max 642783 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 1 1
|
||||
FullNghs: 1.29777e+06 ave 1.30656e+06 max 1.29086e+06 min
|
||||
Histogram: 1 1 0 0 0 0 1 0 0 1
|
||||
Nlocal: 5400 ave 5416 max 5393 min
|
||||
Histogram: 2 1 0 0 0 0 0 0 0 1
|
||||
Nghost: 12902.8 ave 12911 max 12888 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 0 2
|
||||
Neighs: 654322 ave 655602 max 650912 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 3
|
||||
FullNghs: 1.31002e+06 ave 1.31507e+06 max 1.30683e+06 min
|
||||
Histogram: 1 1 0 1 0 0 0 0 0 1
|
||||
|
||||
Total # of neighbors = 5191084
|
||||
Ave neighs/atom = 240.32796
|
||||
Neighbor list builds = 5
|
||||
Total # of neighbors = 5240065
|
||||
Ave neighs/atom = 242.5956
|
||||
Neighbor list builds = 4
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:11
|
||||
|
||||
@ -38,20 +38,25 @@
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
static const char cite_compute_slcsa_atom_c[] =
|
||||
"compute slcsa/atom command: doi:10.1088/0965-0393/21/5/055020\n\n"
|
||||
"@Article{Lafourcade2023,\n"
|
||||
" author = {P. Lafourcade and J.-B. Maillet and C. Denoual and E. Duval and A. Allera and A. M. Goryaeva and M.-C. Marinica},\n"
|
||||
" title = {Robust crystal structure identification at extreme conditions using a density-independent spectral descriptor and supervised learning},\n"
|
||||
" journal = {Computational Materials Science},\n"
|
||||
" year = 2023,\n"
|
||||
" volume = XX,\n"
|
||||
" pages = {XXXXXX}\n"
|
||||
"}\n\n";
|
||||
"compute slcsa/atom command: doi:10.1088/0965-0393/21/5/055020\n\n"
|
||||
"@Article{Lafourcade2023,\n"
|
||||
" author = {P. Lafourcade and J.-B. Maillet and C. Denoual and E. Duval and A. Allera and A. "
|
||||
"M. Goryaeva and M.-C. Marinica},\n"
|
||||
" title = {Robust crystal structure identification at extreme conditions using a "
|
||||
"density-independent spectral descriptor and supervised learning},\n"
|
||||
" journal = {Computational Materials Science},\n"
|
||||
" year = 2023,\n"
|
||||
" volume = XX,\n"
|
||||
" pages = {XXXXXX}\n"
|
||||
"}\n\n";
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), list(nullptr), classification(nullptr)
|
||||
Compute(lmp, narg, arg), list(nullptr), lda_scalings(nullptr),
|
||||
database_mean_descriptor(nullptr), lr_bias(nullptr), lr_decision(nullptr), icov_list(nullptr),
|
||||
mean_projected_descriptors(nullptr), maha_thresholds(nullptr), full_descriptor(nullptr),
|
||||
projected_descriptor(nullptr), scores(nullptr), probas(nullptr), prodright(nullptr),
|
||||
dmaha(nullptr), classification(nullptr)
|
||||
{
|
||||
// command : compute c1 all slcsa/atom jmax nclasses parameters_file.dat
|
||||
// example : compute c1 all slcsa/atom 8 4 slcsa_parameters.dat
|
||||
@ -74,14 +79,16 @@ ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
// # LR bias vector
|
||||
// vector with 1 row x nclasses cols
|
||||
|
||||
if (narg != 11) error->all(FLERR, "Illegal compute slcsa/atom command; wrong numer of arguments");
|
||||
if (narg != 11) utils::missing_cmd_args(FLERR, "compute slcsa/atom", error);
|
||||
|
||||
int twojmax = utils::inumeric(FLERR, arg[3], false, lmp);
|
||||
if (twojmax < 0) error->all(FLERR, "Illegal compute slcsa/atom command: twojmax should be a non-negative integer");
|
||||
if (twojmax < 0)
|
||||
error->all(FLERR, "Illegal compute slcsa/atom command: twojmax must be a non-negative integer");
|
||||
ncomps = compute_ncomps(twojmax);
|
||||
|
||||
nclasses = utils::inumeric(FLERR, arg[4], false, lmp);
|
||||
if (nclasses < 2) error->all(FLERR, "Illegal compute slcsa/atom command: nclasses should be strictly greater than 1");
|
||||
if (nclasses < 2)
|
||||
error->all(FLERR, "Illegal compute slcsa/atom command: nclasses must be greater than 1");
|
||||
|
||||
database_mean_descriptor_file = arg[5];
|
||||
lda_scalings_file = arg[6];
|
||||
@ -89,11 +96,14 @@ ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
lr_bias_file = arg[8];
|
||||
maha_file = arg[9];
|
||||
|
||||
std::cout << " file database mean descriptor : " << database_mean_descriptor_file << std::endl;
|
||||
std::cout << " file lda scalings : " << lda_scalings_file << std::endl;
|
||||
std::cout << " file lr decision : " << lr_decision_file << std::endl;
|
||||
std::cout << " file lr bias : " << lr_bias_file << std::endl;
|
||||
std::cout << " file maha stats : " << maha_file << std::endl;
|
||||
if (comm->me == 0) {
|
||||
auto mesg = fmt::format(
|
||||
"Files used:\n {:24}: {}\n {:24}: {}\n {:24}: {}\n {:24}: {}\n {:24}: {}\n",
|
||||
"database mean descriptor", database_mean_descriptor_file, "lda scalings",
|
||||
lda_scalings_file, "lr decision", lr_decision_file, "lr bias", lr_bias_file, "maha stats",
|
||||
maha_file);
|
||||
utils::logmesg(lmp, mesg);
|
||||
}
|
||||
|
||||
int expand = 0;
|
||||
char **earg;
|
||||
@ -108,7 +118,8 @@ ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
val.which = argi.get_type();
|
||||
val.argindex = argi.get_index1();
|
||||
val.id = argi.get_name();
|
||||
if ((val.which == ArgInfo::FIX) || (val.which == ArgInfo::VARIABLE) || (val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
|
||||
if ((val.which == ArgInfo::FIX) || (val.which == ArgInfo::VARIABLE) ||
|
||||
(val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
|
||||
error->all(FLERR, "Invalid compute slcsa/atom argument: {}", arg[0]);
|
||||
|
||||
// if wildcard expansion occurred, free earg memory from exapnd_args()
|
||||
@ -119,123 +130,126 @@ ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
}
|
||||
|
||||
val.val.c = modify->get_compute_by_id(val.id);
|
||||
if (!val.val.c) error->all(FLERR,"Compute ID {} for fix slcsa/atom does not exist", val.id);
|
||||
if (!val.val.c) error->all(FLERR, "Compute ID {} for fix slcsa/atom does not exist", val.id);
|
||||
if (val.val.c->peratom_flag == 0)
|
||||
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate per-atom values", val.id);
|
||||
if (val.argindex == 0 && val.val.c->size_peratom_cols != 0)
|
||||
error->all(FLERR,"Compute slcsa/atom compute {} does not calculate a per-atom vector", val.id);
|
||||
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate a per-atom vector", val.id);
|
||||
if (val.argindex && val.val.c->size_peratom_cols == 0)
|
||||
error->all(FLERR,"Compute slcsa/atom compute {} does not calculate a per-atom array", val.id);
|
||||
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate a per-atom array", val.id);
|
||||
if (val.argindex && val.argindex > val.val.c->size_peratom_cols)
|
||||
error->all(FLERR,"Compute slcsa/atom compute {} array is accessed out-of-range", val.id);
|
||||
error->all(FLERR, "Compute slcsa/atom compute {} array is accessed out-of-range", val.id);
|
||||
descriptorval = val;
|
||||
memory->create(database_mean_descriptor, ncomps, "slcsa/atom:database_mean_descriptor");
|
||||
memory->create(lda_scalings, ncomps, nclasses-1, "slcsa/atom:lda_scalings");
|
||||
memory->create(lr_decision, nclasses, nclasses-1, "slcsa/atom:lr_decision");
|
||||
memory->create(lda_scalings, ncomps, nclasses - 1, "slcsa/atom:lda_scalings");
|
||||
memory->create(lr_decision, nclasses, nclasses - 1, "slcsa/atom:lr_decision");
|
||||
memory->create(lr_bias, nclasses, "slcsa/atom:lr_bias");
|
||||
memory->create(maha_thresholds, nclasses, "slcsa/atom:maha_thresholds");
|
||||
memory->create(icov_list, nclasses, nclasses-1, nclasses-1, "slcsa/atom:icov_list");
|
||||
memory->create(mean_projected_descriptors, nclasses, nclasses-1, "slcsa/atom:mean_projected_descriptors");
|
||||
memory->create(icov_list, nclasses, nclasses - 1, nclasses - 1, "slcsa/atom:icov_list");
|
||||
memory->create(mean_projected_descriptors, nclasses, nclasses - 1,
|
||||
"slcsa/atom:mean_projected_descriptors");
|
||||
|
||||
if (comm->me == 0) {
|
||||
|
||||
if (strcmp(database_mean_descriptor_file,"NULL") == 0) {
|
||||
error->one(FLERR,"Cannot open database mean descriptor file {}: ", database_mean_descriptor_file, utils::getsyserror());
|
||||
if (strcmp(database_mean_descriptor_file, "NULL") == 0) {
|
||||
error->one(FLERR,
|
||||
"Cannot open database mean descriptor file {}: ", database_mean_descriptor_file,
|
||||
utils::getsyserror());
|
||||
} else {
|
||||
PotentialFileReader reader(lmp,database_mean_descriptor_file, "database mean descriptor file");
|
||||
int nread=0;
|
||||
while (nread<ncomps) {
|
||||
PotentialFileReader reader(lmp, database_mean_descriptor_file,
|
||||
"database mean descriptor file");
|
||||
int nread = 0;
|
||||
while (nread < ncomps) {
|
||||
auto values = reader.next_values(0);
|
||||
database_mean_descriptor[nread]=values.next_double();
|
||||
database_mean_descriptor[nread] = values.next_double();
|
||||
nread++;
|
||||
}
|
||||
}
|
||||
|
||||
if (strcmp(lda_scalings_file,"NULL") == 0) {
|
||||
error->one(FLERR,"Cannot open database linear discriminant analysis scalings file {}: ", lda_scalings_file, utils::getsyserror());
|
||||
if (strcmp(lda_scalings_file, "NULL") == 0) {
|
||||
error->one(FLERR, "Cannot open database linear discriminant analysis scalings file {}: ",
|
||||
lda_scalings_file, utils::getsyserror());
|
||||
} else {
|
||||
PotentialFileReader reader(lmp,lda_scalings_file, "lda scalings file");
|
||||
int nread=0;
|
||||
while (nread<ncomps) {
|
||||
auto values = reader.next_values(nclasses-1);
|
||||
lda_scalings[nread][0]=values.next_double();
|
||||
lda_scalings[nread][1]=values.next_double();
|
||||
lda_scalings[nread][2]=values.next_double();
|
||||
PotentialFileReader reader(lmp, lda_scalings_file, "lda scalings file");
|
||||
int nread = 0;
|
||||
while (nread < ncomps) {
|
||||
auto values = reader.next_values(nclasses - 1);
|
||||
lda_scalings[nread][0] = values.next_double();
|
||||
lda_scalings[nread][1] = values.next_double();
|
||||
lda_scalings[nread][2] = values.next_double();
|
||||
nread++;
|
||||
}
|
||||
}
|
||||
|
||||
if (strcmp(lr_decision_file,"NULL") == 0) {
|
||||
error->one(FLERR,"Cannot open logistic regression decision file {}: ", lr_decision_file, utils::getsyserror());
|
||||
if (strcmp(lr_decision_file, "NULL") == 0) {
|
||||
error->one(FLERR, "Cannot open logistic regression decision file {}: ", lr_decision_file,
|
||||
utils::getsyserror());
|
||||
} else {
|
||||
PotentialFileReader reader(lmp,lr_decision_file, "lr decision file");
|
||||
int nread=0;
|
||||
while (nread<nclasses) {
|
||||
auto values = reader.next_values(nclasses-1);
|
||||
lr_decision[nread][0]=values.next_double();
|
||||
lr_decision[nread][1]=values.next_double();
|
||||
lr_decision[nread][2]=values.next_double();
|
||||
PotentialFileReader reader(lmp, lr_decision_file, "lr decision file");
|
||||
int nread = 0;
|
||||
while (nread < nclasses) {
|
||||
auto values = reader.next_values(nclasses - 1);
|
||||
lr_decision[nread][0] = values.next_double();
|
||||
lr_decision[nread][1] = values.next_double();
|
||||
lr_decision[nread][2] = values.next_double();
|
||||
nread++;
|
||||
}
|
||||
}
|
||||
|
||||
if (strcmp(lr_bias_file,"NULL") == 0) {
|
||||
error->one(FLERR,"Cannot open logistic regression bias file {}: ", lr_bias_file, utils::getsyserror());
|
||||
if (strcmp(lr_bias_file, "NULL") == 0) {
|
||||
error->one(FLERR, "Cannot open logistic regression bias file {}: ", lr_bias_file,
|
||||
utils::getsyserror());
|
||||
} else {
|
||||
PotentialFileReader reader(lmp,lr_bias_file, "lr bias file");
|
||||
PotentialFileReader reader(lmp, lr_bias_file, "lr bias file");
|
||||
auto values = reader.next_values(nclasses);
|
||||
lr_bias[0]=values.next_double();
|
||||
lr_bias[1]=values.next_double();
|
||||
lr_bias[2]=values.next_double();
|
||||
lr_bias[3]=values.next_double();
|
||||
lr_bias[0] = values.next_double();
|
||||
lr_bias[1] = values.next_double();
|
||||
lr_bias[2] = values.next_double();
|
||||
lr_bias[3] = values.next_double();
|
||||
}
|
||||
|
||||
if (strcmp(maha_file,"NULL") == 0) {
|
||||
error->one(FLERR,"Cannot open mahalanobis stats file {}: ", maha_file, utils::getsyserror());
|
||||
if (strcmp(maha_file, "NULL") == 0) {
|
||||
error->one(FLERR, "Cannot open mahalanobis stats file {}: ", maha_file, utils::getsyserror());
|
||||
} else {
|
||||
PotentialFileReader reader(lmp,maha_file, "mahalanobis stats file");
|
||||
int nvalues = nclasses*((nclasses-1)*(nclasses-1)+nclasses);
|
||||
PotentialFileReader reader(lmp, maha_file, "mahalanobis stats file");
|
||||
int nvalues = nclasses * ((nclasses - 1) * (nclasses - 1) + nclasses);
|
||||
auto values = reader.next_values(nvalues);
|
||||
|
||||
for(int i=0;i<nclasses;i++){
|
||||
maha_thresholds[i]=values.next_double();
|
||||
for(int j=0;j<nclasses-1;j++) {
|
||||
mean_projected_descriptors[i][j]=values.next_double();
|
||||
}
|
||||
for(int k=0;k<nclasses-1;k++) {
|
||||
for(int l=0;l<nclasses-1;l++) {
|
||||
icov_list[i][k][l]=values.next_double();
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < nclasses; i++) {
|
||||
maha_thresholds[i] = values.next_double();
|
||||
for (int j = 0; j < nclasses - 1; j++)
|
||||
mean_projected_descriptors[i][j] = values.next_double();
|
||||
for (int k = 0; k < nclasses - 1; k++)
|
||||
for (int l = 0; l < nclasses - 1; l++) icov_list[i][k][l] = values.next_double();
|
||||
}
|
||||
|
||||
for(int i=0;i<nclasses;i++){
|
||||
std::cout << "For class " << i << " maha threshold = " << maha_thresholds[i] << std::endl;
|
||||
std::cout << "\tmean B: " << std::endl;
|
||||
for(int j=0;j<nclasses-1;j++) {
|
||||
std::cout << "\t\t" << mean_projected_descriptors[i][j] << std::endl;
|
||||
}
|
||||
std::cout << "\ticov: " << std::endl;
|
||||
for(int j=0;j<nclasses-1;j++) {
|
||||
std::cout << "\t\t" << icov_list[i][j][0] << " " << icov_list[i][j][1] << " " << icov_list[i][j][2] << std::endl;
|
||||
for (int i = 0; i < nclasses; i++) {
|
||||
auto mesg = fmt::format("For class {} maha threshold = {:.6}\n", i, maha_thresholds[i]);
|
||||
mesg += " mean B:\n";
|
||||
for (int j = 0; j < nclasses - 1; j++)
|
||||
mesg += fmt::format(" {:11.6}\n", mean_projected_descriptors[i][j]);
|
||||
mesg += " icov:\n";
|
||||
for (int j = 0; j < nclasses - 1; j++) {
|
||||
mesg += fmt::format(" {:11.6} {:11.6} {:11.6}\n", icov_list[i][j][0],
|
||||
icov_list[i][j][1], icov_list[i][j][2]);
|
||||
}
|
||||
utils::logmesg(lmp, mesg);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
MPI_Bcast(&database_mean_descriptor[0],ncomps,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lda_scalings[0][0],ncomps*(nclasses-1),MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lr_decision[0][0],nclasses*(nclasses-1),MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lr_bias[0],nclasses,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&maha_thresholds[0],nclasses,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&mean_projected_descriptors[0][0],nclasses*(nclasses-1),MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&icov_list[0][0][0],nclasses*(nclasses-1)*(nclasses-1),MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&database_mean_descriptor[0], ncomps, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&lda_scalings[0][0], ncomps * (nclasses - 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&lr_decision[0][0], nclasses * (nclasses - 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&lr_bias[0], nclasses, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&maha_thresholds[0], nclasses, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&mean_projected_descriptors[0][0], nclasses * (nclasses - 1), MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&icov_list[0][0][0], nclasses * (nclasses - 1) * (nclasses - 1), MPI_DOUBLE, 0, world);
|
||||
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = nclasses+1;
|
||||
ncols=nclasses+1;
|
||||
nmax=0;
|
||||
size_peratom_cols = nclasses + 1;
|
||||
ncols = nclasses + 1;
|
||||
nmax = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -254,6 +268,8 @@ ComputeSLCSAAtom::~ComputeSLCSAAtom()
|
||||
memory->destroy(projected_descriptor);
|
||||
memory->destroy(scores);
|
||||
memory->destroy(probas);
|
||||
memory->destroy(prodright);
|
||||
memory->destroy(dmaha);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -296,76 +312,69 @@ void ComputeSLCSAAtom::compute_peratom()
|
||||
descriptorval.val.c->compute_peratom();
|
||||
descriptorval.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
|
||||
}
|
||||
compute_array = descriptorval.val.c->array_atom;
|
||||
compute_array = descriptorval.val.c->array_atom;
|
||||
}
|
||||
|
||||
memory->create(full_descriptor, ncomps, "slcsa/atom:local descriptor");
|
||||
memory->create(projected_descriptor, nclasses-1, "slcsa/atom:reduced descriptor");
|
||||
memory->create(scores, nclasses, "slcsa/atom:probas");
|
||||
memory->create(projected_descriptor, nclasses - 1, "slcsa/atom:reduced descriptor");
|
||||
memory->create(scores, nclasses, "slcsa/atom:scores");
|
||||
memory->create(probas, nclasses, "slcsa/atom:probas");
|
||||
memory->create(prodright, nclasses-1, "slcsa/atom:prodright");
|
||||
memory->create(prodright, nclasses - 1, "slcsa/atom:prodright");
|
||||
memory->create(dmaha, nclasses, "slcsa/atom:prodright");
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
for (int j = 0; j < ncomps; j++) {
|
||||
full_descriptor[j] = compute_array[i][j];
|
||||
}
|
||||
for (int j = 0; j < ncomps; j++) { full_descriptor[j] = compute_array[i][j]; }
|
||||
// Here comes the LDA + LR process
|
||||
// 1st step : Retrieve mean database descriptor
|
||||
for (int j = 0; j < ncomps; j++) {
|
||||
full_descriptor[j] -= database_mean_descriptor[j];
|
||||
}
|
||||
for (int j = 0; j < ncomps; j++) { full_descriptor[j] -= database_mean_descriptor[j]; }
|
||||
// 2nd step : Matrix multiplication to go from ncompsx1 -> (nclasses-1)*1
|
||||
for (int j = 0; j < nclasses-1; j++) {
|
||||
for (int j = 0; j < nclasses - 1; j++) {
|
||||
projected_descriptor[j] = 0.;
|
||||
for (int k = 0; k < ncomps; k++) {
|
||||
projected_descriptor[j] += full_descriptor[k]*lda_scalings[k][j];
|
||||
projected_descriptor[j] += full_descriptor[k] * lda_scalings[k][j];
|
||||
}
|
||||
}
|
||||
// 3rd step : Matrix multiplication
|
||||
for (int j = 0; j < nclasses; j++) {
|
||||
scores[j] = lr_bias[j];
|
||||
for (int k = 0; k < nclasses-1; k++) {
|
||||
for (int k = 0; k < nclasses - 1; k++) {
|
||||
scores[j] += lr_decision[j][k] * projected_descriptor[k];
|
||||
}
|
||||
}
|
||||
// 4th step : Matrix multiplication
|
||||
double sumexpscores=0.;
|
||||
double sumexpscores = 0.;
|
||||
for (int j = 0; j < nclasses; j++) sumexpscores += exp(scores[j]);
|
||||
for (int j = 0; j < nclasses; j++) {
|
||||
probas[j] = exp(scores[j])/sumexpscores;
|
||||
}
|
||||
for (int j = 0; j < nclasses; j++) { probas[j] = exp(scores[j]) / sumexpscores; }
|
||||
|
||||
classification[i][nclasses]=argmax(probas,nclasses);
|
||||
classification[i][nclasses] = argmax(probas, nclasses);
|
||||
|
||||
// 5th step : Mahalanobis distance
|
||||
for (int j = 0; j < nclasses; j++) {
|
||||
prodright[0]=0.;
|
||||
prodright[1]=0.;
|
||||
prodright[2]=0.;
|
||||
for (int k = 0; k < nclasses-1; k++) {
|
||||
for (int l = 0; l < nclasses-1; l++) {
|
||||
prodright[k] += (icov_list[j][k][l] * (projected_descriptor[k]-mean_projected_descriptors[j][k]));
|
||||
prodright[0] = 0.;
|
||||
prodright[1] = 0.;
|
||||
prodright[2] = 0.;
|
||||
for (int k = 0; k < nclasses - 1; k++) {
|
||||
for (int l = 0; l < nclasses - 1; l++) {
|
||||
prodright[k] +=
|
||||
(icov_list[j][k][l] * (projected_descriptor[k] - mean_projected_descriptors[j][k]));
|
||||
}
|
||||
}
|
||||
double prodleft=0.;
|
||||
for (int k = 0; k < nclasses-1; k++) {
|
||||
prodleft += (prodright[k]*(projected_descriptor[k]-mean_projected_descriptors[j][k]));
|
||||
double prodleft = 0.;
|
||||
for (int k = 0; k < nclasses - 1; k++) {
|
||||
prodleft += (prodright[k] * (projected_descriptor[k] - mean_projected_descriptors[j][k]));
|
||||
}
|
||||
classification[i][j]=sqrt(prodleft);
|
||||
classification[i][j] = sqrt(prodleft);
|
||||
}
|
||||
// 6th step : Sanity check
|
||||
int locclass=classification[i][nclasses];
|
||||
int locclass = classification[i][nclasses];
|
||||
|
||||
if (classification[i][locclass]>maha_thresholds[locclass]){
|
||||
classification[i][nclasses]=-1.;
|
||||
if (classification[i][locclass] > maha_thresholds[locclass]) {
|
||||
classification[i][nclasses] = -1.;
|
||||
}
|
||||
|
||||
} else {
|
||||
for (int j = 0; j < ncols; j++) {
|
||||
classification[i][j] = -1.;
|
||||
}
|
||||
for (int j = 0; j < ncols; j++) { classification[i][j] = -1.; }
|
||||
}
|
||||
}
|
||||
memory->destroy(full_descriptor);
|
||||
@ -373,7 +382,6 @@ void ComputeSLCSAAtom::compute_peratom()
|
||||
memory->destroy(scores);
|
||||
memory->destroy(probas);
|
||||
memory->destroy(prodright);
|
||||
|
||||
}
|
||||
|
||||
int ComputeSLCSAAtom::compute_ncomps(int twojmax)
|
||||
@ -384,24 +392,24 @@ int ComputeSLCSAAtom::compute_ncomps(int twojmax)
|
||||
|
||||
for (int j1 = 0; j1 <= twojmax; j1++)
|
||||
for (int j2 = 0; j2 <= j1; j2++)
|
||||
for (int j = j1 - j2;
|
||||
j <= MIN(twojmax, j1 + j2); j += 2)
|
||||
for (int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
|
||||
if (j >= j1) ncount++;
|
||||
|
||||
return ncount;
|
||||
}
|
||||
|
||||
int ComputeSLCSAAtom::argmax(double arr[], int size) {
|
||||
int maxIndex = 0; // Initialize the index of the maximum value to the first element.
|
||||
double maxValue = arr[0]; // Initialize the maximum value to the first element.
|
||||
int ComputeSLCSAAtom::argmax(double arr[], int size)
|
||||
{
|
||||
int maxIndex = 0; // Initialize the index of the maximum value to the first element.
|
||||
double maxValue = arr[0]; // Initialize the maximum value to the first element.
|
||||
|
||||
for (int i = 1; i < size; ++i) {
|
||||
if (arr[i] > maxValue) {
|
||||
// If a greater value is found, update the maxIndex and maxValue.
|
||||
maxIndex = i;
|
||||
maxValue = arr[i];
|
||||
}
|
||||
for (int i = 1; i < size; ++i) {
|
||||
if (arr[i] > maxValue) {
|
||||
// If a greater value is found, update the maxIndex and maxValue.
|
||||
maxIndex = i;
|
||||
maxValue = arr[i];
|
||||
}
|
||||
}
|
||||
|
||||
return maxIndex;
|
||||
return maxIndex;
|
||||
}
|
||||
|
||||
@ -38,6 +38,7 @@ class ComputeSLCSAAtom : public Compute {
|
||||
// double memory_usage() override;
|
||||
int compute_ncomps(int);
|
||||
int argmax(double *, int);
|
||||
|
||||
private:
|
||||
struct value_t {
|
||||
int which; // type of data: COMPUTE, FIX, VARIABLE
|
||||
@ -86,7 +87,6 @@ class ComputeSLCSAAtom : public Compute {
|
||||
|
||||
// Output array
|
||||
double **classification;
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
Reference in New Issue
Block a user