text:
/* Copyright 2017-2018 PaGMO development team

This file is part of the PaGMO library.

The PaGMO library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 3 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The PaGMO library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the PaGMO library.  If not,
see https://www.gnu.org/licenses/. */

#ifndef PAGMO_NOT_POPULATION_BASED_HPP
#define PAGMO_NOT_POPULATION_BASED_HPP

#include <boost/any.hpp>
#include <cassert>
#include <random>
#include <stdexcept>
#include <string>
#include <utility>

#include <pagmo/exceptions.hpp>
#include <pagmo/population.hpp>
#include <pagmo/rng.hpp>
#include <pagmo/serialization.hpp>
#include <pagmo/types.hpp>

namespace pagmo
{

/// Base class for non population based solvers
/**
 * This class contains common methods useful in the implementation of user-defined algorithms
 * that are not population based and establishes a common interface to a population
 *
 * Currently, this class implements policies for the selection and replacement of a single individual
 * in a population, which are meant to be used in the implementation of the <tt>evolve()</tt> method of the
 * user-defined algorithm (see, e.g., pagmo::nlopt::evolve()).
 */
class not_population_based
{
public:
    /// Default constructor.
    /**
     * The default constructor sets the selection and replacement policies to <tt>"best"</tt>,
     * and it initialises the RNG of the <tt>"random"</tt> selection/replacements policies via
     * random_device::next().
     */
    not_population_based()
        : m_select(std::string("best")), m_replace(std::string("best")), m_rselect_seed(random_device::next()),
          m_e(static_cast<std::mt19937::result_type>(m_rselect_seed))
    {
    }
    /// Set the seed for the ``"random"`` selection/replacement policies.
    /**
     * @param seed the value that will be used to seed the random number generator used by the ``"random"``
     * selection/replacement policies.
     */
    void set_random_sr_seed(unsigned seed)
    {
        m_rselect_seed = seed;
        m_e.seed(static_cast<std::mt19937::result_type>(m_rselect_seed));
    }
    /// Set the individual selection policy.
    /**
     * This method will set the policy that is used to select the individual
     * that will be optimised when calling the ``evolve()`` method of the algorithm.
     *
     * The input string must be one of ``"best"``, ``"worst"`` and ``"random"``:
     * - ``"best"`` will select the best individual in the population,
     * - ``"worst"`` will select the worst individual in the population,
     * - ``"random"`` will randomly choose one individual in the population.
     *
     * set_random_sr_seed() can be used to seed the random number generator used by the ``"random"`` policy.
     *
     * Instead of a selection policy, a specific individual in the population can be selected via
     * set_selection(population::size_type).
     *
     * @param select the selection policy.
     *
     * @throws std::invalid_argument if \p select is not one of ``"best"``, ``"worst"`` or ``"random"``.
     */
    void set_selection(const std::string &select)
    {
        if (select != "best" && select != "worst" && select != "random") {
            pagmo_throw(std::invalid_argument,
                        "the individual selection policy must be one of ['best', 'worst', 'random'], but '" + select
                            + "' was provided instead");
        }
        m_select = select;
    }
    /// Set the individual selection index.
    /**
     * This method will set the index of the individual that is selected for optimisation
     * in the ``evolve()`` method of the algorithm.
     *
     * @param n the index in the population of the individual to be selected for optimisation.
     */
    void set_selection(population::size_type n)
    {
        m_select = n;
    }
    /// Get the individual selection policy or index.
    /**
     * This method will return a \p boost::any containing either the individual selection policy (as an \p std::string)
     * or the individual selection index (as a population::size_type). The selection policy or index is set via
     * set_selection(const std::string &) and set_selection(population::size_type).
     *
     * @return the individual selection policy or index.
     */
    boost::any get_selection() const
    {
        return m_select;
    }
    /// Set the individual replacement policy.
    /**
     * This method will set the policy that is used in the ``evolve()`` method of the algorithm to select the individual
     * that will be replaced by the optimised individual.
     *
     * The input string must be one of ``"best"``, ``"worst"`` and ``"random"``:
     * - ``"best"`` will select the best individual in the population,
     * - ``"worst"`` will select the worst individual in the population,
     * - ``"random"`` will randomly choose one individual in the population.
     *
     * set_random_sr_seed() can be used to seed the random number generator used by the ``"random"`` policy.
     *
     * Instead of a replacement policy, a specific individual in the population can be selected via
     * set_replacement(population::size_type).
     *
     * @param replace the replacement policy.
     *
     * @throws std::invalid_argument if \p replace is not one of ``"best"``, ``"worst"`` or ``"random"``.
     */
    void set_replacement(const std::string &replace)
    {
        if (replace != "best" && replace != "worst" && replace != "random") {
            pagmo_throw(std::invalid_argument,
                        "the individual replacement policy must be one of ['best', 'worst', 'random'], but '" + replace
                            + "' was provided instead");
        }
        m_replace = replace;
    }
    /// Set the individual replacement index.
    /**
     * This method will set the index of the individual that is replaced after the optimisation
     * in the ``evolve()`` method of the algorithm.
     *
     * @param n the index in the population of the individual to be replaced after the optimisation.
     */
    void set_replacement(population::size_type n)
    {
        m_replace = n;
    }
    /// Get the individual replacement policy or index.
    /**
     * This method will return a \p boost::any containing either the individual replacement policy (as an \p
     * std::string) or the individual replacement index (as a population::size_type). The replacement policy or index is
     * set via set_replacement(const std::string &) and set_replacement(population::size_type).
     *
     * @return the individual replacement policy or index.
     */
    boost::any get_replacement() const
    {
        return m_replace;
    }
    /// Save to archive.
    /**
     * @param ar the target archive.
     *
     * @throws unspecified any exception thrown by the serialization of primitive types.
     */
    template <typename Archive>
    void save(Archive &ar) const
    {
        if (boost::any_cast<std::string>(&m_select)) {
            // NOTE: true -> string, false -> idx.
            ar(true);
            ar(boost::any_cast<std::string>(m_select));
        } else {
            ar(false);
            ar(boost::any_cast<population::size_type>(m_select));
        }
        if (boost::any_cast<std::string>(&m_replace)) {
            // NOTE: true -> string, false -> idx.
            ar(true);
            ar(boost::any_cast<std::string>(m_replace));
        } else {
            ar(false);
            ar(boost::any_cast<population::size_type>(m_replace));
        }
        ar(m_rselect_seed, m_e);
    }
    /// Load from archive.
    /**
     * In case of exceptions, \p this will be unaffected.
     *
     * @param ar the source archive.
     *
     * @throws unspecified any exception thrown by the deserialization of primitive types.
     */
    template <typename Archive>
    void load(Archive &ar)
    {
        not_population_based tmp;
        bool flag;
        std::string str;
        population::size_type idx;
        ar(flag);
        if (flag) {
            ar(str);
            tmp.m_select = str;
        } else {
            ar(idx);
            tmp.m_select = idx;
        }
        ar(flag);
        if (flag) {
            ar(str);
            tmp.m_replace = str;
        } else {
            ar(idx);
            tmp.m_replace = idx;
        }
        ar(tmp.m_rselect_seed, tmp.m_e);
        *this = std::move(tmp);
    }

protected:
    /// Select individual.
    /**
     * This method will select a single individual from the input population \p pop, returning
     * its decision vector and fitness as a pair. The selection is done according to the currently
     * active selection policy:
     * - if not_population_based::m_select is <tt>"best"</tt>, then the decision and fitness vectors of the
     *   best individual are returned,
     * - if not_population_based::m_select is <tt>"worst"</tt>, then the decision and fitness vectors of the
     *   worst individual are returned,
     * - if not_population_based::m_select is <tt>"random"</tt>, then the decision and fitness vectors of a
     *   randomly-selected individual are returned,
     * - if not_population_based::m_select is an index, then the decision and fitness vectors of the individual
     *   at the corresponding position in the population are returned.
     *
     * Note that selecting a best or worst individual is meaningful only in single-objective
     * problems.
     *
     * @param pop the input population.
     *
     * @return a pair containing the decision and fitness vectors of the selected individual.
     *
     * @throws std::invalid_argument if not_population_based::m_select is an index and the index is not smaller than the
     * size of \p pop.
     * @throws unspecified any exception thrown by population::best_idx() or population::worst_idx().
     */
    std::pair<vector_double, vector_double> select_individual(const population &pop) const
    {
        vector_double x, f;
        if (boost::any_cast<std::string>(&m_select)) {
            const auto &s_select = boost::any_cast<const std::string &>(m_select);
            if (s_select == "best") {
                x = pop.get_x()[pop.best_idx()];
                f = pop.get_f()[pop.best_idx()];
            } else if (s_select == "worst") {
                x = pop.get_x()[pop.worst_idx()];
                f = pop.get_f()[pop.worst_idx()];
            } else {
                assert(s_select == "random");
                std::uniform_int_distribution<population::size_type> dist(0, pop.size() - 1u);
                const auto idx = dist(m_e);
                x = pop.get_x()[idx];
                f = pop.get_f()[idx];
            }
        } else {
            const auto idx = boost::any_cast<population::size_type>(m_select);
            if (idx >= pop.size()) {
                pagmo_throw(std::invalid_argument, "cannot select the individual at index " + std::to_string(idx)
                                                       + ": the population has a size of only "
                                                       + std::to_string(pop.size()));
            }
            x = pop.get_x()[idx];
            f = pop.get_f()[idx];
        }
        return std::make_pair(std::move(x), std::move(f));
    }
    /// Replace individual.
    /**
     * This method will replace a single individual in the input population \p pop, setting its decision
     * vector to \p x and its fitness vector to \p f. The selection is done according to the currently
     * active selection policy:
     * - if not_population_based::m_replace is <tt>"best"</tt>, then the best individual will be replaced,
     * - if not_population_based::m_replace is <tt>"worst"</tt>, then the worst individual will be replaced,
     * - if not_population_based::m_replace is <tt>"random"</tt>, then a randomly-selected individual will be replaced,
     * - if not_population_based::m_replace is an index, then the individual at the corresponding position in the
     *   population will be replaced.
     *
     * Note that selecting a best or worst individual is meaningful only in single-objective
     * problems.
     *
     * @param pop the input population.
     * @param x the decision vector of the new individual.
     * @param f the fitness vector of the new individual.
     *
     * @throws std::invalid_argument if not_population_based::m_replace is an index and the index is not smaller than
     * the size of \p pop.
     * @throws unspecified any exception thrown by population::best_idx(), population::worst_idx(),
     * or population::set_xf().
     */
    void replace_individual(population &pop, const vector_double &x, const vector_double &f) const
    {
        if (boost::any_cast<std::string>(&m_replace)) {
            const auto &s_replace = boost::any_cast<const std::string &>(m_replace);
            if (s_replace == "best") {
                pop.set_xf(pop.best_idx(), x, f);
            } else if (s_replace == "worst") {
                pop.set_xf(pop.worst_idx(), x, f);
            } else {
                assert(s_replace == "random");
                std::uniform_int_distribution<population::size_type> dist(0, pop.size() - 1u);
                pop.set_xf(dist(m_e), x, f);
            }
        } else {
            const auto idx = boost::any_cast<population::size_type>(m_replace);
            if (idx >= pop.size()) {
                pagmo_throw(std::invalid_argument, "cannot replace the individual at index " + std::to_string(idx)
                                                       + ": the population has a size of only "
                                                       + std::to_string(pop.size()));
            }
            pop.set_xf(idx, x, f);
        }
    }

protected:
    /// Individual selection policy.
    /**
     * This \p boost::any instance must contain either an \p std::string or a population::size_type,
     * otherwise the behaviour will be undefined.
     */
    boost::any m_select;
    /// Individual replacement policy.
    /**
     * This \p boost::any instance must contain either an \p std::string or a population::size_type,
     * otherwise the behaviour will be undefined.
     */
    boost::any m_replace;
    /// Seed for the <tt>"random"</tt> selection/replacement policies.
    unsigned m_rselect_seed;
    /// Random engine for the <tt>"random"</tt> selection/replacement policies.
    mutable detail::random_engine_type m_e;
};
} // namespace pagmo

#endif

max_stars_repo_name: haisenzhao/CarpentryCompiler

max_issues_repo_name: haisenzhao/CarpentryCompiler

max_forks_repo_name: haisenzhao/CarpentryCompiler


text:
#ifdef CXX_USE_BOOST
#define CXX_USE_STD_CHRONO
#include <boost/chrono.hpp>
#define CXX_USE_STD_THREAD
#include <boost/thread.hpp>
#endif
#include <ext/async_iterator>
#include <gtest/gtest.h>

