AdaptiveAdolescentrProject.Rmd

---
title: "The Adaptive Adolescent"
author: "Simy"
date: "04/06/2020"
output: html_document
---

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

pacman::p_load(tidyverse,cowplot,matrixcalc,gganimate,ggformula,pracma)
```

# Adolescence and Ecological Rationality -> UNCERTAINTY NEUTRAL
Often times adolescents are described as incredible risk seekers. This is backed up with evidence from neuroimaging, indicating that adolescents reward sensitivity matures faster than their ability to control their impulses. However,
there is another sociological feature that defines the adolescent period. Adolescents have to find independence. To become independent decision-makers, they need to expierience and explore what it means to make good decisions under uncertainty. Exploration however is often times costly as sometimes it might happen that a risky leap into the unkown results in a bad outcome.
We propose that adolescents prevalent increase in risky behavior can be emergent from their explorative behavior, which is situated in a world that provides an increasing amount of options and possibilites to them. When indivudals enter adoelscence they have expierienced the world within the constraints which their parents set to them. These constraints are then gradually loosened by caregivers. The rich world which lies ahead of them is full of pitfalls but learning how to navigate the worlds uncertainties can only be sucessful on the base of individual expierence.
Whe therefore beleive that an understanding of adolescent risk taking does not necessate the assumption of increased reward sensitivty but can emerge from an exploration policy that seeks out uncertainty in combination with increasing possibilities of how to behave.

# Setup
In our model of risk taking in adolescence, risk taking becomes emergent from adolescnets propensity to explore their environemt, which offers an increasing amount of behaviorual options. We concieve thiese options as probabilisitc patches in a growing environment. 
The difference between adolescnece and childhood is the amount of possible rewards that can be explored. While childhood offers a lot of relativiely safe opportunities, the behaviors that become aware to the adolescent also carry a big detrimental potential which however is unkown to them. 

We setup two environemtns where we understand exploration/exploitation as a multi armed bandid problem. One for childhood, another for adolescence. The child environemtn is a 10x10 grid with each cell in the grid being a bandit that has relatively low variance and losses are not possible.
The adolescent environemnt presents itself as a 100x100 grid where there is more variance and losses are possible.
We start with letting one agent solve this explore exploit dilemma.

#Problem!
We need to do some more thinking on the environment because as it is now lower and upper corners dont differ in how good they are on average. 
That is why agents end up on the right but it doesnt seem to make a difference where exactly on the vertical. Maybe the upper quarter is even more advantageous? 
I now implemented a slightly different environment in which all values can be better than the average but not better than some threshold.
Specificly i re sample the reward if it is over 120. Now, as you may see below, they become "risk-averse".

I am not really happy with that because now, the normative solution is to seek out variance and the story about risk taking is that there is no such normative solution isnt it? 

```{r fig.width=12,fig.height=6}
library("viridis")
envirionmentMeanKids=seq(-50,50,length.out=5)#define the range of kid envirionemt
envirionmentVarianceKids=seq(1,10,length.out=5)
EnvirionemntKids=expand.grid(Mean=envirionmentMeanKids,Variance=envirionmentVarianceKids)# we sample from this later

#define the range of adolescent environment
envirionmentMeanAdol=seq(-100,100,length.out=12)
envirionmentVarianceAdol=seq(1,80,length.out=12)
EnvirionemntAdol=expand.grid(Mean=envirionmentMeanAdol,Variance=envirionmentVarianceAdol)
EnvirionemntAdol%>%ungroup()

EnvirionemntAdol%>%dplyr::rowwise()%>%mutate(value=rnorm(1,mean = Mean,sd=Variance))%>%
  ggplot(aes(x=Mean,y=Variance,fill=value))+geom_tile()+
  scale_fill_viridis(name="Outcome",option="plasma",direction = -1)+
  geom_hline(aes(yintercept=40),size=2)+
  geom_vline(aes(xintercept=0),size=2)+
  annotate("text",x=50,y=38,label=c("Low Risk / Reward"))+
  annotate("text",x=-50,y=38,label=c("Low Risk / Loss"))+
  annotate("text",x=50,y=42,label=c("High Risk / Reward"))+
  annotate("text",x=-50,y=42,label=c("High Risk / Loss"))+
  ggtitle("Adolescent Environment")+theme_minimal()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())->Adolescent


EnvirionemntAdol%>%dplyr::rowwise()%>%mutate(value=rnorm(1,mean = Mean,sd=Variance))%>%
  ggplot(aes(x=Mean,y=Variance,fill=value,alpha=Mean>-46&Variance<40&Mean<46))+geom_tile()+
  scale_fill_viridis(name="Outcome",option="plasma",direction = -1)+
  scale_alpha_discrete(range = c(1, 0.6))+
  geom_hline(aes(yintercept=40),size=2)+
  geom_vline(aes(xintercept=0),size=2)+
  guides(alpha=F)+theme_minimal()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())->Kid

legend=cowplot::get_legend(Kid)
cowplot::plot_grid(Kid+theme(legend.position="none"),
                   Adolescent+theme(legend.position="none",
                                    axis.title.y = element_text(colour="white")),
                   legend,ncol =3,rel_widths = c(1,1,0.2))

ggsave(plot=Kid,filename = "/X_Figures/Environment_new2.png")

ggsave(filename = "../X_Figures/Environments.png")
```

```{r, fig.width=5,fig.height=5}
x=-300:300

  tibble(y=dnorm(x,mean=100,sd=80),
                 x=x)%>%
  ggplot(aes(y=y,x=x,color=x>0))+geom_line(size=3)+
    scale_y_continuous(name="probability")+
    scale_x_continuous(name="Outcome value")+
    scale_color_discrete(name="Outcome",breaks=c("TRUE","FALSE"),labels=c("Loss","Win"))+
    theme_cowplot()+theme(axis.text.y = element_blank(),
                          panel.border = element_rect(colour = "black", fill=NA, size=2),
                          legend.position = c(0.75,0.8))
  ggsave("../X_Figures/HighRiskHighRew.png",dpi=300,height = 5,width = 5)
  
    tibble(y=dnorm(x,mean=100,sd=10),
                 x=x)%>%
