FittedParamsVisExploration.Rmd

---
title: "Fitted PBE parameters exploration"
author: "Dmitry A. Grechka"
date: "July 6, 2016"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

```{r converting data}
library(ggplot2)


prepareVisData <- function(id,lglk,y0,y2,slope_k,bg) {
  x2<-20.0
  m <- matrix(c( 0.0, 0.0, 1.0,
                 x2*x2, x2, 1.0,
                 2.0*x2, 1.0, 0.0),
                nrow=3,byrow = T)
  r <- c(y0, y2, 0.0)
  s <- solve(m,r)
  slope <- s[2]*slope_k
  print(paste0("maximum slope is ",s[2]))
  print(paste0("effective slope is ",slope))
  
  m2 <- matrix(c( 0.0, 0.0, 1.0,
                 x2*x2, x2, 1.0,
                 0.0, 1.0, 0.0),
                nrow=3,byrow = T)
  r2 <- c(y0, y2, slope)
  s2 <- solve(m2,r2)
  
  func <- function(x) {
    return(s2[1]*x*x + s2[2]*x+s2[3])
  }
  
return(
  list(pred=func,bg=bg,lglk_str=paste0('seed ',id,' (lglk ',lglk,')'))
);
}

to_km_per_s <- function(v) {return(v/(1000.0 * 60.0 / 10.0**8.0));}

#docker 0.3 seed 3
#2016-07-11T16:08:11.935964300Z 
#iteration:1105 log-likelihood:11923.2 improvement:0.0278582
#Vfunc nodes: (0;0.180612) (20;0.597612) slope: 0.326881
#Dfunc nodes: (0;0.870576) (20;0.944965) slope: 0.697326
#background wind V:0.259682 D:4.32253
#Noise sigma: 0.0514491
#Earth distance: 1285.718414

vis3d <- prepareVisData(
  3,
  11923.2,
  0.870576,
  0.944965,
  0.697326,
  4.32253)

vis3v <- prepareVisData(
  3,
  11923.2,
  to_km_per_s(0.180612),
  to_km_per_s(0.597612),
  0.326881,
  to_km_per_s(0.259682))

#docker 0.3 seed 8
#2016-07-14T14:55:12.994253100Z 
#iteration:924 log-likelihood:11915 improvement:3.75523
#Vfunc nodes: (0;0.181729) (20;0.867255) slope: 0.0855468
#Dfunc nodes: (0;0.823733) (20;0.834859) slope: 0.55678
#background wind V:0.265428 D:5.82408
#Noise sigma: 0.0520493
#Earth distance: 1230.243211

vis8d <- prepareVisData(
  8,
  11915,
  0.823733,
  0.834859,
  0.55678,
  5.82408)

vis8v <- prepareVisData(
  8,
  11915,
  to_km_per_s(0.181729),
  to_km_per_s(0.867255),
  0.0855468,
  to_km_per_s(0.265428))

#docker 0.3 seed 5
#2016-07-12T05:31:34.962763100Z 
#iteration:583 log-likelihood:11870.8 improvement:88.4558
#Vfunc nodes: (0;0.243192) (20;0.429677) slope: 0.0448303
#Dfunc nodes: (0;0.294873) (20;0.310647) slope: 0.214628
#background wind V:5.20376e-124 D:5.38888e-170
#Noise sigma: 0.0535385
#Earth distance: 1292.090769

vis5d <- prepareVisData(
  5,
  11870.8,
  0.294873,
  0.310647,
  0.214628,
  5.38888e-170)

vis5v <- prepareVisData(
  5,
  11870.8,
  to_km_per_s(0.243192),
  to_km_per_s(0.429677),
  0.0448303,
  to_km_per_s(5.20376e-124))

#docker 0.3 seed 4
#2016-07-10T21:53:49.075621400Z 
#iteration:821 log-likelihood:11855.5 improvement:3.63798e-12
#Vfunc nodes: (0;0.224727) (20;0.40079) slope: 0.35187
#Dfunc nodes: (0;0.0146691) (20;0.0234736) slope: 0.175416
#background wind V:7.88056e-248 D:8.64625e-39
#Noise sigma: 0.0483257
#SEarth distance: 1386.194369

vis4d <- prepareVisData(
  4,
  11855.5,
  0.0146691,
  0.0234736,
  0.175416,
  8.64625e-39)

vis4v <- prepareVisData(
  4,
  11855.5,
  to_km_per_s(0.224727),
  to_km_per_s(0.40079),
  0.35187,
  to_km_per_s(7.88056e-248))

#docker 0.3 seed 1
#2016-07-12T16:21:59.274401800Z 
#iteration:1510 log-likelihood:11851.6 improvement:8.50073
#Vfunc nodes: (0;0.273536) (20;0.430973) slope: 0.735489
#Dfunc nodes: (0;0.00915377) (20;0.984944) slope: 0.897594
#background wind V:0.23322 D:3.04009
#Noise sigma: 0.0542054
#Earth distance: 1613.963631

vis1d <- prepareVisData(
  1,
  11851.6,
  0.00915377,
  0.984944,
  0.897594,
  3.04009)

vis1v <- prepareVisData(
  1,
  11851.6,
  to_km_per_s(0.273536),
  to_km_per_s(0.430973),
  0.735489,
  to_km_per_s(0.23322))

```