#ifdef _EXT_ASYNC_ITERATOR_
#ifdef __cpp_lambdas
TEST(async_iterator_test, int_test) {
  typedef ext::async_result<int> int_result;
  int_result res([](int_result::context &ctx) {
    ctx.begin(4);
    ctx.push(1);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
    ctx.push(2);
    ctx.push(3);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
    ctx.push(4);
    ctx.end();
    std::this_thread::sleep_for(std::chrono::seconds(1));
  });

  int j = 1;
  CXX_FOR(auto &i, res) {
    EXPECT_EQ(j++, i);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
  }

  j = 1;
  CXX_FOR(auto &i, res) {
    EXPECT_EQ(j++, i);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
  }

  CXX_FOR(auto &i, res) { EXPECT_TRUE(false); }

  res = int_result([](int_result::context &ctx) {
    ctx.begin(4);
    ctx.push(5);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
    ctx.push(7);
    ctx.push(9);
    std::this_thread::sleep_for(std::chrono::milliseconds(200));
    ctx.push(11);
    ctx.push(13);
    std::this_thread::sleep_for(std::chrono::seconds(1));
  });
  EXPECT_EQ(res.size(), 5);

  j = 3;
  CXX_FOR(auto &i, res) {
    EXPECT_EQ(j += 2, i);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
  }
}

TEST(async_iterator_test, int_thread_test) {
  typedef ext::async_result<int> int_result;
  const int count = 1000;
  int_result res([count](int_result::context &ctx) {
    ctx.begin(count);
    for (int i = 0; i < count; ++i)
      ctx.push(i);
  });

  std::thread t([&res]() {
    CXX_FOR(auto &i, res) { EXPECT_LT(i, (int)res.size()); }
  });
  std::thread t2([&res]() {
    CXX_FOR(auto &i, res) { EXPECT_LT(i, (int)res.size()); }
  });
  if (t.joinable())
    t.join();
  if (t2.joinable())
    t2.join();
}

TEST(async_iterator_test, int_thread_cancelable_test) {
  typedef ext::async_result<int> int_result;
  const int count = 1000000;
  int processed = 0;
  int_result res([count, &processed](int_result::context &ctx) {
    ctx.begin(count);
    for (int i = 0; i < count; ++i) {
      if (ctx.cancel_requested())
        break;
      processed = i;
      ctx.push(i);
    }
  });

  std::thread t0([&res, &processed]() {
    std::this_thread::sleep_for(std::chrono::microseconds(500));
    res.cancel();
  });
  std::thread t1([&res, &processed]() {
    CXX_FOR(auto &i, res) {
      EXPECT_LE(i, processed);
      EXPECT_LT(i, (int)res.size());
    }
  });
  std::thread t2([&res, &processed]() {
    CXX_FOR(auto &i, res) {
      EXPECT_LE(i, processed);
      EXPECT_LT(i, (int)res.size());
    }
  });

  if (t0.joinable())
    t0.join();
  if (t1.joinable())
    t1.join();
  if (t2.joinable())
    t2.join();
}

TEST(async_iterator_test, string_test) {
  typedef ext::async_result<std::string> str_result;
  str_result res([](str_result::context &ctx) {
    ctx.push("str 1");
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
    ctx.push("str 2");
    ctx.push("str 3");
    std::this_thread::sleep_for(std::chrono::milliseconds(500));
    ctx.push("str 4");
  });

  int j = 1;
  CXX_FOR(auto &s, res) {
    EXPECT_EQ(s, "str " + std::to_string((long long)j++));
  }
}

TEST(async_iterator_test, pair_test) {
  typedef ext::async_result<std::pair<std::string, size_t>> pair_result;
  pair_result res([](pair_result::context &ctx) {
    typedef ext::async_result<std::pair<std::string, size_t>> pair_result;
#ifdef __cpp_variadic_templates
    ctx.emplace("str 1", 1);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
#endif
#ifdef __cpp_aggregate_nsdmi
    ctx.push({"str 2", 2});
#else
    ctx.push(pair_result::type("str 2", 2));
#endif
#ifdef __cpp_variadic_templates
    ctx.emplace("str 3", 3);
    std::this_thread::sleep_for(std::chrono::milliseconds(500));
#endif
#ifdef __cpp_aggregate_nsdmi
    ctx.push({"str 4", 4});
#else
    ctx.push(pair_result::type("str 4", 4));
#endif
  });

  int j = 0;
  CXX_FOR(auto &s, res) {
    EXPECT_EQ(s.first, "str " + std::to_string((long long)s.second));
    EXPECT_EQ(s.second, s.second);
  }
}
#else
typedef ext::async_result<int> int_result;
void int_test_phase0(int_result::context &ctx) {
  ctx.begin(4);
  ctx.push(1);
  std::this_thread::sleep_for(std::chrono::milliseconds(500));
  ctx.push(2);
  ctx.push(3);
  std::this_thread::sleep_for(std::chrono::milliseconds(500));
  ctx.push(4);
  ctx.end();
  std::this_thread::sleep_for(std::chrono::seconds(3));
}
void int_test_phase1(int_result::context &ctx) {
  ctx.begin(4);
  ctx.push(5);
  std::this_thread::sleep_for(std::chrono::milliseconds(100));
  ctx.push(6);
  ctx.push(7);
  std::this_thread::sleep_for(std::chrono::milliseconds(200));
  ctx.push(8);
  ctx.push(9);
  std::this_thread::sleep_for(std::chrono::seconds(3));
}

TEST(async_iterator_test, int_test) {
  typedef ext::async_result<int> int_result;
  int_result res = int_test_phase0;

  int j = 1;
  CXX_FOR(int i, res) {
    EXPECT_EQ(j++, i);
    std::this_thread::sleep_for(std::chrono::seconds(1));
  }

  j = 1;
  CXX_FOR(int i, res) {
    EXPECT_EQ(j++, i);
    std::this_thread::sleep_for(std::chrono::seconds(1));
  }

  CXX_FOR(int &i, res) { EXPECT_TRUE(false); }

  res = int_result(int_test_phase1);
  EXPECT_EQ((int)res.size(), 5);

  j = 5;
  CXX_FOR(int &i, res) {
    EXPECT_EQ(j++, i);
    std::this_thread::sleep_for(std::chrono::milliseconds(100));
  }
}
#endif // __cpp_lambdas
#endif // _EXT_ASYNC_ITERATOR_

ext: cpp


text:
#include <boost/python.hpp>
using namespace boost::python;

BOOST_PYTHON_MODULE(getting_started1)
{

    // Add regular functions to the module.
    def("greet", greet);
    def("square", square);
}

max_stars_repo_path: 1. Hello-World/getting_started1.cpp

max_stars_repo_name: adityagupta-by/boost-examples

max_issues_repo_path: 1. Hello-World/getting_started1.cpp

max_issues_repo_name: adityagupta-by/boost-examples

max_forks_repo_path: 1. Hello-World/getting_started1.cpp

max_forks_repo_name: adityagupta-by/boost-examples


text:
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
// This file is part of the Matrix Template Library
//
// See also license.mtl.txt in the distribution.

#ifndef MTL_STRING_TO_ENUM_INCLUDE
#define MTL_STRING_TO_ENUM_INCLUDE

#include <string>
#include <boost/numeric/mtl/operation/size.hpp>

namespace mtl {

/** Searches string \p s in list \p l of strings and returns enum

    List \p l is given as array of const char*, which is the easiest to
    initialize.  The search is case sensitive, thus (de)-capitalize
    your string upfront, e.g., with boost::to_lower().
    If the string is not found then a runtime_error is thrown.
**/
template <typename EnumType, typename Array>
EnumType inline string_to_enum(const std::string& s, const Array& l, EnumType)
{
    std::size_t i;
    for (i= 0; i < size(l) && std::string(l[i]) != s; i++) {}
    MTL_THROW_IF(i == size(l), runtime_error("Search string not found"));
    return EnumType(i);
}

} // namespace mtl

#endif // MTL_STRING_TO_ENUM_INCLUDE

ext: hpp


text:
/*
 * This file is part of a C++ interface to the Radiative Transfer for Energetics (RTE)
 * and Rapid Radiative Transfer Model for GCM applications Parallel (RRTMGP).
 *
 * The original code is found at https://github.com/RobertPincus/rte-rrtmgp.
 *
 * Contacts: Robert Pincus and Eli Mlawer
 * email: rrtmgp@aer.com
 *
 * Copyright 2015-2020,  Atmospheric and Environmental Research and
 * Regents of the University of Colorado.  All right reserved.
 *
 * This C++ interface can be downloaded from https://github.com/Chiil/rrtmgp_cpp
 *
 * Contact: Chiel van Heerwaarden
 * email: chiel.vanheerwaarden@wur.nl
 *
 * Copyright 2020, Wageningen University & Research.
 *
 * Use and duplication is permitted under the terms of the
 * BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
 *
 */

#include <cmath>
#include <numeric>
#include <boost/algorithm/string.hpp>
// #include <xtensor/xarray.hpp>
// #include <xtensor/xadapt.hpp>
// #include "xtensor/xnoalias.hpp"

#include "Gas_concs.h"
#include "Gas_optics.h"
#include "Array.h"
#include "Optical_props.h"
#include "Source_functions.h"

#include "rrtmgp_kernels.h"
#define restrict __restrict__

namespace
{
    int find_index(
            const Array<std::string,1>& data, const std::string& value)
    {
        auto it = std::find(data.v().begin(), data.v().end(), value);
        if (it == data.v().end())
            return -1;
        else
            return it - data.v().begin() + 1;
    }

    template<typename TF>
    void reduce_minor_arrays(
                const Gas_concs<TF>& available_gases,
                const Array<std::string,1>& gas_names,
                const Array<std::string,1>& gas_minor,
                const Array<std::string,1>& identifier_minor,
                const Array<TF,3>& kminor_atm,
                const Array<std::string,1>& minor_gases_atm,
                const Array<int,2>& minor_limits_gpt_atm,
                const Array<int,1>& minor_scales_with_density_atm, // CvH: logical bool or int?
                const Array<std::string,1>& scaling_gas_atm,
                const Array<int,1>& scale_by_complement_atm, // CvH: bool or int
                const Array<int,1>& kminor_start_atm,

                Array<TF,3>& kminor_atm_red,
                Array<std::string,1>& minor_gases_atm_red,
                Array<int,2>& minor_limits_gpt_atm_red,
                Array<int,1>& minor_scales_with_density_atm_red, // CvH bool or int
                Array<std::string,1>& scaling_gas_atm_red,
                Array<int,1>& scale_by_complement_atm_red,
                Array<int,1>& kminor_start_atm_red)
    {
        int nm = minor_gases_atm.dim(1);
        int tot_g = 0;

        Array<int,1> gas_is_present({nm});

        for (int i=1; i<=nm; ++i)
        {
            const int idx_mnr = find_index(identifier_minor, minor_gases_atm({i}));

            // Search for
            std::string gas_minor_trimmed = gas_minor({idx_mnr});
            boost::trim(gas_minor_trimmed);

            gas_is_present({i}) = available_gases.exists(gas_minor_trimmed);
            if (gas_is_present({i}))
                tot_g += minor_limits_gpt_atm({2,i}) - minor_limits_gpt_atm({1,i}) + 1;
        }

        const int red_nm = std::accumulate(gas_is_present.v().begin(), gas_is_present.v().end(), 0);

        if (red_nm == nm)
        {
            kminor_atm_red = kminor_atm;
            minor_gases_atm_red = minor_gases_atm;
            minor_limits_gpt_atm_red = minor_limits_gpt_atm;
            minor_scales_with_density_atm_red = minor_scales_with_density_atm;
            scaling_gas_atm_red = scaling_gas_atm;
            scale_by_complement_atm_red = scale_by_complement_atm;
            kminor_start_atm_red = kminor_start_atm;
        }
        else
        {
            // Use a lambda function as the operation has to be repeated many times.
            auto resize_and_set = [&](auto& a_red, const auto& a)
            {
                a_red.set_dims({red_nm});
                int counter = 1;
                for (int i=1; i<=gas_is_present.dim(1); ++i)
                {
                    if (gas_is_present({i}))
                    {
                       a_red({counter}) = a({i});
                       ++counter;
                    }
                }
            };

            resize_and_set(minor_gases_atm_red, minor_gases_atm);
            resize_and_set(minor_scales_with_density_atm_red, minor_scales_with_density_atm);
            resize_and_set(scaling_gas_atm_red, scaling_gas_atm);
            resize_and_set(scale_by_complement_atm_red, scale_by_complement_atm);
            resize_and_set(kminor_start_atm_red, kminor_start_atm);

            minor_limits_gpt_atm_red.set_dims({2, red_nm});
            kminor_atm_red.set_dims({tot_g, kminor_atm.dim(2), kminor_atm.dim(3)});

            int icnt = 0;
            int n_elim = 0;
            for (int i=1; i<=nm; ++i)
            {
                int ng = minor_limits_gpt_atm({2,i}) - minor_limits_gpt_atm({1,i}) + 1;
                if (gas_is_present({i}))
                {
                    ++icnt;
                    minor_limits_gpt_atm_red({1,icnt}) = minor_limits_gpt_atm({1,i});
                    minor_limits_gpt_atm_red({2,icnt}) = minor_limits_gpt_atm({2,i});
                    kminor_start_atm_red({icnt}) = kminor_start_atm({i}) - n_elim;

                    for (int j=1; j<=ng; ++j)
                        for (int i2=1; i2<=kminor_atm.dim(2); ++i2)
                            for (int i3=1; i3<=kminor_atm.dim(3); ++i3)
                                kminor_atm_red({kminor_start_atm_red({icnt})+j-1,i2,i3}) =
                                        kminor_atm({kminor_start_atm({i})+j-1,i2,i3});
                }
                else
                    n_elim += ng;
            }
        }
    }

