fastscapelib_cpp_design_proto

fastscapelib_cpp_design_proto.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <cstddef>\n",
    "#include <type_traits>\n",
    "\n",
    "#import \"xtensor/xtensor.hpp\"\n",
    "#import \"xtensor/xarray.hpp\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## xtensor container tags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct xtensor_tag {};\n",
    "struct xarray_tag {};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct pytensor_tag {};\n",
    "struct pyarray_tag {};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class Tag, class T, std::size_t N = 0>\n",
    "struct xt_container;\n",
    "\n",
    "template <class T, std::size_t N>\n",
    "struct xt_container<xtensor_tag, T, N>\n",
    "{\n",
    "    using type = xt::xtensor<T, N>;\n",
    "};\n",
    "\n",
    "template <class T>\n",
    "struct xt_container<xarray_tag, T>\n",
    "{\n",
    "    using type = xt::xarray<T>;\n",
    "};\n",
    "\n",
    "template <class Tag, class T, std::size_t N = 0>\n",
    "using xt_container_t = typename xt_container<Tag, T, N>::type;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid\n",
    "\n",
    "Grid classes hold grid geometry, topology and boundary conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "enum class node_status : std::uint8_t\n",
    "{\n",
    "    core = 0,\n",
    "    fixed_value_boundary = 1,\n",
    "    fixed_gradient_boundary = 2,\n",
    "    looped_boundary = 3\n",
    "};\n",
    "\n",
    "\n",
    "struct border_status\n",
    "{\n",
    "    node_status top = node_status::core;\n",
    "    node_status bottom = node_status::core;\n",
    "    node_status left = node_status::core;\n",
    "    node_status right = node_status::core;\n",
    "    \n",
    "    border_status(node_status border)\n",
    "        : top(border), right(border), bottom(border), left(border)\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    border_status(std::initializer_list<node_status> borders)\n",
    "    {\n",
    "        auto iter = borders.begin();\n",
    "        top = *iter++;\n",
    "        bottom = *iter++;\n",
    "        left = *iter++;\n",
    "        right = *iter++;\n",
    "    }\n",
    "};\n",
    "\n",
    "struct raster_node\n",
    "{\n",
    "    size_t row;\n",
    "    size_t col;\n",
    "    node_status status;\n",
    "};\n",
    "\n",
    "struct node\n",
    "{\n",
    "    size_t idx;\n",
    "    node_status status;\n",
    "};\n",
    "\n",
    "struct raster_neighbor\n",
    "{\n",
    "    size_t flatten_idx;\n",
    "    size_t row;\n",
    "    size_t col;\n",
    "    node_status status;\n",
    "    double distance;\n",
    "};\n",
    "\n",
    "struct neighbor\n",
    "{\n",
    "    size_t idx;\n",
    "    node_status status;\n",
    "    double distance;\n",
    "};\n",
    "\n",
    "enum class raster_connectivity\n",
    "{\n",
    "    king = 4,\n",
    "    queen = 8\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class container_tag> \n",
    "class raster_grid_gen\n",
    "{\n",
    "public:\n",
    "    static const std::uint8_t ndims = 2;\n",
    "    using ctag = container_tag;\n",
    "    using data_type = xt_container_t<container_tag, node_status, ndims>;\n",
    "\n",
    "    template <raster_connectivity C = raster_connectivity::queen>\n",
    "    using neighbor_list = std::array<raster_neighbor, static_cast<std::size_t>(C)>;\n",
    "    \n",
    "    template <class S>\n",
    "    raster_grid_gen(S&& shape,\n",
    "                    const std::array<double, 2>& spacing,\n",
    "                    const border_status& status_at_borders,\n",
    "                    const std::vector<raster_node>& status_at_nodes = {})\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    raster_grid_gen(std::initializer_list<std::size_t> shape,\n",
    "                    const std::array<double, 2>& spacing,\n",
    "                    const border_status& status_at_borders,\n",
    "                    const std::vector<raster_node>& status_at_nodes = {})\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    template <raster_connectivity C,\n",
    "              std::enable_if_t<C == raster_connectivity::queen, int> = 0>\n",
    "    neighbor_list<C> neighbors(size_t row, size_t col) const\n",
    "    {\n",
    "        \n",
    "    }\n",
    "    \n",
    "    template <raster_connectivity C,\n",
    "              std::enable_if_t<C == raster_connectivity::king, int> = 0> \n",
    "    neighbor_list<C> neighbors(size_t row, size_t col) const\n",
    "    {\n",
    "        \n",
    "    }\n",
    "    \n",
    "    const data_type& get_node_status() const {\n",
    "        return m_node_status;\n",
    "    }\n",
    "\n",
    "private:\n",
    "    data_type m_node_status;\n",
    "};\n",
    "\n",
    "\n",
    "using raster_grid = raster_grid_gen<xtensor_tag>;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class container_tag> \n",
    "class channel_grid_gen\n",
    "{\n",
    "public:\n",
    "    static const std::uint8_t ndims = 1;\n",
    "    using ctag = container_tag;\n",
    "    using data_type = xt_container_t<container_tag, node_status, ndims>;\n",
    "    using neighbor_list = std::array<neighbor, 2>;\n",
    "    \n",
    "    channel_grid_gen(std::size_t size,\n",
    "                     const double spacing,\n",
    "                     const std::vector<node>& status_at_nodes)\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    neighbor_list neighbors(std::size_t idx) const\n",
    "    {\n",
    "        \n",
    "    }\n",
    "    \n",
    "    const data_type& get_node_status() const {\n",
    "        return m_node_status;\n",
    "    }\n",
    "\n",
    "private:\n",
    "    data_type m_node_status;\n",
    "};\n",
    "\n",
    "\n",
    "using channel_grid = channel_grid_gen<xtensor_tag>;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "border_status borders {node_status::fixed_value_boundary,\n",
    "                       node_status::fixed_value_boundary,\n",
    "                       node_status::fixed_value_boundary,\n",
    "                       node_status::fixed_value_boundary};\n",
    "\n",
    "std::array<size_t, 2> shape {10, 10};\n",
    "\n",
    "auto rgrid = raster_grid({10, 5}, {100, 100}, borders, {});\n",
    "\n",
    "// border_status implicit constructor used here\n",
    "auto rgrid2 = raster_grid(shape, {100, 100}, node_status::fixed_value_boundary);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "auto cgrid = channel_grid(100, 10, {{0, node_status::fixed_value_boundary}});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "// const auto& nb = g.neighbors(row, col);\n",
    "// for (const auto& k : nb) {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Flow routing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sink resolvers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "enum class sink_resolve_algo\n",
    "{\n",
    "    none = 0,\n",
    "    fill_pflood,\n",
    "    fill_mst_kruskal,\n",
    "    fill_mst_boruvka,\n",
    "    fill_auto,         // either pflood or mst depending on number of sinks\n",
    "    carve_mst_kruskal,\n",
    "    carve_mst_boruvka,   \n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG>\n",
    "class sink_resolver\n",
    "{\n",
    "public:\n",
    "    using elev_type = typename FG::elev_type;\n",
    "    \n",
    "    // Entity semantic\n",
    "    virtual ~sink_resolver() = default;\n",
    "    \n",
    "    sink_resolver(const sink_resolver&) = delete;\n",
    "    sink_resolver(sink_resolver&&) = delete;\n",
    "    sink_resolver& operator=(const sink_resolver&) = delete;\n",
    "    sink_resolver& operator=(sink_resolver&&) = delete;\n",
    "    \n",
    "    virtual void resolve_before_route(const elev_type& elevation, FG& flow_graph) = 0;\n",
    "    virtual void resolve_after_route(const elev_type& elevation, FG& flow_graph) = 0;\n",
    "    \n",
    "protected:\n",
    "    sink_resolver() = default;\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG>\n",
    "class dummy_resolver : public sink_resolver<FG>\n",
    "{\n",
    "public:\n",
    "    using base_type = sink_resolver<FG>;\n",
    "    using elev_type = typename base_type::elev_type;\n",
    "    \n",
    "    dummy_resolver() = default;\n",
    "    \n",
    "    virtual ~dummy_resolver() = default;\n",
    "    \n",
    "    void resolve_before_route(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    void resolve_after_route(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG, raster_connectivity C>\n",
    "class pflood_resolver : public sink_resolver<FG>\n",
    "{\n",
    "public:\n",
    "    using base_type = sink_resolver<FG>;\n",
    "    using elev_type = typename base_type::elev_type;\n",
    "    \n",
    "    pflood_resolver() = default;\n",
    "    \n",
    "    virtual ~pflood_resolver() = default;\n",
    "    \n",
    "    void resolve_before_route(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "        // run priority flood\n",
    "    }\n",
    "    \n",
    "    void resolve_after_route(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "template<class FG, raster_connectivity C>\n",
    "sink_resolver<FG>* make_sink_resolver_impl(sink_resolve_algo algo)\n",
    "{\n",
    "    switch(algo)\n",
    "    {\n",
    "        case sink_resolve_algo::none:\n",
    "        default:\n",
    "            return new dummy_resolver<FG>();\n",
    "        \n",
    "        case sink_resolve_algo::fill_pflood:\n",
    "            return new pflood_resolver<FG, C>();\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "template<class FG>\n",
    "sink_resolver<FG>* make_sink_resolver(sink_resolve_algo algo,\n",
    "                                      raster_connectivity conn)\n",
    "{\n",
    "    switch(conn)\n",
    "    {\n",
    "    case raster_connectivity::king:\n",
    "        return make_sink_resolver_impl<FG, raster_connectivity::king>(algo);\n",
    "    case raster_connectivity::queen:\n",
    "    default:\n",
    "        return make_sink_resolver_impl<FG, raster_connectivity::queen>(algo);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow routers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "enum class flow_route_algo\n",
    "{\n",
    "    one_channel = 0,\n",
    "    single,\n",
    "    multiple,  // there are here parameters (floating point exponents)\n",
    "    single_parallel\n",
    "};\n",
    "\n",
    "struct flow_route_method\n",
    "{\n",
    "    flow_route_algo algo;\n",
    "    double param1;\n",
    "    double param2;\n",
    "    \n",
    "    flow_route_method(flow_route_algo algo, double p1 = 0., double p2 = 0.)\n",
    "        : algo(algo), param1(p1), param2(p2)\n",
    "    {\n",
    "    }\n",
    "};\n",
    "\n",
    "// for 1-d case, possible to restrict to flow_routing_algorithm::single and sink_resolving_algorithm::none?\n",
    "// in fact if there is no option here that makes sense, p_strategy will refer to the 1-d case (one class)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG>\n",
    "class flow_router\n",
    "{\n",
    "public:\n",
    "    using elev_type = typename FG::elev_type;\n",
    "    \n",
    "    // Entity semantic\n",
    "    virtual ~flow_router() = default;\n",
    "    \n",
    "    flow_router(const flow_router&) = delete;\n",
    "    flow_router(flow_router&&) = delete;\n",
    "    flow_router& operator=(const flow_router&) = delete;\n",
    "    flow_router& operator=(flow_router&&) = delete;\n",
    "    \n",
    "    void route(const elev_type& elevation,\n",
    "               FG& flow_graph,\n",
    "               sink_resolver<FG>& s_resolver)\n",
    "    {\n",
    "        s_resolver.resolve_before_route(elevation, flow_graph);\n",
    "        route_impl(elevation, flow_graph);\n",
    "        s_resolver.resolve_after_route(elevation, flow_graph);\n",
    "    }\n",
    "\n",
    "private:\n",
    "    virtual void route_impl(const elev_type& elevation, FG& flow_graph) = 0;\n",
    "    \n",
    "protected:\n",
    "    flow_router() = default;\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG>\n",
    "class one_channel_router : public flow_router<FG>\n",
    "{\n",
    "public:\n",
    "    using base_type = flow_router<FG>;\n",
    "    using elev_type = typename base_type::elev_type;\n",
    "    \n",
    "    one_channel_router() = default;\n",
    "    \n",
    "    virtual ~one_channel_router() = default;\n",
    "    \n",
    "private:\n",
    "    void route_impl(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG, raster_connectivity C>\n",
    "class single_flow_router : public flow_router<FG>\n",
    "{\n",
    "public:\n",
    "    using base_type = flow_router<FG>;\n",
    "    using elev_type = typename base_type::elev_type;\n",
    "    \n",
    "    single_flow_router() = default;\n",
    "    \n",
    "    virtual ~single_flow_router() = default;\n",
    "\n",
    "private:\n",
    "    void route_impl(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class FG, raster_connectivity C>\n",
    "class multiple_flow_router : public flow_router<FG>\n",
    "{\n",
    "public:\n",
    "    using base_type = flow_router<FG>;\n",
    "    using elev_type = typename base_type::elev_type;\n",
    "    \n",
    "    multiple_flow_router(double slope_exp, double dist_exp)\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    virtual ~multiple_flow_router() = default;\n",
    "    \n",
    "private:\n",
    "    void route_impl(const elev_type& elevation, FG& flow_graph) override\n",
    "    {\n",
    "    }\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "template<class FG, raster_connectivity C>\n",
    "flow_router<FG>* make_flow_router_impl(flow_route_method method)\n",
    "{\n",
    "    switch(method.algo)\n",
    "    {\n",
    "        case flow_route_algo::one_channel:\n",
    "            return new one_channel_router<FG>();\n",
    "        \n",
    "        case flow_route_algo::single:\n",
    "        default:\n",
    "            return new single_flow_router<FG, C>();\n",
    "        \n",
    "        case flow_route_algo::multiple:\n",
    "            return new multiple_flow_router<FG, C>(method.param1, method.param2);\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "template<class FG>\n",
    "flow_router<FG>* make_flow_router(flow_route_method method,\n",
    "                                  raster_connectivity conn)\n",
    "{\n",
    "    switch(conn)\n",
    "    {\n",
    "    case raster_connectivity::king:\n",
    "        return make_flow_router_impl<FG, raster_connectivity::king>(method);\n",
    "    \n",
    "    case raster_connectivity::queen:\n",
    "    default:\n",
    "        return make_flow_router_impl<FG, raster_connectivity::queen>(method);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "template <class G, class elev_t>\n",
    "class flow_graph\n",
    "{\n",
    "public:\n",
    "    using self_type = flow_graph<G, elev_t>;\n",
    "    using elev_type = xt_container_t<typename G::ctag, elev_t, G::ndims>;\n",
    "    \n",
    "    using flow_router_ptr = std::unique_ptr<flow_router<self_type>>;\n",
    "    using sink_resolver_ptr = std::unique_ptr<sink_resolver<self_type>>;\n",
    "    \n",
    "    // use user-defined flow_router and sink_resolver\n",
    "    flow_graph(const G& grid,\n",
    "               flow_router_ptr f_router,\n",
    "               sink_resolver_ptr s_resolver)\n",
    "        : p_flow_router(f_router), p_sink_resolver(s_resolver)\n",
    "    {\n",
    "    }\n",
    "    \n",
    "    // use built-in flow_router and sink_resolver\n",
    "    flow_graph(const G& grid,\n",
    "               flow_route_method route_method,\n",
    "               sink_resolve_algo sink_algo = sink_resolve_algo::none,\n",
    "               raster_connectivity conn = raster_connectivity::queen)\n",
    "        : p_flow_router(make_flow_router<self_type>(route_method, conn)),\n",
    "          p_sink_resolver(make_sink_resolver<self_type>(sink_algo, conn))\n",
    "    {    \n",
    "    }\n",
    "    \n",
    "    // single channel profile constructor\n",
    "    template<class PG, std::enable_if_t<std::is_same<PG, channel_grid>::value, int> = 0>\n",
    "    flow_graph(const PG& channel_grid)\n",
    "        : flow_graph(channel_grid,\n",
    "                     flow_route_algo::one_channel,\n",
    "                     sink_resolve_algo::none)\n",
    "    {\n",
    "    }\n",
    "     \n",
    "    void update(const elev_type& elevation)\n",
    "    {\n",
    "        p_flow_router->route(elevation, *this, *p_sink_resolver);\n",
    "    }\n",
    "\n",
    "private:\n",
    "    flow_router_ptr p_flow_router;\n",
    "    sink_resolver_ptr p_sink_resolver;\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "// flow_route_algo implicit constructor used here\n",
    "auto fg = flow_graph<raster_grid, double>(rgrid, flow_route_algo::single);\n",
    "\n",
    "auto mdf_method = flow_route_method(flow_route_algo::multiple, 2., 3.);\n",
    "\n",
    "auto fg3 = flow_graph<raster_grid, double>(rgrid, mdf_method);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "auto fg4 = flow_graph<channel_grid, double>(cgrid);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Erosion\n",
    "\n",
    "### Bedrock channel (stream power)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct stream_power_params\n",
    "{\n",
    "    double k_coef = 1e-7;\n",
    "    double m_exp = 0.5;\n",
    "    double n_exp = 1.;\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "template<class Er, class El, class FG>\n",
    "void compute_erosion(Er& erosion,\n",
    "                     const El& elevation,\n",
    "                     const stream_power_params& params,\n",
    "                     const FG& flow_graph);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hillslope (diffusion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct diffusion_params\n",
    "{\n",
    "    double k_coef = 1e-3;\n",
    "};"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "template<class Er, class El, class G>\n",
    "void compute_erosion(Er& erosion,\n",
    "                     const El& elevation,\n",
    "                     const diffusion_params& params,\n",
    "                     const G& grid);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "C++14 [conda env:fastscapelib]",
   "language": "C++14",
   "name": "conda-env-fastscapelib-xeus-cling-cpp14"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".cpp",
   "mimetype": "text/x-c++src",
   "name": "c++",
   "version": "-std=c++14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}