ggplot(aes(y=y,x=x,color=x>0))+geom_line(size=3)+
    scale_y_continuous(name="probability")+
    scale_x_continuous(name="Outcome value")+
    scale_color_discrete(name="Outcome",breaks=c("TRUE","FALSE"),labels=c("Loss","Win"))+
    theme_cowplot()+theme(axis.text.y = element_blank(),
                          panel.border = element_rect(colour = "black", fill=NA, size=2),
                          legend.position = c(0.75,0.8))
  ggsave("../X_Figures/LowRiskHighRew.png",dpi=300,height = 5,width = 5)

      tibble(y=dnorm(x,mean=-100,sd=80),
                 x=x)%>%
  ggplot(aes(y=y,x=x,color=x>0))+geom_line(size=3)+
    scale_y_continuous(name="probability")+
    scale_x_continuous(name="Outcome value")+
    scale_color_discrete(name="Outcome",breaks=c("TRUE","FALSE"),labels=c("Loss","Win"))+
    theme_cowplot()+theme(axis.text.y = element_blank(),
                          panel.border = element_rect(colour = "black", fill=NA, size=2),
                          legend.position = c(0.75,0.8))
  ggsave("../X_Figures/HighRiskLowRew.png",dpi=300,height = 5,width = 5)

    tibble(y=dnorm(x,mean=-100,sd=10),
                 x=x)%>%
ggplot(aes(y=y,x=x,color=x>0))+geom_line(size=3)+
    scale_y_continuous(name="probability")+
    scale_x_continuous(name="Outcome value")+
    scale_color_discrete(name="Outcome",breaks=c("TRUE","FALSE"),labels=c("Loss","Win"))+
    theme_cowplot()+theme(axis.text.y = element_blank(),
                          panel.border = element_rect(colour = "black", fill=NA, size=2),
                          legend.position = c(0.75,0.8))
  ggsave("../X_Figures/LowRiskLowRew.png",dpi=300,height = 5,width = 5)

```


# Bayes update
```{r}
library(viridis)
x=-20:50
NormalSamples=tibble(values=rnorm(20,mean=25,sd=6))

rbind(
tibble(y=dnorm(x,mean=20,sd=10),x=x)%>%mutate(type="1"),
tibble(y=dnorm(x,mean=25,sd=6),x=x)%>%mutate(type="2")
)%>%
ggplot()+
  geom_line(aes(y=y,x=x,color=type),size=3)+
  scale_color_viridis_d(name="",option="plasma")+
  #scale_alpha_discrete(range=c(0.4,1))+
  geom_dotplot(data = NormalSamples,aes(x=values,color="3"),fill="#F1F592",binwidth = 1,dotsize = 1.5)+
  #scale_fill_viridis()+
  coord_cartesian(ylim=c(0,0.15))+
  theme_cowplot(20)+
  theme(axis.line.y=element_blank(),
        axis.title.y =element_blank(),
        axis.text.y =element_blank(),
        axis.ticks.y=element_blank()
        )->Posterior

rbind(
tibble(y=dnorm(x,mean=20,sd=10),x=x)%>%mutate(type="1")
#tibble(y=dnorm(x,mean=25,sd=3),x=x)%>%mutate(type="Posterior")
)%>%
ggplot()+
  geom_line(aes(y=y,x=x,color=type),size=3)+
  scale_color_viridis_d(name="",option="plasma")+
  #geom_dotplot(data = NormalSamples,aes(x=values,color="Observations (Likelihood)"),fill="grey",binwidth = 1,dotsize = 1.5)+
  #scale_fill_viridis()+
  theme_cowplot(20)+
  coord_cartesian(ylim=c(0,0.15))+
  theme(axis.line.y=element_blank(),
        axis.title.y =element_blank(),
        axis.text.y =element_blank(),
        axis.ticks.y=element_blank()
        )->Prior

rbind(
tibble(y=dnorm(x,mean=20,sd=10),x=x)%>%mutate(type="1")
#tibble(y=dnorm(x,mean=25,sd=3),x=x)%>%mutate(type="Posterior")
)%>%
ggplot()+
  geom_line(aes(y=y,x=x,color=type),size=3)+
  scale_color_viridis_d(name="",option="plasma")+
  geom_dotplot(data = NormalSamples,aes(x=values,color="2"),fill="#F1F592",binwidth = 1,dotsize = 1.5)+
  #scale_fill_viridis()+
  coord_cartesian(ylim=c(0,0.15))+theme_cowplot(20)+
  theme(axis.line.y=element_blank(),
        axis.title.y =element_blank(),
        axis.text.y =element_blank(),
        axis.ticks.y=element_blank()
        )->Outcomes

ggsave(plot=Posterior,filename="../X_Figures/Posterior.png",dpi=330)
ggsave(plot=Outcomes,filename="../X_Figures/Outcomes.png",dpi=330)
ggsave(plot=Prior,filename="../X_Figures/Prior.png",dpi=330)

```

# load the kalman filter and the UCB rule
Here i define the kalman filter and the choice rule.
The Kalman Agent explores the bandits by representing their mean reward rate and its uncertainty about that mean. Both, beliefs about the mean reward rate and uncertainty about that belief are updated every iteration after observing a new outcome. The Agents behavior can be gouverend by 3 parameters: 

* mu0, or Optimism: its prior assumption about how rewarding the envirionment will be on average. This Paramter will also gouvern the extend of exploration

* var0, or speed of learning: its prior uncertainty about that mean. the higher prior uncertainty is, the more extremely will the agent update its beliefs in the beginning after observing a new sample

* vare, or imprecision: Imprecision will impact the speed of updating as well. The higher the imprecision parameter, the more uncertain the agent will remain about one options mean, even after observing a lot of samples.

New samples are generated using an upper (lower) confidence bound rule where the uncertainty about an option is added to (subtracted from) the beleif about the mean of this option. This policy regulates exploration as well where we can vary uncertaintyphile or phobe behavioral ploicies. This proces is subject to paramter

* beta, or uncertainty aversion (seeking)
```{r}
bayesianMeanTracker <- function(x, y, theta, prevPost=NULL,mu0Par,var0Par){ 
  #Updates the previous posterior based on a single observation
  #parameters
  mu0 <- mu0Par #prior mean
  var0 <- var0Par #prior variance
  vare <- theta[1] #error varriance
  if (is.null(prevPost)){#if no posterior prior, assume it is the first observation
    predictions <- data.frame(mu=rep(mu0,144), sig=rep(var0,144))
  }else{#if previous posterior is provided, update
    predictions <- prevPost
  }
  #Which of the 121 options were chosen at time?
  allopts<-expand.grid(1:12, 1:12)
  chosen <- which(allopts$Var1==x[1] & allopts$Var2==x[2])
  #Kalman gain
  kGain <- predictions$sig[chosen] / (predictions$sig[chosen] + 360)#EnvirionemntAdol$Variance[chosen]^2)# feed the uncertainty in here.
  #update mean
  predictions$mu[chosen] <- predictions$mu[chosen] + (kGain * (y-predictions$mu[chosen]))
  #update variance for observed arm
  predictions$sig[chosen] <- predictions$sig[chosen] * (1 - kGain)
  #return output
  return(predictions)
}
class(bayesianMeanTracker)<- c(class(bayesianMeanTracker), "KalmanFilter")