```{r plotting velocities}

linetype_scale <- c('1'='solid','2'='dashed')
linetype_scale_labels <- c('1'='burst func','2'='background level')
color_scale <- c(
        '3'='red',
        '8'='green',
        '5'='blue',
        '1'='purple',
        '4'='black')
color_scale_labels <- c(
        '3'=vis3v$lglk_str,
        '8'=vis8v$lglk_str,
        '5'=vis5v$lglk_str,
        '4'=vis4v$lglk_str,
        '1'=vis1v$lglk_str
        )

a <- ggplot(data.frame(x = c(0, 20.0)))+
  theme_bw()+
  ggtitle("Particle burst velocity as a function of coronal hole area")+
  xlab("coronal holes total area (%)")+
  ylab("particle burst velocity (km/s)")+

  scale_linetype_manual(values= linetype_scale,labels=linetype_scale_labels)+
  scale_color_manual(values= color_scale,labels=color_scale_labels)+
  
  stat_function(aes(x,linetype='1',color='3'),fun=vis3v$pred) +
  geom_hline(aes(yintercept=vis3v$bg,linetype='2',color='3'))+
  
  stat_function(aes(x,linetype='1',color='8'),fun=vis8v$pred) +
  geom_hline(aes(yintercept=vis8v$bg,linetype='2',color='8')) +
  
  stat_function(aes(x,linetype='1',color='5'),fun=vis5v$pred) +
  geom_hline(aes(yintercept=vis5v$bg,linetype='2',color='5')) +
  
  stat_function(aes(x,linetype='1',color='4'),fun=vis4v$pred) +
  geom_hline(aes(yintercept=vis4v$bg,linetype='2',color='4')) +
  
  stat_function(aes(x,linetype='1',color='1'),fun=vis1v$pred) +
  geom_hline(aes(yintercept=vis1v$bg,linetype='2',color='1'))
  
print(a)
```

```{r burst densities}

a = ggplot(data.frame(x = c(0, 20.0)))+
  theme_bw()+
  ggtitle("Particle burst density as a function of coronal hole area")+
  xlab("coronal holes total area (%)")+
  ylab("particle burst density (portion)")+

  scale_linetype_manual(values= linetype_scale,labels=linetype_scale_labels)+
  scale_color_manual(values= color_scale,labels=color_scale_labels) +
  
  stat_function(aes(x,linetype='1',color='3'),fun=vis3d$pred) +
  #geom_hline(aes(yintercept=vis3d$bg,linetype='2',color='3'))+
  
  stat_function(aes(x,linetype='1',color='8'),fun=vis8d$pred) +
  #geom_hline(aes(yintercept=vis8d$bg,linetype='2',color='8')) +
  
  stat_function(aes(x,linetype='1',color='5'),fun=vis5d$pred) +
  #geom_hline(aes(yintercept=vis5d$bg,linetype='2',color='5')) +
  
  stat_function(aes(x,linetype='1',color='4'),fun=vis4d$pred) +
  #geom_hline(aes(yintercept=vis4d$bg,linetype='2',color='4')) +
  
  stat_function(aes(x,linetype='1',color='1'),fun=vis1d$pred)
  #geom_hline(aes(yintercept=vis1d$bg,linetype='2',color='1'))

print(a)
```


Plotting simulations