    void create_idx_minor(
            const Array<std::string,1>& gas_names,
            const Array<std::string,1>& gas_minor,
            const Array<std::string,1>& identifier_minor,
            const Array<std::string,1>& minor_gases_atm,
            Array<int,1>& idx_minor_atm)
    {
        Array<int,1> idx_minor_atm_out({minor_gases_atm.dim(1)});

        for (int imnr=1; imnr<=minor_gases_atm.dim(1); ++imnr)
        {
            // Find identifying string for minor species in list of possible identifiers (e.g. h2o_slf)
            const int idx_mnr = find_index(identifier_minor, minor_gases_atm({imnr}));

            // Find name of gas associated with minor species identifier (e.g. h2o)
            idx_minor_atm_out({imnr}) = find_index(gas_names, gas_minor({idx_mnr}));
        }

        idx_minor_atm = idx_minor_atm_out;
    }

    void create_idx_minor_scaling(
            const Array<std::string,1>& gas_names,
            const Array<std::string,1>& scaling_gas_atm,
            Array<int,1>& idx_minor_scaling_atm)
    {
        Array<int,1> idx_minor_scaling_atm_out({scaling_gas_atm.dim(1)});

        for (int imnr=1; imnr<=scaling_gas_atm.dim(1); ++imnr)
            idx_minor_scaling_atm_out({imnr}) = find_index(gas_names, scaling_gas_atm({imnr}));

        idx_minor_scaling_atm = idx_minor_scaling_atm_out;
    }

    void create_key_species_reduce(
            const Array<std::string,1>& gas_names,
            const Array<std::string,1>& gas_names_red,
            const Array<int,3>& key_species,
            Array<int,3>& key_species_red,
            Array<int,1>& key_species_present_init)
    {
        const int np = key_species.dim(1);
        const int na = key_species.dim(2);
        const int nt = key_species.dim(3);

        key_species_red.set_dims({key_species.dim(1), key_species.dim(2), key_species.dim(3)});
        key_species_present_init.set_dims({gas_names.dim(1)});

        for (int i=1; i<=key_species_present_init.dim(1); ++i)
            key_species_present_init({i}) = 1;

        for (int ip=1; ip<=np; ++ip)
            for (int ia=1; ia<=na; ++ia)
                for (int it=1; it<=nt; ++it)
                {
                    const int ks = key_species({ip,ia,it});
                    if (ks != 0)
                    {
                        const int ksr = find_index(gas_names_red, gas_names({ks}));
                        key_species_red({ip,ia,it}) = ksr;
                        if (ksr == -1)
                            key_species_present_init({ks}) = 0;
                    }
                    else
                        key_species_red({ip,ia,it}) = ks;
                }
    }

    void check_key_species_present_init(
            const Array<std::string,1>& gas_names,
            const Array<int,1>& key_species_present_init
            )
    {
        for (int i=1; i<=key_species_present_init.dim(1); ++i)
        {
            if (key_species_present_init({i}) == 0)
            {
                std::string error_message = "Gas optics: required gas " + gas_names({i}) + " is missing";
                throw std::runtime_error(error_message);
            }
        }
    }

    void create_flavor(
            const Array<int,3>& key_species,
            Array<int,2>& flavor)
    {
        Array<int,2> key_species_list({2, key_species.dim(3)*2});

        // Prepare list of key species.
        int i = 1;
        for (int ibnd=1; ibnd<=key_species.dim(3); ++ibnd)
            for (int iatm=1; iatm<=key_species.dim(1); ++iatm)
            {
                key_species_list({1,i}) = key_species({1,iatm,ibnd});
                key_species_list({2,i}) = key_species({2,iatm,ibnd});
                ++i;
            }

        // Rewrite single key_species pairs.
        for (int i=1; i<=key_species_list.dim(2); ++i)
        {
            if ( key_species_list({1,i}) == 0 && key_species_list({2,i}) == 0 )
            {
                key_species_list({1,i}) = 2;
                key_species_list({2,i}) = 2;
            }
        }

        // Count unique key species pairs.
        int iflavor = 0;
        for (int i=1; i<=key_species_list.dim(2); ++i)
        {
            bool pair_exists = false;
            for (int ii=1; ii<=i-1; ++ii)
            {
                if ( (key_species_list({1,i}) == key_species_list({1,ii})) &&
                     (key_species_list({2,i}) == key_species_list({2,ii})) )
                {
                    pair_exists = true;
                    break;
                }
            }
            if (!pair_exists)
                ++iflavor;
        }

        // Fill flavors.
        flavor.set_dims({2,iflavor});
        iflavor = 0;
        for (int i=1; i<=key_species_list.dim(2); ++i)
        {
            bool pair_exists = false;
            for (int ii=1; ii<=i-1; ++ii)
            {
                if ( (key_species_list({1,i}) == key_species_list({1,ii})) &&
                     (key_species_list({2,i}) == key_species_list({2,ii})) )
                {
                    pair_exists = true;
                    break;
                }
            }
            if (!pair_exists)
            {
                ++iflavor;
                flavor({1,iflavor}) = key_species_list({1,i});
                flavor({2,iflavor}) = key_species_list({2,i});
            }
        }
    }

    int key_species_pair2flavor(
            const Array<int,2>& flavor,
            const Array<int,1>& key_species_pair)
    {
        // Search for match.
        for (int iflav=1; iflav<=flavor.dim(2); ++iflav)
        {
            if ( key_species_pair({1}) == flavor({1, iflav}) &&
                 key_species_pair({2}) == flavor({2, iflav}) )
                return iflav;
        }

        // No match found.
        return -1;
    }

    void create_gpoint_flavor(
            const Array<int,3>& key_species,
            const Array<int,1>& gpt2band,
            const Array<int,2>& flavor,
            Array<int,2>& gpoint_flavor)
    {
        const int ngpt = gpt2band.dim(1);
        gpoint_flavor.set_dims({2,ngpt});

        for (int igpt=1; igpt<=ngpt; ++igpt)
            for (int iatm=1; iatm<=2; ++iatm)
            {
                int pair_1 = key_species( {1, iatm, gpt2band({igpt})} );
                int pair_2 = key_species( {2, iatm, gpt2band({igpt})} );

                // Rewrite species pair.
                Array<int,1> rewritten_pair({2});
                if (pair_1 == 0 && pair_2 == 0)
                {
                    rewritten_pair({1}) = 2;
                    rewritten_pair({2}) = 2;
                }
                else
                {
                    rewritten_pair({1}) = pair_1;
                    rewritten_pair({2}) = pair_2;
                }

                // Write the output.
                gpoint_flavor({iatm,igpt}) = key_species_pair2flavor(
                        flavor, rewritten_pair);
            }
    }

    /*
    template<typename TF>
    inline void reorder123x321_test(
            TF* restrict out, const TF* restrict in,
            const int d1, const int d2, const int d3)
    {
        const int jj_in = d1;
        const int kk_in = d1*d2;

        const int ii_out = d3*d2;
        const int jj_out = d3;

        for (int i=0; i<d1; ++i)
            for (int j=0; j<d2; ++j)
                #pragma GCC ivdep
                for (int k=0; k<d3; ++k)
                {
                    const int ijk_in  = i + j*jj_in  + k*kk_in ;
                    const int ijk_out = k + j*jj_out + i*ii_out;
                    out[ijk_out] = in[ijk_in];
                }
    }

    template<typename TF>
    inline void reorder123x321_test(
            TF* out, const TF* in,
            const size_t d1, const size_t d2, const size_t d3)
    {
        const size_t size = d1*d2*d3;
        const std::array<size_t, 3> in_shape = { d1, d2, d3 };
        const std::array<size_t, 3> out_shape = { d3, d2, d1 };

        const auto a_in = xt::adapt<xt::layout_type::column_major>(in, size, xt::no_ownership(), in_shape);
        auto a_out = xt::adapt<xt::layout_type::column_major>(out, size, xt::no_ownership(), out_shape);

        xt::noalias(a_out) = xt::transpose(a_out);
    }
    */
}

// IMPLEMENTATION OF CLASS FUNCTIONS.
// Constructor of longwave variant.
template<typename TF>
Gas_optics<TF>::Gas_optics(
        const Gas_concs<TF>& available_gases,
        const Array<std::string,1>& gas_names,
        const Array<int,3>& key_species,
        const Array<int,2>& band2gpt,
        const Array<TF,2>& band_lims_wavenum,
        const Array<TF,1>& press_ref,
        const TF press_ref_trop,
        const Array<TF,1>& temp_ref,
        const TF temp_ref_p,
        const TF temp_ref_t,
        const Array<TF,3>& vmr_ref,
        const Array<TF,4>& kmajor,
        const Array<TF,3>& kminor_lower,
        const Array<TF,3>& kminor_upper,
        const Array<std::string,1>& gas_minor,
        const Array<std::string,1>& identifier_minor,
        const Array<std::string,1>& minor_gases_lower,
        const Array<std::string,1>& minor_gases_upper,
        const Array<int,2>& minor_limits_gpt_lower,
        const Array<int,2>& minor_limits_gpt_upper,
        const Array<int,1>& minor_scales_with_density_lower,
        const Array<int,1>& minor_scales_with_density_upper,
        const Array<std::string,1>& scaling_gas_lower,
        const Array<std::string,1>& scaling_gas_upper,
        const Array<int,1>& scale_by_complement_lower,
        const Array<int,1>& scale_by_complement_upper,
        const Array<int,1>& kminor_start_lower,
        const Array<int,1>& kminor_start_upper,
        const Array<TF,2>& totplnk,
        const Array<TF,4>& planck_frac,
        const Array<TF,3>& rayl_lower,
        const Array<TF,3>& rayl_upper) :
            Optical_props<TF>(band_lims_wavenum, band2gpt),
            totplnk(totplnk),
            planck_frac(planck_frac)
{
    // Initialize the absorption coefficient array, including Rayleigh scattering
    // tables if provided.
    init_abs_coeffs(
            available_gases,
            gas_names, key_species,
            band2gpt, band_lims_wavenum,
            press_ref, temp_ref,
            press_ref_trop, temp_ref_p, temp_ref_t,
            vmr_ref,
            kmajor, kminor_lower, kminor_upper,
            gas_minor,identifier_minor,
            minor_gases_lower, minor_gases_upper,
            minor_limits_gpt_lower,
            minor_limits_gpt_upper,
            minor_scales_with_density_lower,
            minor_scales_with_density_upper,
            scaling_gas_lower, scaling_gas_upper,
            scale_by_complement_lower,
            scale_by_complement_upper,
            kminor_start_lower,
            kminor_start_upper,
            rayl_lower, rayl_upper);

    // Temperature steps for Planck function interpolation.
    // Assumes that temperature minimum and max are the same for the absorption coefficient grid and the
    // Planck grid and the Planck grid is equally spaced.
    totplnk_delta = (temp_ref_max - temp_ref_min) / (totplnk.dim(1)-1);
}

