-
Notifications
You must be signed in to change notification settings - Fork 0
/
edgeParticleSimulator_id007.R
152 lines (129 loc) · 6.67 KB
/
edgeParticleSimulator_id007.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
# Copyright 2019 A. London, B. Jenkins
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
# Edge particle simulation
StartTime <- Sys.time()
# Generate a set of non-colliding particles (possibly with a core-shell structure)
setwd(dir = "../posgen")
source('../AlphaEdge/particlePlace.R')
source('../AlphaEdge/edgeParticleSimulator_xmlGen_all.R')
source('../AlphaEdge/findEdgeClustersConvex.R')
source('../AlphaEdge/read.pos.sampled.R')
source('../AlphaEdge/read.pos.R')
AtomicDensity <- 20
boxSize <- 40
runs <- 2000
expId <- "id007"
# Generate pos file using posgen xml script
# Read template, write particles
summaryList <- vector("list", length = runs)
for (fileId in 1:runs) {
print(fileId)
particleCompositionCore <- t(matrix(c(1, 1), nrow = 2))
particleCompositionShell <- t(matrix(c(2, 1), nrow = 2))
particles <- edgeParticleSim_coreShell(fileName = sprintf("%03d_gen.xml",fileId),
posFileName = sprintf("%03d_full.pos",fileId),
ionDensity = AtomicDensity,
numParticles = 40,
boxSize = boxSize,
boxShape = "sphere",
rCore = 1.587,
rShell = 2,
compShell = particleCompositionShell,
compCore = particleCompositionCore,
matrixComp = data.frame(mass=c(1,2,3),count=c(.1,.1,99.8)),
clusterstatsFileName = sprintf("%03d_stats.txt",fileId),
unclusterstatsFileName = sprintf("%03d_stats_matrix.txt",fileId),
sizedistFileName = sprintf("%03d_sizeDist.txt",fileId),
clusteredposFileName = sprintf("%03d_cluster.pos",fileId),
unclusteredposFileName = sprintf("%03d_matrix.pos",fileId),
clusteridposFileName = sprintf("%03d_cluster.index.pos",fileId))
# Run posgen
# Perform cluster selection as part of posgen
# This generates the files listed in the input XML file
posgenOut <- system(sprintf("posgen.exe %03d_gen.xml",fileId), wait = TRUE, intern=TRUE)
# Find edge clusters
edgeClustersAlpha <- findEdgeClustersConvex(posFileName = sprintf("%03d_full.pos",fileId),
clusterStatsFile = sprintf("%03d_stats.txt",fileId),
clusterIndexPosFile = sprintf("%03d_cluster.index.pos",fileId),
AtomicDensity = AtomicDensity,
DetectionEfficiency = 1,
SamplingFraction = 0.01,
AlphaValue = 12)
# Analyse results
# Make data frame of: testNo, clusterNo, simSize, detected?, alphaEdge?, actualEdge, composition, measuredRg, soluteCount, matrixCount,
# Need to analyse composion
# Which clusters were identified, and which ones were edge clusters (make truth table)
# match generated to detected clusters (by nearest centre distance)
df1 <- select(particles,"x","y","z")
df2 <- select(edgeClustersAlpha$ClusterImportedData, x=xpos, y=ypos, z=zpos)
distances <- as.matrix(dist(bind_rows(df1, df2)))
row.start <- nrow(df1)+1
row.end <- nrow(df1) + nrow(df2)
col.start <- 1
col.end <- nrow(df1)
# This \/ returns a list as long as the detected clusters, with which number it matches in the original list
distanceIndex<-apply(distances[row.start:row.end, col.start:col.end], 1, which.min)
# This is the distance between those points (should be within original cluster radius)
#d<-apply(distances[row.start:row.end, col.start:col.end], 1, min)
# Construct data frame of results
summaryResults <- particles
colnames(summaryResults) <- paste("sim", colnames(summaryResults), sep = "_")
temp <- bind_cols(summaryResults[distanceIndex,], edgeClustersAlpha$ClusterImportedData)
summaryResults <- bind_rows(temp, summaryResults[-distanceIndex,])
# Define if they were edge or not based on the simulated positions
# If any of (x,y,z) +/- size exceeds the upper or lower bounding box limit then they are edge clusters
summaryResults$sim_edge <- ((summaryResults$sim_R+summaryResults$sim_r)> boxSize/2)
summaryResults$detected <- !is.na(summaryResults$Y)
summaryResults$edge <- FALSE
summaryResults$edge[edgeClustersAlpha$TotalEdgeClusters] <- TRUE
summaryResults$ID <- fileId
summaryList[[fileId]] <- summaryResults
# clean up
# copy generated file to separate folder "results/"
if (!dir.exists("results")){
# create results sub-directory
dir.create("results")
}
if (!dir.exists(paste0("results/results_", expId))){
# create results sub-directory
dir.create(paste0("results/results_", expId))
}
file.copy(from = c(
sprintf("%03d_stats.txt",fileId),
sprintf("%03d_sizeDist.txt",fileId),
sprintf("%03d_gen.xml",fileId)
),
to = paste0("results/results_", expId),
overwrite = TRUE)
# delete files to free up space
if (TRUE) {
file.remove(
sprintf("%03d_full.pos",fileId),
sprintf("%03d_cluster.index.pos",fileId),
sprintf("%03d_stats.txt",fileId),
sprintf("%03d_sizeDist.txt",fileId),
sprintf("%03d_gen.xml",fileId))
}
## Timings ##
iterAvgTime <- (Sys.time() - StartTime)/fileId
# Time remaining = iters remaining * iterAvgTime
# End time = sys.time + timeRemaining
print(paste(c("finish time:", as.character(Sys.time()+iterAvgTime*(runs-fileId)))))
# save the data as an RData format
assign(expId, bind_rows(summaryList))
save("id007", file="id007_data")
}
EndTime <- Sys.time()
TotalTime <- EndTime - StartTime
print(TotalTime)