ucb<-function(out, pars, refactor=F){
  if (refactor==TRUE){
    gamma <- pars[1]
    beta_star<-pars[2]
    #calulate all the upper confidence bounds
    outtotal<-(gamma*out$mu)#+(beta_star*sqrt(out$sig)) #refactored parameters in combination with softmax tau, where gamma = 1/tau and beta_star = beta/tau
    #avoid borderline cases
    #outtotal[outtotal<=0]<-0.0001
    #outtotal[outtotal>100]<-100
    outtotal<-matrix(outtotal, ncol=1, byrow=TRUE)
  }else{
    beta <- pars[1]
    #calulate all the upper confidence bounds
    outtotal<-out$mu+(beta*sqrt(out$sig)) #refactored parameters in combination with softmax tau, where gamma = 1/tau and beta_star = beta/tau
    #avoid borderline cases
    #outtotal[outtotal<=0]<-0.0001
    #outtotal[outtotal>99]<-99
    outtotal<-matrix(outtotal, ncol=1, byrow=TRUE)
  }
  #return them
  return(outtotal)
}
```

### some values to set.
```{r setupObser}
# get sd of whole environemt for normalizing model input
set.seed(as.numeric(Sys.time()))

#get lambda
lambda=0.8
#get beta
beta<-0# this scales risk attitude.
#get tau
tau<-0.8
mu0<-100#exploration bonus
var0<-80
#create a parameter vector
parVec <- c(lambda, lambda, 1, .0001) 
#
ExploreBonus=0
#kernel is RBF
#k<-rbf
#loop through trials
out=NULL
AllChoices=NULL
dummy=NULL
overallCnt=1
dat=expand.grid(x1=1:12,x2=1:12)
```

# Make observations
Here i let the agent learn about the envirionment.
There are three phases of the simulation. Phase one is "childhood". During childhood Agents can explore only the lower middle quadrant. Here all decisions have low risk and there are some low rewards possible. 
Phase two, that occurs after some learning experience (400 samples), can be understood as the onset of adolescence. Here the whole environment becomes availible but the agent does not know about it so they have to learn.
Then after another 400 samples, the agent transitions into "adulthood" where the same environment is still present but the learning expierience now lead to greater exploitation of presumably advantengeous options.

```{r}
for (nround in 1:3){
  #get parameters for participant on that round
  if (nround==1){
    # define vectors that are used by the kalman filter
    lowestx=4
    highestx=9
    sampleVec=as.numeric(rownames(dat[dat$x1>=lowestx & dat$x1<=highestx  & dat$x2<7,]))# here you define where a child should sample from
    ind<-sample(sampleVec,1)
    nTrials=400
  }else {
    ind<-sample(1:144,1)
    nTrials=400
  }
  #random initialization as observation t=0
  #y matrix
  if (nround==1 & overallCnt==1){
    X<-as.matrix(dat[ind,1:2])# generate a new vector of Xs
    y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
  }else if(overallCnt==1) {
    print("Youre an adolescent now")
    X<-as.matrix(dat[ind,1:2])# generate a new vector of Xs
    y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
  }
  #X-start, i.e. all possible observations
  Xstar<-as.matrix(dat[,1:2])
  
  for (trial in 1:nTrials){
    #output by GP with particular parameter settings
    #don't forget mean centering and standardization.... mean is already 0 :)
    if (overallCnt>1){
      out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt], theta=lambda, prevPost=out,mu0Par=mu0,var0Par = var0)
    }else{
      out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt],theta=lambda, prevPost=NULL,mu0Par=mu0,var0Par=var0)
    }
    #utility vector. transpose if you use greedyMean
    #where is everybody?
    #here i need a function that calls bayesianMeanTracker. n times and returns the values X for each n. Also, i need some kind of list, where i save the prior for each instance....
    #
    utilityVec<-ucb(out,beta)
    utilities <- utilityVec - max(utilityVec)
    #softmaximization
    p <- exp(utilities/tau)
    #probabilities
    p <- p/colSums(p)
    #numerical overflow
    p <- (pmax(p, 0.00001))
    p <- (pmin(p, 0.99999))
    #index is sampled proprotionally to softmaxed utitily vector
    if (nround==1){# subset the probability vector so that it corresponds to the right tiles.
      ind<-sample(sampleVec,1,prob=p[dat$x1>=lowestx & dat$x1<=highestx & dat$x2<7,])# sample from a childhood environemnt
      #this monster just scales exploration boni
    }else {
      ind<-sample(1:144, 1, prob=p)# sample from an adolescent environemnt
      # print(ind)
    }
    X<-rbind(X, as.matrix(dat[ind,1:2]))
    #bind y-observations
    y<-rbind(y, as.matrix(rnorm(n=1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance)))# change this into a sample.
    #if(y[overallCnt]<0){
    #  y[overallCnt]=-1*y[overallCnt]^2# make losses more severe. 
    #}
    
    dummy<-data.frame(trial=overallCnt, x=as.numeric(X[overallCnt,1]), y=as.numeric(X[overallCnt,2]),
                      z=as.numeric(y[overallCnt]),round=nround)
    
    AllChoices<-rbind(AllChoices,dummy)
    overallCnt=overallCnt+1
  }
  #dummy data frame
}
#}
#This Here is for Plotting
Plot_dat=expand.grid(x=1:12,y=1:12,trials=0:max(dummy$trial))
Plot_dat$sample=0
Plot_dat$out=0


for (i in 1:length(AllChoices$x)){
  AllChoices$y[i]
  AllChoices$x[i]
  Plot_dat[Plot_dat$x==AllChoices$x[i] & Plot_dat$y==AllChoices$y[i] & Plot_dat$trials==AllChoices$trial[i],]$sample=1
  Plot_dat[Plot_dat$trials==AllChoices$trial[i],]$out=AllChoices$z[i]
}
```

# Animate the exploration.
Here, i plot how often each bandit gets sampled. In each timestep i show the sum of decisions in favor of each option in the color code where brighter colors indicate more samples.

```{r}
library(gganimate)# you can plot the trajectory with gganimate. 
Plot_dat%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%
  ggplot(aes(x=x,y=y,fill=cumDens))+
  geom_tile()+
  scale_fill_viridis_c()+
  scale_y_continuous(name="Variance in Outcomes",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(1,80,length.out=6))+
  scale_x_continuous(name="Mean Outcome",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(-100,100,length.out=6))+
  ggtitle("Sampling Density trial {frame_time}")+
  transition_time(trials)+theme_cowplot()->anim#plot(anim)
