#' Compile results and generate graphs for one species.
#'
#' Creates graphs by taking X and Y series from slots in a Species object, and
#' combines some of those data into a data frame that is returned with the
#' plots. Reads instructions from the plotlist.csv file to decide which plots
#' and series to produce (there are different sets of plots for prey and
#' predators). Note that since a dataframe cannot be a ragged array, the number
#' of rows in the data frame that is returned by the function is either nYears or
#' nTimes, but certain items (catch, recruits) are *plotted* by year, even when
#' the species is tracked by timestep.
#'
#' @param species A Species object
#' @param slotlist A list of slots, from which graphs are to be made.
#' @param nYears Number of years in the simulation.
#' @param nSteps Number of timesteps per year (ignored for some species)
#'
#' @return A list with the values of certain vectors assembled into a data frame, and a list of plots.
#' @export
#' @family results
#' @keywords internal
compile_results<-function(species,plot_master,nYears,nSteps,preynames){
if(!is(species,"Species")) stop("getresults() takes a Species object")
pm<-plot_master
pm[pm==""]<-NA #Don't skip this step!
#Check to make sure we're not trying to plot a non-existent slot
for (i in 1:nrow(pm)){
if(is(species,pm$obj_type[i])){
if(!is.na(pm$slot.x[i]) && !pm$slot.x[i] %in% slotNames(species)){
stop(paste("Cannot plot",pm$slot.x[i],": slot not found in the",species@name,"object."))
}
if(!is.na(pm$slot.y[i]) && !pm$slot.y[i] %in% slotNames(species)){
stop(paste("Cannot plot",pm$slot.y[i],": slot not found in the",species@name,"object."))
}
}
}#for 1 to nrow(pm)
speciesname<-species@name
res<-list()
plotlist<-list()
#Set up two alternative X vectors (nYears, and nTimes) for the output dataframe
#(NOTE: this x is only used for plotting if the xtype column of pm is "default".)
x_ts<-1:(nSteps * nYears)
x_yr<-1:nYears
#Initialize the results dataframe (to which Y columns will be added).
#For convenience, add a intra-year timestep to by-timestep dataframes
#This dataframe is for output, not plotting.
if(species@uses_timestep){
data<-data.frame(x_ts,rep(1:nSteps,nYears))
names(data)<-c("Timestep","Intr_year_step")
} else{
data<-data.frame(x_yr)
names(data)<-"Year"
}
#Step through the plotmaster dataframe; make plots and add data to the result
#dataframe as instructed.
for(i in 1:nrow(pm)){
yname<-NA
pdata<-NA
ptitle<-NA
#Get x and y labels, and plot title if specified
if(!is.na(pm$xlab[i])){
xname<-pm$xlab[i]
}
if(!is.na(pm$title[i])){
ptitle<-pm$title[i]
}
#Get the y range option
y_range<-pm$y_range[i]
#Get the predicted option
predicted<-pm$predicted[i]
series_type<-pm$series_type[i]
#Make one set of plots for pm$obj_type=="Prey", a different set for Predators
if(is(species,pm$obj_type[i])){
if(!pm$plottype[i]=="availability"){
#Get y axis label
if(!is.na(pm$ylab[i])){
yname<-pm$ylab[i]
}else{
yname<-pm$slot.y[i] #default; may be overridden
}
#Get y values
vals<-slot(species,pm$slot.y[i])
#Skip this row in pm if there are no data in the y slot.
if(length(vals)==0){
next
#print(paste("There is no data for",yname,"in",speciesname)) #next
}
#Get x values if the x slot is specified, or create an x axis (time or years)
if(!is.na(pm$slot.x[i])){
px<-slot(species,pm$slot.x[i])
xname<-pm$slot.x[i]
} else{
#Specify an x-axis for timesteps or years
if(length(vals)==length(x_ts)){
px<-x_ts
xname<-"Timestep"
} else {
px<-x_yr
xname<-"Year"
}
}
#Make a small dataframe (pdata) for plotting
if(length(vals) != length(px)){
stop(paste("Cannot plot",yname,"; The length of",yname,"does not match the length of x (",length(px),")"))
} else{
pdata<-data.frame(px,vals)
}
#Subset it, in the case that the data are in timestep units but we want the plot in years,
#and replace the x-axis with years.
if(pm$xtype[i]=="step1only"){
if(species@uses_timestep && nrow(pdata)==length(x_ts)){
#First, subset the data
pdata<-get_every_nth(nSteps,pdata,xname,yname) #subset the df
#Don't change the x-axis or label unless the name is Timestep
#because xy plots may have a different x-axis
if(xname == "Timestep"){
pdata<-data.frame(x_yr,pdata[,match(yname,names(pdata))])#replace the x column with year
xname<-"Year"
}
names(pdata)<-c(xname,yname)
}
}
#Add the (unshortened) Y timeseries to the results dataframe if in_dataframe=1 in plotmaster
if(pm$in_dataframe[i]){
data<-data.frame(data,vals)
names(data)[ncol(data)]<-pm$slot.y[i] #name the column
}
}#plot type is not availability
#Several plotting options:
#Plot a time series
if(pm$plottype[i]=="timeseries"){
newplot<-makeplot(pdata,xname,yname,speciesname,ptitle,y_range,series_type)
}
#Plot a time series, summed by year
if(pm$plottype[i]=="summed_by_year"){
if(pm$slot.y[i]=="Ct" && species@uses_timestep && species@fished){
xname<-"Year"
newplot<-plot_summed_by_year(nYears,nSteps,pdata,xname,yname,speciesname,ptitle,y_range,series_type)
}
}
#Plot two variables X and Y against each other
if(pm$plottype[i]=="xy"){
if(!is.na(predicted)){
#get the predicted data
fitted_data<-calc_predicted(species,predicted)
names(fitted_data)<-names(pdata)
#Send the observed and predicted data separately to ggplot
newplot<-plot_fitted(fitted_data,pdata,xname,yname,speciesname,ptitle,y_range)
} else{
#pdata<-na.exclude(pdata) #remove rows with missing values
pdata<-pdata[order(pdata[,1]),] #sort the data by x values.
newplot<-makeplot(pdata,xname,yname,speciesname,ptitle,y_range,series_type)
}
}#plottype=xy
#Plot just a function (no observed data)
if(pm$plottype[i]=="function_only"){
if(pm$predicted[i]=="availability"){
fndata<-calc_predicted(species,predicted)
newplot<-plot_function_only(fndata,xname,yname,speciesname,ptitle)
}
}
#A special n-panel (1 per prey) plot of prey availability per predator
if(pm$plottype[i]=="availability"){
observed<-data.frame()
predicted<-data.frame()
for(k in 1:length(preynames)){
#Put observed data in a dataframe
x<-species@prey_availability[[preynames[[k]]]]$B_agepref
y<-species@prey_availability[[preynames[[k]]]]$B_avail
preysp<-rep(preynames[k],length(x))
obs<-data.frame(preysp,x,y,rep("obs",length(x)))
observed<-rbind(observed,obs)
#Get predicted values
pred<-calc_predicted(species,"availability",preynames[k])
preysp<-rep(preynames[k],nrow(pred))
pred<-data.frame(preysp,pred,rep("pred",nrow(pred))) #add the preyname column
predicted<-rbind(predicted,pred)
}#for prey species
names(observed)<-c("species","x","y","dtype")
names(predicted)<-c("species","x","y","dtype")
pdata<-rbind(observed,predicted)
yname<-pm$ylab[i]
newplot<-plot_fitted_multipanel(pdata,xname,yname,speciesname,ptitle)
}#plottype="availability
if(pm$plottype[i]=="compare2xy" && pm$slot.y[i]=="R_Punt"){
observed<-data.frame()
predicted<-data.frame()
#Put observed data in a dataframe
xobs<-species@Nm
yobs<-species@R
xobs_Punt<-species@Nm
yobs_Punt<-species@R_Punt
method<-c(rep("Hilborn",length(x)),rep("Punt",length(x)))
observed<-data.frame(method,c(xobs,xobs_Punt),c(yobs,yobs_Punt),rep("obs",length(x)))
names(observed)<-c("method","x","y","dtype")
#Get predicted values
pred_R<-calc_predicted(species,"recruitment")
pred_R<-data.frame(rep("Hilborn",nrow(pred_R)),pred_R,rep("pred",nrow(pred_R)))
names(pred_R)<-c("method","x","y","dtype")
pred_Punt<-calc_predicted(species,"R_Punt")
pred_Punt<-data.frame(rep("Punt",nrow(pred_Punt)),pred_Punt,rep("pred",nrow(pred_Punt)))
names(pred_Punt)<-c("method","x","y","dtype")
predicted<-rbind(pred_R,pred_Punt)
pdata<-rbind(observed,predicted)
newplot<-plot_compare2fitted_multipanel(pdata,xname,yname,speciesname,ptitle)
}#plottype=compare2xy
#put the new plot in the growing plotlist
plotlist<-c(plotlist,list(newplot))
#name the list element
names(plotlist)[length(plotlist)]<-pm$plotname[i]
}#If the species is this obj_type
}#for i in 1:nrow(pm)
#Add the plot and dataframe to a list, and name them
reslist<-list(data,plotlist)
names(reslist)<-c("data","plots")
reslist
} #compile_results()
#' Plot a time series (a utility function)
#'
#' @param data A data frame in which the X series to be plotted is in column 1.
#' @param xname A display name for the X-axis
#' @param yname A display name for the Y-axis
#' @param speciesname A display name for the species
#' @param ptitle A partial title
#' @param y_range Character; "include_0" to include zero, or "auto"
#'
#' @return A plot
#' @export
#' @keywords internal
#' @family results
makeplot<-function(pdata,xname,yname,speciesname,ptitle,y_range,series_type){
if(y_range=="include_0"){
if(series_type=="line"){
tmp<-ggplot2::ggplot(pdata,aes(x=pdata[,1],y=pdata[,2])) +
geom_line() + ylim(c(0,NA)) + xlab(xname) +
ylab(yname) + ggtitle(paste(speciesname,ptitle))
}
if(series_type=="point"){
tmp<-ggplot2::ggplot(pdata,aes(x=pdata[,1],y=pdata[,2])) +
geom_point() + ylim(c(0,NA)) + xlab(xname) +
ylab(yname) + ggtitle(paste(speciesname,ptitle))
}
if(series_type=="both" || series_type=="auto"){
tmp<-ggplot2::ggplot(pdata,aes(x=pdata[,1],y=pdata[,2])) +
geom_point() + geom_line() + ylim(c(0,NA)) + xlab(xname) +
ylab(yname) + ggtitle(paste(speciesname,ptitle))
}
} else{
tmp<-ggplot2::ggplot(pdata,aes(x=pdata[,1],y=pdata[,2])) +
geom_point() + geom_line() + xlab(xname) + ylab(yname) +
ggtitle(paste(speciesname,ptitle))
}
newname<-paste0("p.",yname)
assign(newname,tmp)
tmp
}
#' Sum a by-timestep series by years and plot it.
#'
#' A utility function for plotting.
#'
#' @param data A small data frame with the by-timestep series in it
#' @param item The name of the time series (matches a colname in data)
#' @param xname The name of the x series
#' @param speciesname The name of the species
#' @param nYears The number of years in the time series
#' @param nSteps The number of timesteps per year
#' @return A plot of values of a vector, summed by year.
#' @export
#' @keywords internal
#' @family results
plot_summed_by_year<-function(nYears,nSteps,pdata,xname,yname,...){
#For species that use timesteps but where we want to plot values by years, make
#a vector of the year number, e.g.,: 1 1 1 1 2 2 2 2 3 3 3 3 (here, nSteps=4).
ts<-1:(nYears * nSteps)
yrindex<-(ts + (nSteps-1)) %/% nSteps
years<-1:nYears
#Sum the y-values by year
yrvals<-tapply(pdata[,2],yrindex,sum,na.rm=T)
yrsumdata<-data.frame(years,yrvals)
#Plot the new dataframe
newplot<-makeplot(yrsumdata,xname,yname,...)
newplot
}
#' Get every nth value of a time series
#'
#' A utility function that returns every nth intra-annual step from a time
#' series that includes intra-annual steps.
#'
#' @param nSteps The number of timesteps per year
#' @param pdata A small data frame with the by-timestep series in it
#' @param xname The name of the x series
#' @param yname The name of the y series
#' @return A data frame (x,y) of every nth value of a time series
#' @export
#' @keywords internal
#' @family results
get_every_nth<-function(nSteps,pdata,xname,yname){
#Take every nth value of the item to plot (here, the 1st timestep of each year)
x_nth<-pdata[seq(1,nrow(pdata),nSteps),1]
y_nth<-pdata[seq(1,nrow(pdata),nSteps),2]
#put it into a dataframe and return it
nth_data<-data.frame(x_nth,y_nth)
names(nth_data)<-c(xname,yname)
nth_data
}
#' Plot a fitted function (observed vs. expected)
#'
#' @param fitted_data A dataframe of expected data
#' @param pdata A dataframe of observed data
#' @param xname X-axis label
#' @param yname Y-axis label
#' @param speciesname Species name (for title)
#' @param ptitle Partial title
#' @return A plot
#' @export
#' @keywords internal
#' @family results
plot_fitted<-function(fitted_data,pdata,xname,yname,speciesname,ptitle,...){
tmp<-ggplot2::ggplot(pdata,aes(x=pdata[,1],y=pdata[,2])) +
geom_point() + geom_line(data=fitted_data,aes(x=fitted_data[,1],y=fitted_data[,2])) +
xlab(xname) + ylab(yname) + ggtitle(paste(speciesname,ptitle))
tmp
}
#' Plot a multi-panel fitted function (observed vs. expected)
#'
#' @param observed.df A dataframe of observed data
#' @param predicted.df A dataframe of predicted data
#' @param xname X-axis label
#' @param yname Y-axis label
#' @param speciesname Species name (for title)
#' @param ptitle Partial title
#' @return A plot
#' @export
#' @keywords internal
#' @family results
plot_fitted_multipanel<-function(pdata,xname,yname,speciesname,ptitle,...){
tmp<-ggplot2::ggplot(pdata,aes(x=x,y=y,group=dtype)) +
geom_point(data=pdata[pdata$dtype=="obs",]) +
geom_line(data=pdata[pdata$dtype=="pred",]) +
facet_wrap(~species,scales="free") +
xlab(xname) + ylab(yname) + ggtitle(paste(speciesname,ptitle))
tmp
}
plot_compare2fitted_multipanel<-function(pdata,xname,yname,speciesname,ptitle,...){
tmp<-ggplot2::ggplot(pdata,aes(x=x,y=y,colour=method)) +
geom_point(data=pdata[pdata$dtype=="obs",]) +
geom_line(data=pdata[pdata$dtype=="pred",]) +
xlab(xname) + ylab(yname) + ggtitle(paste(speciesname,ptitle))
tmp
}
#' Plot a function (expected values only)
#'
#' @param fndata A dataframe of expected data
#' @param xname X-axis label
#' @param yname Y-axis label
#' @param speciesname Species name (for title)
#' @param ptitle Partial title
#' @return A plot
#' @export
#' @keywords internal
#' @family results
plot_function_only<-function(fndata,xname,yname,speciesname,ptitle){
tmp<-ggplot2::ggplot(fndata,aes(x=fndata[,1],y=fndata[,2])) +
geom_line() + xlab(xname) + ylab(yname) + ggtitle(paste(speciesname,ptitle))
tmp
}
#' Calculate predicted values
#'
#' A utility function for plotting, which calculates various kinds of predicted curves
#' for overlay in plots.
#' @param species A species object (whatever is being plotted)
#' @param predicted Character; the name of a predicted series
#' @param preyname Character; the name of a prey species (only used for availabililty)
#'
#' @return A data frame of predicted data (x and y; 100 rows)
#' @export
#' @keywords internal
#' @family results
calc_predicted<-function(species,predicted,preyname){
if(predicted=="recruitment" && is(species,"Prey") && !is(species,"Other_forage")){
alpha<-species@recruitment$alpha
beta<-species@recruitment$beta
bmax<-max(species@B_sp,na.rm=T)
x<-seq(0,bmax,bmax/100)
y<-sapply(x,function(b) calc_Ricker(alpha,beta,b))
}
if(predicted=="recruitment" && is(species,"Other_forage")){
maxb<-max(species@B_sp,na.rm=T)
x<-seq(0,maxb,maxb/100)
y<-rep(species@R_baseyear,length(x))
}
if(predicted=="recruitment" && is(species,"Predator")){
p<-species@recruitment$BHp
c<-species@recruitment$BHc
maxnm<-max(species@Nm,na.rm=T)
x<-seq(0,maxnm,maxnm/100)
y<-sapply(x,function(n) calc_Bev_Holt(n,p,c))
}
if(predicted=="prey_effect" && is(species,"Predator")){
maxd<-max(species@Dy,na.rm=T)
x<-seq(0,maxd,maxd/100)
y<-sapply(x,function(d) calc_p(species,d))
}
#For plotting R_Punt, phi=1 (no effect of prey)
if(predicted=="R_Punt" && is(species,"Predator")){
maxd<-max(species@Nm,na.rm=T)
x<-seq(0,maxd,maxd/100)
y<-numeric(length(x))
K1P<-species@K1P
#Estimate N1P from Nm, via ratios of survival
am<-species@am #age at maturity
nAges<-species@P + 1
S1<-exp(-(species@Mv[2])) #age-1 survival
Sx<-sapply(1:(nAges-1),function(x) S1^(x-1))
for(i in 1:length(x)){
Nm<-x[i]
Stot<-0
Sam<-0
for (j in 1:(nAges - 1)){
Stot<-Stot + Sx[j]
if(j >= am){
Sam<-Sam + Sx[j]
}
}#for j
N1P<-Nm * (Stot/Sam)
y[i]<-calc_predRecruits_Punt(species,Nm,1,N1P,K1P)
}
}
#For availability, create a data frame for the prey in question
#(availability parameters vary by predator/prey combinations)
if(predicted=="availability" && is(species,"Predator")){
ap<-species@prey_availability[[preyname]]$ap
ac<-species@prey_availability[[preyname]]$ac
maxb<-max(species@prey_availability[[preyname]]$B_agepref,na.rm=T)
x<-seq(0,maxb,maxb/100)
y<-sapply(x,function(b) calc_availability(ap,ac,b))
}
#put the results in a data frame and add a column called "id" labeling the values as predicted
df<-data.frame(x,y)
df
}
#' Make a 4-panel plot
#'
#' Make a 4-panel plot that shows the abundance of the key species (spawning
#' biomass for prey; mature adults for predators) .
#'
#' @param prey A list containing the Prey objects
#' @param predator A list containing the Predator objects
#'
#' @return A single 4-panel plot showing the abundance of 4 key species
#' @export
#'
#' @keywords internal
make_4panel<-function(prey,predator,scenario,nSteps){
pl_list<-c("Sardine","Anchovy","Other_forage","Sealion","Pelican")
data4p<-data.frame()
#Create a data frame of B_sp for prey, and Nm for predators
for(preyx in prey){
if(preyx@name %in% pl_list){
names4p<-preyx@name
#pick out every timestep 1 if necessary
if(preyx@uses_timestep){
xdata<-1:length(preyx@B_sp)
ydata<-preyx@B_sp
temp<-data.frame(xdata,ydata)
temp<-get_every_nth(nSteps,temp,"Year","B_sp")
newd<-data.frame(rep(preyx@name,nrow(temp)),temp)
}else{
newd<-data.frame(rep(preyx@name,length(preyx@B1PperB0)),1:length(preyx@B1PperB0),preyx@B1PperB0)
}
names(newd)<-c("species","year","spawners")
if(nrow(data4p)==0){
data4p<-newd
} else{
data4p<-rbind(data4p,newd)
}
}#preyname in pl_list
}#preyx in prey
for(predx in predator){
if(predx@name %in% pl_list){
names4p<-predx@name
newd<-data.frame(rep(predx@name,length(predx@NperK)),1:length(predx@NperK),predx@NperK)
names(newd)<-c("species","year","spawners")
if(nrow(data4p)==0){
data4p<-newd
} else{
data4p<-rbind(data4p,newd)
}
}#predname in pl_list
}#predx in predator
#Plot the data and return the plot
p4<-ggplot(data4p,aes(x=year,y=spawners)) +
geom_line() + facet_wrap(~species,scales="free",ncol=2) + ylim(c(0,NA)) +
xlab("Year") + ylab("B1+, B1+/B0, N, or N/K") + ggtitle(paste("Scenario:",scenario))
p4
}
#show all plots from a result object
#' Show all plots in a result object
#'
#' @param res A result object
#'
#' @return all the plots (in a list, but they display by default)
#' @export
#'
allplots<-function(res){
species<-names(res$species)
plist<-list()
for(i in 1:length(species)){
plist<-c(plist,list(res[[species[i]]]$plots))
}
plist
}
#' Show one type of plot for all species
#'
#' @param res Results object
#' @param plotname The name of a plot type (e.g., "B1P" or "R")
#'
#' @return Returns a list of plots (also displays them by default in RStudio)
#' @export
#'
#' @examples
#' \dontrun{
#' justplot(res,"R") #plots the "R" (recruitment) plots for each species
#' #To find plotnames, look at the plotnames for the species you're interested
#' in, for example;
#' names(res$Anchovy$plots)
#' #or, to see what a particular plot looks like:
#' res$Anchovy$plots$R
#' }
justplot<-function(res,plotname){
species<-names(res$species)
plist<-list()
for(i in 1:length(species)){
if(exists(plotname,where=res[[species[i]]]$plots)){
plist<-c(plist,list(res[[species[i]]]$plots[[plotname]]))
}
}
plist
}
make_panel_plot<-function(res){
scenario<-res$input_params@scenario
p1<-res$Anchovy$plots$B1P
p1$labels$title<-"Anchovy"
p2<-res$Sardine$plots$B1P
p2$labels$title<-"Sardine"
p3<-res$Other_forage$plots$B1PperB0
p3$labels$title<-"Other_forage"
p4<-res$Sealion$plots$NperK
p4$labels$title<-"Sea Lion"
p5<-res$Pelican$plots$NperK
p5$labels$title<-"Pelican"
grid.arrange(p1,p2,p3,p4,p5,ncol=2,top=paste("Scenario =",scenario))
}#splot
getStats<-function(res){
scenario<-res$input_params@scenario
nSteps<-res$input_params@nSteps
C<-res$objects$Sardine@Cy
#Get the B1P estimates for sardine and anchovy (they are by year)
BsarperB0<-res$objects$Sardine@B1PperB0
BanchperB0<-res$objects$Anchovy@B1PperB0
#Get the actual sardine biomass, B1P, which is by step, so we select only step 1
Bsar<-res$objects$Sardine@B1P
Bsar<-Bsar[seq(1,length(Bsar),nSteps)]
#Remove missing values (better to do this here, because you have to remove them from both numerator and denominator later)
Bsar<-Bsar[!is.na(Bsar)]
#And the data for predators
NperK_pel<-res$objects$Pelican@NperK
NperK_seal<-res$objects$Sealion@NperK
#remove missing values
NperK_pel<-NperK_pel[!is.na(NperK_pel)]
NperK_seal<-NperK_seal[!is.na(NperK_seal)]
#Assemble the pieces for the table
meanC<-mean(C)
pClt50<-length(C[C<50])/length(C)
meanBsarperB0<-mean(BsarperB0,na.rm=T)
meanBanchperB0<-mean(BanchperB0,na.rm=T)
pBsargt400<-length(Bsar[Bsar > 400])/length(Bsar)
pBsarlt150<-length(Bsar[Bsar < 150])/length(Bsar)
meanPelNperK<-mean(NperK_pel)
plt5_pel<-length(NperK_pel[NperK_pel < .5])/length(NperK_pel)
plt1_pel<-length(NperK_pel[NperK_pel < .1])/length(NperK_pel)
meanSealNperK<-mean(NperK_seal)
plt5_seal<-length(NperK_seal[NperK_seal < .5])/length(NperK_seal)
plt1_seal<-length(NperK_seal[NperK_seal < .1])/length(NperK_seal)
#Stick them into a 1-row data frame
nlabels<-c("C_mean","p(Catch < 50kt)","mean B/B0_sard","mean B/B0_anch","p(B_sard > 400kt)",
"p(B_sard < 150kt)","mean N/K_pel","p(N_pel < .5K)","p(N_pel < .1K)",
"mean N/K_seal","p(N_seal < .5K)","p(N_seal < .1K)")
numbers<-data.frame(meanC,pClt50,meanBsarperB0,meanBanchperB0,pBsargt400,pBsarlt150,meanPelNperK,
plt5_pel,plt1_pel,meanSealNperK,plt5_seal,plt1_seal)
stats<-data.frame(rbind(nlabels,numbers))
names(stats)<-c("C","C50","B_sar","B_anch","Bsar400","Bsar150","Npel",
"Npel5","Npel1","Nseal","Nseal5","Nseal1")
#return it
stats
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.