fix deletion of seedlot by deleting from stock_owner first closes #1732
[sgn.git] / R / MADII_layout_function.R
blob5b125d784953fe529598a6bbfd81005c5d63797a
1 ### Function for generating a randomized, MADII experimental design
2 ## Tyler Tiede. University of Minnesota. March 7, 2014
3 #cat (" \n** Before using MADIIdgn be sure 
4 #that there are no files and/or folders in your Working Directory named\nthe same as 
5 #the character string you plan to designate for 'designID'; they will be overwritten.\n ")
7 #List of checks starting with the primary; assumes you have a primary check
9 MADIIdgn <- function(enviro="Eretz", 
10                      entries= NULL, 
11                      num.entries= NULL, 
12                      chk.names= NULL, 
13                      num.sec.chk= NULL, 
14                      num.rows= NULL, 
15                      num.cols= NULL,
16                      
17                      plot.start = 1001, 
18                      designID=NULL, 
19                      annoy=F)
20   {
22   ## Ensure that no folders/files in WD will be overwritten and then create folder for results
23   #JL "is.null(designID)==F"  is equivalent to "!is.null(designID)"
24   if(annoy==T & is.null(designID)==F){
25   YorN <- NA    
26     while(YorN %in% c("y","n","Y","N","yes","no","YES","NO","Yes","No") == F){
27       YorN <- readline("Confirm there is Not a folder and/or file in your Working Directory with\nthe same name as 'designID' (type 'y' to proceed, or 'n' to stop): ")
28     }  
29     #JL It doesn't look like typing 'n' will actually stop the function...
30     # I would put in something like if(YorN %in% c("n", "N", "no", "No", "NO")) stop("User requested stop")
31     if(YorN %in% c("y","Y","yes","YES","Yes")==T){ 
32       unlink(designID, recursive=T)
33       dir.create(path=designID, recursive=F) 
34     }
35   } else{
36     #if(is.null(designID)==F){
37       unlink(designID, recursive=T)
38       dir.create(path=designID, recursive=F) 
39     #}
40   }
41    
42   
43   ## Load necessary packag#es
44   require("grid", quietly=T) ; library(grid)
45   
46   ## QC of function inputs
47   #JL why are these only warnings and not errors that cause the function to stop?
48   if(is.null(entries) & is.null(num.entries)){
49     warning("Must provide an entry list (entries=) OR the number of entries desired (num.entries=).")
50   }
51   
52   if(is.null(chk.names) & is.null(num.sec.chk)){
53     warning("Must provide a list of check names (chk.names=) with the primary check listed first\n OR the number of SECONDARY checks desired (num.sec.chk=).")
54   }
55   
56   if(is.null(num.rows)){
57     warning("Provide the number of rows (sometimes called ranges or beds) (num.rows=).")
58   }
59   
60   if(num.rows %% 3 != 0){
61     warning("The MADII design requires that the number of rows be a multiple of 3.")
62   }
63   
64   
65   ## Develop other non-input function parameters
66   if(is.null(entries)){
67     entries <- as.matrix(paste("entry", 1:num.entries, sep="_"))
68   }
69   
70   if(is.null(num.entries)){
71     entries <- as.matrix(entries) # If user input of entries was a list, will convert it to matrix
72     num.entries <- nrow(entries)
73   }
74   
76   ## This warning is dependent on the number of entries
77   if(is.null(num.cols)==F){
78     if(((num.cols / 5) * (num.rows / 3) + num.entries) > num.cols*num.rows){
79       warning("The minimum number of plots has not been met by the given row/column dimensions.")
80     }
81   }
82   
83   if(is.null(chk.names)){
84     sec.chks <- as.character(2:(num.sec.chk+1)) ## Do need as seperate object for later on in function
85     chk.names <- paste("chk", c(1,sec.chks), sep="") ## All generic check names
86   }
87   
88   if(is.null(num.sec.chk)){
89     sec.chks <- chk.names[-1] ## The primary check must be listed first in the function input
90     num.sec.chk <- length(sec.chks)
91   }
92   
93   blk.rows <- num.rows / 3 # This is the standard for the MADII design
94   
95   ## If the number of columns is provided then it is straight forward, otherwise a set of algorithms will develop and optimize the paramaters
96   if(is.null(num.cols)==F){
97     
98     if(num.cols%%5 != 0){
99       warning("The MADII design requires that the number of columns be a multiple of 5.")
100     }
101     
102     blk.cols <- num.cols / 5 # This is the standard for the MADII design
103     
104     num.blks <- blk.rows * blk.cols
105     exp.size <- num.blks * 15
106     num.chks <- exp.size - num.entries
107     #JL is there a rule that if a block has secondary checks, then it must have _all_ secondary checks?
108     # That would not seem like a good rule.  Maybe it is required by some MADII analyses but it would not be
109     # required by the moving average analysis or by an AR1 analysis, and my guess is those are better.
110     num.sec.chk.blks <- floor((num.chks - num.blks) / num.sec.chk) # one primary check in each block
111     # Number of checks to total number of plots: I think per.chks gets awfully high depending on the design.  I think there
112     # should be no more than two check plots per block of 15.  It's just too much effort otherwise.
113     per.chks <- (num.blks + (num.sec.chk.blks*num.sec.chk)) / (num.entries + num.blks + (num.sec.chk.blks*num.sec.chk))
114     num.fill <- exp.size - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
115     #JL do Fill plots actually get planted to something?  I am thinking that they would be.  If they are, it would seem better to 
116     #put secondary checks in them and distribute them randomly around the experiment
117     
118     if(is.null(designID)==F){
119       write.table(per.chks, paste(designID, "/", "%checks_in_", designID, ".txt", sep=""))
120     }else{
121       write.table(per.chks, "%checks.txt")
122     }
123     
124     
125   }else{
126     ## If the number of columns is not specified, below algorithms will develop the necessary design
127     ## Calculate starting (non-optimized paramaters)
128     per.chks <- 0.10 #JL yes, this is a good number
129     
130     ## Number of total checks in experiment (primary + secondary) ; calculated as percent of entries
131     num.chks <- ceiling((per.chks * num.entries) / (1-per.chks)) 
132     entries.plus.chks <- num.entries + num.chks
133     
134     num.cols <- ceiling(entries.plus.chks / num.rows)
135     #JL num.cols <- ceiling(num.cols / 5) * 5
136     # so the statement above could just be "num.cols <- ceiling(entries.plus.chks / num.rows / 5) * 5"
137     while(num.cols %% 5 != 0){
138       num.cols <- num.cols + 1
139     }
140     
141     blk.cols <- num.cols / 5 # This is the standard for the MADII design
142     num.blks <- blk.rows * blk.cols
143     exp.size <- num.blks*15 # 15 plots per block
144     num.sec.chk.blks <- ceiling((num.chks - num.blks) / num.sec.chk) # one primary check in each block
145     
146     ## If the ratio of blk.cols to num.sec.chk.blks does not allow each blk.col to have a sec.chk blk in it, then optimize per.chks
147     while((blk.cols > num.sec.chk.blks) & (num.sec.chk.blks <  num.blks)){
148       per.chks <- per.chks + 0.0001
149       num.chks <- ceiling((per.chks * num.entries) / (1-per.chks))
150       entries.plus.chks <- num.entries + num.chks
151       
152       num.cols <- ceiling(entries.plus.chks / num.rows)
153       
154       while(num.cols %% 5 != 0){
155         num.cols <- num.cols + 1
156       }
157       
158       blk.cols <- num.cols / 5 # This is the standard for the MADII design
159       num.blks <- blk.rows * blk.cols
160       exp.size <- num.blks*15 # 15 plots per block
161       num.sec.chk.blks <- ceiling((num.chks - num.blks) / num.sec.chk) # one primary check in each block
162       
163     }
164     
165     num.fill <- num.blks*15 - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
166   
167     ## Increase number of checks to minimize number of Fill plots
168     while(num.fill >= num.sec.chk & (num.sec.chk.blks <  num.blks)){ 
169       per.chks <- per.chks + 0.0001
170       num.chks <- ceiling((per.chks * num.entries) / (1-per.chks))
171       entries.plus.chks <- num.entries + num.chks
172       
173       num.sec.chk.blks <- floor((num.chks - num.blks) / num.sec.chk) # one primary check in each block
174       
175       num.fill <-num.blks*15 - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
176       
177     }
178     
179     if(is.null(designID)==F){
180       write.table(per.chks, paste(designID, "/", "%checks_in_", designID, ".txt", sep=""))
181     }else{
182       write.table(per.chks, "%checks.txt")
183     }
184   
185   }
186   
187   
188   #################################################################
189   ############### Build Field File ################################
190   
191   ## Put together field design columns; plot, row, col, blk, row.blk, col.blk
192   fld.dgn <- as.data.frame(matrix(data=c((plot.start:(plot.start-1+exp.size)), 
193                                           rep(1:num.rows, each=num.cols), rep(c(1:num.cols, num.cols:1), length.out=exp.size), 
194                                           rep(NA, times=exp.size), rep(1:blk.rows, each=(exp.size/blk.rows)), 
195                                           rep(c(1:blk.cols, blk.cols:1), each=5, length.out=exp.size), 
196                                           rep(NA, times=2*exp.size)), nrow=exp.size, byrow=F))
197   
198   colnames(fld.dgn) <- c("Plot", "Row", "Col", "Blk", "Row.Blk", "Col.Blk", "Line.Code", "Entry")
199   
200   if(num.fill>0){
201     fld.dgn[(exp.size-num.fill+1):exp.size, 7] <- "F"
202   }
203   
204   blk.list <- 1:num.blks
205   for(b in 1:blk.rows){
206     if((b %% 2 == 0)==T){
207       blk.list[(1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols)] <- rev((1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols)) 
208     } else{
209       blk.list[(1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols)] <- ((1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols))
210     }
211   }
212   
213   ## Assign plots to blocks
214   count <- 1
215   for(b in 1:blk.rows){
216     for(c in 1:blk.cols){
217       blk <- blk.list[count]
218       r1 <- (1+seq(0,1000,by=3))[b]
219       r2 <- (seq(0,1000,by=3))[b+1]
220       
221       c1 <- (1+seq(0,250000,by=5))[c]
222       c2 <- (seq(0,250000,by=5))[c+1]
223       
224       fld.dgn[which(fld.dgn$Row %in% r1:r2 & fld.dgn$Col %in% c1:c2), 4] <- blk
225       
226       count <- count+1
227     } 
228   }
229   
230   
231   ## Selecting secondary blocks randomly, but ensuring that each blk.row and each blk.col is represented at least once
233   if((num.sec.chk.blks > blk.cols)==T){
234     satisfied <- F
235     col.blk.list <- c(1:blk.cols)
236     row.blk.list <- c(1:blk.rows)
237     while(satisfied != T){
238       blk.cols.rep <- NULL
239       blk.rows.rep <- NULL
240       sample.kept <- c()
241       new.length <- c()
242       while((length(blk.cols.rep) < blk.cols)==T){
243         length <- new.length
244         sample <- sample(1:(num.blks), 1, replace=F)
245         blk.cols.rep <- c(blk.cols.rep, unique(fld.dgn[which(fld.dgn$Blk==sample), 6]))
246         blk.cols.rep <- unique(blk.cols.rep)
247         new.length <- length(blk.cols.rep)
248         if(new.length == 1){
249           sample.kept <- c(sample.kept, sample)
250         }else{
251           if(new.length == (length+1))
252             sample.kept <- c(sample.kept, sample)
253         }
254       }
255       
256       for(s in sample.kept){
257         blk.rows.rep <- c(blk.rows.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 5]))
258         blk.rows.rep <- unique(blk.rows.rep)
259       } 
260     
261       if(all(col.blk.list %in% blk.cols.rep) & all(row.blk.list %in% blk.rows.rep)){
262         satisfied <- T
263       }   
264       
265     }
266     
267     ## Need to bring "sample" back up to the number of secondary check blocks with some randomly chosen blocks
268     cat((num.sec.chk.blks - length(sample.kept)), "\n")
269     sample2 <- sample(setdiff(1:num.blks, sample.kept), (num.sec.chk.blks - length(sample.kept)))
270     sample <- c(sample.kept, sample2)
271     
272   }else{
273     ## Selecting secondary blocks randomly, but ensuring that if sec.chk.blk < blk.cols then >1 sec.chk.block does not end up in a column
274     satisfied <- F
275     col.blk.list <- c(1:blk.cols)
276     row.blk.list <- c(1:blk.rows)
277     while(satisfied != T){
278       blk.cols.rep <- NULL
279       blk.rows.rep <- NULL
280       #sample <- sample(1:(num.blks-(ceiling(num.fill/5))), num.sec.chk.blks, replace=F) 
281       # Use this line instead of below line if secondary checks are not wanted in partial plots due to Filler
282       sample <- sample(1:(num.blks), num.sec.chk.blks, replace=F) # Using this line will allow secondary blocks to contain fill plots
283       for(s in sample){
284         blk.cols.rep <- c(blk.cols.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 6]))
285         blk.rows.rep <- c(blk.rows.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 5]))
286       }
287       if(length(unique(blk.cols.rep))==length(blk.cols.rep) & all(row.blk.list %in% blk.rows.rep)){
288         satisfied <- T
289       }   
290     }
291   }
292   
293   ## Assign primary checks to field design
294   for(b in 1:num.blks){
295     blk <- fld.dgn[which(fld.dgn$Blk==b), ]
296     row <- which(fld.dgn$Row == mean(blk$Row))
297     col <- which(fld.dgn$Col == mean(blk$Col))
298     fld.dgn[row[which(row  %in% col ==T)], 7] <- 1
299   }
300   
301   ## Assign secondary checks to field
302   for(s in sample){
303     sec.blk <- fld.dgn[which((fld.dgn$Blk==s) & is.na(fld.dgn$Line.Code)),][,1]
304     sec.plots <- sec.blk[sample(1:length(sec.blk), num.sec.chk)]
305     for(i in 1:length(sec.plots)){
306       plot <- sec.plots[i]
307       fld.dgn[which(fld.dgn$Plot==plot), 7] <- (i+1)
308     }
309   }
311   
312   fld.dgn[which(is.na(fld.dgn$Line.Code)), 7] <- 0
313   
314   options(warn=-1) ## If no fill plots then the lines below will throw an error
315   ## Assign entry names and check names
316   fld.dgn[which(fld.dgn$Line.Code == 0), 8] <- entries[order(sample(1:length(entries), length(entries)))]
317   if(num.fill>0){
318     fld.dgn[which(fld.dgn$Line.Code == "F"), 8] <- "Fill"
319   }
320   fld.dgn[which(fld.dgn$Line.Code == "F"), 7] <- 0 # Change the fill lines back to experimental entries for purposes of MADIIadj
321   options(warn=0) # Turn warnings back on
323   
324   for(c in 1:length(chk.names)){
325     chk <- chk.names[c]
326     fld.dgn[which(fld.dgn$Line.Code==c), 8] <- chk
327   }
328   
329   # Instead of re-working code, just re-order fld.dgn df before reading out
330   
331   
332   fld.dgn.mod <- cbind(matrix(enviro, nrow=num.rows, ncol=1), fld.dgn) ; colnames(fld.dgn.mod)[c(1,8)] <- c("Enviro", "Check")
333   fld.dgn.mod <- fld.dgn.mod[,c(1,2,9,8, 3:7)]
334   
335   if(is.null(designID)==F){
336     assign(designID, fld.dgn.mod, envir=.GlobalEnv)
337     cat("\nYour MADII design is available in your working directory and R environment; named per ID provided.")
338     write.csv(get(designID), paste(designID,"/",designID, ".csv", sep=""), row.names=F)
339   } else{
340     assign("YourMADIIdgn", fld.dgn.mod, envir=.GlobalEnv)
341     cat("\nYour MADII design is available in your working directory and R environment; named 'YourMADIIdgn'.")
342     write.csv(YourMADIIdgn, "YourMADIIdgn.csv", row.names=F)
343   }
344   
345   ########## Generating Visual Field Map ##########
346   
347   ## First will generate a series of .csv files for excel
348   ### Plots
349   for.fld.plot <- fld.dgn[order(fld.dgn$Row, fld.dgn$Col),1:3]
350   fld.plot <- apply(t(matrix(for.fld.plot$Plot, nrow=num.rows, ncol=num.cols, byrow=T)), 1, rev) 
351   # transposes matrix then flips all rows horizontally
352   fld.plot <- cbind(matrix(num.rows:1, nrow=num.rows, ncol=1), fld.plot)
353   fld.plot <- rbind(fld.plot, matrix(c(" ", 1:num.cols), nrow=1, ncol=ncol(fld.plot)))
354   
355   ### Line codes
356   for.line.code <- fld.dgn[order(fld.dgn$Row, fld.dgn$Col), c(2,3,7)]
357   line.code <- apply(t(matrix(for.line.code$Line.Code, nrow=num.rows, ncol=num.cols, byrow=T)), 1, rev) 
358   #transposes matrix then flips all rows horizontally
359   line.code <- cbind(matrix(num.rows:1, nrow=num.rows, ncol=1), line.code)
360   line.code <- rbind(line.code, matrix(c(" ", 1:num.cols), nrow=1, ncol=ncol(line.code)))
361   
362   ### Write out the design .csv files
363   if(is.null(designID)==F){
364     write.table(fld.plot, paste(designID,"/",designID, "_plot.lay.csv", sep=""), row.names=F, col.names=F, sep=",")
365     write.table(line.code, paste(designID, "/", designID, "_line.codes.lay.csv", sep=""), row.names=F, col.names=F, sep=",")
366   }else{
367     write.table(fld.plot, paste(designID, "plot.lay.csv", sep="."), row.names=F, col.names=F, sep=",")
368     write.table(line.code, paste(designID, "line.codes.lay.csv", sep="."), row.names=F, col.names=F, sep=",")
369   }
370   
371   
372   ###### Now a .png image ######
373   options(warn=-1)
374   plot.new()
375   
376   if(is.null(designID)==F){
377     dev.copy(png, paste(designID, "/", designID, "_plot.png", sep=""))
378   }else{
379     dev.copy(png, "YourMADIIdgn.plog.png")
380   }
381   
382   
383   field.plot <- matrix(nrow=num.rows, ncol=num.cols)
384   field.NA.plot <- t(matrix(NA, nrow=num.rows, ncol=num.cols))
385   for(r in 1:num.rows){
386     for(c in 1:num.cols){
387       field.plot[r,c] <- fld.dgn[which(fld.dgn$Row==r & fld.dgn$Col==c), 1]
388     }
389   }
390   field.plot <- t(field.plot)
391   
392   for(i in 1:(num.sec.chk+1)){
393     chk.plots <- as.list(fld.dgn$Plot[which(fld.dgn$Line.Code==i)])
394     field.NA.plot[which(field.plot %in% chk.plots)] <- i
395   }
396   
397   fill.plots <- as.list(fld.dgn$Plot[which(fld.dgn$Entry=="Fill")])
398   field.NA.plot[which(field.plot %in% fill.plots)] <- (i+1)
399   
400   colors <- c("red", "blue", "green", "orange", "gray", "brown", "yellow", "purple", "cyan")
401   par(new=T,mar=c(4.5,2,.5,.5), srt=0, xpd=T, font=2, bg="white")
402   s
403   image(field.plot, add=F, col=0, xaxt='n', yaxt='n' ,xlab="Columns", ylab="Rows", line=0.5)
404   image(field.NA.plot, add=T, col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])
405   
406   box(lwd=3)
407   
408   if(num.fill>0){
409     legend("bottom",legend=c(chk.names, "Fill"), bty="n", horiz=T, cex=.75, inset=c(0,-.125), 
410            xjust=0.5, text.col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])}
411   else{
412     legend("bottom",legend=c(chk.names), bty="n", horiz=T, cex=.75, inset=c(0,-.125), 
413            xjust=0.5, text.col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])
414   }
415   
416   
417   par(xpd=F)
418   grid(ny=num.rows, nx=num.cols, col=1, lty=1 )
419   grid(ny=blk.rows, nx=blk.cols, col=1, lty=1, lwd=3) 
420   
421   dev.off()
422   
423   cat("\n\nA graphic of the design can also be found in the project's folder.\n ")
424   
425   options(warn=0)
426   
427 } # End of MADIIdgn loop
429 test <- MADIIdgn(num.entries=330, num.rows=9, num.cols=NULL, num.sec.chk=3, designID="tester1", annoy=T)
431 #test2 <- MADIIdgn(num.entries=100, num.rows=6, num.cols=5, num.sec.chk=3, designID="tester1", annoy=T)
433 #layoutParameters <- function(num.entries=330, num.rows=9, num.cols=NULL, num.sec.chk=3, designID="tester1", annoy=T){
434 #  minBlks <- ceiling(num.entries / 14)
435 #