// Constructor of the shortwave variant.
template<typename TF>
Gas_optics<TF>::Gas_optics(
        const Gas_concs<TF>& available_gases,
        const Array<std::string,1>& gas_names,
        const Array<int,3>& key_species,
        const Array<int,2>& band2gpt,
        const Array<TF,2>& band_lims_wavenum,
        const Array<TF,1>& press_ref,
        const TF press_ref_trop,
        const Array<TF,1>& temp_ref,
        const TF temp_ref_p,
        const TF temp_ref_t,
        const Array<TF,3>& vmr_ref,
        const Array<TF,4>& kmajor,
        const Array<TF,3>& kminor_lower,
        const Array<TF,3>& kminor_upper,
        const Array<std::string,1>& gas_minor,
        const Array<std::string,1>& identifier_minor,
        const Array<std::string,1>& minor_gases_lower,
        const Array<std::string,1>& minor_gases_upper,
        const Array<int,2>& minor_limits_gpt_lower,
        const Array<int,2>& minor_limits_gpt_upper,
        const Array<int,1>& minor_scales_with_density_lower,
        const Array<int,1>& minor_scales_with_density_upper,
        const Array<std::string,1>& scaling_gas_lower,
        const Array<std::string,1>& scaling_gas_upper,
        const Array<int,1>& scale_by_complement_lower,
        const Array<int,1>& scale_by_complement_upper,
        const Array<int,1>& kminor_start_lower,
        const Array<int,1>& kminor_start_upper,
        const Array<TF,1>& solar_source_quiet,
        const Array<TF,1>& solar_source_facular,
        const Array<TF,1>& solar_source_sunspot,
        const TF tsi_default,
        const TF mg_default,
        const TF sb_default,
        const Array<TF,3>& rayl_lower,
        const Array<TF,3>& rayl_upper) :
            Optical_props<TF>(band_lims_wavenum, band2gpt)
{
    // Initialize the absorption coefficient array, including Rayleigh scattering
    // tables if provided.
    init_abs_coeffs(
            available_gases,
            gas_names, key_species,
            band2gpt, band_lims_wavenum,
            press_ref, temp_ref,
            press_ref_trop, temp_ref_p, temp_ref_t,
            vmr_ref,
            kmajor, kminor_lower, kminor_upper,
            gas_minor,identifier_minor,
            minor_gases_lower, minor_gases_upper,
            minor_limits_gpt_lower,
            minor_limits_gpt_upper,
            minor_scales_with_density_lower,
            minor_scales_with_density_upper,
            scaling_gas_lower, scaling_gas_upper,
            scale_by_complement_lower,
            scale_by_complement_upper,
            kminor_start_lower,
            kminor_start_upper,
            rayl_lower, rayl_upper);

    // Compute the solar source.
    this->solar_source_quiet = solar_source_quiet;
    this->solar_source_facular = solar_source_facular;
    this->solar_source_sunspot = solar_source_sunspot;

    this->solar_source.set_dims(solar_source_quiet.get_dims());

    set_solar_variability(mg_default, sb_default);
}

template<typename TF>
void Gas_optics<TF>::init_abs_coeffs(
        const Gas_concs<TF>& available_gases,
        const Array<std::string,1>& gas_names,
        const Array<int,3>& key_species,
        const Array<int,2>& band2gpt,
        const Array<TF,2>& band_lims_wavenum,
        const Array<TF,1>& press_ref,
        const Array<TF,1>& temp_ref,
        const TF press_ref_trop,
        const TF temp_ref_p,
        const TF temp_ref_t,
        const Array<TF,3>& vmr_ref,
        const Array<TF,4>& kmajor,
        const Array<TF,3>& kminor_lower,
        const Array<TF,3>& kminor_upper,
        const Array<std::string,1>& gas_minor,
        const Array<std::string,1>& identifier_minor,
        const Array<std::string,1>& minor_gases_lower,
        const Array<std::string,1>& minor_gases_upper,
        const Array<int,2>& minor_limits_gpt_lower,
        const Array<int,2>& minor_limits_gpt_upper,
        const Array<int,1>& minor_scales_with_density_lower,
        const Array<int,1>& minor_scales_with_density_upper,
        const Array<std::string,1>& scaling_gas_lower,
        const Array<std::string,1>& scaling_gas_upper,
        const Array<int,1>& scale_by_complement_lower,
        const Array<int,1>& scale_by_complement_upper,
        const Array<int,1>& kminor_start_lower,
        const Array<int,1>& kminor_start_upper,
        const Array<TF,3>& rayl_lower,
        const Array<TF,3>& rayl_upper)
{
    // Which gases known to the gas optics are present in the host model (available_gases)?
    std::vector<std::string> gas_names_to_use;

    for (const std::string &s : gas_names.v())
    {
        if (available_gases.exists(s))
            gas_names_to_use.push_back(s);
    }

    // Now the number of gases is the union of those known to the k-distribution and provided
    // by the host model.
    const int n_gas = gas_names_to_use.size();
    Array<std::string,1> gas_names_this(std::move(gas_names_to_use), {n_gas});
    this->gas_names = gas_names_this;

    // Initialize the gas optics object, keeping only those gases known to the
    // gas optics and also present in the host model.
    // Add an offset to the indexing to interface the negative ranging of fortran.
    Array<TF,3> vmr_ref_red({vmr_ref.dim(1), n_gas + 1, vmr_ref.dim(3)});
    vmr_ref_red.set_offsets({0, -1, 0});

    // Gas 0 is used in single-key species method, set to 1.0 (col_dry)
    for (int i1=1; i1<=vmr_ref_red.dim(1); ++i1)
        for (int i3=1; i3<=vmr_ref_red.dim(3); ++i3)
            vmr_ref_red({i1, 0, i3}) = vmr_ref({i1, 1, i3});

    for (int i=1; i<=n_gas; ++i)
    {
        int idx = find_index(gas_names, this->gas_names({i}));
        for (int i1=1; i1<=vmr_ref_red.dim(1); ++i1)
            for (int i3=1; i3<=vmr_ref_red.dim(3); ++i3)
                vmr_ref_red({i1, i, i3}) = vmr_ref({i1, idx+1, i3}); // CvH: why +1?
    }

    this->vmr_ref = std::move(vmr_ref_red);

    // Reduce minor arrays so variables only contain minor gases that are available.
    // Reduce size of minor Arrays.
    Array<std::string, 1> minor_gases_lower_red;
    Array<std::string, 1> scaling_gas_lower_red;
    Array<std::string, 1> minor_gases_upper_red;
    Array<std::string, 1> scaling_gas_upper_red;

    reduce_minor_arrays(
            available_gases,
            gas_names,
            gas_minor, identifier_minor,
            kminor_lower,
            minor_gases_lower,
            minor_limits_gpt_lower,
            minor_scales_with_density_lower,
            scaling_gas_lower,
            scale_by_complement_lower,
            kminor_start_lower,
            this->kminor_lower,
            minor_gases_lower_red,
            this->minor_limits_gpt_lower,
            this->minor_scales_with_density_lower,
            scaling_gas_lower_red,
            this->scale_by_complement_lower,
            this->kminor_start_lower);

    reduce_minor_arrays(
            available_gases,
            gas_names,
            gas_minor,
            identifier_minor,
            kminor_upper,
            minor_gases_upper,
            minor_limits_gpt_upper,
            minor_scales_with_density_upper,
            scaling_gas_upper,
            scale_by_complement_upper,
            kminor_start_upper,
            this->kminor_upper,
            minor_gases_upper_red,
            this->minor_limits_gpt_upper,
            this->minor_scales_with_density_upper,
            scaling_gas_upper_red,
            this->scale_by_complement_upper,
            this->kminor_start_upper);

    // Arrays not reduced by the presence, or lack thereof, of a gas
    this->press_ref = press_ref;
    this->temp_ref = temp_ref;
    this->kmajor = kmajor;

    // Create a new vector that consists of rayl_lower and rayl_upper stored in one variable.
    if (rayl_lower.size() > 0)
    {
        this->krayl.set_dims({rayl_lower.dim(1), rayl_lower.dim(2), rayl_lower.dim(3), 2});
        for (int i=0; i<rayl_lower.size(); ++i)
        {
            this->krayl.v()[i                    ] = rayl_lower.v()[i];
            this->krayl.v()[i + rayl_lower.size()] = rayl_upper.v()[i];
        }
    }

    // ---- post processing ----
    //  creates log reference pressure
    this->press_ref_log = this->press_ref;
    for (int i1=1; i1<=this->press_ref_log.dim(1); ++i1)
        this->press_ref_log({i1}) = std::log(this->press_ref_log({i1}));

    // log scale of reference pressure
    this->press_ref_trop_log = std::log(press_ref_trop);

    // Get index of gas (if present) for determining col_gas
    create_idx_minor(
            this->gas_names, gas_minor, identifier_minor, minor_gases_lower_red, this->idx_minor_lower);
    create_idx_minor(
            this->gas_names, gas_minor, identifier_minor, minor_gases_upper_red, this->idx_minor_upper);

    // Get index of gas (if present) that has special treatment in density scaling
    create_idx_minor_scaling(
            this->gas_names, scaling_gas_lower_red, this->idx_minor_scaling_lower);
    create_idx_minor_scaling(
            this->gas_names, scaling_gas_upper_red, this->idx_minor_scaling_upper);

    // Create flavor list.
    // Reduce (remap) key_species list; checks that all key gases are present in incoming
    Array<int, 3> key_species_red;
    Array<int, 1> key_species_present_init; // CvH bool or int?

    create_key_species_reduce(
            gas_names, this->gas_names, key_species, key_species_red, key_species_present_init);

    check_key_species_present_init(gas_names, key_species_present_init);

    // create flavor list
    create_flavor(key_species_red, this->flavor);

    // create gpoint flavor list
    create_gpoint_flavor(
            key_species_red, this->get_gpoint_bands(), this->flavor, this->gpoint_flavor);

    // minimum, maximum reference temperature, pressure -- assumes low-to-high ordering
    // for T, high-to-low ordering for p
    this->temp_ref_min = this->temp_ref({1});
    this->temp_ref_max = this->temp_ref({temp_ref.dim(1)});
    this->press_ref_min = this->press_ref({press_ref.dim(1)});
    this->press_ref_max = this->press_ref({1});

    // creates press_ref_log, temp_ref_delta
    this->press_ref_log_delta =
            (std::log(this->press_ref_min) - std::log(this->press_ref_max)) / (this->press_ref.dim(1) - 1);
    this->temp_ref_delta = (this->temp_ref_max - this->temp_ref_min) / (this->temp_ref.dim(1) - 1);

    // Which species are key in one or more bands?
    // this->flavor is an index into this->gas_names
    // if (allocated(this%is_key)) deallocate(this%is_key) ! Shouldn't ever happen...
    Array<int,1> is_key({get_ngas()}); // CvH bool, defaults to 0.?

    for (int j=1; j<=this->flavor.dim(2); ++j)
        for (int i=1; i<=this->flavor.dim(1); ++i)
        {
            if (this->flavor({i, j}) != 0)
                is_key({this->flavor({i, j})}) = true;
        }

    this->is_key = is_key;
}

template<typename TF>
void Gas_optics<TF>::set_solar_variability(
        const TF mg_index, const TF sb_index)
{
    constexpr TF a_offset = TF(0.1495954);
    constexpr TF b_offset = TF(0.00066696);

    for (int igpt=1; igpt<=this->solar_source_quiet.dim(1); ++igpt)
    {
        this->solar_source({igpt}) = this->solar_source_quiet({igpt})
                + (mg_index - a_offset) * this->solar_source_facular({igpt})
                + (sb_index - b_offset) * this->solar_source_sunspot({igpt});
    }

    // Scale solar source to input TSI value
    // if (present(tsi)) error_msg = this%set_tsi(tsi)
}

// Calculate the molecules of dry air.
template<typename TF>
void Gas_optics<TF>::get_col_dry(
        Array<TF,2>& col_dry, const Array<TF,2>& vmr_h2o,
        const Array<TF,2>& plev)
{
    // CvH: RRTMGP uses more accurate method based on latitude.
    constexpr TF g0 = 9.80665;

    constexpr TF avogad = 6.02214076e23;
    constexpr TF m_dry = 0.028964;
    constexpr TF m_h2o = 0.018016;

    Array<double,2> delta_plev({col_dry.dim(1), col_dry.dim(2)});
    Array<double,2> m_air     ({col_dry.dim(1), col_dry.dim(2)});

    for (int ilay=1; ilay<=col_dry.dim(2); ++ilay)
        for (int icol=1; icol<=col_dry.dim(1); ++icol)
            delta_plev({icol, ilay}) = std::abs(plev({icol, ilay}) - plev({icol, ilay+1}));

    for (int ilay=1; ilay<=col_dry.dim(2); ++ilay)
        for (int icol=1; icol<=col_dry.dim(1); ++icol)
            m_air({icol, ilay}) = (m_dry + m_h2o * vmr_h2o({icol, ilay})) / (1. + vmr_h2o({icol, ilay}));

    for (int ilay=1; ilay<=col_dry.dim(2); ++ilay)
        for (int icol=1; icol<=col_dry.dim(1); ++icol)
        {
            col_dry({icol, ilay}) = TF(10.) * delta_plev({icol, ilay}) * avogad / (TF(1000.)*m_air({icol, ilay})*TF(100.)*g0);
            col_dry({icol, ilay}) /= (TF(1.) + vmr_h2o({icol, ilay}));
        }
}

