library(readr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
bw_data <- read.csv("~/Desktop/FRESCA Data Ofav/Ofav BW ORGANIZED.csv")
sa_data <- read.csv("~/Desktop/FRESCA Data Ofav/Ofav_SA.csv")
master_data <- read.csv("~/Desktop/FRESCA Data Ofav/OFAV MASTER.csv")
DATA <- left_join(bw_data, master_data, by="ID")    # keep all data in bw_data but add columns from master_data
DATA <- left_join(DATA, sa_data, by="ID")           # keep all data in DATA but add columns from sa_data

DATA$pH <- as.factor(DATA$pH)
DATA$O2_umol_L <- as.factor(DATA$O2_umol_L)
DATA$Date <- as.Date(DATA$Date)
# 1 # calculate density (long ass equation) #
DATA$calcDensity <- (1000*(1-(DATA$Temp_C+288.9414)/(508929.2*(DATA$Temp_C+68.12963))*(DATA$Temp_C-3.9863)^2)) + (0.824493 - 0.0040899*DATA$Temp_C + 0.000076438*DATA$Temp_C^2 -0.00000082467*DATA$Temp_C^3 + 0.0000000053675*DATA$Temp_C^4)*DATA$Salinity +(-0.005724 + 0.00010227*DATA$Temp_C - 0.0000016546*DATA$Temp_C^2)*DATA$Salinity^(3/2) + 0.00048314*DATA$Salinity^2    

# 2 # calculate density of seawater in g/cm3 (density/1000)
DATA$Density_seawater <- DATA$calcDensity/1000

# 3 # add density of aragonite (g/cm3) column for all (2.93 for all)
DATA$Density_arag <- 2.93

# 4 # calculate (density of seawater)/(density of aragonite)
DATA$Density_sea_arag <- DATA$Density_seawater / DATA$Density_arag

# 5 # calculate 1-[Density (SW/A)]
DATA$Density_diff <- 1-DATA$Density_sea_arag

# 6 # calculate dry weight (wet weight/line above)
DATA$Dry_weight_g <- DATA$Wet_weight_g/DATA$Density_diff
DATA <- DATA %>%
  arrange(ID, Date) %>%    # arrange by coral ID and date
  group_by(ID) %>%         # group by coral ID
  mutate(
    prev_date = lag(Date),
    weight_diff = Dry_weight_g - lag(Dry_weight_g),      # calculate dry weight difference
    time_diff = as.numeric(Date - lag(Date)),            # calculate time difference
    growth_rate = weight_diff / (time_diff)*(SA_cm2*1000),    # calculate GR 
    GR_label = paste0("GR", row_number() - 1)            # label GR based on time points 
  ) %>%
  filter(!is.na(growth_rate))    

DATA <- DATA %>%
  mutate(O2_umol_L = as.numeric(as.character(O2_umol_L)),
         O2_umol_L = ifelse(GR_label %in% c("GR1", "GR2"), 200, O2_umol_L),
         O2_umol_L = as.factor(O2_umol_L))
plot_ly(data = DATA,
        x = ~Date,
        y = ~Dry_weight_g,
        color = ~ID,
        type = 'scatter',
        mode = 'lines+markers',
        hoverinfo = 'text',
        text = ~paste("ID:", ID,
                      "<br>Date:", Date,
                      "<br>Dry Weight:", Dry_weight_g)) %>%
  layout(
    title = "Coral Dry Weight Over Time",
    xaxis = list(title = "Date"),
    yaxis = list(title = "Dry Weight (g)"),
    showlegend = FALSE  # <-- removes the big legend
  )
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
# Filter for corals with weight_diff > 1 or weight_diff < 0
filtered <- DATA[DATA$weight_diff > 1 | DATA$weight_diff < 0, ]

# Optional: Show only the coral_id and weight_diff columns
result <- filtered[, c("ID", "weight_diff", "Date")]      # add any other columns you want to see

# View the result
print(result)
## # A tibble: 93 × 3
## # Groups:   ID [67]
##    ID         weight_diff Date      
##    <chr>            <dbl> <date>    
##  1 Ofav34-804    8.46     2025-06-03
##  2 Ofav34-804   -8.31     2025-06-16
##  3 Ofav34-809   -0.0349   2025-06-16
##  4 Ofav34-812   -0.494    2025-06-03
##  5 Ofav34-814   -0.152    2025-06-03
##  6 Ofav34-818   -0.495    2025-06-03
##  7 Ofav34-819   -0.312    2025-06-16
##  8 Ofav34-820   -0.000162 2025-06-16
##  9 Ofav34-821   -0.132    2025-06-03
## 10 Ofav34-824   -0.00606  2025-06-16
## # ℹ 83 more rows
# Create facet variable
DATA <- DATA %>%
  mutate(facet_group =  as.character(O2_umol_L))

# Reorder facet levels
DATA$facet_group <- factor(DATA$facet_group, 
                                 levels = c("200", "125", "62"))

# Custom labels
facet_labels <- c(
  "200" = "200 μmol L⁻¹",
  "125" = "125 μmol L⁻¹",
  "62"  = "62 μmol L⁻¹")

# Plot - growth rate by individual coral
g <- ggplot(DATA, aes(x = GR_label, y = growth_rate, fill=Genotype,
        text = paste("ID:", ID, 
        "<br>Growth Rate:", growth_rate, 
        "<br>Label:", GR_label))) +
  geom_jitter()+
   facet_wrap(~facet_group, scales = "free_x", nrow = 1,
             labeller = as_labeller(facet_labels))

ggplotly(g, tooltip = "text") # make it interactive