This note contains technical details to the vignette Studying curvature measures of symmetric cones. We will provide the code that generates those plots, which are only included as image files and not created in the markdown file.

Packages used:

library(symconivol)
library(conivol)
library(tidyverse)
library(knitr)
library(Rmisc)
library(rstan)
library(png)
library(latex2exp)

Sampling index constrained eigenvalues: In the following lines we construct Stan models for index constrained eigenvalues, run the models to obtain ten thousand sample points for each model, and save these sample points in correspondingly named files.

warmup <- 1e3
num_samp <- 1e5
fname_model <- "tmp.stan"

# define the number of negative, free, positive eigenvalues
nn_nf_np_list <- tribble(~nn, ~nf, ~np,
                           5,  10,  25,
                           0,  15,  25,
                           2,   0,   8,
                           4,   0,   6,
                           8,   0,  32,
                          16,   0,  24,
                          35,   5,   0,
                          20,   0,   0,
                           0,  10,   0,
                           0,  40,   0,
                           0,   0,  40)

for (i in 1:length(nn_nf_np_list)) {
    nn <- nn_nf_np_list$nn[i]
    nf <- nn_nf_np_list$nf[i]
    np <- nn_nf_np_list$np[i]
    
    # construct the Stan model file
    constr_eigval(pos=(np>0), free=(nf>0), neg=(nn>0), filename=fname_model)
    for (beta in c(1,2,4)) {
        fname_evals <- str_c("evals_beta_nn_nf_np_N=",beta,"_",nn,"_",nf,
                             "_",np,"_",num_samp,"_wmup=",warmup,".RData")
        data = list(beta=beta)
        if (np>0) data$np <- np
        if (nf>0) data$nf <- nf
        if (nn>0) data$nn <- nn
        
        # run the Stan model
        stan_samp <- stan( file=fname_model, data=data, chains=1, warmup=warmup,
                           iter=num_samp+warmup, cores=2, refresh=1e4 )
        
        # extract the sample points
        samp <- list()
        if (np>0) samp$ep <- rstan::extract(stan_samp)$ep
        if (nf>0) samp$ef <- rstan::extract(stan_samp)$ef
        if (nn>0) samp$en <- rstan::extract(stan_samp)$en
        
        save(samp, file=fname_evals)
    }
    file.remove(fname_model)
}

Plots of empirical distribution of index constrained eigenvalues: The following lines produce the plots of the empirical distribution of the index constrained eigenvalues sampled above.

colors <- hcl(h = seq(15, 375, length=4), l = 65, c = 100)[c(2,3,1)]
ymax <- c(0.18,0.16,0.16,0.24,0.6,0.5,0.083,0,rep(0.9,2),rep(1.3,2))
I <- c(1,2,6,7,8,10,11)

for (i in I) {
    nn <- nn_nf_np_list[[i,1]]
    nf <- nn_nf_np_list[[i,2]]
    np <- nn_nf_np_list[[i,3]]
    n <- nn+nf+np
    plotsPB <- list()
    j <- 0
    for (beta in c(1,2,4)) {
        j <- j+1
        fname_evals <- str_c("evals_beta_nn_nf_np_N=",beta,"_",nn,"_",nf,
                             "_",np,"_",num_samp,"_wmup=",warmup,".RData")
        load(file=fname_evals)
        tmp <- tibble()
        if (np>0) tmp <- bind_rows(tmp,samp$ep %>% as_tibble() %>%
                                       gather(factor_key=TRUE) %>% add_column(type="0"))
        if (nf>0) tmp <- bind_rows(tmp,samp$ef %>% as_tibble() %>%
                                       gather(factor_key=TRUE) %>% add_column(type="1"))
        if (nn>0) tmp <- bind_rows(tmp,samp$en %>% as_tibble() %>%
                                       gather(factor_key=TRUE) %>% add_column(type="2"))
        this_colors <- colors[c(np>0,nf>0,nn>0)]
        plotsPB[[j]] <- ggplot() +
            geom_density(data=tmp, aes(x=value, y=..count.., group=type, color=type,
                                       fill=type), alpha=0.5, bw=0.03) +
            theme_bw() + xlab(TeX(sprintf("$\\beta = %d$", beta))) +
            theme(axis.title.y=element_blank(),
                  axis.text.y=element_blank(), axis.ticks.y=element_blank(),
                  panel.border = element_blank(), panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  axis.line = element_line(colour = "black")) +
            coord_cartesian(ylim = c(0, ymax[i]/sqrt(beta)*num_samp*n)) +
            theme(legend.position="none") +
            scale_fill_manual(values=this_colors) + scale_color_manual(values=this_colors)
    }
    plotFileName <-str_c("nn_nf_np=",nn,"_",nf,"_",np,".png")
    ggsave(plotFileName, plot=multiplot(plotlist = plotsPB, cols = 3), device="png",
           width=178, height=80, units="mm", dpi=300)
}

