knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "curv-meas_figures/" )
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.