// Gas optics solver longwave variant.
template<typename TF>
void Gas_optics<TF>::gas_optics(
        const Array<TF,2>& play,
        const Array<TF,2>& plev,
        const Array<TF,2>& tlay,
        const Array<TF,1>& tsfc,
        const Gas_concs<TF>& gas_desc,
        std::unique_ptr<Optical_props_arry<TF>>& optical_props,
        Source_func_lw<TF>& sources,
        const Array<TF,2>& col_dry,
        const Array<TF,2>& tlev) const
{
    const int ncol = play.dim(1);
    const int nlay = play.dim(2);
    const int ngpt = this->get_ngpt();
    const int nband = this->get_nband();

    Array<int,2> jtemp({play.dim(1), play.dim(2)});
    Array<int,2> jpress({play.dim(1), play.dim(2)});
    Array<int,2> tropo({play.dim(1), play.dim(2)});
    Array<TF,6> fmajor({2, 2, 2, this->get_nflav(), play.dim(1), play.dim(2)});
    Array<int,4> jeta({2, this->get_nflav(), play.dim(1), play.dim(2)});

    // Gas optics.
    compute_gas_taus(
            ncol, nlay, ngpt, nband,
            play, plev, tlay, gas_desc,
            optical_props,
            jtemp, jpress, jeta, tropo, fmajor,
            col_dry);

    // External sources.
    source(
            ncol, nlay, nband, ngpt,
            play, plev, tlay, tsfc,
            jtemp, jpress, jeta, tropo, fmajor,
            sources, tlev);
}

// Gas optics solver shortwave variant.
template<typename TF>
void Gas_optics<TF>::gas_optics(
        const Array<TF,2>& play,
        const Array<TF,2>& plev,
        const Array<TF,2>& tlay,
        const Gas_concs<TF>& gas_desc,
        std::unique_ptr<Optical_props_arry<TF>>& optical_props,
        Array<TF,2>& toa_src,
        const Array<TF,2>& col_dry) const
{
    const int ncol = play.dim(1);
    const int nlay = play.dim(2);
    const int ngpt = this->get_ngpt();
    const int nband = this->get_nband();

    Array<int,2> jtemp({play.dim(1), play.dim(2)});
    Array<int,2> jpress({play.dim(1), play.dim(2)});
    Array<int,2> tropo({play.dim(1), play.dim(2)});
    Array<TF,6> fmajor({2, 2, 2, this->get_nflav(), play.dim(1), play.dim(2)});
    Array<int,4> jeta({2, this->get_nflav(), play.dim(1), play.dim(2)});

    // Gas optics.
    compute_gas_taus(
            ncol, nlay, ngpt, nband,
            play, plev, tlay, gas_desc,
            optical_props,
            jtemp, jpress, jeta, tropo, fmajor,
            col_dry);

    // External source function is constant.
    for (int igpt=1; igpt<=ngpt; ++igpt)
        for (int icol=1; icol<=ncol; ++icol)
            toa_src({icol, igpt}) = this->solar_source({igpt});
}

namespace rrtmgp_kernel_launcher
{
    template<typename TF> void zero_array(
            int ni, int nj, int nk, Array<TF,3>& array)
    {
        rrtmgp_kernels::zero_array_3D(&ni, &nj, &nk, array.ptr());
    }

    template<typename TF> void zero_array(
            int ni, int nj, int nk, int nl, Array<TF,4>& array)
    {
        rrtmgp_kernels::zero_array_4D(&ni, &nj, &nk, &nl, array.ptr());
    }

    template<typename TF>
    void interpolation(
            int ncol, int nlay,
            int ngas, int nflav, int neta, int npres, int ntemp,
            const Array<int,2>& flavor,
            const Array<TF,1>& press_ref_log,
            const Array<TF,1>& temp_ref,
            TF press_ref_log_delta,
            TF temp_ref_min,
            TF temp_ref_delta,
            TF press_ref_trop_log,
            const Array<TF,3>& vmr_ref,
            const Array<TF,2>& play,
            const Array<TF,2>& tlay,
            Array<TF,3>& col_gas,
            Array<int,2>& jtemp,
            Array<TF,6>& fmajor, Array<TF,5>& fminor,
            Array<TF,4>& col_mix,
            Array<int,2>& tropo,
            Array<int,4>& jeta,
            Array<int,2>& jpress)
    {
        rrtmgp_kernels::interpolation(
                &ncol, &nlay,
                &ngas, &nflav, &neta, &npres, &ntemp,
                const_cast<int*>(flavor.ptr()),
                const_cast<TF*>(press_ref_log.ptr()),
                const_cast<TF*>(temp_ref.ptr()),
                &press_ref_log_delta,
                &temp_ref_min,
                &temp_ref_delta,
                &press_ref_trop_log,
                const_cast<TF*>(vmr_ref.ptr()),
                const_cast<TF*>(play.ptr()),
                const_cast<TF*>(tlay.ptr()),
                col_gas.ptr(),
                jtemp.ptr(),
                fmajor.ptr(), fminor.ptr(),
                col_mix.ptr(),
                tropo.ptr(),
                jeta.ptr(),
                jpress.ptr());
    }

    template<typename TF>
    void compute_tau_absorption(
            int ncol, int nlay, int nband, int ngpt,
            int ngas, int nflav, int neta, int npres, int ntemp,
            int nminorlower, int nminorklower,
            int nminorupper, int nminorkupper,
            int idx_h2o,
            const Array<int,2>& gpoint_flavor,
            const Array<int,2>& band_lims_gpt,
            const Array<TF,4>& kmajor,
            const Array<TF,3>& kminor_lower,
            const Array<TF,3>& kminor_upper,
            const Array<int,2>& minor_limits_gpt_lower,
            const Array<int,2>& minor_limits_gpt_upper,
            const Array<int,1>& minor_scales_with_density_lower,
            const Array<int,1>& minor_scales_with_density_upper,
            const Array<int,1>& scale_by_complement_lower,
            const Array<int,1>& scale_by_complement_upper,
            const Array<int,1>& idx_minor_lower,
            const Array<int,1>& idx_minor_upper,
            const Array<int,1>& idx_minor_scaling_lower,
            const Array<int,1>& idx_minor_scaling_upper,
            const Array<int,1>& kminor_start_lower,
            const Array<int,1>& kminor_start_upper,
            Array<int,2>& tropo,
            Array<TF,4>& col_mix, Array<TF,6>& fmajor, Array<TF,5>& fminor,
            const Array<TF,2>& play, const Array<TF,2>& tlay, Array<TF,3>& col_gas,
            Array<int,4>& jeta, Array<int,2>& jtemp, Array<int,2>& jpress,
            Array<TF,3>& tau)
    {
        rrtmgp_kernels::compute_tau_absorption(
            &ncol, &nlay, &nband, &ngpt,
            &ngas, &nflav, &neta, &npres, &ntemp,
            &nminorlower, &nminorklower,
            &nminorupper, &nminorkupper,
            &idx_h2o,
            const_cast<int*>(gpoint_flavor.ptr()),
            const_cast<int*>(band_lims_gpt.ptr()),
            const_cast<TF*>(kmajor.ptr()),
            const_cast<TF*>(kminor_lower.ptr()),
            const_cast<TF*>(kminor_upper.ptr()),
            const_cast<int*>(minor_limits_gpt_lower.ptr()),
            const_cast<int*>(minor_limits_gpt_upper.ptr()),
            const_cast<int*>(minor_scales_with_density_lower.ptr()),
            const_cast<int*>(minor_scales_with_density_upper.ptr()),
            const_cast<int*>(scale_by_complement_lower.ptr()),
            const_cast<int*>(scale_by_complement_upper.ptr()),
            const_cast<int*>(idx_minor_lower.ptr()),
            const_cast<int*>(idx_minor_upper.ptr()),
            const_cast<int*>(idx_minor_scaling_lower.ptr()),
            const_cast<int*>(idx_minor_scaling_upper.ptr()),
            const_cast<int*>(kminor_start_lower.ptr()),
            const_cast<int*>(kminor_start_upper.ptr()),
            tropo.ptr(),
            col_mix.ptr(), fmajor.ptr(), fminor.ptr(),
            const_cast<TF*>(play.ptr()), const_cast<TF*>(tlay.ptr()), col_gas.ptr(),
            jeta.ptr(), jtemp.ptr(), jpress.ptr(),
            tau.ptr());
    }

    template<typename TF>
    void compute_tau_rayleigh(
            int ncol, int nlay, int nband, int ngpt,
            int ngas, int nflav, int neta, int npres, int ntemp,
            const Array<int,2>& gpoint_flavor,
            const Array<int,2>& band_lims_gpt,
            const Array<TF,4>& krayl,
            int idx_h2o, const Array<TF,2>& col_dry, const Array<TF,3>& col_gas,
            const Array<TF,5>& fminor, const Array<int,4>& jeta,
            const Array<int,2>& tropo, const Array<int,2>& jtemp,
            Array<TF,3>& tau_rayleigh)
    {
        rrtmgp_kernels::compute_tau_rayleigh(
                &ncol, &nlay, &nband, &ngpt,
                &ngas, &nflav, &neta, &npres, &ntemp,
                const_cast<int*>(gpoint_flavor.ptr()),
                const_cast<int*>(band_lims_gpt.ptr()),
                const_cast<TF*>(krayl.ptr()),
                &idx_h2o,
                const_cast<TF*>(col_dry.ptr()), const_cast<TF*>(col_gas.ptr()),
                const_cast<TF*>(fminor.ptr()), const_cast<int*>(jeta.ptr()),
                const_cast<int*>(tropo.ptr()), const_cast<int*>(jtemp.ptr()),
                tau_rayleigh.ptr());
    }

    template<typename TF>
    void reorder123x321(
            const Array<TF,3>& data,
            Array<TF,3>& data_out)
    {
        int dim1 = data.dim(1);
        int dim2 = data.dim(2);
        int dim3 = data.dim(3);
        rrtmgp_kernels::reorder_123x321_kernel(
                &dim1, &dim2, &dim3,
                const_cast<TF*>(data.ptr()),
                data_out.ptr());
    }

    template<typename TF>
    void combine_and_reorder_2str(
            int ncol, int nlay, int ngpt,
            const Array<TF,3>& tau_local, const Array<TF,3>& tau_rayleigh,
            Array<TF,3>& tau, Array<TF,3>& ssa, Array<TF,3>& g)
    {
        rrtmgp_kernels::combine_and_reorder_2str(
                &ncol, &nlay, &ngpt,
                const_cast<TF*>(tau_local.ptr()), const_cast<TF*>(tau_rayleigh.ptr()),
                tau.ptr(), ssa.ptr(), g.ptr());
    }

    template<typename TF>
    void compute_Planck_source(
            int ncol, int nlay, int nbnd, int ngpt,
            int nflav, int neta, int npres, int ntemp, int nPlanckTemp,
            const Array<TF,2>& tlay, const Array<TF,2>& tlev, const Array<TF,1>& tsfc, int sfc_lay,
            const Array<TF,6>& fmajor, const Array<int,4>& jeta, const Array<int,2>& tropo, const Array<int,2>& jtemp, const Array<int,2>& jpress,
            const Array<int,1>& gpoint_bands, const Array<int,2>& band_lims_gpt, const Array<TF,4>& pfracin, TF temp_ref_min,
            TF totplnk_delta, const Array<TF,2>& totplnk, const Array<int,2>& gpoint_flavor,
            Array<TF,2>& sfc_src, Array<TF,3>& lay_src, Array<TF,3>& lev_src_inc, Array<TF,3>& lev_src_dec,
            Array<TF,2>& sfc_src_jac)
    {
        rrtmgp_kernels::compute_Planck_source(
                &ncol, &nlay, &nbnd, &ngpt,
                &nflav, &neta, &npres, &ntemp, &nPlanckTemp,
                const_cast<TF*>(tlay.ptr()),
                const_cast<TF*>(tlev.ptr()),
                const_cast<TF*>(tsfc.ptr()),
                &sfc_lay,
                const_cast<TF*>(fmajor.ptr()),
                const_cast<int*>(jeta.ptr()),
                const_cast<int*>(tropo.ptr()),
                const_cast<int*>(jtemp.ptr()),
                const_cast<int*>(jpress.ptr()),
                const_cast<int*>(gpoint_bands.ptr()), const_cast<int*>(band_lims_gpt.ptr()), const_cast<TF*>(pfracin.ptr()), &temp_ref_min,
                &totplnk_delta, const_cast<TF*>(totplnk.ptr()), const_cast<int*>(gpoint_flavor.ptr()),
                sfc_src.ptr(), lay_src.ptr(), lev_src_inc.ptr(), lev_src_dec.ptr(),
                sfc_src_jac.ptr());
    }
}