Plots of empirical distribution of free eigenvalues including limit distribution: The following lines produce the plots of the empirical distribution of unconstrained eigenvalues sampled above including its limiting semicircle distribution.

t <- seq(-pi/2,pi/2,length=1e3)
C <- tibble(x=sqrt(2)*sin(t),y=sqrt(2)/pi*cos(t))
for (i in c(9,10)) {
    n <- format[[i,2]]
    plotsPB <- list()
    j <- 0
    for (beta in c(1,2,4)) {
        j <- j+1
        fname_evals <- str_c("beta_nn_nf_np=",beta,"_0_",n,"_0_wmup=",warmup,".RData")
        load(file=fname_evals)
        tmp <- (1/sqrt(beta*n)*samp$ef) %>% as_tibble() %>%
                gather(factor_key=TRUE) %>% add_column(type="1")
        plotsPB[[j]] <- ggplot() +
            geom_density(data=tmp, aes(x=value, group=type, color=type, fill=type),
                         alpha=0.5, bw=0.03) +
            theme_bw() + xlab(TeX(sprintf("$\\beta = %d$", beta))) +
            theme(axis.title.y=element_blank(),
                  axis.text.y=element_blank(), axis.ticks.y=element_blank(),
                  panel.border = element_blank(), panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  axis.line = element_line(colour = "black")) +
            coord_cartesian(ylim = c(0, 1.2*sqrt(2)/pi)) +
            theme(legend.position="none") + scale_color_manual(values=colors[2]) +
            scale_fill_manual(values=colors[2]) +
            geom_line(data=C, aes(x=x,y=y), color="red",size=2,alpha=0.3)
    }
    plotFileName <-str_c("semic_n=",n,".png")
    ggsave(plotFileName, plot=multiplot(plotlist = plotsPB, cols = 3), device="png",
           width=178, height=80, units="mm", dpi=300)
}

Plots of empirical distribution of index constrained eigenvalues including limit distribution: The following lines produce the plots of the empirical distribution of the index constrained eigenvalues sampled above.

c_a <- function(a) {
    integrand <- function(x) return(sqrt((1-x)/x*(x^2+x+(a-1)/a^2)))
    return(2/pi/(1-(a-1)/a^2)*integrate(integrand, lower=0, upper=1)$value)
}
L_a <- function(a) return(a*sqrt(2/(a^2-a+1)))
bnd_a <- function(a) {
    L <- L_a(a)
    return(list(neg_low=-L/a, neg_upp=-L*(1-1/a), pos_low=0, pos_upp=L))
}
rho_a <- function(a,moment=0) {
    L <- L_a(a)
    return( function(x) x^moment/pi*sqrt( (L-x)/x * (x+L/a) * (x+(1-1/a)*L) ) )
}
a <- seq(1,2,length.out=1e3)
c <- sapply(a, c_a)