animate(anim)
anim_save(animation = animate(anim), filename ="../X_Figures/Sampling_Timecourse_indi.gif")
```




# Now lets contrast the relative & the total number of risky decisions.
Here i compare different developmental phases with respect to how many rewards and losses have been encounterd as a consequence of exploring this envirionment. An exploring agent likely samples highly disadventegous options which agents with more expierience may avoid. Therefore I also show how severe the losses were on average in each developmental phase. An interesting divergence occurs here:
While encountered losses decrease and rewards increase across development; the severity of losses is highest for the adolescent agents who still need to explore the environment.

```{r fig.width=10}
Plot_dat%>%group_by(x,y)%>%filter(trials<=400)%>%mutate(cumDens=cumsum(sample),
                                                        Risk=case_when(
                                                          (out<0)~"Loss",
                                                          (out>=0)~"Win"
                                                        )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out))->Kids

Plot_dat%>%group_by(x,y)%>%filter(trials>=400 & trials<800)%>%mutate(cumDens=cumsum(sample),
                                                                     Risk=case_when(
                                                                       (out<0)~"Loss",
                                                                       (out>=0)~"Win"
                                                                     )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out))->Adolescents

Plot_dat%>%group_by(x,y)%>%filter(trials>800 & trials<2400)%>%mutate(cumDens=cumsum(sample),
                                                                     Risk=case_when(
                                                                       (out<0)~"Loss",
                                                                       (out>=0)~"Win"
                                                                     )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out)
)->YoungAdults

Adolescents$Stage="2"
Kids$Stage="1"
YoungAdults$Stage="3"

Adolescents$Soc="0"
Kids$Soc="0"
YoungAdults$Soc="0"

Plot_labels=c("Kids","Adolescents","Adults")

rbind(Kids,Adolescents,YoungAdults)->AllIndi

```

```{r}
#first show distribution of risky decisions as kid
ggplot(Kids,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="#")+
  ggtitle("outcomes children\nsolo")+
  guides(fill=F)+
  theme_minimal(8)->One

#then show distribution of risky decisions as adolescent
ggplot(Adolescents,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="")+
  ggtitle("outcomes adolescents\nsolo")+
  guides(fill=F)+
  theme_minimal(8)->Two

#then Adults distribution of risky decisions as adolescent
ggplot(YoungAdults,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="")+
  guides(fill=F)+
  ggtitle("outcomes adults\nsolo")+
  theme_minimal(8)->Three

#Now, look at the outcome
ggplot(AllIndi[AllIndi$Risk=="Loss",],aes(x=Stage,y=howMuch,fill=Stage))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  scale_y_continuous(name="Cumulative Loss")+
  scale_x_discrete("Developmental Stage",breaks=c("1","2","3"),labels=Plot_labels)+
  ggtitle("Severety of bad\noutcomes")+
  scale_fill_brewer(palette = "Set2")+
  coord_cartesian(ylim=c(-2000000,0))+
  guides(fill=F)+
  theme_minimal(8)->Four

cowplot::plot_grid(One,Two,Three,NULL,Four,ncol=5,rel_widths = c(1,1,1,0.3,1))->Solo
ggsave(plot = Solo,"../X_Figures/EmergentBehavior_summary.png",width = 10,height = 4)
```

# Number of explorative decisions

Not only the outcomes are interesting but also how much exploraiton happens. 
So i now declare an explorative decision as switching a bandit. This is done by observing whether sampling of one xy pair changes from 0 to 1. 

```{r}
#this is why i love dplyr!
#create count for each new decision.
Plot_dat%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%ungroup()%>% group_by(x,y, trials) %>%
  mutate(newC = ifelse(cumDens == 1 & lag(cumDens)==0, 1, 0))->newExp# if i take the derivative of this this might be exaclty what wouter was intersted in.
newExp[is.na(newExp$newC),]$newC=1#
#now i have to "reset" so that each trial there can only be one new decision and because otherwise this will scale up the cumulative sum too much
newExp$newExplore=0
newExp$generalExplore=0
for( i in 1:length(unique(newExp$trials))){
  newExp[newExp$trials==i,]$newExplore=xor(newExp[newExp$trials==i,]$newC,newExp[newExp$trials==i-1,]$newC)#
  #  newExp[newExp$trials==i,]$generalExplore=xor(newExp[newExp$trials==i,]$sample,newExp[newExp$trials==i-1,]$sample)
}

saveRDS(newExp,"Derivatives_TrialWise.rds")
# to compute the derivative
#newExp<-readRDS("Derivatives_TrialWise.rds")

```

# Plot

I then compute the cumulative sum of this count and subtract the cumulative count from 50 trials earlier from it, which gives me a "slope" of new explorative decisions. But i also contrast this to the raw cumulative count, in order for us to inspect how these values correspond.

```{r}
newExp%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,50)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%
  ggplot(aes(x=trials,y=explore2))+
  #geom_smooth(aes(linetype=Stage),color="black",method="loess")+
  geom_line()+
  # geom_point(size=2)+
  #geom_rect(aes(xmin=50,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  #geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  #geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  #geom_spline(size=1)+
  geom_point(aes(shape=Stage),size=2)+
  scale_shape_manual(name="Developmental Stage",values=c(0,1,2))+
  geom_vline(xintercept = 420, linetype="dotted",color="red")+
  geom_vline(xintercept = 410, linetype="dotted",color="red")+
  geom_vline(xintercept = 800, linetype="dotted",color="red")+
  geom_vline(xintercept = 790, linetype="dotted",color="red")+
  annotate("text",x=150,y=100,label=c("Childhood"))+
  annotate("text",x=600,y=100,label=c("Adolescence"))+
  annotate("text",x=950,y=100,label=c("Adulthood"))+
  #scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
  coord_cartesian(ylim=c(0,130))+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="iteration")+
  ggtitle("Exploration -> Solo")+
  guides(linetype=F)+
  theme_cowplot()