template<typename TF>
void Gas_optics<TF>::compute_gas_taus(
        const int ncol, const int nlay, const int ngpt, const int nband,
        const Array<TF,2>& play,
        const Array<TF,2>& plev,
        const Array<TF,2>& tlay,
        const Gas_concs<TF>& gas_desc,
        std::unique_ptr<Optical_props_arry<TF>>& optical_props,
        Array<int,2>& jtemp, Array<int,2>& jpress,
        Array<int,4>& jeta,
        Array<int,2>& tropo,
        Array<TF,6>& fmajor,
        const Array<TF,2>& col_dry) const
{
    Array<TF,3> tau({ngpt, nlay, ncol});
    Array<TF,3> tau_rayleigh({ngpt, nlay, ncol});
    Array<TF,3> vmr({ncol, nlay, this->get_ngas()});
    Array<TF,3> col_gas({ncol, nlay, this->get_ngas()+1});
    col_gas.set_offsets({0, 0, -1});
    Array<TF,4> col_mix({2, this->get_nflav(), ncol, nlay});
    Array<TF,5> fminor({2, 2, this->get_nflav(), ncol, nlay});

    // CvH add all the checking...
    const int ngas = this->get_ngas();
    const int nflav = this->get_nflav();
    const int neta = this->get_neta();
    const int npres = this->get_npres();
    const int ntemp = this->get_ntemp();

    const int nminorlower = this->minor_scales_with_density_lower.dim(1);
    const int nminorklower = this->kminor_lower.dim(1);
    const int nminorupper = this->minor_scales_with_density_upper.dim(1);
    const int nminorkupper = this->kminor_upper.dim(1);

    for (int igas=1; igas<=ngas; ++igas)
    {
        const Array<TF,2>& vmr_2d = gas_desc.get_vmr(this->gas_names({igas}));

        // Fill array with constant value.
        if (vmr_2d.dim(1) == 1 && vmr_2d.dim(2) == 1)
        {
            const TF vmr_c = vmr_2d({1, 1});
            for (int ilay=1; ilay<=nlay; ++ilay)
                for (int icol=1; icol<=ncol; ++icol)
                    vmr({icol, ilay, igas}) = vmr_c;
        }
        // Fill array with constant profile.
        else if (vmr_2d.dim(1) == 1)
        {
            for (int ilay=1; ilay<=nlay; ++ilay)
            {
                const TF vmr_lay = vmr_2d({1, ilay});
                for (int icol=1; icol<=ncol; ++icol)
                    vmr({icol, ilay, igas}) = vmr_lay;
            }
        }
        // Fill array with full 2d data.
        else
        {
            for (int ilay=1; ilay<=nlay; ++ilay)
                for (int icol=1; icol<=ncol; ++icol)
                    vmr({icol, ilay, igas}) = vmr_2d({icol, ilay});
        }
    }

    // CvH: Assume that col_dry is provided.
    for (int ilay=1; ilay<=nlay; ++ilay)
        for (int icol=1; icol<=ncol; ++icol)
            col_gas({icol, ilay, 0}) = col_dry({icol, ilay});

    for (int igas=1; igas<=ngas; ++igas)
        for (int ilay=1; ilay<=nlay; ++ilay)
            for (int icol=1; icol<=ncol; ++icol)
                col_gas({icol, ilay, igas}) = vmr({icol, ilay, igas}) * col_dry({icol, ilay});

    // Call the fortran kernels
    rrtmgp_kernel_launcher::zero_array(ngpt, nlay, ncol, tau);

    rrtmgp_kernel_launcher::interpolation(
            ncol, nlay,
            ngas, nflav, neta, npres, ntemp,
            this->flavor,
            this->press_ref_log,
            this->temp_ref,
            this->press_ref_log_delta,
            this->temp_ref_min,
            this->temp_ref_delta,
            this->press_ref_trop_log,
            this->vmr_ref,
            play,
            tlay,
            col_gas,
            jtemp,
            fmajor, fminor,
            col_mix,
            tropo,
            jeta, jpress);

    int idx_h2o = -1;
    for (int i=1; i<=this->gas_names.dim(1); ++i)
        if (gas_names({i}) == "h2o")
        {
            idx_h2o = i;
            break;
        }

    if (idx_h2o == -1)
        throw std::runtime_error("idx_h2o cannot be found");

    rrtmgp_kernel_launcher::compute_tau_absorption(
            ncol, nlay, nband, ngpt,
            ngas, nflav, neta, npres, ntemp,
            nminorlower, nminorklower,
            nminorupper, nminorkupper,
            idx_h2o,
            this->gpoint_flavor,
            this->get_band_lims_gpoint(),
            this->kmajor,
            this->kminor_lower,
            this->kminor_upper,
            this->minor_limits_gpt_lower,
            this->minor_limits_gpt_upper,
            this->minor_scales_with_density_lower,
            this->minor_scales_with_density_upper,
            this->scale_by_complement_lower,
            this->scale_by_complement_upper,
            this->idx_minor_lower,
            this->idx_minor_upper,
            this->idx_minor_scaling_lower,
            this->idx_minor_scaling_upper,
            this->kminor_start_lower,
            this->kminor_start_upper,
            tropo,
            col_mix, fmajor, fminor,
            play, tlay, col_gas,
            jeta, jtemp, jpress,
            tau);

    bool has_rayleigh = (this->krayl.size() > 0);

    if (has_rayleigh)
    {
        rrtmgp_kernel_launcher::compute_tau_rayleigh(
                ncol, nlay, nband, ngpt,
                ngas, nflav, neta, npres, ntemp,
                this->gpoint_flavor,
                this->get_band_lims_gpoint(),
                this->krayl,
                idx_h2o, col_dry, col_gas,
                fminor, jeta, tropo, jtemp,
                tau_rayleigh);
    }

    combine_and_reorder(tau, tau_rayleigh, has_rayleigh, optical_props);
}

template<typename TF>
void Gas_optics<TF>::combine_and_reorder(
        const Array<TF,3>& tau,
        const Array<TF,3>& tau_rayleigh,
        const bool has_rayleigh,
        std::unique_ptr<Optical_props_arry<TF>>& optical_props) const
{
    int ncol = tau.dim(3);
    int nlay = tau.dim(2);
    int ngpt = tau.dim(1);

    if (!has_rayleigh)
    {
        // CvH for 2 stream and n-stream zero the g and ssa
        rrtmgp_kernel_launcher::reorder123x321(tau, optical_props->get_tau());
        // reorder123x321_test(optical_props->get_tau().ptr(), tau.ptr(), ngpt, nlay, ncol);

        // rrtmgp_kernel_launcher::zero_array(ngpt, nlay, ncol, optical_props->get_ssa());
        // rrtmgp_kernel_launcher::zero_array(ngpt, nlay, ncol, optical_props->get_g  ());
    }
    else
    {
        // In case of 1scl type
        // rrtmgp_kernel_launcher::reorder123x321(tau, optical_props->get_tau());

        // In case of 2str type
        rrtmgp_kernel_launcher::combine_and_reorder_2str(
                ncol, nlay, ngpt,
                tau, tau_rayleigh,
                optical_props->get_tau(), optical_props->get_ssa(), optical_props->get_g());
    }
}

template<typename TF>
void Gas_optics<TF>::source(
        const int ncol, const int nlay, const int nbnd, const int ngpt,
        const Array<TF,2>& play, const Array<TF,2>& plev,
        const Array<TF,2>& tlay, const Array<TF,1>& tsfc,
        const Array<int,2>& jtemp, const Array<int,2>& jpress,
        const Array<int,4>& jeta, const Array<int,2>& tropo,
        const Array<TF,6>& fmajor,
        Source_func_lw<TF>& sources,
        const Array<TF,2>& tlev) const
{
    // CvH Assume tlev is available.
    // Compute internal (Planck) source functions at layers and levels,
    // which depend on mapping from spectral space that creates k-distribution.
    const int nflav = this->get_nflav();
    const int neta = this->get_neta();
    const int npres = this->get_npres();
    const int ntemp = this->get_ntemp();
    const int nPlanckTemp = this->get_nPlanckTemp();
    auto gpoint_bands = this->get_gpoint_bands();
    auto band_lims_gpoint = this->get_band_lims_gpoint();

    Array<TF,3> lay_source_t({ngpt, nlay, ncol});
    Array<TF,3> lev_source_inc_t({ngpt, nlay, ncol});
    Array<TF,3> lev_source_dec_t({ngpt, nlay, ncol});
    Array<TF,2> sfc_source_t({ngpt, ncol});
    Array<TF,2> sfc_source_jac({ngpt, ncol});

    int sfc_lay = play({1, 1}) > play({1, nlay}) ? 1 : nlay;
    rrtmgp_kernel_launcher::compute_Planck_source(
            ncol, nlay, nbnd, ngpt,
            nflav, neta, npres, ntemp, nPlanckTemp,
            tlay, tlev, tsfc, sfc_lay,
            fmajor, jeta, tropo, jtemp, jpress,
            gpoint_bands, band_lims_gpoint, this->planck_frac, this->temp_ref_min,
            this->totplnk_delta, this->totplnk, this->gpoint_flavor,
            sfc_source_t, lay_source_t, lev_source_inc_t, lev_source_dec_t,
            sfc_source_jac);

    // CvH this transpose is super slow.
    for (int j=1; j<=sfc_source_t.dim(2); ++j)
        for (int i=1; i<=sfc_source_t.dim(1); ++i)
        {
            sources.get_sfc_source    ()({j, i}) = sfc_source_t  ({i, j});
            sources.get_sfc_source_jac()({j, i}) = sfc_source_jac({i, j});
        }

    rrtmgp_kernel_launcher::reorder123x321(lay_source_t, sources.get_lay_source());
    rrtmgp_kernel_launcher::reorder123x321(lev_source_inc_t, sources.get_lev_source_inc());
    rrtmgp_kernel_launcher::reorder123x321(lev_source_dec_t, sources.get_lev_source_dec());
    // reorder123x321_test(sources.get_lay_source    ().ptr(), lay_source_t    .ptr(), ngpt, nlay, ncol);
    // reorder123x321_test(sources.get_lev_source_inc().ptr(), lev_source_inc_t.ptr(), ngpt, nlay, ncol);
    // reorder123x321_test(sources.get_lev_source_dec().ptr(), lev_source_dec_t.ptr(), ngpt, nlay, ncol);
}

#ifdef FLOAT_SINGLE_RRTMGP
template class Gas_optics<float>;
#else
template class Gas_optics<double>;
#endif

max_stars_repo_name: Chiil/rte-rrtmgp-cpp

max_issues_repo_name: Chiil/rte-rrtmgp-cpp

max_forks_repo_name: Chiil/rte-rrtmgp-cpp


text: #include <boost/regex/v4/regex_match.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/*! \file
    \brief A JSON writer.

    Copyright (C) 2019-2021 kaoru  https://www.tetengo.org/
*/

#if !defined(TETENGO_PROPERTY_JSONWRITER_HPP_)
#define TETENGO_PROPERTY_JSONWRITER_HPP_

#include <memory>
#include <ostream>

#include <boost/core/noncopyable.hpp>

#include <tetengo/property/storage.hpp>

namespace tetengo::property
{
    class json_writer : private boost::noncopyable
    {
    public:
        // types

        using value_map_type = storage::value_map_type;

        // static functions

        static const json_writer& instance();

        // constructors and destructor

        ~json_writer();

        // functions

        void write(const value_map_type& value_map, std::ostream& stream) const;

    private:
        // types

        class impl;

        // variables

        const std::unique_ptr<impl> m_p_impl;

        // constructors

        json_writer();
    };

}

#endif

max_stars_repo_name: kaorut/tetengo

max_issues_repo_name: kaorut/tetengo

max_forks_repo_name: kaorut/tetengo