```{r simulations}
sim1 <- read.csv('ResultData/docker-0.3-seed-1-iter-1510.csv')
sim1.sigma <- to_km_per_s(0.0542054)
sim1$lb68 <- sim1$predMean-sim1.sigma
sim1$ub68 <- sim1$predMean+sim1.sigma
sim1q1 <- sim1[sim1$obs_time<=95430,] #first 60 days
sim1q2 <- sim1[sim1$obs_time>95430 & sim1$obs_time<=181830,] #next 60 days

sim3 <- read.csv('ResultData/docker-0.3-seed-3-iter-1105.csv')
sim3.sigma <- to_km_per_s(0.0542054)
sim3$lb68 <- sim3$predMean-sim3.sigma
sim3$ub68 <- sim3$predMean+sim3.sigma
sim3q1 <- sim3[sim3$obs_time<=95430,] #first 60 days
sim3q2 <- sim3[sim3$obs_time>95430 & sim3$obs_time<=181830,] #next 60 days

sim4 <- read.csv('ResultData/docker-0.3-seed-4-iter-821.csv')
sim4.sigma <- to_km_per_s(0.0542054)
sim4$lb68 <- sim4$predMean-sim4.sigma
sim4$ub68 <- sim4$predMean+sim4.sigma
sim4q1 <- sim4[sim4$obs_time<=95430,] #first 60 days
sim4q2 <- sim4[sim4$obs_time>95430 & sim4$obs_time<=181830,] #next 60 days

sim5 <- read.csv('ResultData/docker-0.3-seed-5-iter-583.csv')
sim5.sigma <- to_km_per_s(0.0542054)
sim5$lb68 <- sim5$predMean-sim5.sigma
sim5$ub68 <- sim5$predMean+sim5.sigma
sim5q1 <- sim5[sim5$obs_time<=95430,] #first 60 days
sim5q2 <- sim5[sim5$obs_time>95430 & sim5$obs_time<=181830,] #next 60 days

sim8 <- read.csv('ResultData/docker-0.3-seed-8-iter-924.csv')
sim8.sigma <- to_km_per_s(0.0542054)
sim8$lb68 <- sim8$predMean-sim8.sigma
sim8$ub68 <- sim8$predMean+sim8.sigma
sim8q1 <- sim8[sim8$obs_time<=95430,] #first 60 days
sim8q2 <- sim8[sim8$obs_time>95430 & sim8$obs_time<=181830,] #next 60 days

```

```{r}
b <- ggplot() +
    theme_bw()+
    scale_color_manual(values= color_scale,labels=color_scale_labels)+
    scale_fill_manual(values= color_scale,labels=color_scale_labels)+
    xlab("model time ticks")+
    ylab("particle velocity (km/s)")+

  
    geom_point(data=sim1q1,aes(x=obs_time,y=obs),shape=3) +
  
    geom_line(data=sim1q1,aes(x=obs_time,y=predMean,color='1'),lwd=1) +
    geom_ribbon(data=sim1q1,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='1',fill='1'),
                alpha=.06) +
  
    geom_line(data=sim3q1,aes(x=obs_time,y=predMean,color='3'),lwd=1) +
    geom_ribbon(data=sim3q1,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='3',fill='3'),
                alpha=.06) +
  
    geom_line(data=sim4q1,aes(x=obs_time,y=predMean,color='4'),lwd=1) +
    geom_ribbon(data=sim4q1,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='4',fill='4'),
                alpha=.06) +

    geom_line(data=sim5q1,aes(x=obs_time,y=predMean,color='5'),lwd=1) +
    geom_ribbon(data=sim5q1,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='5',fill='5'),
                alpha=.06) +
  
    geom_line(data=sim8q1,aes(x=obs_time,y=predMean,color='8'),lwd=1) +
    geom_ribbon(data=sim8q1,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='8',fill='8'),
                alpha=.06)
print(b)
```

```{r}
d <- ggplot() +
    theme_bw()+
    scale_color_manual(values= color_scale,labels=color_scale_labels)+
    scale_fill_manual(values= color_scale,labels=color_scale_labels)+
  
    geom_point(data=sim1q2,aes(x=obs_time,y=obs),shape=1) +
  
    geom_line(data=sim1q2,aes(x=obs_time,y=predMean,color='1'),lwd=1) +
    geom_ribbon(data=sim1q2,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='1',fill='1'),
                alpha=.06) +
  
    geom_line(data=sim3q2,aes(x=obs_time,y=predMean,color='3'),lwd=1) +
    geom_ribbon(data=sim3q2,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='3',fill='3'),
                alpha=.06) +
  
    geom_line(data=sim4q2,aes(x=obs_time,y=predMean,color='4'),lwd=1) +
    geom_ribbon(data=sim4q2,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='4',fill='4'),
                alpha=.06) +

    geom_line(data=sim5q2,aes(x=obs_time,y=predMean,color='5'),lwd=1) +
    geom_ribbon(data=sim5q2,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='5',fill='5'),
                alpha=.06) +
  
    geom_line(data=sim8q2,aes(x=obs_time,y=predMean,color='8'),lwd=1) +
    geom_ribbon(data=sim8q2,
                aes(x=obs_time,
                    ymin=lb68,
                    ymax=ub68,color='8',fill='8'),
                alpha=.06)
print(d)
```