for (i in 3:6) {
    nn <- format[[i,1]]
    np <- format[[i,3]]
    n <- nn+np
    c0 <- np/n
    ind <- which.min((c-c0)^2)
    a0 <- a[ind]
    bnd <- bnd_a(a0)
    x_neg <- seq(bnd$neg_low,bnd$neg_upp,length.out=1e3)
    x_pos <- seq(bnd$pos_low,bnd$pos_upp,length.out=1e3)
    data_lim_neg <- tibble(x=x_neg, y=sapply(x_neg, rho_a(a0))*n*num_samp)
    data_lim_pos <- tibble(x=x_pos, y=sapply(x_pos, rho_a(a0))*n*num_samp)
    plotsPB <- list()
    j <- 0
    for (beta in c(1,2,4)) {
        j <- j+1
        load(file=str_c("dbeta_nn_nf_np=",beta,"_",nn,"_0_",np,"_wmup=",warmup,".RData"))
        data_pos <- (1/sqrt(beta*n)*samp$ep) %>% as_tibble() %>% gather()
        data_neg <- (1/sqrt(beta*n)*samp$en) %>% as_tibble() %>% gather()
        plotsPB[[j]] <- ggplot() +
            geom_density(data=data_pos, aes(x=value, y=..count..), color=colors[1],
                         fill=colors[1], alpha=0.5, bw=0.01) +
            geom_density(data=data_neg, aes(x=value, y=..count..), color=colors[3],
                         fill=colors[3], alpha=0.5, bw=0.01) +
            geom_line(data=data_lim_neg, aes(x=x,y=y),color=colors[2],size=1,alpha=0.7) +
            geom_line(data=data_lim_pos, aes(x=x,y=y),color=colors[2],size=1,alpha=0.7) +
            theme_bw() + xlab(TeX(sprintf("$\\beta = %d$", beta))) +
            theme(axis.title.y=element_blank(),
                  axis.text.y=element_blank(), axis.ticks.y=element_blank(),
                  panel.border = element_blank(), panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  axis.line = element_line(colour = "black")) +
            coord_cartesian(ylim = c(0, ymax[i]*n*num_samp)) +
            theme(legend.position="none")
    }
    plotFileName <-str_c("indlim_n_cPerc=",n,"_",round(c0*100),".png")
    ggsave(plotFileName, plot=multiplot(plotlist = plotsPB, cols = 3), device="png",
           width=178, height=80, units="mm", dpi=300)
}

Reconstructing index normalized curvature measures: We find the limit curves of the dimension normalized curvature measures for different values of the constant \(C\).

L <- leigh(1e5)
R <- rate(1e5)
Nk <- 1e4
kappa <- (0:Nk)/Nk
Nr <- 1e4
rho <- (0:Nr)/Nr
Lkappa <- L$lkup_kappa(rho)
Rkappa <- R$lkup_R(rho)

C_L <- c(0.15,0.18,0.2,0.22,0.25)
hues = seq(15, 375, length = length(C_L) + 1)
mycolors <- hcl(h = hues, l = 65, c = 100)[1:length(C_L)]
data_rhomin <- tibble()
for (C in C_L) {
    rhomin <- sapply(kappa, function(x)
        return( max(1-sqrt(1-x), min(sqrt(x),rho[which.min(
            (x-Lkappa)^2 / (1-abs(2*rho-1)) + 4*C*Rkappa )] ) ) ) )
    data_rhomin <- bind_rows(data_rhomin,
                             tibble(kappa=kappa,
                                    rho=rhomin,
                                    C=as.factor(C)))
}
y_pat <- seq(0,1,length.out=1e2)
pat_bd_low <- tibble(patX=y_pat^2, patY=y_pat)
pat_bd_upp <- tibble(patX=1-y_pat^2, patY=1-y_pat)
P <- ggplot() + theme_bw() +
    geom_line(data=pat_bd_low, aes(patY, patX), linetype=2, alpha=0.5) +
    geom_line(data=pat_bd_upp, aes(patY, patX), linetype=2, alpha=0.5) +
    geom_line(data=data_rhomin, aes(x=rho,y=kappa,group=C,color=C)) +
    scale_colour_manual(values=mycolors) +
    xlab(TeX(sprintf("$\\rho$"))) + ylab(TeX(sprintf("$\\kappa$")))
ggsave("lim_PhiDim.png", plot=P, device="png",
       width=178, height=120, units="mm", dpi=300)