ggsave("../X_Figures/TimecourseExplore.png",width=7,height=4)
```
# Area under the curve

```{r}
newExp%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,50)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"1",
      (trials>400&trials<800)~"2",
      (trials>800)~"3"
    )
  )%>%filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%group_by(Stage)%>%
  summarize(AUC=trapz(trials,explore2))->AUCSolo

# set limits so i can compare the AUC of solo/social later.
ylim=max(AUCSolo$AUC)+50

  AUCSolo%>%ggplot(aes(x=Stage,y=AUC))+geom_bar(stat="identity",color="black",size=1,width = 0.8,alpha=0.5)+
  scale_x_discrete(name="Developmental Stage",breaks=c(1,2,3),labels=c("Kids","Adolescents","Adults"))+
  ggtitle("Exploration (AUC) - Solo")+
  scale_y_continuous(name=c("AUC Explore"),limits=c(0,ylim))+
  theme_cowplot()->AUCSolo

print(AUCSolo)
```

# Other ways to visualize it.
```{r}
newExp%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,50)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200 & trials>50)%>%
  ggplot(aes(x=trials,y=explore2))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  geom_line(size=2)+
  #geom_path()+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="trials")+
  #transition_time(trials)+
  ggtitle("Exploration - Raw")+
  theme_cowplot()
ggsave("../X_Figures/TimecourseExplore_Raw.png")


newExp%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=300,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=300,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=300,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
  scale_y_continuous(name="Explorative Decision")+
  scale_x_continuous(name="trials")+
  ggtitle("Exploration - Cumulative")+
  theme_cowplot()

ggsave("../X_Figures/TimecourseExplore_Cumulative.png")

```


# Social
One particular feature about adolescents is that they are sensitive to social information. In order to invesitage how such a sensitivity impacts 1) exploration behavior and 2) encountered outcomes.
For this i let _n_ agents perform the task in parallel and introduce a bonus to the choice utilites that depends on how many others have visited this option in the previous trial.

```{r setupsims}
# get sd of whole environemt for normalizing model input
set.seed(as.numeric(Sys.time()))

###
###
###
###
### PARAMETERS
###
###
###
#get lambda
#get beta
beta<-0# this scales risk attitude.
#get tau
tau<-0.8
mu0<-100#exploration bonus
var0<-5
#create a parameter vector
parVec <- c(lambda, lambda, 1, .0001) 
#loop through trials
out=NULL
AllChoices_social=NULL
dummy=NULL
overallCnt=1
dat_social=expand.grid(x1=1:12,x2=1:12)
dat_social$others=0#to get the right indices
otherLoc=dat_social$others# location of others in last turn i need it to pass it to the social updating function to avoid social info to have a cumulative effect relative to trials.

HowManyOthers=18
diminishingSocial=0.77
# info about the agents
others <- vector(mode = "list", length = HowManyOthers)# environment for everybody_needed for learning
X_oth<- vector(mode = "list",length = HowManyOthers)#new sample
y_oth<- vector(length = HowManyOthers)#new outcome


```

# Social function
This here is a function that returns the location of other kalman agent in the space. 
You can modify the agents risk attitudes by changing the prior mu0 in the beginning.

```{r}
###
###
### This here is a function that returns the location of other kalman agent in the space. 
### You can modify their risk attitudes by changing the prior mu0 in the beginning.
###
WhereIsEverybody<-function(HowManyOthers,others,otherLoc,dat_social,X_oth,y_oth){
  mu0=100
  if (HowManyOthers==1){
    return(list(Others=dat_social$others,LastSamples=X_oth,LastReturn=y_oth,OthersUtil=others))
  }
  else{
    ####
    ####Update_socialing process for each individual
    ####
    
    if (nround==1 & overallCnt==1){
      sampleVec=as.numeric(rownames(dat_social[dat_social$x1>=lowestx & dat_social$x1<=highestx  & dat_social$x2<7,]))# here you define where a child should sample from
      ind<-sample(sampleVec,1)
      X_oth[[HowManyOthers]]<-as.matrix(dat_social[ind,1:2])# generate a new vector of Xs
      y_oth[HowManyOthers]<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
    }else if(overallCnt==1) {
      ind<-sample(1:144,1)
      print("Youre an adolescent now")
      X_oth[[HowManyOthers]]<-as.matrix(dat_social[ind,1:2])# generate a new vector of Xs
      y_oth[HowManyOthers]<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
    }
    if (overallCnt>1){
      others[[HowManyOthers]]<-bayesianMeanTracker(X_oth[[HowManyOthers]],y_oth[HowManyOthers], theta=lambda, prevPost=others[[HowManyOthers]],mu0Par=mu0,var0Par = var0)
    }else{
      others[[HowManyOthers]]<-bayesianMeanTracker(X_oth[[HowManyOthers]],y_oth[HowManyOthers],theta=lambda, prevPost=NULL,mu0Par=mu0,var0Par=var0)
    }     
    utilityVec<-ucb(others[[HowManyOthers]],beta)
    #browser()
    if(nround==1){
      # no social impact in "kids environment"
      utilityVec=utilityVec#+otherLoc
    } else{
      #social impact follows a power law starting in adolescnece. 
      # browser()
      utilityVec=utilityVec+otherLoc^diminishingSocial
    }
    utilities <- utilityVec - max(utilityVec)
    # utilities=utilities
    #softmaximization
    p <- exp(utilities/tau)
    #probabilities
    p <- p/colSums(p)
    #numerical overflow
    p <- (pmax(p, 0.00001))
    p <- (pmin(p, 0.99999))
    #index is sampled proprotionally to softmaxed utitily vector
    if (nround==1){# subset the probability vector so that it corresponds to the right tiles.
      ind<-sample(sampleVec,1,prob=p[dat_social$x1>=lowestx & dat_social$x1<=highestx & dat_social$x2<7,])# sample from a childhood environemnt
      #this monster just scales exploration boni
    }else {
      ind<-sample(1:144, 1, prob=p)# sample from an adolescent environemnt
      # print(ind)
    }
    X_oth[[HowManyOthers]]<-as.matrix(dat_social[ind,1:2])
    #bind y-observations
    y_oth[HowManyOthers]<-as.matrix(rnorm(n=1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))# change this into a sample.
    dat_social[ind,]$others=dat_social[ind,]$others+1
    
    ####
    ####Now, let the others play
    ####
    WhereIsEverybody(HowManyOthers-1,others,otherLoc,dat_social,X_oth,y_oth)# recursion
  }
}

