DEC_RevBayes_Orectolobiformes.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Missing Variable:\tVariable os does not exist\r\n",
      "\n"
     ]
    }
   ],
   "source": [
    "range_fn = \"Data_orectolobiformes/Orectolobiformes_range.nex\"\n",
    "tree_fn = \"Data_orectolobiformes/Orectolobiformes_extant_tree.nex\"\n",
    "out_fn = \"Output_Orectolobiformes/out\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_fn = \"Data_orectolobiformes/Orectolobiformes\"\n",
    "times_fn = geo_fn + \".times.txt\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "moves = VectorMoves()\n",
    "monitors = VectorMonitors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Successfully read one character matrix from file 'Data_orectolobiformes/Orectolobiformes_range.nex'\r\n"
     ]
    }
   ],
   "source": [
    "dat_range_01 = readDiscreteCharacterData(range_fn)\n",
    "n_areas <- dat_range_01.nchar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_areas <- 5\n",
    "n_states <- 0\n",
    "for (k in 0:max_areas) n_states += choose(n_areas, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat_range_n = formatDiscreteCharacterData(dat_range_01, \"DEC\", n_states)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_desc = dat_range_n.getStateDescriptions()\n",
    "state_desc_str = \"state,range\\n\"\n",
    "for (i in 1:state_desc.size()){\n",
    "    state_desc_str += (i-1) + \",\" + state_desc[i] + \"\\n\"\n",
    "}\n",
    "write(state_desc_str, file=out_fn+\".state_labels.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Attempting to read the contents of file \"Orectolobiformes_extant_tree.nex\"\r\n",
      "   Successfully read file\r\n"
     ]
    }
   ],
   "source": [
    "tree <- readTrees(tree_fn)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_bounds <- readDataDelimitedFile(file=times_fn, delimiter=\" \")\n",
    "n_epochs <- time_bounds.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:n_epochs) {\n",
    "      epoch_fn[i] = geo_fn + \"_connectivity_\" + i + \".txt\"\n",
    "      connectivity[i] <- readDataDelimitedFile(file=epoch_fn[i], delimiter=\" \")\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_bg ~ dnLoguniform(1E-4,1E2)\n",
    "rate_bg.setValue(1E-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "moves.append( mvSlide(rate_bg, weight=4) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersal_rate <- 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:n_epochs) {\n",
    "    for (j in 1:n_areas) {\n",
    "        for (k in 1:n_areas) {\n",
    "            dr[i][j][k] <- 0.0\n",
    "                if (connectivity[i][j][k] > 0) {\n",
    "                    dr[i][j][k]  := dispersal_rate\n",
    "                }\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_sd <- 0.5\n",
    "log_mean <- ln(1) - 0.5*log_sd^2\n",
    "extirpation_rate ~ dnLognormal(mean=log_mean, sd=log_sd)\n",
    "moves.append( mvScale(extirpation_rate, weight=2) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:n_epochs) {\n",
    "    for (j in 1:n_areas) {\n",
    "        for (k in 1:n_areas) {\n",
    "            er[i][j][k] <- 0.0\n",
    "        }\n",
    "        er[i][j][j] := extirpation_rate\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:n_epochs) {\n",
    "    Q_DEC[i] := fnDECRateMatrix(dispersalRates=dr[i],\n",
    "                extirpationRates=er[i],\n",
    "                maxRangeSize=max_areas)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (i in 1:n_epochs) {\n",
    "    time_max[i] <- time_bounds[i][1]\n",
    "    time_min[i] <- time_bounds[i][2]\n",
    "    if (i != n_epochs) {\n",
    "        epoch_times[i] ~ dnUniform(time_min[i], time_max[i])\n",
    "        moves.append( mvSlide(epoch_times[i], delta=(time_max[i]-time_min[i])/2) )\n",
    "    } else {\n",
    "        epoch_times[i] <- 0.0\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_DEC_epoch := fnEpoch(Q=Q_DEC, times=epoch_times, rates=rep(1,n_epochs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "clado_event_types <- [ \"s\", \"a\" ]\n",
    "p_sympatry ~ dnUniform(0,1)\n",
    "p_allopatry := abs(1.0 - p_sympatry)\n",
    "clado_type_probs := simplex(p_sympatry, p_allopatry)\n",
    "moves.append( mvSlide(p_sympatry, weight=2) )\n",
    "P_DEC := fnDECCladoProbs(eventProbs=clado_type_probs,\n",
    "                         eventTypes=clado_event_types,\n",
    "                         numCharacters=n_areas,\n",
    "                         maxRangeSize=max_areas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_DEC_tmp <- rep(0, n_states)\n",
    "rf_DEC_tmp[2] <- 1\n",
    "rf_DEC <- simplex(rf_DEC_tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_bg ~ dnPhyloCTMCClado(tree=tree,\n",
    "                        Q=Q_DEC_epoch,\n",
    "                        cladoProbs=P_DEC,\n",
    "                        branchRates=rate_bg,\n",
    "                        rootFrequencies=rf_DEC,\n",
    "                        type=\"NaturalNumbers\",\n",
    "                        nSites=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_bg.clamp(dat_range_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "monitors.append( mnScreen(printgen=100, rate_bg, extirpation_rate) )\n",
    "monitors.append( mnModel(file=out_fn+\".model.log\", printgen=10) )\n",
    "monitors.append( mnFile(tree, filename=out_fn+\".tre\", printgen=10) )\n",
    "monitors.append( mnJointConditionalAncestralState(tree=tree,\n",
    "                                                  ctmc=m_bg,\n",
    "                                                  type=\"NaturalNumbers\",\n",
    "                                                  withTips=true,\n",
    "                                                  withStartStates=true,\n",
    "                                                  filename=out_fn+\".states.log\",\n",
    "                                                  printgen=10) )\n",
    "\n",
    "monitors.append( mnStochasticCharacterMap(ctmc=m_bg,\n",
    "                                          filename=out_fn+\".stoch.log\",\n",
    "                                           printgen=100) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "mymodel = model(m_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   \n",
      "   Running MCMC simulation\n",
      "   This simulation runs 1 independent replicate.\n",
      "   The simulator uses 16 different moves in a random move schedule with 21 moves per iteration\n",
      "   \n",
      "\n",
      "Iter        |      Posterior   |     Likelihood   |          Prior   |   extirpatio..   |        rate_bg   |    elapsed   |        ETA   |\n",
      "------------------------------------------------------------------------------------------------------------------------------------------\n",
      "0           |       -121.629   |         -86.75   |       -34.8793   |       1.102313   |           0.01   |   00:00:00   |   --:--:--   |"
     ]
    }
   ],
   "source": [
    "mymcmc = mcmc(mymodel, moves, monitors)\n",
    "mymcmc.run(10000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RevBayes",
   "language": "bash",
   "name": "revbayes_kernel"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "r"
   },
   "file_extension": ".Rev",
   "help_links": [
    {
     "text": "RevBayes",
     "url": "https://revbayes.org"
    },
    {
     "text": "RevBayes Kernel",
     "url": "https://github.com/sdwfrost/revbayes_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-rsrc",
   "name": "RevBayes",
   "pygments_lexer": "R"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}