replicate.patterns.function.only.R

replicate_patterns <- function(Status,method,mark, pattern_to_murge) {
                                        #
    #method == "remove" OR "size500" if we remove 500 size plot or put every plot at 500 size respectively
    library(data.table)
    library(spatstat)
    mainDir <- "/media/tiphaine/73134CEF234D24BB/These/PostDoc/Analyses"#path do the main folder without "/" to the end

    pot500 <- spatstat::disc(radius=12.62,centre=c(0,0),poly=256)
    pot1000 <- spatstat::disc(radius=17.846,centre=c(0,0),poly=256)
    curve_set_pooled_all <- NULL
      #  if(method=="size500") {
      #      pattern_to_murgea <- pattern_to_murge_500
      #      }

    for (i in 1:length(names(pattern_to_murge))) {
       country_name <- names(pattern_to_murge[i]) #by country

       full_name_stand <- gsub("([^_]*)_[^_]*.*", "\\1", names(pattern_to_murge[[i]]))
        unique_stand_name <- unique(full_name_stand)

   #     if(method=="remove") {

            random_labelling_result_country_select <-  pattern_to_murge[[i]]
    #        }

#        if(method=="size500") {
 #           random_labelling_result_country_select <-pattern_to_murge_500[[i]]
 #           }


        data_plot_selected <- data_plot[as.character(data_plot$Site_and_Plot_ID)%in%names(pattern_to_murge[[i]]),] # select the data plot to have the plot size using the names of the selected stands
                                        #   sapply(tapply( data_plot_selected$Plot_Size,as.character(data_plot_selected$Site_ID),unique),length)>1 #compute if there is more than 1 plot size per stands (output TRUE ofr more than 1, else FALSE)
                                        #   tapply( as.factor(data_plot_selected$Plot_Size),as.character(data_plot_selected$Site_ID),summary)[1]#compute the count of stands by plot size 500 and 1000
                                        #
                                        #   tapply( as.factor(data_plot_selected$Plot_Size),as.character(data_plot_selected$Site_ID),summary)%in% sapply(tapply( as.factor(data_plot_selected$Plot_Size),as.character(data_plot_selected$Site_ID),summary),max)#compute the count of stands by plot size 500 and 1000

        table_plot_size_Log <- NULL
        curve_set_pooled_per_country <- NULL





 library(doMC)#only for Linux,used doSNOW for Windows and uncomment the above lines
  #      clusterN <-  detectCores()-1  ### choose number of nodes to add to cluster in doSNOW ; Compute the number of core available from the computer less 1
  #      cl = makeCluster(clusterN, rscript="Rscript.exe", type='SOCK')  ###Cluster for doSNOW
        core <- detectCores()-1#Compute the number of core available from the computer less 1
        if(is.na(core)) {
            core <- 8-1}
        registerDoMC(core)#reg


         master_list <- NULL
         curve_set_pooled_per_country <- foreach(y=1:length( unique_stand_name),.packages=c("spatstat","spptest","marksummary","data.table"),.multicombine=TRUE) %dopar% {


   #     for (y in 1:length( unique_stand_name)) { #by stand



             select_in_rank_result <- which(full_name_stand==unique_stand_name[y])
            rank_result_stand <- random_labelling_result_country_select[select_in_rank_result]
                                        #names(rank_result_stand)
            names_rank_result_selected <- names( rank_result_stand )


            data_tree_selected <- data[as.character(data$Site_and_Plot_ID)%in% names_rank_result_selected,]


              data_tree_selected$distance <- sqrt ( data_tree_selected$x_m^2+  data_tree_selected$y_m^2)





            data_tree_selected_stands_normal <-   data_tree_selected[data_tree_selected$distance<=17.846]#selecte tree of plot size 1000 within a plot size 500 radius of 12.62m from the center

            number_tree_per_for_each_stands_normal<- tapply( data_tree_selected_stands_normal$Global_ID,as.character( data_tree_selected_stands_normal$Site_and_Plot_ID)   ,length)

            number_of_trees_and_plot_size_normal <- rbind( number_tree_per_for_each_stands_normal, data_plot_selected[select_in_rank_result]$plotsize )

            number_of_trees_and_plot_size_normal <- as.data.table(t(number_of_trees_and_plot_size_normal))
            colnames(number_of_trees_and_plot_size_normal) <- c("Number_of_trees","plotsize")
            rownames(number_of_trees_and_plot_size_normal) <-  names( rank_result_stand )
            resume_number_of_tree_per_plot_size_normal <- tapply(number_of_trees_and_plot_size_normal$Number_of_trees,as.character(number_of_trees_and_plot_size_normal$plotsize),sum)#compute the number of trees per plot size for the stand "y"

             if ( length(unique(data_plot_selected$plotsize))>2) {
                 resume_number_of_tree_per_plot_size_500 <- NA
                 } else {

            data_tree_selected_stands_500 <-   data_tree_selected[data_tree_selected$distance<=12.62]#selecte tree of plot size 1000 within a plot size 500 radius of 12.62m from the center
            number_tree_per_for_each_stands_500 <- tapply(data_tree_selected_stands_500$Global_ID,as.character(data_tree_selected_stands_500$Site_and_Plot_ID)   ,length)


                     number_of_trees_and_plot_size_500 <- rbind( number_tree_per_for_each_stands_500, gsub("1000", "1000_resize_500",  data_plot_selected[select_in_rank_result]$plotsize))
            assign("last.warning", NULL, envir = baseenv())#remove warnings() from rbind

            number_of_trees_and_plot_size_500 <- as.data.table(t(number_of_trees_and_plot_size_500))
            colnames(number_of_trees_and_plot_size_500) <- c("Number_of_trees","plotsize")
            rownames(number_of_trees_and_plot_size_500) <-  names( rank_result_stand )
            resume_number_of_tree_per_plot_size_500 <- tapply(as.numeric(number_of_trees_and_plot_size_500$Number_of_trees),as.character(number_of_trees_and_plot_size_500$plotsize),sum)#compute the number of trees per plot size for the stand "y"
            resume_number_of_tree_per_plot_size_500 <- resume_number_of_tree_per_plot_size_500[1]
                 }


            inter <- tapply( factor(data_plot_selected$plotsize,levels= as.character( unique(data_plot$plotsize))),as.character(data_plot_selected$Site_ID),summary)[y]
                                        # names(inter)
                                        # names(inter[[1]])
                                        # as.numeric(inter[[1]])
            data_plot_size_stand <- cbind(country_name,
                                          names(inter),
                                          names(inter[[1]]),
                                          as.numeric(inter[[1]]))
            colnames( data_plot_size_stand ) <- c("Country","Stand_name","plotsize","Number_of_plots")
            data_plot_size_stand <- as.data.table(data_plot_size_stand)

            data_plot_size_stand_normal <- cbind(data_plot_size_stand[ as.character(plotsize) %in% names(resume_number_of_tree_per_plot_size_normal),],Number_of_trees=resume_number_of_tree_per_plot_size_normal)#add the number of trees per plot size for the selected stand 'y' and remove Plot__size with 0 trees
            #compute number of trees for resize 500 plots

            if(dim( data_plot_size_stand_normal )[1] >= 2) { #run the 1000 resize 500 only for stands with more than 1 type size
                plot_1000_resize_500 <-  data_plot_size_stand_normal[plotsize==1000]
                plot_1000_resize_500 <- as.data.table(plot_1000_resize_500)
                plot_1000_resize_500 <- plot_1000_resize_500[plotsize == plot_1000_resize_500[,plotsize],plotsize := names( resume_number_of_tree_per_plot_size_500)]
                plot_1000_resize_500 <- plot_1000_resize_500[Number_of_trees==plot_1000_resize_500[,Number_of_trees],Number_of_trees := resume_number_of_tree_per_plot_size_500[[1]]][]
                                        #end computing resize of 500 plots
                data_plot_size_stand <- rbind( plot_1000_resize_500, data_plot_size_stand_normal )#combine normal and resize information for LOG
            } else {
                data_plot_size_stand <-   data_plot_size_stand_normal
            }

            table_plot_size_Log <- rbind(  table_plot_size_Log, data_plot_size_stand)#use for LOG to know the number of trees per stands per plot size class


            if(sapply(tapply( data_plot_selected$plotsize,as.character(data_plot_selected$Site_ID),unique),length)[[y]]==1) { #if result is TRUE then compute plot size normally otherwise need to choose the method

                r <-  rank_result_stand[[1]]$r
                pooled_sims <- matrix(nrow=nrow(rank_result_stand[[1]]$sim_m), ncol=ncol(rank_result_stand[[1]]$sim_m))
                obs_df <- NULL
                pooled_sim_temp <- NULL
                number_sim <- dim(rank_result_stand[[1]]$sim_m)[2]
                for (z in  names(rank_result_stand)) {

                                        #   if (data_plot[data_plot$Site_and_Plot_ID==z,]$Plot_Size[1]==500) {
                                        #       r <- seq(0,max(pot500$xrange)/2,length.out=513)# CAUTIION HERE. Need a plot with 500m2
                                        #   }
                                        #
                                        #   if(data_plot[data_plot$Site_and_Plot_ID==z,]$Plot_Size[1]==1000) {
                                        #       r <- seq(0,max(pot1000$xrange)/2,length.out=513)# CAUTIION HERE. Need a plot with 1000m2
                                        #   }

                    data_curve_rank <- rank_result_stand[[z]]$obs
                    obs_df <- cbind(obs_df, data_curve_rank)
                                        #  temp_sim_m <-rank_result_stand[[z]]$sim_m
                                        #  colnames( temp_sim_m) <- rep(z,dim(rank_result_stand[[z]]$sim_m)[2])
                                        #  pooled_sim_temp <-  cbind(pooled_sim_temp, temp_sim_m)
                }
                pooled_obs <- apply(obs_df, MARGIN=1, FUN=mean)

                for(w in 1:ncol(pooled_sims)) {
                    table_sim_m <- NULL
                    for (zz in  names(rank_result_stand)) {
                        temp_sim_m <-  rank_result_stand[[zz]]$sim_m[,w]
                        table_sim_m <- cbind( table_sim_m, temp_sim_m)
                    }
                    pooled_sims[,w] <- apply( table_sim_m, MARGIN=1, FUN=mean)
                }
                                        # Collect the above components together to a curve_set object
                                        # (with is_residual=TRUE, we tell that the functions are L_mm(r) - L(r))
                curve_set_pooled <- create_curve_set(list(r = r, obs = pooled_obs, sim_m = pooled_sims, is_residual=TRUE))


            }  else { #end if TRUE with only one plot size for one stands

          #      if(method=="remove") {

                   # stand_names_select <-  data_plot_selected[ data_plot_selected$Site_and_Plot_ID==names(rank_result_stand),]


                    stand_names_select <-  data_plot_selected[ data_plot_selected$Site_and_Plot_ID %in% names(rank_result_stand),]


                    stand_selected_diff_1000 <-  stand_names_select[ stand_names_select$plotsize!=1000]
                    names_stand_selected_diff_1000 <- as.character(  stand_selected_diff_1000[,(Site_and_Plot_ID)])

                    selected_diff_1000_randomlabelling_country <- rank_result_stand[  names( rank_result_stand) != names_stand_selected_diff_1000]

                    r <-   selected_diff_1000_randomlabelling_country[[1]]$r
                    pooled_sims <- matrix(nrow=nrow(selected_diff_1000_randomlabelling_country [[1]]$sim_m), ncol=ncol(selected_diff_1000_randomlabelling_country [[1]]$sim_m))
                    obs_df <- NULL
                    pooled_sim_temp <- NULL
                    number_sim <- dim(selected_diff_1000_randomlabelling_country [[1]]$sim_m)[2]
                    for (z in  names(selected_diff_1000_randomlabelling_country )) {
                        data_curve_rank <- selected_diff_1000_randomlabelling_country [[z]]$obs
                        obs_df <- cbind(obs_df, data_curve_rank)
                                        #  temp_sim_m <-rank_result_stand[[z]]$sim_m
                                        #  colnames( temp_sim_m) <- rep(z,dim(rank_result_stand[[z]]$sim_m)[2])
                                        #  pooled_sim_temp <-  cbind(pooled_sim_temp, temp_sim_m)
                    }
                    pooled_obs <- apply(obs_df, MARGIN=1, FUN=mean)

                    for(w in 1:ncol(pooled_sims)) {
                        table_sim_m <- NULL
                        for (zz in  names(selected_diff_1000_randomlabelling_country)) {
                            temp_sim_m <-  selected_diff_1000_randomlabelling_country[[zz]]$sim_m[,w]
                            table_sim_m <- cbind( table_sim_m, temp_sim_m)
                        }
                        pooled_sims[,w] <- apply( table_sim_m, MARGIN=1, FUN=mean)
                    }
                                        # Collect the above components together to a curve_set object
                                        # (with is_residual=TRUE, we tell that the functions are L_mm(r) - L(r))
                    curve_set_pooled <- create_curve_set(list(r = r, obs = pooled_obs, sim_m = pooled_sims, is_residual=TRUE))
  #              }#stop  IF with REMOVE method

        }#stop ELSE when there is 500 and 1000 plots size
            master_list <-  list(curve_set_pooled,  data.table(table_plot_size_Log))
             return(  master_list)
         }# end of y    #end loop by stand END METHOD

        # for (t in 1:length( unique_stand_name)) {


      #  curve_set_pooled_per_country



       # write.table(  table_plot_size_Log, paste0("Replicated Logs",".txt"),sep="\t",row.names = FALSE) #write a file with the plot name where there is point out of the window

        names(curve_set_pooled_per_country) <-  unique_stand_name

        curve_set_pooled_all <-  append(curve_set_pooled_all,list(lapply(curve_set_pooled_per_country, `[[`, 1))) #Select all first object in the list which are the pooled function


        #Print log table per country


       if (file.exists(file.path(mainDir,Status,mark,method,"replicate", country_name ))){
                setwd(file.path(mainDir,Status,mark,method,"replicate", country_name ))
            } else {
                dir.create(file.path(mainDir,Status,mark,method,"replicate", country_name ), showWarnings = FALSE,recursive=TRUE)
                setwd(file.path(mainDir,Status,mark,method,"replicate", country_name ))
            }


       table_plot_size_Log <- do.call(rbind.data.frame, lapply(curve_set_pooled_per_country, `[[`, 2)) #select all second object of the list which are the Log file per country

        write.table( table_plot_size_Log, paste0("Decision Logs ",method, " ", country_name,".txt"),sep="\t",row.names = FALSE)
       setwd(mainDir)

    } #end  loop country

    names(curve_set_pooled_all) <- names(pattern_to_murge)

    return(curve_set_pooled_all)

}#end function
despresT/replicate.patterns documentation built on May 15, 2019, 5:05 a.m.