text:
/*
 * Copyright 2019 OmniSci, Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "QueryEngine/TableFunctions/TableFunctionCompilationContext.h"

#include <llvm/IR/Verifier.h>
#include <llvm/Support/raw_os_ostream.h>
#include <algorithm>
#include <boost/algorithm/string.hpp>

#include "QueryEngine/CodeGenerator.h"

extern std::unique_ptr<llvm::Module> g_rt_module;
extern std::unique_ptr<llvm::Module> rt_udf_cpu_module;
extern std::unique_ptr<llvm::Module> rt_udf_gpu_module;

namespace {

llvm::Function* generate_entry_point(const CgenState* cgen_state) {
  auto& ctx = cgen_state->context_;
  const auto pi8_type = llvm::PointerType::get(get_int_type(8, ctx), 0);
  const auto ppi8_type = llvm::PointerType::get(pi8_type, 0);
  const auto pi64_type = llvm::PointerType::get(get_int_type(64, ctx), 0);
  const auto ppi64_type = llvm::PointerType::get(pi64_type, 0);
  const auto i32_type = get_int_type(32, ctx);

  const auto func_type = llvm::FunctionType::get(
      i32_type, {ppi8_type, pi64_type, ppi64_type, pi64_type}, false);

  auto func = llvm::Function::Create(func_type,
                                     llvm::Function::ExternalLinkage,
                                     "call_table_function",
                                     cgen_state->module_);
  auto arg_it = func->arg_begin();
  const auto input_cols_arg = &*arg_it;
  input_cols_arg->setName("input_col_buffers");
  const auto input_row_count = &*(++arg_it);
  input_row_count->setName("input_row_count");
  const auto output_buffers = &*(++arg_it);
  output_buffers->setName("output_buffers");
  const auto output_row_count = &*(++arg_it);
  output_row_count->setName("output_row_count");
  return func;
}

inline llvm::Type* get_llvm_type_from_sql_column_type(const SQLTypeInfo elem_ti,
                                                      llvm::LLVMContext& ctx) {
  if (elem_ti.is_fp()) {
    switch (elem_ti.get_size()) {
      case 4:
        return llvm::Type::getFloatPtrTy(ctx);
      case 8:
        return llvm::Type::getDoublePtrTy(ctx);
    }
  }
  if (elem_ti.is_boolean()) {
    return llvm::Type::getInt8PtrTy(ctx);
  }
  CHECK(elem_ti.is_integer());
  switch (elem_ti.get_size()) {
    case 1:
      return llvm::Type::getInt8PtrTy(ctx);
    case 2:
      return llvm::Type::getInt16PtrTy(ctx);
    case 4:
      return llvm::Type::getInt32PtrTy(ctx);
    case 8:
      return llvm::Type::getInt64PtrTy(ctx);
  }
  LOG(FATAL) << "get_llvm_type_from_sql_column_type: not implemented for "
             << ::toString(elem_ti);
  return nullptr;
}

llvm::Value* alloc_column(std::string col_name,
                          const SQLTypeInfo& data_target_info,
                          llvm::Value* data_ptr,
                          llvm::Value* data_size,
                          llvm::LLVMContext& ctx,
                          llvm::IRBuilder<>& ir_builder,
                          bool byval) {
  /*
    Creates a new Column instance of given element type and initialize
    its data ptr and sz members. If data ptr or sz are unspecified
    (have nullptr values) then the corresponding members are
    initialized with NULL and -1, respectively.
   */
  llvm::Type* data_ptr_llvm_type =
      get_llvm_type_from_sql_column_type(data_target_info, ctx);
  llvm::StructType* col_struct_type =
      llvm::StructType::get(ctx,
                            {
                                data_ptr_llvm_type,         /* T* ptr */
                                llvm::Type::getInt64Ty(ctx) /* int64_t sz */
                            });
  auto col = ir_builder.CreateAlloca(col_struct_type);
  col->setName(col_name);
  auto col_ptr_ptr = ir_builder.CreateStructGEP(col_struct_type, col, 0);
  auto col_sz_ptr = ir_builder.CreateStructGEP(col_struct_type, col, 1);
  col_ptr_ptr->setName(col_name + ".ptr");
  col_sz_ptr->setName(col_name + ".sz");

  if (data_ptr != nullptr) {
    if (data_ptr->getType() == data_ptr_llvm_type->getPointerElementType()) {
      ir_builder.CreateStore(data_ptr, col_ptr_ptr);
    } else {
      auto tmp = ir_builder.CreateBitCast(data_ptr, data_ptr_llvm_type);
      ir_builder.CreateStore(tmp, col_ptr_ptr);
    }
  } else {
    ir_builder.CreateStore(llvm::Constant::getNullValue(data_ptr_llvm_type), col_ptr_ptr);
  }
  if (data_size != nullptr) {
    auto data_size_type = data_size->getType();
    if (data_size_type->isPointerTy()) {
      CHECK(data_size_type->getPointerElementType()->isIntegerTy(64));
      auto val = ir_builder.CreateLoad(data_size);
      ir_builder.CreateStore(val, col_sz_ptr);
    } else {
      CHECK(data_size_type->isIntegerTy(64));
      ir_builder.CreateStore(data_size, col_sz_ptr);
    }
  } else {
    auto const_minus1 = llvm::ConstantInt::get(llvm::Type::getInt64Ty(ctx), -1, true);
    ir_builder.CreateStore(const_minus1, col_sz_ptr);
  }

  if (byval) {
    return ir_builder.CreateLoad(col);
  } else {
    auto col_ptr = ir_builder.CreatePointerCast(
        col_ptr_ptr, llvm::PointerType::get(llvm::Type::getInt8Ty(ctx), 0));
    col_ptr->setName(col_name + "_ptr");
    return col_ptr;
  }
}

}  // namespace

TableFunctionCompilationContext::TableFunctionCompilationContext()
    : cgen_state_(std::make_unique<CgenState>(/*num_query_infos=*/0,
                                              /*contains_left_deep_outer_join=*/false)) {
  auto cgen_state = cgen_state_.get();
  CHECK(cgen_state);

  std::unique_ptr<llvm::Module> module(runtime_module_shallow_copy(cgen_state));
  cgen_state->module_ = module.get();

  entry_point_func_ = generate_entry_point(cgen_state);
  module_ = std::move(module);
}

void TableFunctionCompilationContext::compile(const TableFunctionExecutionUnit& exe_unit,
                                              const CompilationOptions& co,
                                              Executor* executor) {
  generateEntryPoint(exe_unit);
  if (co.device_type == ExecutorDeviceType::GPU) {
    generateGpuKernel();
  }
  finalize(co, executor);
}

void TableFunctionCompilationContext::generateEntryPoint(
    const TableFunctionExecutionUnit& exe_unit) {
  CHECK(entry_point_func_);
  auto arg_it = entry_point_func_->arg_begin();
  const auto input_cols_arg = &*arg_it;
  const auto input_row_count = &*(++arg_it);
  const auto output_buffers_arg = &*(++arg_it);
  const auto output_row_count_ptr = &*(++arg_it);

  auto cgen_state = cgen_state_.get();
  CHECK(cgen_state);
  auto& ctx = cgen_state->context_;

  const auto bb_entry = llvm::BasicBlock::Create(ctx, ".entry", entry_point_func_, 0);
  cgen_state->ir_builder_.SetInsertPoint(bb_entry);

  const auto bb_exit = llvm::BasicBlock::Create(ctx, ".exit", entry_point_func_);

  const auto func_body_bb = llvm::BasicBlock::Create(
      ctx, ".func_body", cgen_state->ir_builder_.GetInsertBlock()->getParent());
  cgen_state->ir_builder_.SetInsertPoint(func_body_bb);

  auto col_heads = generate_column_heads_load(
      exe_unit.input_exprs.size(), input_cols_arg, cgen_state->ir_builder_, ctx);
  CHECK_EQ(exe_unit.input_exprs.size(), col_heads.size());

  // The column arguments of C++ UDTFs processed by clang must be
  // passed by reference, see rbc issue 200.
  auto pass_column_by_value = exe_unit.table_func.isRuntime();
  std::vector<llvm::Value*> func_args;
  for (size_t i = 0; i < exe_unit.input_exprs.size(); i++) {
    const auto& expr = exe_unit.input_exprs[i];
    const auto& ti = expr->get_type_info();
    if (ti.is_fp()) {
      auto r = cgen_state->ir_builder_.CreateBitCast(
          col_heads[i], llvm::PointerType::get(get_fp_type(get_bit_width(ti), ctx), 0));
      func_args.push_back(cgen_state->ir_builder_.CreateLoad(r));
    } else if (ti.is_integer()) {
      auto r = cgen_state->ir_builder_.CreateBitCast(
          col_heads[i], llvm::PointerType::get(get_int_type(get_bit_width(ti), ctx), 0));
      func_args.push_back(cgen_state->ir_builder_.CreateLoad(r));
    } else if (ti.is_column()) {
      auto col = alloc_column(std::string("input_col.") + std::to_string(i),
                              ti.get_elem_type(),
                              col_heads[i],
                              input_row_count,
                              ctx,
                              cgen_state_->ir_builder_,
                              pass_column_by_value);
      func_args.push_back(col);
    } else {
      throw std::runtime_error(
          "Only integer and floating point columns or scalars are supported as inputs to "
          "table "
          "functions, got " +
          ti.get_type_name());
    }
  }
  std::vector<llvm::Value*> output_col_args;
  for (size_t i = 0; i < exe_unit.target_exprs.size(); i++) {
    auto output_load = cgen_state->ir_builder_.CreateLoad(
        cgen_state->ir_builder_.CreateGEP(output_buffers_arg, cgen_state_->llInt(i)));
    const auto& expr = exe_unit.target_exprs[i];
    const auto& ti = expr->get_type_info();
    CHECK(!ti.is_column());  // UDTF output column type is its data type
    auto col = alloc_column(std::string("output_col.") + std::to_string(i),
                            ti,
                            output_load,
                            output_row_count_ptr,
                            ctx,
                            cgen_state_->ir_builder_,
                            pass_column_by_value);
    func_args.push_back(col);
  }
  auto func_name = exe_unit.table_func.getName();
  boost::algorithm::to_lower(func_name);
  const auto table_func_return =
      cgen_state->emitExternalCall(func_name, get_int_type(32, ctx), func_args);
  table_func_return->setName("table_func_ret");

  // If table_func_return is non-negative then store the value in
  // output_row_count and return zero. Otherwise, return
  // table_func_return that negative value contains the error code.
  const auto bb_exit_0 = llvm::BasicBlock::Create(ctx, ".exit0", entry_point_func_);

  auto const_zero = llvm::ConstantInt::get(table_func_return->getType(), 0, true);
  auto is_ok = cgen_state_->ir_builder_.CreateICmpSGE(table_func_return, const_zero);
  cgen_state_->ir_builder_.CreateCondBr(is_ok, bb_exit_0, bb_exit);

  cgen_state_->ir_builder_.SetInsertPoint(bb_exit_0);
  auto r = cgen_state->ir_builder_.CreateIntCast(
      table_func_return, get_int_type(64, ctx), true);
  cgen_state->ir_builder_.CreateStore(r, output_row_count_ptr);
  cgen_state->ir_builder_.CreateRet(const_zero);

  cgen_state->ir_builder_.SetInsertPoint(bb_exit);
  cgen_state->ir_builder_.CreateRet(table_func_return);

  cgen_state->ir_builder_.SetInsertPoint(bb_entry);
  cgen_state->ir_builder_.CreateBr(func_body_bb);

  /*
  std::cout << "=================================" << std::endl;
  entry_point_func_->print(llvm::outs());
  std::cout << "=================================" << std::endl;
  */

  verify_function_ir(entry_point_func_);
}

void TableFunctionCompilationContext::generateGpuKernel() {
  CHECK(entry_point_func_);
  std::vector<llvm::Type*> arg_types;
  arg_types.reserve(entry_point_func_->arg_size());
  std::for_each(entry_point_func_->arg_begin(),
                entry_point_func_->arg_end(),
                [&arg_types](const auto& arg) { arg_types.push_back(arg.getType()); });
  CHECK_EQ(arg_types.size(), entry_point_func_->arg_size());

  auto cgen_state = cgen_state_.get();
  CHECK(cgen_state);
  auto& ctx = cgen_state->context_;

  std::vector<llvm::Type*> wrapper_arg_types(arg_types.size() + 1);
  wrapper_arg_types[0] = llvm::PointerType::get(get_int_type(32, ctx), 0);
  wrapper_arg_types[1] = arg_types[0];

  for (size_t i = 1; i < arg_types.size(); ++i) {
    wrapper_arg_types[i + 1] = arg_types[i];
  }

  auto wrapper_ft =
      llvm::FunctionType::get(llvm::Type::getVoidTy(ctx), wrapper_arg_types, false);
  kernel_func_ = llvm::Function::Create(wrapper_ft,
                                        llvm::Function::ExternalLinkage,
                                        "table_func_kernel",
                                        cgen_state->module_);

  auto wrapper_bb_entry = llvm::BasicBlock::Create(ctx, ".entry", kernel_func_, 0);
  llvm::IRBuilder<> b(ctx);
  b.SetInsertPoint(wrapper_bb_entry);
  std::vector<llvm::Value*> loaded_args = {kernel_func_->arg_begin() + 1};
  for (size_t i = 2; i < wrapper_arg_types.size(); ++i) {
    loaded_args.push_back(kernel_func_->arg_begin() + i);
  }
  auto error_lv = b.CreateCall(entry_point_func_, loaded_args);
  b.CreateStore(error_lv, kernel_func_->arg_begin());
  b.CreateRetVoid();
}