```



```{r}
#otherLoc=0
overallCnt=1
AllOthers=NULL
for (nround in 1:3){
  #get parameters for participant on that round
  if (nround==1){
    # define vectors that are used by the kalman filter
    lowestx=4
    highestx=9
    sampleVec=as.numeric(rownames(dat_social[dat_social$x1>=lowestx & dat_social$x1<=highestx  & dat_social$x2<7,]))# here you define where a child should sample from
    ind<-sample(sampleVec,1)
    nTrials=400
  }else {
    ind<-sample(1:144,1)
    nTrials=400
  }
  #random initialization as observation t=0
  #y matrix
  if (nround==1 & overallCnt==1){
    X<-as.matrix(dat_social[ind,1:2])# generate a new vector of Xs
    y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
  }else if(overallCnt==1) {
    print("Youre an adolescent now")
    X<-as.matrix(dat_social[ind,1:2])# generate a new vector of Xs
    y<-as.matrix(rnorm(1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance))
  }
  #X-start, i.e. all possible observations
  Xstar<-as.matrix(dat_social[,1:2])
  
  for (trial in 1:nTrials){
    #dat_social$others=0
    #output by GP with particular parameter settings
    #don't forget mean centering and standardization.... mean is already 0 :)
    if (overallCnt>1){
      out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt], theta=lambda, prevPost=out,mu0Par=mu0,var0Par = var0)
    }else{
      out<-bayesianMeanTracker(X[overallCnt,1:2],y[overallCnt],theta=lambda, prevPost=NULL,mu0Par=mu0,var0Par=var0)
    }
    #utility vector. transpose if you use greedyMean
    #where is everybody?
    #here i need a function that calls bayesianMeanTracker. n times and returns the values X for each n. Also, i need some kind of list, where i save the prior for each instance....
    ####
    ####
    Out_Others=WhereIsEverybody(HowManyOthers,others,otherLoc,dat_social,X_oth,y_oth)
    otherLoc=Out_Others$Others
    X_oth=Out_Others$LastSamples
    y_oth=Out_Others$LastReturn
    others=Out_Others$OthersUtil
    
    #add social info
    #browser()
    ####
    ####
    #print(otherLoc)
    utilityVec<-ucb(out,beta)
    if(nround==1){
      # no social impact in "kids environment"
      utilityVec=utilityVec#+otherLoc
    } else{
      #social impact follows a power law starting in adolescnece. 
      utilityVec=utilityVec+otherLoc^diminishingSocial
    }
    #utilityVec=utilityVec+otherLoc^diminishingSocial#add social info
    utilities <- utilityVec - max(utilityVec)
    #softmaximization
    p <- exp(utilities/tau)
    #probabilities
    p <- p/colSums(p)
    #numerical overflow
    p <- (pmax(p, 0.00001))
    p <- (pmin(p, 0.99999))
    #index is sampled proprotionally to softmaxed utitily vector
    if (nround==1){# subset the probability vector so that it corresponds to the right tiles.
      ind<-sample(sampleVec,1,prob=p[dat_social$x1>=lowestx & dat_social$x1<=highestx & dat_social$x2<7,])# sample from a childhood environemnt
      #this monster just scales exploration boni
    }else {
      ind<-sample(1:144, 1, prob=p)# sample from an adolescent environemnt
      # print(ind)
    }
    X<-rbind(X, as.matrix(dat_social[ind,1:2]))
    #bind y-observations
    y<-rbind(y, as.matrix(rnorm(n=1,mean = EnvirionemntAdol[ind,]$Mean,sd=EnvirionemntAdol[ind,]$Variance)))# change this into a sample.
    #if(y[overallCnt]<0){
    #  y[overallCnt]=-1*y[overallCnt]^2# make losses more severe. 
    #}
    AllOthers=rbind(AllOthers,data.frame(Loc=otherLoc,trial=overallCnt))
    dummy<-data.frame(trial=overallCnt, x=as.numeric(X[overallCnt,1]), y=as.numeric(X[overallCnt,2]),
                      z=as.numeric(y[overallCnt]),round=nround)
    
    AllChoices_social<-rbind(AllChoices_social,dummy)
    overallCnt=overallCnt+1
  }
  #dummy data frame
}
#}
#This Here is for Plotting
Plot_dat_social=expand.grid(x=1:12,y=1:12,trials=0:max(dummy$trial))
Plot_dat_social$sample=0
Plot_dat_social$out=0


for (i in 1:length(AllChoices_social$x)){
  AllChoices_social$y[i]
  AllChoices_social$x[i]
  Plot_dat_social[Plot_dat_social$x==AllChoices_social$x[i] & Plot_dat_social$y==AllChoices_social$y[i] & Plot_dat_social$trials==AllChoices_social$trial[i],]$sample=1
  Plot_dat_social[Plot_dat_social$trials==AllChoices_social$trial[i],]$out=AllChoices_social$z[i]
}

## here add the tally of "others".
Plot_dat_social$Others=0
for(k in unique(Plot_dat_social$trials)){
  if (k>0){
    Plot_dat_social[Plot_dat_social$trials==k,]$Others=AllOthers[AllOthers$trial==k,]$Loc
  }
}


```


# Animate

This image is analogous to the above.

```{r}
library(gganimate)# you can plot the trajectory with gganimate. 
Plot_dat_social%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%
  ggplot(aes(x=x,y=y,fill=cumDens))+
  geom_tile()+
  scale_fill_viridis_c(option="inferno")+
  scale_y_continuous(name="Variance in Outcomes",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(1,80,length.out=6))+
  scale_x_continuous(name="Mean Outcome",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(-100,100,length.out=6))+
  ggtitle("Sampling Density Social {frame_time}")+
  transition_time(trials)+theme_cowplot()->animInd#plot(anim)
animate(animInd)
```

# Where are the others in that space

Now i show you the same stuff in the figure but in this image, brighter colors signal that more agents visited this patch on a given trial, as opposed to darker colors.

```{r}
Plot_dat_social%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%
  ggplot(aes(x=x,y=y,fill=Others))+
  geom_tile()+
  scale_fill_distiller(palette = "RdPu")+
  scale_y_continuous(name="Variance in Outcomes",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(1,80,length.out=6))+
  scale_x_continuous(name="Mean Outcome",breaks=c(0,2.5,5.0,7.5,10.0,12.5),labels=seq(-100,100,length.out=6))+
  ggtitle("Everybody Else {frame_time}")+
  transition_time(trials)+theme_cowplot()->animSoc#plot(anim)
animate(animSoc)

anim_save(animation = animate(animSoc), filename ="../X_Figures/Sampling_TimecourseSoc_Oth.gif")
anim_save(animation = animate(animInd), filename ="../X_Figures/Sampling_TimecourseSoc_Ind.gif")

