Map of thickness and current intensity

Bin the data and fit the shape does not work. Using 2d gaussian kernel density estimation instead.

def bin_df(df: pl.DataFrame, col_to_bin, bins=10):
    binned_col = f"{col_to_bin}_bin"
    
    return (
        df.with_columns(
            pl.col(col_to_bin).qcut(bins).alias(binned_col),
        )
        .group_by(binned_col)
        .agg(cs.numeric().median(), pl.count().alias("bin_count"))
        .drop(binned_col)
    )

col_to_bin="L_mn_norm_log"
# col_to_bin="j0_norm_log"

all_events_l1_L_binned = pl.concat(
    [
        data.pipe(bin_df, col_to_bin=col_to_bin, bins=64).with_columns(sat= pl.lit(name))
        for name, data in all_events_l1.group_by("sat")
    ]
)

jno_events_l1_L_binned = pl.concat(
    [
        data.pipe(bin_df, col_to_bin=col_to_bin, bins=64).with_columns(sat= pl.lit(name))
        for name, data in jno_candidates_l1.group_by("r_bin")
    ]
)
%R -i all_events_l1_L_binned -c conv_pl

df = all_events_l1.filter(pl.col('L_mn_norm_log').is_not_nan())
# df = jno_candidates_l1.filter(pl.col('L_mn_norm_log').is_not_nan())
%R -i df -c conv_pl

import lmfit
import numpy as np

def gaussian2d(x, y, amplitude=1., centerx=0., centery=0., sigmax=1., sigmay=1.,
                 rotation=0, A0=0.):
    """Return a two dimensional lorentzian.

    The maximum of the peak occurs at ``centerx`` and ``centery``
    with widths ``sigmax`` and ``sigmay`` in the x and y directions
    respectively. The peak can be rotated by choosing the value of ``rotation``
    in radians.
    """
    xp = (x - centerx)*np.cos(rotation) - (y - centery)*np.sin(rotation)
    yp = (x - centerx)*np.sin(rotation) + (y - centery)*np.cos(rotation)
    R = (xp/sigmax)**2 + (yp/sigmay)**2

    return A0 + amplitude * np.exp(-R/2)

model = lmfit.Model(gaussian2d, independent_vars=['x', 'y'])
# params = model.make_params(amplitude=10, centerx=x[np.argmax(z)], centery=y[np.argmax(z)])

%%R
fit_gaussian_2D_pdf(df)

%%R
library(purrr)
library(tidyr)
fit_gaussian_2D_pdf <- function(data) {
  kde_result <- MASS::kde2d(data$L_mn_norm_log, data$j0_norm_log)
  x_values <- rep(kde_result$x, each = length(kde_result$y))
  y_values <- rep(kde_result$y, length(kde_result$x))
  response <- as.vector(kde_result$z)

  density <- data.frame(X_values = x_values, Y_values = y_values, response = response)
  model <- fit_gaussian_2D(density)
  return(model)
}

results <- df %>% 
  # group_by(r_bin) %>% 
  nest() %>% 
  mutate(fitted = map(data, fit_gaussian_2D_pdf))

%%R -w 1000 -h 500
# Creating a list of layers for the binned data
# model <- lm(j0_norm_log ~ L_mn_norm_log, data = all_events_l1_L_binned)
# slope <- coef(model)[2]

binned_layer <- list(
  geom_line(data = all_events_l1_L_binned, color = 'blue'),
  geom_point(data = all_events_l1_L_binned, color = 'blue'), 
  geom_smooth(data = all_events_l1_L_binned, method = "glm", color = 'red')
)


# Plot creation
p <- ggplot(mapping = aes(x = L_mn_norm_log, y = j0_norm_log)) +
  geom_density_2d(data = all_events_l1) +
  # stat_density_2d(data = all_events_l1, aes(fill = after_stat(density)), geom = "raster", contour = FALSE) +
  binned_layer +
  facet_wrap(~ sat, scales = "free")

  
# Print the plot
print(p)

%%R
p <- ggplot() +
  geom_point(data = all_events_l1, aes(x = L_mn_norm_log, y = j0_norm_log)) +
  binned_layer +
  facet_wrap(~ sat, scales = "free")

print(p)

%%R
# Fit a linear model to the log-transformed data
lm_fit <- lm(j0_norm_log ~ L_mn_norm_log, data = all_events_l1)

# Extract the coefficients
intercept <- coef(lm_fit)[1]
slope <- coef(lm_fit)[2]

# Create a scatter plot with the log-log transformation
p <- ggplot(all_events_l1, aes(x = L_mn_norm_log, y = j0_norm_log)) +
  geom_point() + # Add the scatter points
  geom_abline(intercept = intercept, slope = slope, color = 'blue', size = 1) + # Add the fitted line
  facet_wrap(~ sat, scales = "free") + # Facet by 'sat'
  labs(x = "Log10(L_mn_norm)", y = "Log10(j0_norm)") # Label axes

print(p)

%%R
# Plot creation
p <- ggplot(all_events_l1_L_binned, aes(x = L_mn_norm_log, y = j0_norm_log)) +
    geom_line(color = 'blue') +
    geom_point(color = 'blue') +
    geom_smooth(method = "glm", color = 'red') +
    facet_wrap(~ sat, scales = "free") +
    stat_regline_equation()

  
# Print the plot
print(p)

%%R -i jno_events_l1_L_binned -c conv_pl
# Plot creation
p <- ggplot(jno_events_l1_L_binned, aes(x = L_mn_norm_log, y = j0_norm_log)) +
    geom_line(color = 'blue') +
    geom_point(color = 'blue') +
    geom_smooth(method = "glm", color = 'red') +
    facet_wrap(~ r_bin, scales = "free") +
    stat_regline_equation()

  
# Print the plot
print(p)

%%R -i jno_candidates_l1 -c conv_pl

p <- ggplot(jno_candidates_l1, aes(x = L_mn_norm, y = j0_norm)) +
  stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
  facet_wrap(~ r_bin, nrow = length(unique(jno_candidates_l1$r_bin))) +
  scale_x_log10() + 
  scale_y_log10() +
  labs(fill = "Density")


print(p)