void TableFunctionCompilationContext::finalize(const CompilationOptions& co,
                                               Executor* executor) {
  /*
    TODO 1: eliminate need for OverrideFromSrc
    TODO 2: detect and link only the udf's that are needed
  */
  if (co.device_type == ExecutorDeviceType::GPU && rt_udf_gpu_module != nullptr) {
    CodeGenerator::link_udf_module(rt_udf_gpu_module,
                                   *module_,
                                   cgen_state_.get(),
                                   llvm::Linker::Flags::OverrideFromSrc);
  }
  if (co.device_type == ExecutorDeviceType::CPU && rt_udf_cpu_module != nullptr) {
    CodeGenerator::link_udf_module(rt_udf_cpu_module,
                                   *module_,
                                   cgen_state_.get(),
                                   llvm::Linker::Flags::OverrideFromSrc);
  }

  module_.release();
  // Add code to cache?

  LOG(IR) << "Table Function Entry Point IR\n"
          << serialize_llvm_object(entry_point_func_);

  if (co.device_type == ExecutorDeviceType::GPU) {
    LOG(IR) << "Table Function Kernel IR\n" << serialize_llvm_object(kernel_func_);

    CHECK(executor);
    executor->initializeNVPTXBackend();
    const auto cuda_mgr = executor->catalog_->getDataMgr().getCudaMgr();
    CHECK(cuda_mgr);

    CodeGenerator::GPUTarget gpu_target{executor->nvptx_target_machine_.get(),
                                        cuda_mgr,
                                        executor->blockSize(),
                                        cgen_state_.get(),
                                        false};
    gpu_code_ = CodeGenerator::generateNativeGPUCode(entry_point_func_,
                                                     kernel_func_,
                                                     {entry_point_func_, kernel_func_},
                                                     co,
                                                     gpu_target);
  } else {
    auto ee =
        CodeGenerator::generateNativeCPUCode(entry_point_func_, {entry_point_func_}, co);
    func_ptr = reinterpret_cast<FuncPtr>(ee->getPointerToFunction(entry_point_func_));
    own_execution_engine_ = std::move(ee);
  }

  LOG(IR) << "End of IR";
}

max_stars_repo_name: ngocbd/omniscidb

max_issues_repo_name: ngocbd/omniscidb

max_forks_repo_name: ngocbd/omniscidb


text:
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh/IO.h>

#include <boost/property_map/property_map.hpp>

#include <iostream>
#include <fstream>
#include <set>
#include <unordered_map>

#include <CGAL/boost/graph/selection.h>
#include <CGAL/boost/graph/helpers.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> SM;
typedef boost::graph_traits<SM>::face_descriptor face_descriptor;

void border_cases()
{
  SM sm_ref;
  std::ifstream input("data/nm_selection_removal.off");
  input >> sm_ref;

  {
    SM sm = sm_ref;
    std::vector<SM::Face_index> faces_to_remove;
    faces_to_remove.push_back(SM::Face_index(5));
    faces_to_remove.push_back(SM::Face_index(6));
    SM::Property_map<SM::Face_index, bool> is_selected = sm.add_property_map<SM::Face_index, bool>("f:is_selected", false).first;
    is_selected[SM::Face_index(5)]=true;
    is_selected[SM::Face_index(6)]=true;
    CGAL::expand_face_selection_for_removal(faces_to_remove,
                                            sm,
                                            is_selected);
    int i=0;
    for(face_descriptor fh : sm.faces())
      if(!is_selected[fh]) ++i;
    assert(i==4);
  }

  {
    SM sm = sm_ref;
    std::vector<SM::Face_index> faces_to_remove;
    faces_to_remove.push_back(SM::Face_index(4));
    faces_to_remove.push_back(SM::Face_index(6));
    SM::Property_map<SM::Face_index, bool> is_selected = sm.add_property_map<SM::Face_index, bool>("f:is_selected", false).first;
    is_selected[SM::Face_index(4)]=true;
    is_selected[SM::Face_index(6)]=true;
    CGAL::expand_face_selection_for_removal(faces_to_remove,
                                            sm,
                                            is_selected);
    int i=0;
    for(face_descriptor fh : sm.faces())
      if(!is_selected[fh]) ++i;

    assert(i==1 || i==4); // depends on the start point
  }

  {
    SM sm = sm_ref;
    std::vector<SM::Face_index> faces_to_remove;
    faces_to_remove.push_back(SM::Face_index(4));
    faces_to_remove.push_back(SM::Face_index(5));
    faces_to_remove.push_back(SM::Face_index(6));
    SM::Property_map<SM::Face_index, bool> is_selected = sm.add_property_map<SM::Face_index, bool>("f:is_selected", false).first;
    is_selected[SM::Face_index(4)]=true;
    is_selected[SM::Face_index(5)]=true;
    is_selected[SM::Face_index(6)]=true;
    CGAL::expand_face_selection_for_removal(faces_to_remove,
                                            sm,
                                            is_selected);
    int i=0;
    for(face_descriptor fh : sm.faces())
      if(!is_selected[fh]) ++i;
    assert(i==4); // depends on the start point
  }

}

int main()
{
  SM sm;
  std::ifstream input(CGAL::data_file_path("meshes/head.off"));
  input >> sm;

// define my selection of faces to remove
  std::unordered_map<face_descriptor, bool> is_selected_map;

  const int selection_indices[30] = {652,18,328,698,322,212,808,353,706,869,646,352,788,696,714,796,937,2892,374,697,227,501,786,794,345,16,21,581,347,723};
  std::set<int> index_set(&selection_indices[0], &selection_indices[0]+30);

  std::vector<face_descriptor> faces_to_remove;
  int index = 0;
  for(face_descriptor fh : faces(sm))
  {
    if(index_set.count(index)==0)
      is_selected_map[fh]=false;
    else
    {
      faces_to_remove.push_back(fh);
      is_selected_map[fh]=true;
    }
    ++index;
  }

  std::size_t nb_input_faces = sm.number_of_faces();

  expand_face_selection_for_removal(faces_to_remove,
                                    sm,
                                    boost::make_assoc_property_map(is_selected_map));

  index=0;
  for(face_descriptor fh : faces(sm))
  {
    if (is_selected_map[fh])
    {
      CGAL::Euler::remove_face(halfedge(fh, sm), sm);
      ++index;
    }
  }

  CGAL_USE(nb_input_faces);
  assert( sm.number_of_faces()+30 < nb_input_faces);
  assert(is_valid_polygon_mesh(sm));

  border_cases();

  return 0;
}

max_stars_repo_name: ffteja/cgal

max_issues_repo_name: ffteja/cgal

max_forks_repo_name: ffteja/cgal


text:
//==============================================================================
//         Copyright 2003 - 2011 LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011 LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef NT2_SWAR_FUNCTIONS_COMPLEX_SIMD_COMMON_SPLIT_HPP_INCLUDED
#define NT2_SWAR_FUNCTIONS_COMPLEX_SIMD_COMMON_SPLIT_HPP_INCLUDED
#include <nt2/swar/functions/split.hpp>
#include <nt2/include/functions/load.hpp>
#include <nt2/include/functions/store.hpp>
#include <boost/simd/sdk/memory/aligned_type.hpp>
#include <boost/simd/sdk/meta/cardinal_of.hpp>
#include <boost/simd/sdk/meta/scalar_of.hpp>
#include <boost/fusion/tuple.hpp>
#include <nt2/sdk/complex/complex.hpp>
#include <nt2/sdk/complex/meta/as_complex.hpp>
#include <nt2/sdk/complex/meta/as_real.hpp>
#include <nt2/sdk/complex/meta/as_dry.hpp>
#include <nt2/include/functions/imag.hpp>
#include <nt2/include/functions/real.hpp>
// TODO pb simd_ is not recognized here !?
namespace nt2 { namespace ext
{
//   NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(A1)(X),
//                               ((simd_<complex_<arithmetic_<A0> >,X>))
//                               ((simd_<complex_<arithmetic_<A1> >,X>))
//                               ((simd_<complex_<arithmetic_<A1> >,X>))
//                             )
//   {
//     typedef int result_type;
//     inline result_type operator()(A0 const& a0,A1 & a1, A1 & a2) const
//     {
//       typedef typename meta::as_real<A1>::type rA1;
//       rA1 ar0, ar1, ai0, ai1;
//       split(nt2::real(a0), ar0, ar1);
//       split(nt2::imag(a0), ai0, ai1);
//       a1 = A1(ar0, ai0);
//       a2 = A1(ar1, ai1);
//       return 0;
//     }
//   };

//   NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(X),
//                               ((simd_<complex_<arithmetic_<A0> >,X>))
//                             )
//   {
//     typedef typename meta::as_real<A0>::type rA0;
//     typedef typename boost::dispatch::meta::upgrade<rA0>::type rU;
//     typedef typename meta::as_complex<rU>::type crU;
//     typedef boost::fusion::tuple<crU, crU> result_type;

//     NT2_FUNCTOR_CALL_REPEAT(1)
//     {
//       result_type res;
//       nt2::split(a0,  boost::fusion::at_c<0>(res),  boost::fusion::at_c<1>(res));
//       return res;
//     }
//   };

//    NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(A1)(X),
//                               ((simd_<imaginary_<arithmetic_<A0> >,X>))
//                               ((simd_<imaginary_<arithmetic_<A1> >,X>))
//                               ((simd_<imaginary_<arithmetic_<A1> >,X>))
//                             )
//   {
//     typedef int result_type;
//     inline result_type operator()(A0 const& a0,A1 & a1, A1 & a2) const
//     {
//       typedef typename meta::as_real<A1>::type rA1;
//       rA1 ai0, ai1;
//       split(nt2::imag(a0), ai0, ai1);
//       a1 = A1(ai0);
//       a2 = A1(ai1);
//       return 0;
//     }
//   };

//   NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(X),
//                               ((simd_<imaginary_<arithmetic_<A0> >,X>))
//                             )
//   {
//     typedef typename meta::as_real<A0>::type rA0;
//     typedef typename boost::dispatch::meta::upgrade<rA0>::type rU;
//     typedef typename meta::as_imaginary<rU>::type crU;
//     typedef boost::fusion::tuple<crU, crU> result_type;

//     NT2_FUNCTOR_CALL_REPEAT(1)
//     {
//       result_type res;
//       nt2::split(a0,  boost::fusion::at_c<0>(res),  boost::fusion::at_c<1>(res));
//       return res;
//     }
//   };

//   NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(A1)(X),
//                               ((simd_<dry_<arithmetic_<A0> >,X>))
//                               ((simd_<dry_<arithmetic_<A1> >,X>))
//                               ((simd_<dry_<arithmetic_<A1> >,X>))
//                             )
//   {
//     typedef int result_type;
//     inline result_type operator()(A0 const& a0,A1 & a1, A1 & a2) const
//     {
//       typedef typename meta::as_real<A1>::type rA1;
//       rA1 ai0, ai1;
//       split(nt2::imag(a0), ai0, ai1);
//       a1 = A1(ai0);
//       a2 = A1(ai1);
//       return 0;
//     }
//   };

//   NT2_FUNCTOR_IMPLEMENTATION( nt2::tag::split_, tag::cpu_,
//                               (A0)(X),
//                               ((simd_<dry_<arithmetic_<A0> >,X>))
//                             )
//   {
//     typedef typename meta::as_real<A0>::type rA0;
//     typedef typename boost::dispatch::meta::upgrade<rA0>::type rU;
//     typedef typename meta::as_dry<rU>::type crU;
//     typedef boost::fusion::tuple<crU, crU> result_type;

//     NT2_FUNCTOR_CALL_REPEAT(1)
//     {
//       result_type res;
//       nt2::split(a0,  boost::fusion::at_c<0>(res),  boost::fusion::at_c<1>(res));
//       return res;
//     }
//   };
} }

#endif