// Modelo de Trait-based competition model siguiendo a Kunstler et al 2016 
data {
  int<lower=0> N; // Número total de observaciones
  vector[N] Biomasa; // Biomasa observada
  vector[N] tf; // Rasgo del especie focal
  vector[N] tc; // Rasgos del especie competidor
  vector[N] Bc; // Biomasa del competidor
  int<lower=1> N_species; // Número de especies diferentes
  int<lower=1, upper=N_species> species[N]; // Índice de especie para cada observación
  int<lower=0, upper=1> C[N]; // Indicador de competencia intra (1) vs. interespecífica (0) para cada observación
  int<lower=0, upper=1> water_status[N]; // Estado hídrico para cada observación
}

parameters {
  // Parámetros para condiciones secas
  real m0_dry;
  real m1_dry;
  vector[N_species] epsilon_Bmax_f_dry;
  vector[N_species] epsilon_alpha0intra_dry;
  vector[N_species] epsilon_alpha0inter_dry;
  real alpha0intra_species_dry;
  real alpha0inter_species_dry;
  real alpha_t_dry; // Tolerancia a la competencia por el rasgo tf en condiciones secas
  real alpha_e_dry; // Efecto competitivo por el rasgo tc en condiciones secas
  real alpha_d_dry; // Efecto de la disimilitud de rasgos en condiciones secas

  // Parámetros para condiciones húmedas
  real m0_wet;
  real m1_wet;
  vector[N_species] epsilon_Bmax_f_wet;
  vector[N_species] epsilon_alpha0intra_wet;
  vector[N_species] epsilon_alpha0inter_wet;
  real alpha0intra_species_wet;
  real alpha0inter_species_wet;
  real alpha_t_wet; // Tolerancia a la competencia por el rasgo tf en condiciones húmedas
  real alpha_e_wet; // Efecto competitivo por el rasgo tc en condiciones húmedas
  real alpha_d_wet; // Efecto de la disimilitud de rasgos en condiciones húmedas

  // Parámetros generales (que no varían con el water_status)
  real<lower=0> sigma_species_Bmax;
  real<lower=0> sigma_alpha_intra;
  real<lower=0> sigma_alpha_inter;
  real<lower=0> sigma;
}


model {
  // Distribuciones previas para los parámetros
  m0_dry ~ normal(0, 1);
  m1_dry ~ normal(0, 1);
  m0_wet ~ normal(0, 1);
  m1_wet ~ normal(0, 1);
  sigma_species_Bmax ~ normal(0, 0.5);
  sigma_alpha_intra ~ normal(0, 0.5);
  sigma_alpha_inter ~ normal(0, 0.5);
  sigma ~ normal(0, 1);
  
  // Efectos aleatorios y parámetros específicos de condiciones secas y húmedas
  epsilon_Bmax_f_dry ~ normal(0, sigma_species_Bmax);
  epsilon_Bmax_f_wet ~ normal(0, sigma_species_Bmax);
  epsilon_alpha0intra_dry ~ normal(0, sigma_alpha_intra);
  epsilon_alpha0intra_wet ~ normal(0, sigma_alpha_intra);
  epsilon_alpha0inter_dry ~ normal(0, sigma_alpha_inter);
  epsilon_alpha0inter_wet ~ normal(0, sigma_alpha_inter);
  alpha0intra_species_dry ~ normal(0, 1); 
  alpha0inter_species_dry ~ normal(0, 1);
  alpha0intra_species_wet ~ normal(0, 1);
  alpha0inter_species_wet ~ normal(0, 1);
  alpha_t_dry ~ normal(0, 1);
  alpha_e_dry ~ normal(0, 1);
  alpha_d_dry ~ normal(0, 1);
  alpha_t_wet ~ normal(0, 1);
  alpha_e_wet ~ normal(0, 1);
  alpha_d_wet ~ normal(0, 1);

  // Modelo log-linear de Biomasa ajustado para cada observación
  for (i in 1:N) {
    // Selección de parámetros basados en water_status
    real m0 = water_status[i] == 1 ? m0_wet : m0_dry;
    real m1 = water_status[i] == 1 ? m1_wet : m1_dry;
    real alpha_t = water_status[i] == 1 ? alpha_t_wet : alpha_t_dry;
    real alpha_e = water_status[i] == 1 ? alpha_e_wet : alpha_e_dry;
    real alpha_d = water_status[i] == 1 ? alpha_d_wet : alpha_d_dry;
    vector[N_species] epsilon_Bmax_f = water_status[i] == 1 ? epsilon_Bmax_f_wet : epsilon_Bmax_f_dry;
    vector[N_species] epsilon_alpha0intra = water_status[i] == 1 ? epsilon_alpha0intra_wet : epsilon_alpha0intra_dry;
    vector[N_species] epsilon_alpha0inter = water_status[i] == 1 ? epsilon_alpha0inter_wet : epsilon_alpha0inter_dry;
    real alpha0intra_species = water_status[i] == 1 ? alpha0intra_species_wet : alpha0intra_species_dry;
    real alpha0inter_species = water_status[i] == 1 ? alpha0inter_species_wet : alpha0inter_species_dry;

    // Cálculo de Bmax_i y competition_effect
    real Bmax_i = exp(m0 + m1 * tf[i] + epsilon_Bmax_f[species[i]]);
    real competition_effect = 0;

    if (C[i] == 1) {
      competition_effect += alpha0intra_species + epsilon_alpha0intra[species[i]];
    } else {
      competition_effect += alpha0inter_species + epsilon_alpha0inter[species[i]];
    }

    competition_effect += -alpha_t * tf[i] + alpha_e * tc[i] + alpha_d * fabs(tc[i] - tf[i]);

    // Aplicación del modelo final para cada observación
    Biomasa[i] ~ normal(Bmax_i * exp(-competition_effect * Bc[i]), sigma);
  }
}