DEC_RevBayes_Orectolobiformes_Replicated.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/posterior_distribution_100_tree_orecto.nex\"\n",
    "out_fn = \"Output_Orectolobiformes/out_tree\" + args[1]"
   ]
  },
  {
   "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": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Attempting to read the contents of file \"posterior_distribution_100_tree_orecto.nex\"\r\n",
      "   Successfully read file\r\n"
     ]
    }
   ],
   "source": [
    "tree <- readTrees(tree_fn)[args[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_bounds <- readDataDelimitedFile(file=times_fn, delimiter=\" \")\n",
    "n_epochs <- time_bounds.size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_bg ~ dnLoguniform(1E-4,1E2)\n",
    "rate_bg.setValue(1E-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "moves.append( mvSlide(rate_bg, weight=4) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispersal_rate <- 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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": 30,
   "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": 31,
   "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": 32,
   "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": 33,
   "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": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_DEC_epoch := fnEpoch(Q=Q_DEC, times=epoch_times, rates=rep(1,n_epochs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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": 36,
   "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": 37,
   "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": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_bg.clamp(dat_range_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "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": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "mymodel = model(m_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "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           |       -120.236   |       -84.7184   |       -35.5177   |       1.435415   |           0.01   |   00:00:00   |   --:--:--   |\n",
      "100         |       -115.782   |       -80.9564   |       -34.8256   |      0.6980479   |     0.01480196   |   00:12:11   |   --:--:--   |\n",
      "200         |        -115.81   |        -80.662   |       -35.1479   |      0.6190812   |     0.01999972   |   00:25:39   |   10:15:36   |\n",
      "300         |       -116.884   |       -78.8146   |       -38.0696   |       1.522726   |      0.1070878   |   00:37:52   |   09:53:14   |\n",
      "400         |       -116.907   |       -81.8359   |       -35.0715   |      0.8241903   |     0.01772768   |   00:50:24   |   09:39:36   |\n",
      "500         |       -118.168   |       -80.6789   |       -37.4893   |       1.141099   |      0.1270821   |   01:03:07   |   09:28:03   |\n",
      "600         |       -114.797   |       -78.6695   |       -36.1274   |       1.355528   |     0.02163756   |   01:15:39   |   09:14:46   |\n",
      "700         |       -114.938   |       -78.8394   |       -36.0984   |       1.280196   |     0.02439165   |   01:28:09   |   09:01:29   |\n",
      "800         |       -116.357   |       -80.7619   |       -35.5951   |       1.017304   |     0.02350344   |   01:40:45   |   08:48:56   |\n",
      "900         |       -116.698   |       -81.7308   |        -34.967   |       0.873331   |      0.0152094   |   01:53:14   |   08:35:50   |\n",
      "1000        |        -118.92   |       -83.1335   |       -35.7867   |      0.9142904   |     0.03289816   |   02:05:47   |   08:23:08   |\n",
      "1100        |       -117.868   |       -80.9975   |       -36.8704   |       1.382779   |     0.04306166   |   02:18:10   |   08:09:51   |\n",
      "1200        |       -117.095   |       -80.8123   |       -36.2825   |       1.207362   |     0.03369188   |   02:29:59   |   07:54:56   |\n",
      "1300        |        -115.73   |       -79.9949   |       -35.7352   |       1.115305   |     0.02301342   |   02:42:24   |   07:42:12   |\n",
      "1400        |        -115.88   |       -79.9857   |       -35.8946   |       1.139133   |     0.02588555   |   02:54:20   |   07:28:17   |\n",
      "1500        |       -115.784   |       -80.1686   |       -35.6156   |      0.9157721   |     0.02767574   |   03:06:22   |   07:14:51   |\n",
      "1600        |       -117.575   |       -79.9093   |       -37.6655   |       1.205132   |      0.1348873   |   03:18:45   |   07:02:20   |\n",
      "1700        |       -116.247   |       -78.6971   |       -37.5502   |       1.470559   |     0.07100472   |   03:31:02   |   06:49:39   |\n",
      "1800        |       -118.378   |       -80.3147   |       -38.0636   |       1.438128   |      0.1268477   |   03:43:42   |   06:37:41   |\n",
      "1900        |       -119.945   |       -83.6548   |         -36.29   |       0.922622   |     0.05385098   |   03:55:45   |   06:24:38   |\n",
      "\n",
      "\n",
      "Iter        |      Posterior   |     Likelihood   |          Prior   |   extirpatio..   |        rate_bg   |    elapsed   |        ETA   |\n",
      "------------------------------------------------------------------------------------------------------------------------------------------\n",
      "2000        |       -115.618   |        -80.066   |       -35.5519   |       1.281258   |     0.01409354   |   04:08:08   |   06:12:12   |\n",
      "2100        |       -115.919   |       -80.4193   |       -35.4997   |       1.036943   |     0.02071876   |   04:20:23   |   05:59:34   |\n",
      "2200        |       -117.886   |       -79.5297   |       -38.3564   |       1.813921   |      0.0768903   |   04:32:57   |   05:47:23   |\n",
      "2300        |       -118.497   |       -79.8352   |       -38.6614   |       1.702065   |      0.1324725   |   04:45:31   |   05:35:10   |\n",
      "2400        |       -117.634   |         -80.15   |       -37.4842   |       1.379316   |     0.08011332   |   04:57:41   |   05:22:29   |\n",
      "2500        |       -115.343   |       -80.6343   |       -34.7084   |      0.7630646   |     0.01288633   |   05:10:15   |   05:10:15   |\n",
      "2600        |       -117.875   |       -79.9411   |       -37.9336   |       1.352466   |      0.1325264   |   05:22:24   |   04:57:36   |\n",
      "2700        |       -116.014   |       -80.7796   |       -35.2346   |       1.026162   |     0.01616624   |   05:34:34   |   04:45:00   |\n",
      "2800        |       -116.527   |       -80.2993   |       -36.2276   |       1.437561   |     0.02025017   |   05:46:50   |   04:32:30   |\n",
      "2900        |        -120.88   |       -84.7337   |       -36.1465   |         0.9514   |     0.04490899   |   05:59:09   |   04:20:04   |\n",
      "3000        |       -119.775   |       -83.2267   |        -36.548   |      0.7449901   |     0.08183141   |   06:11:42   |   04:07:48   |\n",
      "3100        |       -114.008   |       -79.0382   |       -34.9693   |      0.8977614   |     0.01482405   |   06:24:07   |   03:55:25   |\n",
      "3200        |       -120.989   |       -82.7225   |       -38.2668   |       1.475217   |      0.1439828   |   06:36:25   |   03:42:59   |\n",
      "3300        |       -117.218   |       -79.8919   |       -37.3257   |       1.137779   |      0.1085497   |   06:48:46   |   03:30:34   |\n",
      "3400        |       -118.935   |       -82.7669   |       -36.1678   |       1.273301   |     0.02649805   |   07:01:04   |   03:18:08   |\n",
      "3500        |       -117.109   |       -81.1065   |       -36.0025   |       1.218838   |     0.02492207   |   07:13:25   |   03:05:45   |\n",
      "3600        |       -114.606   |       -79.4413   |       -35.1643   |       1.011345   |     0.01541797   |   07:25:28   |   02:53:14   |\n",
      "3700        |       -118.669   |       -81.0778   |       -37.5907   |       1.408064   |     0.08406389   |   07:37:52   |   02:40:52   |\n",
      "3800        |       -116.521   |       -81.3514   |       -35.1698   |      0.7976907   |     0.01998694   |   07:50:26   |   02:28:33   |\n",
      "3900        |       -115.432   |       -79.4483   |       -35.9839   |       1.408772   |      0.0168327   |   08:02:41   |   02:16:08   |\n",
      "\n",
      "\n",
      "Iter        |      Posterior   |     Likelihood   |          Prior   |   extirpatio..   |        rate_bg   |    elapsed   |        ETA   |\n",
      "------------------------------------------------------------------------------------------------------------------------------------------\n",
      "4000        |       -116.033   |       -80.7475   |       -35.2857   |        0.85637   |     0.02129868   |   08:14:39   |   02:03:39   |\n",
      "4100        |       -116.604   |       -79.6477   |       -36.9561   |       1.137067   |     0.07510445   |   08:26:47   |   01:51:14   |\n",
      "4200        |       -115.155   |       -80.0052   |       -35.1496   |      0.9635445   |     0.01629701   |   08:39:10   |   01:38:53   |\n",
      "4300        |       -117.726   |       -80.0493   |       -37.6766   |       1.429169   |     0.08773695   |   08:51:43   |   01:26:33   |\n",
      "4400        |       -118.979   |       -83.9447   |       -35.0343   |      0.6795147   |     0.01824222   |   09:03:45   |   01:14:08   |\n",
      "4500        |       -116.731   |       -79.3044   |       -37.4268   |       1.187637   |      0.1097387   |   09:16:14   |   01:01:48   |\n",
      "4600        |       -116.169   |       -79.9293   |       -36.2394   |       1.144103   |     0.03621999   |   09:28:27   |   00:49:25   |\n",
      "4700        |       -116.237   |       -80.2855   |       -35.9514   |       1.337082   |     0.01882797   |   09:40:37   |   00:37:03   |\n",
      "4800        |       -119.791   |       -81.2211   |       -38.5698   |        1.89215   |     0.08050248   |   09:53:17   |   00:24:43   |\n",
      "4900        |        -115.93   |       -78.5801   |       -37.3503   |       1.424759   |     0.06388524   |   10:05:43   |   00:12:21   |\n",
      "5000        |       -119.236   |       -81.4652   |       -37.7708   |       1.449389   |     0.09248204   |   10:18:06   |   00:00:00   |\n"
     ]
    }
   ],
   "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
}