#cowplot::plot_grid(animInd,animSoc)
```



# Now lets contrast the relative & the total number of risky decisions.
This plot shows the same as above. Counts of Wins and Losses and Severity in Losses.
Generally having others in this patch seems to help making better decisions, also during adolescence!

```{r fig.width=10}
Plot_dat_social[length(Plot_dat_social$x),]$out=-0.1# dummy to make the x axis for the adult plot right
Plot_dat_social%>%group_by(x,y)%>%filter(trials<=400)%>%mutate(cumDens=cumsum(sample),
                                                               Risk=case_when(
                                                                 (out<0)~"Loss",
                                                                 (out>=0)~"Win"
                                                               )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out))->Kids

Plot_dat_social%>%group_by(x,y)%>%filter(trials>=400 & trials<800)%>%mutate(cumDens=cumsum(sample),
                                                                            Risk=case_when(
                                                                              (out<0)~"Loss",
                                                                              (out>=0)~"Win"
                                                                            )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out))->Adolescents

Plot_dat_social%>%group_by(x,y)%>%filter(trials>800 & trials<1200)%>%mutate(cumDens=cumsum(sample),
                                                                            Risk=case_when(
                                                                              (out<0)~"Loss",
                                                                              (out>=0)~"Win"
                                                                            )
)%>%group_by(Risk)%>%summarise(howMany=n()/144,
                               howMuch=sum(out)
)->YoungAdults

Adolescents$Stage="2"
Kids$Stage="1"
YoungAdults$Stage="3"

Adolescents$Soc="1"
Kids$Soc="1"
YoungAdults$Soc="1"

Plot_labels=c("Kids","Adolescents","Adults")

rbind(Kids,Adolescents,YoungAdults)->AllSoc
```

```{r}

#first show distribution of risky decisions as kid
ggplot(Kids,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="#")+
  ggtitle("outcomes children\nsocial")+
  guides(fill=F)+
  theme_minimal(8)->One

#then show distribution of risky decisions as adolescent
ggplot(Adolescents,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="")+
  ggtitle("outcomes adolescence\nsocial")+
  guides(fill=F)+
  theme_minimal(8)->Two

#then Adults distribution of risky decisions as adolescent
ggplot(YoungAdults,aes(x=Risk,y=howMany,fill=Risk))+
  stat_summary(geom="bar",color="black",size=1,width=0.9)+
  coord_cartesian(ylim=c(0,700))+
  scale_y_continuous(name="")+
  guides(fill=F)+
  ggtitle("outcomes adults\nsocial")+
  theme_minimal(8)->Three

#Now, look at the outcome
ggplot(AllSoc[AllSoc$Risk=="Loss",],aes(x=Stage,y=howMuch,fill=Stage))+
  stat_summary(geom="bar",color="black")+
  scale_y_continuous(name="Cumulative Loss")+
  scale_x_discrete("Developmental Stage",breaks=c("1","2","3"),labels=Plot_labels)+
  ggtitle("Severety of bad\noutcomes")+
  scale_fill_brewer(palette = "Set2")+ 
  coord_cartesian(ylim=c(-2000000,0))+
  guides(fill=F)+
  theme_minimal(8)->Four

cowplot::plot_grid(One,Two,Three,NULL,Four,ncol=5,rel_widths = c(1,1,1,0.3,1))->Social
ggsave(plot = Social,filename="../X_Figures/EmergentBehavior_summary_Social.png",width = 10,height = 4)


```

# joint plot
```{r}
cowplot::plot_grid(Solo,NULL,Social,nrow=3,rel_heights = c(1,0.2,1))->both
ggsave(filename="../X_Figures/Solo_Social_Outcomes.png",plot = both, width = 9, height=6)
```



# Number of explorative decisions
This here is equivalent to the above, but now for social agents. In comparison to solo agents we can see that introducing social information and a social following rule, 

```{r}
#this is why i love dplyr!
#create count for each new decision.
Plot_dat_social%>%group_by(x,y)%>%mutate(cumDens=cumsum(sample))%>%ungroup()%>% group_by(x,y, trials) %>%
  mutate(newC = ifelse(cumDens == 1 & lag(cumDens)==0, 1, 0))->newExpSoc# if i take the derivative of this this might be exaclty what wouter was intersted in.
newExpSoc[is.na(newExpSoc$newC),]$newC=1#
#now i have to "reset" so that each trial there can only be one new decision and because otherwise this will scale up the cumulative sum too much
newExpSoc$newExplore=0
newExpSoc$generalExplore=0
for( i in 1:length(unique(newExpSoc$trials))){
  newExpSoc[newExpSoc$trials==i,]$newExplore=xor(newExpSoc[newExpSoc$trials==i,]$newC,newExpSoc[newExpSoc$trials==i-1,]$newC)#
  #  newExp[newExp$trials==i,]$generalExplore=xor(newExp[newExp$trials==i,]$sample,newExp[newExp$trials==i-1,]$sample)
}

saveRDS(newExpSoc,"Derivatives_TrialWise_Social.rds")
# to compute the derivative
```

```{r}
# newExp%>%ungroup()%>%
#   mutate(newCum=cumsum(newExplore))%>%
#   group_by(trials)%>%
#   summarise(explore=max(newCum))%>%mutate(
#     explore2=explore-lag(explore,50)
#   )%>%mutate(
#     Stage=case_when(
#       (trials<401)~"Kids",
#       (trials>400&trials<800)~"Adolescents",
#       (trials>800)~"Adults"
#     )
#   )%>%filter(trials<1200)%>%
#   ggplot(aes(x=trials,y=explore2))+
#   # geom_point(size=2)+
#   geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
#   geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
#   geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
#   geom_spline(size=2,)+
#   geom_vline(xintercept = 400, linetype="dotted",size=2)+
#   geom_vline(xintercept = 800, linetype="dotted",size=2)+
#   scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
#   scale_y_continuous(name="Explorative Decisions")+
#   scale_x_continuous(name="trials")+
#   ggtitle("Exploration - Spline")+
#   theme_cowplot()
# ggsave("../X_Figures/TimecourseExplore_Spline.png")


newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,50)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%
  ggplot(aes(x=trials,y=explore2))+
  geom_line()+
  #geom_rect(aes(xmin=50,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  #geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  #geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  #geom_spline(size=1)+
  geom_point(aes(shape=Stage),size=2)+
  scale_shape_manual(name="Developmental Stage",values=c(0,1,2))+
  geom_vline(xintercept = 420, linetype="dotted",color="royalblue")+
  geom_vline(xintercept = 410, linetype="dotted",color="royalblue")+
  geom_vline(xintercept = 800, linetype="dotted",color="royalblue")+
  geom_vline(xintercept = 790, linetype="dotted",color="royalblue")+
  annotate("text",x=150,y=100,label=c("Childhood"))+
  annotate("text",x=600,y=100,label=c("Adolescence"))+
  annotate("text",x=950,y=100,label=c("Adulthood"))+
  #scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
  coord_cartesian(ylim=c(0,130))+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="iteration")+
  ggtitle("Exploration -> Social")+
  guides(linetype=F)+
  theme_cowplot()
ggsave("../X_Figures/TimecourseExplore_Raw_Social.png",width = 7,height = 4)

```

#Area under the curve

```{r}
# Area under the curve
newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,50)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"1",
      (trials>400&trials<800)~"2",
      (trials>800)~"3"
    )
  )%>%filter(!is.na(Stage) & trials %in% round(seq(50,1200,length.out = 24)))%>%group_by(Stage)%>%
  summarize(AUC=trapz(trials,explore2))%>%ggplot(aes(x=Stage,y=AUC))+geom_bar(stat="identity",color="black",size=1,width = 0.8,alpha=0.5)+
  scale_x_discrete(name="Developmental Stage",breaks=c(1,2,3),labels=c("Kids","Adolescents","Adults"))+
  ggtitle("Exploration (AUC) - Social")+
  scale_y_continuous(name=c("AUC Explore"),limits=c(0,ylim))+
  theme_cowplot()->AUCSocial
```

#Compare Solo & Social AUCs
```{r}
cowplot::plot_grid(AUCSolo,AUCSocial+theme(axis.text.y=element_text(color="white")))
```

```{r}
newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(newExplore))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=300,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=300,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=300,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "YlGnBu")+
  scale_y_continuous(name="Explorative Decision")+
  scale_x_continuous(name="trials")+
  ggtitle("Exploration - Cumulative")+
  theme_cowplot()

ggsave("../X_Figures/TimecourseExplore_Cumulative_Social.png")

```



# 4 types of Social decisions.
But The social dimension also allows for examining additional dimensions of behavior.
Specifically a four fold pattern of exploration emerges in respect to social content.

1) Individuals could switch to a patch that is favoured by SI 
2) Individuals could switch to a patch that is _not_ favoured by SI
3) Individuals could stay at a patch that has been favoured by SI 
4) Individuals could stay at a patch that has been _not_ favoured by SI 

## Going Where Others have been before.

```{r}
newExpSoc$GoToOthers=0
newExpSoc$GoToNoone=0
newExpSoc$StayWithOthers=0
newExpSoc$StayWithNoone=0

for( i in 1:length(unique(newExpSoc$trials))){
  #If you explored a new option that was sampled my more than 0 others before.
  newExpSoc[newExpSoc$trials==i,]$GoToOthers=as.numeric(newExpSoc[newExpSoc$trials==i,]$newExplore & (newExpSoc[newExpSoc$trials==i-1,]$Others!=0))#
  
  #If you explored a new option that was sampled my 0 others before.
  newExpSoc[newExpSoc$trials==i,]$GoToNoone=as.numeric(newExpSoc[newExpSoc$trials==i,]$newExplore & (newExpSoc[newExpSoc$trials==i-1,]$Others==0))#
  #  newExp[newExp$trials==i,]$generalExplore=xor(newExp[newExp$trials==i,]$sample,newExp[newExp$trials==i-1,]$sample)
  
  # If you sample an option, where others have been and still are and you dont explore again.
  newExpSoc[newExpSoc$trials==i,]$StayWithOthers=as.numeric(newExpSoc[newExpSoc$trials==i,]$sample & (newExpSoc[newExpSoc$trials==i-1,]$Others!=0) & (newExpSoc[newExpSoc$trials==i,]$Others!=0) & (newExpSoc[newExpSoc$trials==i,]$newExplore==0))#
  
  # If you sample an option, where no others have been and still none are and you dont explore again.
  newExpSoc[newExpSoc$trials==i,]$StayWithNoone=as.numeric(newExpSoc[newExpSoc$trials==i,]$sample & (newExpSoc[newExpSoc$trials==i-1,]$Others==0) & (newExpSoc[newExpSoc$trials==i,]$Others==0) &(newExpSoc[newExpSoc$trials==i,]$newExplore==0))#
}




newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(GoToNoone))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,100)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore2))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "Reds")+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="trials")+
  ggtitle("Switch Where Noone Is")+
  theme_cowplot()
ggsave("../X_Figures/TimecourseExplore_Raw_Social.png")





newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(GoToOthers))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,100)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore2))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "Reds")+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="trials")+
  ggtitle("Switch to Where Others Are")+
  theme_cowplot()





newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(StayWithOthers))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,100)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore2))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "Reds")+
  scale_y_continuous(name="Explorative Decisions")+
  scale_x_continuous(name="trials")+
  ggtitle("Stay With Others")+
  theme_cowplot()




newExpSoc%>%ungroup()%>%
  mutate(newCum=cumsum(StayWithNoone))%>%
  group_by(trials)%>%
  summarise(explore=max(newCum))%>%mutate(
    explore2=explore-lag(explore,100)
  )%>%mutate(
    Stage=case_when(
      (trials<401)~"Kids",
      (trials>400&trials<800)~"Adolescents",
      (trials>800)~"Adults"
    )
  )%>%filter(trials<1200)%>%
  ggplot(aes(x=trials,y=explore2))+
  # geom_point(size=2)+
  geom_rect(aes(xmin=0,xmax=400,ymin=-5,ymax=100,fill="1"),alpha=0.1)+
  geom_rect(aes(xmin=400,xmax=800,ymin=-5,ymax=100,fill="2"),alpha=0.1)+
  geom_rect(aes(xmin=801,xmax=1200,ymin=-5,ymax=100,fill="3"),alpha=0.1)+
  geom_spline(size=2)+
  geom_vline(xintercept = 400, linetype="dotted",size=2)+
  geom_vline(xintercept = 800, linetype="dotted",size=2)+
  scale_fill_brewer(name="Developmental Stage",labels=c("Kids","Adolescents","Adults"),palette = "Reds")+
  scale_y_continuous(name="Sum Previous Decisions")+
  scale_x_continuous(name="trials")+
  ggtitle("Stay Alone")+
  theme_cowplot()

```