macrotaskq.h

Go to the documentation of this file.

/**
 \file macrotaskq.h
 \brief Declares the \c macrotaskq and MacroTaskBase classes
 \ingroup mra
 
 A MacroTaskq executes tasks on World objects, e.g. differentiation of a function or other
 arithmetic. Complex algorithms can be implemented.
 
 The universe world is split into subworlds, each of them executing macrotasks of the task queue.
 This improves locality and speedups for large number of compute nodes, by reducing communications
 within worlds.
 
 The user defines a task (examples are found in test_vectormacrotask.cc), the tasks are
 lightweight and carry only bookkeeping information, actual input and output are stored in a
 cloud (see cloud.h)
 
 The user-defined macrotask is derived from MacroTaskIntermediate and must implement the run()
 method. A heterogeneous task queue is possible.
 
 The result of a macrotask is an object that lives in the universe and that is accessible from all
 subworlds. The result is accumulated in the universe and must therefore be a WorldObject. Currently we
 have implemented
  - Function<T,NDIM>
  - std::vector<Function<T,NDIM>> (a vector of Function<T,NDIM>)
  - ScalarResultImpl<T> (a scalar value)
  - std::vector<std::shared_ptr<ScalarResultImpl<T>>> (a vector of scalar values), shared_ptr for technical reasons
 
  - std::tuple<std::vector<XXX>, std::vector<YYY>> (a tuple of n vectors of WorldObjects: XXX, YYY, .. = {Function, ScalarResultImpl, ...})
 
 
 TODO: task submission from inside task (serialize task instead of replicate)
 TODO: update documentation
 TODO: consider serializing task member variables
 
*/
 
 
 
#ifndef SRC_MADNESS_MRA_MACROTASKQ_H_
#define SRC_MADNESS_MRA_MACROTASKQ_H_
 
#include <madness/world/cloud.h>
#include <madness/world/world.h>
#include <madness/mra/macrotaskpartitioner.h>
#include <madness/mra/memory_measurement.h>
 
namespace madness {
 
/// helper class for returning the result of a task, which is not a madness Function, but a simple scalar
 
/// the result value is accumulated via gaxpy in universe rank=0, after completion of the taskq the final
/// value can be obtained via get(), which includes a broadcast of the final value to all processes
template<typename T>
class ScalarResultImpl : public WorldObject<ScalarResultImpl<T>> {
public:
    typedef T value_type;
    ScalarResultImpl(World &world) : WorldObject<ScalarResultImpl<T>>(world) {
        this->process_pending();
    }
 
    /// Disable the default copy constructor
    ScalarResultImpl(const ScalarResultImpl& other) = delete;
    // ScalarResultImpl<T>(ScalarResultImpl<T>&& ) = default;
 
    /// disable assignment operator
    ScalarResultImpl<T>& operator=(const ScalarResultImpl<T>& other) = delete;
 
    ~ScalarResultImpl() {
    }
 
    /// simple assignment of the scalar value
    ScalarResultImpl<T>& operator=(const T& x) {
        value = x;
        return *this;
    }
 
    ScalarResultImpl<T>& operator+= (const T& x) {
        gaxpy(1.0, x, 1.0,true);
        return *this;
    }
 
    /// accumulate, optional fence
    void gaxpy(const double a, const T& right, double b, const bool fence=true) {
        if (this->get_world().rank()==0) {
            value =a*value + b * right;
        }
        else this->send(0, &ScalarResultImpl<T>::gaxpy, a, right, b, fence);
    }
 
    template<typename Archive>
    void serialize(Archive &ar) {
        ar & value;
    }
 
    /// after completion of the taskq get the final value
    T get() {
        this->get_world().gop.broadcast_serializable(*this, 0);
        return value;
    }
 
    /// get the local value of this rank, which might differ for different ranks
    /// for the final value use get()
    T get_local() const {
        return value;
    }
 
private:
    /// the scalar value
    T value=T();
};
 
template<typename T=double>
class ScalarResult {
public:
    typedef ScalarResultImpl<T> implT;
    std::shared_ptr<implT> impl;
 
    ScalarResult() = default;
    ScalarResult(World &world) : impl(new implT(world)) {}
    ScalarResult(const std::shared_ptr<implT>& impl) : impl(impl) {}
    ScalarResult& operator=(const T& x) {
        *(this->impl) = x;
        return *this;
    }
 
    std::shared_ptr<implT> get_impl() const {
        return impl;
    }
 
    void set_impl(const std::shared_ptr<implT>& newimpl) {
        impl=newimpl;
    }
 
    uniqueidT id() const {
        return impl->id();
    }
 
    /// accumulate, optional fence
    void gaxpy(const double a, const T& right, double b, const bool fence=true) {
        impl->gaxpy(a,right,b,fence);
    }
 
    template<typename Archive>
    void serialize(Archive &ar) {
        ar & impl;
    }
 
    /// after completion of the taskq get the final value
    T get() {
        return impl->get();
    }
 
    /// after completion of the taskq get the final value
    T get_local() const {
        return impl->get_local();
    }
 
};
 
/// helper function to create a vector of ScalarResultImpl, circumventing problems with the constructors
template<typename T>
std::vector<ScalarResult<T>> scalar_result_vector(World& world, std::size_t n)  {
    std::vector<ScalarResult<T>> v;
    for (std::size_t i=0; i<n; ++i) v.emplace_back(ScalarResult<T>(world));
    return v;
}
 
 
// type traits to check if a template parameter is a WorldContainer
template<typename>
struct is_scalar_result_ptr : std::false_type {};
 
template <typename T>
struct is_scalar_result_ptr<std::shared_ptr<madness::ScalarResultImpl<T>>> : std::true_type {};
 
template<typename>
struct is_scalar_result_ptr_vector : std::false_type {
};
 
template<typename T>
struct is_scalar_result_ptr_vector<std::vector<std::shared_ptr<typename madness::ScalarResultImpl<T>>>> : std::true_type {
};
 
// type traits to check if a template parameter is a WorldContainer
template<typename>
struct is_scalar_result : std::false_type {};
 
template <typename T>
struct is_scalar_result<madness::ScalarResult<T>> : std::true_type {};
 
template<typename>
struct is_scalar_result_impl : std::false_type {};
 
template <typename T>
struct is_scalar_result<madness::ScalarResultImpl<T>> : std::true_type {};
 
template<typename>
struct is_scalar_result_vector : std::false_type {
};
 
template<typename T>
struct is_scalar_result_vector<std::vector<typename madness::ScalarResult<T>>> : std::true_type {
};
 
/// check if type is a valid task result: it must be a WorldObject and must implement gaxpy
template <typename T>
inline constexpr bool is_valid_task_result_v =
                      is_madness_function<T>::value             // Function<T,NDIM>
                      || is_madness_function_vector<T>::value   // std::vector<Function<T,NDIM>>
                      || is_scalar_result<T>::value         // ScalarResultImpl<T>
                      || is_scalar_result_vector<T>::value  // std::vector<std::shared_ptr<ScalarResultImpl<T>>>
                      || is_scalar_result_ptr<T>::value         // ScalarResultImpl<T>
                      || is_scalar_result_ptr_vector<T>::value; // std::vector<std::shared_ptr<ScalarResultImpl<T>>>
 
 
template<typename> struct is_tuple : std::false_type { };
template<typename ...T> struct is_tuple<std::tuple<T...>> : std::true_type { };
 
/// given a tuple check recursively if all elements are valid task results
template<typename tupleT, std::size_t I>
constexpr bool check_tuple_is_valid_task_result() {
 
    typedef decay_tuple <tupleT> argtupleT;   // removes const, &, etc
 
    if constexpr(I >= std::tuple_size_v<tupleT>) {
        // Last case, if nothing is left to iterate, then exit the function
        return true;
    } else {
        using typeT = typename std::tuple_element<I, argtupleT>::type;// use decay types for determining a vector
        if constexpr (not is_valid_task_result_v<typeT>) {
            return false;
        } else {
            // Going for next element.
            return check_tuple_is_valid_task_result<tupleT,I+1>();
        }
    }
}
 
 
 
/// the result type of a macrotask must implement gaxpy
template<typename T>
void gaxpy(const double a, ScalarResult<T>& left, const double b, const T& right, const bool fence=true) {
    left.gaxpy(a, right, b, fence);
}
 
template <class Archive, typename T>
struct madness::archive::ArchiveStoreImpl<Archive, std::shared_ptr<ScalarResultImpl<T>>> {
    static void store(const Archive& ar, const std::shared_ptr<ScalarResultImpl<T>>& ptr) {
        bool exists=(ptr) ? true : false;
        ar & exists;
        if (exists) ar & ptr->id();
    }
};
 
 
template <class Archive, typename T>
struct madness::archive::ArchiveLoadImpl<Archive, std::shared_ptr<ScalarResultImpl<T>>> {
    static void load(const Archive& ar, std::shared_ptr<ScalarResultImpl<T>>& ptr) {
        bool exists=false;
        ar & exists;
        if (exists) {
            uniqueidT id;
            ar & id;
            World* world = World::world_from_id(id.get_world_id());
            MADNESS_ASSERT(world);
            auto ptr_opt = (world->ptr_from_id<  ScalarResultImpl<T> >(id));
            if (!ptr_opt)
            MADNESS_EXCEPTION("ScalarResultImpl: remote operation attempting to use a locally uninitialized object",0);
            ptr.reset(ptr_opt.value(), [] (ScalarResultImpl<T> *p_) -> void {}); // disable destruction
            if (!ptr)
            MADNESS_EXCEPTION("ScalarResultImpl<T> operation attempting to use an unregistered object",0);
        } else {
            ptr=nullptr;
        }
    }
};
 
struct MacroTaskInfo {
    enum StoragePolicy {
        StoreFunction,            ///< store a madness function in the cloud  -- can have a large memory impact
        StorePointerToFunction,   ///< store the pointer to the function in the cloud, the actual function lives in the universe and
        ///< its coefficients can be copied to the subworlds (e.g. by macrotaskq) when needed.
        ///< The task itself is responsible for handling data movement
        StoreFunctionViaPointer   ///< store a pointer to the function in the cloud, but macrotaskq will move the
        ///< coefficients to the subworlds when the task is started. This is the default policy.
    };
 
    friend std::ostream& operator<<(std::ostream& os, const StoragePolicy sp) {
        if (sp==StoreFunction) os << "StoreFunction";
        if (sp==StorePointerToFunction) os << "StorePointerToFunction";
        if (sp==StoreFunctionViaPointer) os << "StoreFunctionViaPointer";
        return os;
    }
 
    /// given the MacroTask's storage policy return the corresponding Cloud storage policy
    static Cloud::StoragePolicy to_cloud_storage_policy(MacroTaskInfo::StoragePolicy policy) {
        switch (policy) {
        case MacroTaskInfo::StoreFunction: return Cloud::StoreFunction;
        case MacroTaskInfo::StorePointerToFunction: return Cloud::StoreFunctionPointer;
        case MacroTaskInfo::StoreFunctionViaPointer: return Cloud::StoreFunctionPointer;
        default: MADNESS_EXCEPTION("unknown storage policy",0);
        }
    }
 
    /// return some preset policies
    /// - "default": StoreFunctionViaPointer, cloud rank-replicated, initial functions node-replicated
    /// - "small_memory": StoreFunctionViaPointer, cloud rank-replicated, initial functions distributed
    /// - "large_memory": StoreFunction, cloud rank-replicated, initial functions distributed
    /// the user can also set the policies manually
    /// note: the policies are checked for consistency when the MacroTaskQ is created
    static MacroTaskInfo preset(const std::string name) {
        MacroTaskInfo info;
        if (name=="default") {
            info.storage_policy=MacroTaskInfo::StoreFunction;
            info.cloud_distribution_policy=DistributionType::RankReplicated;
            info.ptr_target_distribution_policy=DistributionType::Distributed;
        } else if (name=="node_replicated_target") {
            info.storage_policy=MacroTaskInfo::StoreFunctionViaPointer;
            info.cloud_distribution_policy=DistributionType::RankReplicated;
            info.ptr_target_distribution_policy=DistributionType::NodeReplicated;
        } else if (name=="small_memory") {
            info.storage_policy=MacroTaskInfo::StoreFunctionViaPointer;
            info.cloud_distribution_policy=DistributionType::RankReplicated;
            info.ptr_target_distribution_policy=DistributionType::Distributed;
        } else if (name=="small_memory_owner") {
            // for tasks that fetch their own operands: the cloud holds pointers and the
            // queue does not copy coefficients into the subworld, see
            // handles_own_data_movement
            info.storage_policy=MacroTaskInfo::StorePointerToFunction;
            info.cloud_distribution_policy=DistributionType::RankReplicated;
            info.ptr_target_distribution_policy=DistributionType::Distributed;
        } else if (name=="large_memory") {
            info.storage_policy=MacroTaskInfo::StoreFunction;
            info.cloud_distribution_policy=DistributionType::RankReplicated;
            info.ptr_target_distribution_policy=DistributionType::Distributed;
        } else {
            std::string msg="MacroTaskQFactory::preset: unknown preset "+name;
            MADNESS_EXCEPTION(msg.c_str(),0);
        }
        return info;
    }
 
    static std::vector<std::string> get_all_preset_names() {
        return {"default","node_replicated_target","small_memory","small_memory_owner","large_memory"};
    }
 
    /// helper function to return all presets
    static std::vector<MacroTaskInfo> get_all_presets() {
        std::vector<MacroTaskInfo> result;
        for (const auto& name : get_all_preset_names()) {
            result.push_back(preset(name));
        }
        return result;
    }
 
    /// make sure the policies are consistent
    bool check_consistency() const {
        bool store_pointer_in_cloud = (storage_policy==MacroTaskInfo::StorePointerToFunction
            or storage_policy==MacroTaskInfo::StoreFunctionViaPointer);
        bool good=true;
 
        if (storage_policy==MacroTaskInfo::StoreFunction) {
            // if functions are stored in the cloud, the initial functions should be distributed
            good=ptr_target_distribution_policy==DistributionType::Distributed;
 
        } else if (store_pointer_in_cloud) {
            // if pointers are stored in the cloud, the initial functions can be distributed or replicated,
            // the cloud should be rank-replicated
            good=(cloud_distribution_policy==DistributionType::RankReplicated);
 
        }
        if (not good) std::cout << *this ;
 
        return good;
    }
 
    /// set policy from a vector of strings, assuming the order is storage policy, cloud distribution policy, ptr target distribution policy
    void from_vector_of_strings(const std::vector<std::string>& vec) {
        if (vec.size()!=3) {
            std::string msg="expected 3 policies, got "+std::to_string(vec.size());
            MADNESS_EXCEPTION(msg.c_str(),0);
        }
        auto remove_quotes = [](const std::string& s) {
            std::string result=s;
            if (s.size()>=2 and s.front()=='"' and s.back()=='"') {
                result=s.substr(1,s.size()-2);
            }
            return result;
        };
 
        std::string sstorage=remove_quotes(vec[0]);
        if (sstorage=="storefunction") storage_policy=MacroTaskInfo::StoreFunction;
        else if (sstorage=="storepointertofunction") storage_policy=MacroTaskInfo::StorePointerToFunction;
        else if (sstorage=="storefunctionviapointer") storage_policy=MacroTaskInfo::StoreFunctionViaPointer;
        else {
            std::string msg="unknown storage policy: "+sstorage;
            print("msg",msg);
            MADNESS_CHECK_THROW(0, "Oh no1");
        }
 
        std::string scloud=remove_quotes(vec[1]);
        if (scloud=="rankreplicated") cloud_distribution_policy=DistributionType::RankReplicated;
        else if (scloud=="nodereplicated") cloud_distribution_policy=DistributionType::NodeReplicated;
        else if (scloud=="distributed") cloud_distribution_policy=DistributionType::Distributed;
        else {
            std::string msg="unknown cloud distribution policy: "+scloud;
            print("msg",msg);
            MADNESS_CHECK_THROW(0, "Oh no2");
        }
 
        std::string sptrtarget=remove_quotes(vec[2]);
        if (sptrtarget=="rankreplicated") ptr_target_distribution_policy=DistributionType::RankReplicated;
        else if (sptrtarget=="nodereplicated") ptr_target_distribution_policy=DistributionType::NodeReplicated;
        else if (sptrtarget=="distributed") ptr_target_distribution_policy=DistributionType::Distributed;
        else {
            std::string msg="unknown ptr target distribution policy: "+sptrtarget;
            print("msg",msg);
            MADNESS_CHECK_THROW(0, "Oh no3");
        }
        // NB: parsing does not log — callers print the resolved policy
        // gated on verbosity (see ExchangeImpl::set_macro_task_info).
    }
 
    StoragePolicy storage_policy=StoreFunctionViaPointer;
    DistributionType cloud_distribution_policy=DistributionType::RankReplicated;
    DistributionType ptr_target_distribution_policy=DistributionType::NodeReplicated;
 
    friend std::ostream& operator<<(std::ostream& os, const MacroTaskInfo policy) {
        os << "StoragePolicy:                  " << policy.storage_policy << std::endl;
        os << "cloud_storage_policy:           " << to_cloud_storage_policy(policy.storage_policy) << std::endl;
        os << "cloud_distribution_policy:      " << policy.cloud_distribution_policy << std::endl;
        os << "ptr_target_distribution_policy: " << policy.ptr_target_distribution_policy << std::endl;
        return os;
    }
 
    friend std::string to_string(const MacroTaskInfo::StoragePolicy sp) {
        std::ostringstream os;
        os << sp;
        return os.str();
    }
 
    static StoragePolicy policy_to_string(const std::string policy) {
        std::string policy_lc=commandlineparser::tolower(policy);
        if (policy=="function") return MacroTaskInfo::StoreFunction;
        if (policy=="pointer") return MacroTaskInfo::StorePointerToFunction;
        if (policy=="functionviapointer") return MacroTaskInfo::StoreFunctionViaPointer;
        std::string msg="unknown policy: "+policy;
        MADNESS_EXCEPTION(msg.c_str(),1);
        return MacroTaskInfo::StorePointerToFunction;
    }
 
    nlohmann::json to_json() const {
        nlohmann::json j;
        j["storage_policy"]=to_string(storage_policy);
        j["cloud_distribution_policy"]=to_string(cloud_distribution_policy);
        j["ptr_target_distribution_policy"]=to_string(ptr_target_distribution_policy);
        return j;
    }
 
 
};
 
template<typename T=double>
std::ostream& operator<<(std::ostream& os, const typename MacroTaskInfo::StoragePolicy sp) {
    if (sp==MacroTaskInfo::StoreFunction) os << "Function";
    if (sp==MacroTaskInfo::StorePointerToFunction) os << "PointerToFunction";
    if (sp==MacroTaskInfo::StoreFunctionViaPointer) os << "FunctionViaPointer";
    return os;
}
 
/// base class
class MacroTaskBase {
public:
 
    typedef std::vector<std::shared_ptr<MacroTaskBase> > taskqT;
 
    MacroTaskBase() {}
    virtual ~MacroTaskBase() {};
 
    double priority=1.0;
    long owner_slot=-1;   ///< the subworld that should run this task; -1 means any
    enum Status {Running, Waiting, Complete, Unknown} stat=Unknown;
 
    void set_complete() {stat=Complete;}
    void set_running() {stat=Running;}
    void set_waiting() {stat=Waiting;}
 
    bool is_complete() const {return stat==Complete;}
    bool is_running() const {return stat==Running;}
    bool is_waiting() const {return stat==Waiting;}
 
    virtual void run(World& world, Cloud& cloud, taskqT& taskq, const long element,
        const bool debug, const MacroTaskInfo policy) = 0;
    virtual void cleanup() = 0;     // clear static data (presumably persistent input data)
 
    /// Finalize, for tasks that accumulate their own output.
 
    /// run_all holds only base pointers, so these have to be virtual. Stage 1 reduces the
    /// per-subworld buffers within a node, stage 2 reduces across nodes into the universe
    /// result; a nodeworld of nullptr means one node, so stage 1 has nothing to do.
    virtual bool wants_node_local_reduction() const { return false; }
    virtual void finalize_stage1(World& /*subworld*/, World* /*nodeworld*/, Cloud& /*cloud*/) {}
    virtual void finalize_stage2(World& /*subworld*/, World* /*nodeworld*/, Cloud& /*cloud*/) {}
 
    virtual void print_me(std::string s="") const {
        printf("this is task with priority %4.1f\n",priority);
    }
    virtual void print_me_as_table(std::string s="") const {
        print("nothing to print");
    }
    std::string print_priority_and_status_to_string() const {
        std::stringstream ss;
        ss << std::setw(5) << this->get_priority() << "  " <<this->stat;
        return ss.str();
    }
 
    double get_priority() const {return priority;}
    void set_priority(const double p) {priority=p;}
    long get_owner_slot() const {return owner_slot;}
    void set_owner_slot(const long owner) {owner_slot=owner;}
    /// an unowned task may be run by any subworld, an owned one only by its own
    bool is_owned_by(const long slot) const {return (owner_slot<0) or (owner_slot==slot);}
 
    friend std::ostream& operator<<(std::ostream& os, const MacroTaskBase::Status s) {
        if (s==MacroTaskBase::Status::Running) os << "Running";
        if (s==MacroTaskBase::Status::Waiting) os << "Waiting";
        if (s==MacroTaskBase::Status::Complete) os << "Complete";
        if (s==MacroTaskBase::Status::Unknown) os << "Unknown";
        return os;
    }
};
 
 
template<typename macrotaskT>
class MacroTaskIntermediate : public MacroTaskBase {
 
public:
 
    MacroTaskIntermediate() {}
 
    ~MacroTaskIntermediate() {}
 
    void cleanup() {};
};
 
 
    /// Factory for the MacroTaskQ
    class MacroTaskQFactory {
    public:
        long printlevel=0;
        World& world;
        long nworld=1;
 
        MacroTaskInfo policy=MacroTaskInfo::preset("default");
 
        MacroTaskQFactory(World& universe) : world(universe), nworld(universe.size()) {}
 
        MacroTaskQFactory& set_nworld(const long n) {
            nworld = n;
            return *this;
        }
 
        MacroTaskQFactory& set_printlevel(const long p) {
            printlevel = p;
            return *this;
        }
 
        MacroTaskQFactory& preset(const std::string name) {
            return *this;
        }
 
        MacroTaskQFactory& set_policy(const MacroTaskInfo p) {
            policy = p;
            return *this;
        }
 
        MacroTaskQFactory& set_storage_policy(const MacroTaskInfo::StoragePolicy sp) {
            policy.storage_policy = sp;
            return *this;
        }
 
        MacroTaskQFactory& set_cloud_distribution_policy(const DistributionType dp) {
            policy.cloud_distribution_policy = dp;
            return *this;
        }
 
        MacroTaskQFactory& set_ptr_target_distribution_policy(const DistributionType dp) {
            policy.ptr_target_distribution_policy = dp;
            return *this;
        }
 
    };
 
 
 
class MacroTaskQ : public WorldObject< MacroTaskQ> {
 
    World& universe;
    std::shared_ptr<World> subworld_ptr;
    /// node-scoped World for the two-stage finalize; created on first use, reused after
    std::shared_ptr<World> nodeworld_ptr;
    MacroTaskBase::taskqT taskq;
    std::mutex taskq_mutex;
    long printlevel=0;
    long nsubworld=1;
    nlohmann::json cloud_statistics;    ///< save cloud statistics after run_all()
    nlohmann::json taskq_statistics;    ///< save taskq statistics after run_all()
 
    const MacroTaskInfo policy;         ///< storage and distribution policy
 
    /// set the process map for the subworld
    std::shared_ptr< WorldDCPmapInterface< Key<1> > > pmap1;
    std::shared_ptr< WorldDCPmapInterface< Key<2> > > pmap2;
    std::shared_ptr< WorldDCPmapInterface< Key<3> > > pmap3;
    std::shared_ptr< WorldDCPmapInterface< Key<4> > > pmap4;
    std::shared_ptr< WorldDCPmapInterface< Key<5> > > pmap5;
    std::shared_ptr< WorldDCPmapInterface< Key<6> > > pmap6;
 
    bool printdebug() const {return printlevel>=10;}
    bool printprogress() const {return (printlevel>=4) and (not (printdebug()));}
    bool printtimings() const {return universe.rank()==0 and printlevel>=3;}
    bool printtimings_detail() const {return universe.rank()==0 and printlevel>=5;}
 
public:
 
    madness::Cloud cloud;
    World& get_subworld() {return *subworld_ptr;}
    long get_nsubworld() const {return nsubworld;}
    void set_printlevel(const long p) {printlevel=p;}
 
    MacroTaskInfo get_policy() const {
        return policy;
    }
 
    nlohmann::json get_cloud_statistics() const {
        return cloud_statistics;
    }
 
    nlohmann::json get_taskq_statistics() const {
        return taskq_statistics;
    }
 
    /// create an empty taskq and initialize the subworlds
    explicit MacroTaskQ(const MacroTaskQFactory factory)
      : WorldObject<MacroTaskQ>(factory.world)
      , universe(factory.world)
      , taskq()
      , printlevel(factory.printlevel)
      , nsubworld(factory.nworld)
      , policy(factory.policy)
      , cloud(factory.world)
    {
        subworld_ptr=create_worlds(universe,nsubworld);
        MADNESS_CHECK_THROW(policy.check_consistency(),"MacroTaskQ: inconsistent storage policy");
        cloud.set_replication_policy(policy.cloud_distribution_policy);
        cloud.set_storing_policy(MacroTaskInfo::to_cloud_storage_policy(policy.storage_policy));
 
        if (printdebug()) print(policy);
        this->process_pending();
    }
 
    ~MacroTaskQ() {}
 
    /// for each process create a world using a communicator shared with other processes by round-robin
    /// copy-paste from test_world.cc
    static std::shared_ptr<World> create_worlds(World& universe, const std::size_t nsubworld) {
 
        int color = universe.rank() % nsubworld;
        SafeMPI::Intracomm comm = universe.mpi.comm().Split(color, universe.rank() / nsubworld);
 
        std::shared_ptr<World> all_worlds;
        all_worlds.reset(new World(comm));
 
        universe.gop.fence();
        return all_worlds;
    }
 
    /// a World spanning the ranks that share memory with this one
 
    /// Collective, and worth creating only once: the two-stage finalize reduces within a
    /// node first, so it needs a World with node scope.
    static std::shared_ptr<World> create_node_world(World& universe) {
        SafeMPI::Intracomm comm = universe.mpi.comm().Split_type(
            SafeMPI::Intracomm::SHARED_SPLIT_TYPE, universe.rank());
        std::shared_ptr<World> node_world(new World(comm));
        universe.gop.fence();
        return node_world;
    }
 
    /// run all tasks
 
    /// Put the cloud where the tasks will read it from.
    
    /// Rank- or node-replication of the cloud, and of the targets that stored pointers refer to.
    /// Both are no-ops under a distributed policy.
    void replicate_inputs() {
        auto replication_policy = cloud.get_replication_policy();
        if (replication_policy!=DistributionType::Distributed) {
            double cpu0=cpu_time();
            if (replication_policy==DistributionType::RankReplicated) cloud.replicate();
            if (replication_policy==DistributionType::NodeReplicated) cloud.replicate_per_node(); // replicate to all hosts
            universe.gop.fence();
            double cpu1=cpu_time();
            if (printtimings_detail()) print("cloud replication wall time",cpu1-cpu0);
        }
 
        // replicate the targets (not the cloud) if needed
        {
            double cpu0=cpu_time();
            const bool need_replication_of_target=(policy.ptr_target_distribution_policy!=DistributionType::Distributed)
                    and (policy.storage_policy!=MacroTaskInfo::StoreFunction);
            Cloud::cloudtimer timer(universe,cloud.target_replication_time);
 
 
            if (need_replication_of_target) {
                auto dt=policy.ptr_target_distribution_policy;
                for (auto wo : cloud.world_object_base_list) {
                    loop_types<Cloud::DistributeFunctor, double, float, double_complex, float_complex>(std::tuple<DistributionType>(dt),wo);
                }
            }
 
            // if (need_replication_of_target) cloud.distribute_targets(policy.ptr_target_distribution_policy);
            double cpu1=cpu_time();
            if (printtimings_detail()) print("target replication wall time to ",policy.ptr_target_distribution_policy,cpu1-cpu0);
        }
    }
 
    /// Run this rank's share of the queue.
    
    /// Two ways in: if every task names an owner and there is one subworld per rank, each rank walks
    /// the queue and runs what it owns, since the assignment is already decided and the pull
    /// scheduler's round trip per task would buy nothing. Otherwise task numbers are handed out
    /// dynamically. Both converge on the same accounting.
    /// \return the cpu time this rank spent inside task bodies
    double execute_tasks() {
        World& subworld=get_subworld();
        double tasktime=0.0;
 
        // If every task names an owner and there is one subworld per rank, each rank can walk
        // the queue and run its own tasks: the assignment is already decided, so the pull
        // scheduler's round trip per task buys nothing. Falls through to the dynamic loop
        // whenever any task is unowned, which is the case for every task that does not use
        // the owner_hint hook.
        const long requester_slot = universe.rank() % std::max<long>(1, nsubworld);
        const bool all_tasks_owned = (not taskq.empty()) and std::all_of(taskq.begin(), taskq.end(),
                [](const std::shared_ptr<MacroTaskBase>& t) { return t->get_owner_slot() >= 0; });
 
        if (all_tasks_owned and (nsubworld==universe.size())) {
            for (long element=0; element<long(taskq.size()); ++element) {
                std::shared_ptr<MacroTaskBase> task=taskq[element];
                if (not task->is_owned_by(requester_slot)) continue;
                double cpu0=cpu_time();
                if (printdebug()) print("starting task no",element, "in subworld",subworld.id(),"at time",wall_time());
                task->run(subworld,cloud, taskq, element, printdebug(), policy);
                double cpu1=cpu_time();
                tasktime+=(cpu1-cpu0);
                if (printdebug()) printf("completed task %3ld after %6.1fs at time %6.1fs\n",element,cpu1-cpu0,wall_time());
            }
            if (printdebug() and subworld.rank()==0) {
                printf("rank %3d (subworld %3lu) finished its task queue at time %6.1fs (own tasktime %4.1fs)\n",
                       universe.rank(), static_cast<unsigned long>(subworld.id()), wall_time(), tasktime);
            }
        } else {
 
        if (printprogress() and universe.rank()==0) std::cout << "progress in percent: " << std::flush;
        while (true) {
            long element=get_scheduled_task_number(subworld);
            double cpu0=cpu_time();
            if (element<0) break;
            std::shared_ptr<MacroTaskBase> task=taskq[element];
            if (printdebug()) print("starting task no",element, "in subworld",subworld.id(),"at time",wall_time(), "no node", get_hostname());
 
            task->run(subworld,cloud, taskq, element, printdebug(), policy);
 
            double cpu1=cpu_time();
            set_complete(element);
            tasktime+=(cpu1-cpu0);
            if (printdebug()) printf("completed task %3ld after %6.1fs at time %6.1fs\n",element,cpu1-cpu0,wall_time());
 
            // print progress
            const std::size_t ntask=taskq.size();
            // return percentile of ntask for element
            auto in_percentile = [&ntask](const long element) {
                return std::floor(element/(0.1*(ntask+1)));
            };
            auto is_first_in_percentile = [&](const long element) {
                return (in_percentile(element)!=in_percentile(element-1));
            };
            if (printprogress() and is_first_in_percentile(element)) {
                std::cout << int(in_percentile(element)*10) << " " << std::flush;
            }
        }
        }   // both the owned-walk and the pull loop converge here
        return tasktime;
    }
 
    /// Move the results of tasks that accumulated their own output into the universe result.
    
    /// The ordering here is load-bearing and the comments inside say why. Two stages, so a node-local
    /// reduction can sit between them and only one rank per node scatters between nodes.
    void drain_own_output_buffers() {
        World& subworld=get_subworld();
        universe.gop.fence();
        // A task that accumulates its own output does so without fencing, once per batch, so
        // its buffer still has operations in flight here. The drains below read that buffer,
        // so the subworld has to be quiesced first -- once, not per batch, which is the whole
        // point of accumulating locally. Without this the drain reads a partly-written buffer
        // and the result is wrong in a way that varies from run to run.
        subworld.gop.fence();
 
        // Drain the buffers of tasks that accumulated their own output. Two stages: reduce
        // within a node, then across nodes into the universe result.
        //
        // The accumulator buffers MUST outlive the fences below, so a hook must not release
        // or clear them: a remote gaxpy that arrives after the buffer is gone finds an impl
        // that was released but is still registered, and throws bad_weak_ptr intermittently.
        // Release them in cleanup(), after the universe fence. If this ever fences per
        // subworld instead, re-audit every task that accumulates its own output.
        const bool want_node_reduction = (not taskq.empty())
            and taskq.front()->wants_node_local_reduction();
        if (want_node_reduction and not nodeworld_ptr)
            nodeworld_ptr = create_node_world(universe);
        World* nodeworld = want_node_reduction ? nodeworld_ptr.get() : nullptr;
        // one node means the node World is the universe, so stage 1 would only add a fence
        if (nodeworld and nodeworld->size() == universe.size()) nodeworld = nullptr;
 
        for (auto& t : taskq) t->finalize_stage1(subworld, nodeworld, cloud);
        if (nodeworld) nodeworld->gop.fence();
        for (auto& t : taskq) t->finalize_stage2(subworld, nodeworld, cloud);
        universe.gop.fence();
    }
    void run_all() {
 
        sort_tasks_by_priority();
 
        if (printdebug()) print_taskq();
        if (printtimings_detail()) {
            if (universe.rank()==0) {
                print("number of tasks in taskq",taskq.size());
                print("redirecting output to files task.#####");
            }
        }
        taskq_statistics["number_tasks"]=taskq.size();
 
        // replicate the cloud (not necessarily the target if pointers are stored)
        replicate_inputs();
 
        if (printdebug()) cloud.print_size(universe);
        cloud_statistics=cloud.get_statistics(universe);    // get stats before clearing the cloud
        universe.gop.fence();
        universe.gop.set_forbid_fence(true); // make sure there are no hidden universe fences
        pmap1=FunctionDefaults<1>::get_pmap();
        pmap2=FunctionDefaults<2>::get_pmap();
        pmap3=FunctionDefaults<3>::get_pmap();
        pmap4=FunctionDefaults<4>::get_pmap();
        pmap5=FunctionDefaults<5>::get_pmap();
        pmap6=FunctionDefaults<6>::get_pmap();
        set_pmap(get_subworld());
 
        double cpu00=cpu_time();
 
        World& subworld=get_subworld();
//      if (printdebug()) print("I am subworld",subworld.id());
        double tasktime=execute_tasks();   // not const: gop.sum reduces it in place below
        universe.gop.set_forbid_fence(false);
        drain_own_output_buffers();
 
        universe.gop.sum(tasktime);
        if (printprogress() and universe.rank()==0) std::cout << std::endl;
        cloud_statistics.update(cloud.gather_timings(universe));    // get stats before clearing the cloud
        double cpu11=cpu_time();
        if (printlevel>=4) {
            if (universe.rank()==0) {
                Cloud::print_timings(cloud_statistics);
                Cloud::print_memory_statistics(cloud_statistics);
                print("all tasks complete");
            }
            MemoryMeasurer::measure_and_print(universe);
        }
        if (printtimings_detail()) {
            printf("completed taskqueue after    %4.1fs at time %4.1fs\n", cpu11 - cpu00, wall_time());
            printf(" total cpu time / per world  %4.1fs %4.1fs\n", tasktime, tasktime / universe.size());
        }
        taskq_statistics["elapsed_time"]=cpu11-cpu00;
        taskq_statistics["cpu_time_per_world"]=tasktime/universe.size();
        taskq_statistics["total_cpu_time"]=tasktime;
 
        // cleanup task-persistent input data
        for (auto& task : taskq) task->cleanup();
        cloud.clear_cache(subworld);
        subworld.gop.fence();
        subworld.gop.fence();
        universe.gop.fence();
        FunctionDefaults<1>::set_pmap(pmap1);
        FunctionDefaults<2>::set_pmap(pmap2);
        FunctionDefaults<3>::set_pmap(pmap3);
        FunctionDefaults<4>::set_pmap(pmap4);
        FunctionDefaults<5>::set_pmap(pmap5);
        FunctionDefaults<6>::set_pmap(pmap6);
        universe.gop.fence();
        // restore targets to their original state
        if (policy.ptr_target_distribution_policy!=Distributed) cloud.distribute_targets();
        universe.gop.fence();
        subworld.gop.fence();
        cloud.clear();
        universe.gop.fence();
        subworld.gop.fence();
        universe.gop.fence();
    }
 
    void add_tasks(MacroTaskBase::taskqT& vtask) {
        for (const auto& t : vtask) {
            if (universe.rank()==0) t->set_waiting();
            add_replicated_task(t);
        }
    }
 
    void print_taskq() const {
        universe.gop.fence();
        if (universe.rank()==0) {
            print("\ntaskq on universe rank",universe.rank());
            print("total number of tasks: ",taskq.size());
            print(" task                                   batch                 priority  status");
            for (const auto& t : taskq) t->print_me_as_table();
        }
        universe.gop.fence();
    }
 
 
    void sort_tasks_by_priority() {
        std::stable_sort(taskq.begin(),taskq.end(),[](const std::shared_ptr<MacroTaskBase>& a, const std::shared_ptr<MacroTaskBase>& b) {
            return a->get_priority() > b->get_priority();
        });
    }
private:
 
    void add_replicated_task(const std::shared_ptr<MacroTaskBase>& task) {
        taskq.push_back(task);
    }
 
    /// scheduler is located on universe.rank==0
    long get_scheduled_task_number(World& subworld) {
        long number=0;
        if (subworld.rank()==0) {
          Future<long> r = this->send(ProcessID(0), &MacroTaskQ::get_scheduled_task_number_local);
          number=r.get();
        }
        subworld.gop.broadcast_serializable(number, 0);
        subworld.gop.fence();
        return number;
 
    }
 
    long get_scheduled_task_number_local() {
        MADNESS_ASSERT(universe.rank()==0);
        std::lock_guard<std::mutex> lock(taskq_mutex);
 
        auto is_Waiting = [](const std::shared_ptr<MacroTaskBase>& mtb_ptr) {return mtb_ptr->is_waiting();};
        auto it=std::find_if(taskq.begin(),taskq.end(),is_Waiting);
        if (it!=taskq.end()) {
            it->get()->set_running();
            long element=it-taskq.begin();
            return element;
        }
//      print("could not find task to schedule");
        return -1;
    }
 
    /// scheduler is located on rank==0
    void set_complete(const long task_number) const {
        this->task(ProcessID(0), &MacroTaskQ::set_complete_local, task_number);
    }
 
    /// scheduler is located on rank==0
    void set_complete_local(const long task_number) const {
        MADNESS_ASSERT(universe.rank()==0);
        taskq[task_number]->set_complete();
    }
 
public:
    void static set_pmap(World& world) {
        FunctionDefaults<1>::set_default_pmap(world);
        FunctionDefaults<2>::set_default_pmap(world);
        FunctionDefaults<3>::set_default_pmap(world);
        FunctionDefaults<4>::set_default_pmap(world);
        FunctionDefaults<5>::set_default_pmap(world);
        FunctionDefaults<6>::set_default_pmap(world);
    }
private:
    std::size_t size() const {
        return taskq.size();
    }
 
};
 
 
template<typename taskT>
class MacroTask {
    using partitionT = MacroTaskPartitioner::partitionT;
 
    typedef typename taskT::resultT resultT;
    typedef typename taskT::argtupleT argtupleT;
    typedef Cloud::recordlistT recordlistT;
 
 
    template<typename Q>
    struct is_vector : std::false_type {
    };
    template<typename Q>
    struct is_vector<std::vector<Q>> : std::true_type {
    };
 
    // the next two type traits figure out if a taskT has a member function: compute_priority(Batch&, argtupleT&)
    template <typename T,typename argtupleT>
    using member_compute_priority_t = decltype(std::declval<T>().compute_priority(std::declval<Batch&>(), std::declval<argtupleT&>()));
 
    template <typename T, typename argtupleT>
    using has_member_compute_priority = madness::meta::is_detected<member_compute_priority_t, T, argtupleT>;
 
    // ---- optional task hooks ----
    //
    // A task opts into a piece of framework machinery by declaring the corresponding member. Each
    // hook is a named detection trait plus an `if constexpr`, the same shape madness::archive uses
    // to ask whether a type can serialize itself (see has_member_serialize_v in type_traits.h): the
    // trait names the expression that has to compile, and the dispatch reads as an ordinary branch.
    // A task that declares nothing takes the empty branch and behaves exactly as before.
    //
    // Only the hooks whose signature depends on the task's own argument tuple or result type live
    // here; the ones with a fixed signature are plain virtuals on MacroTaskOperationBase
    // (owner_hint, accumulates_own_output, handles_own_data_movement, cleanup) or on MacroTaskBase
    // (wants_node_local_reduction, finalize_stage1, finalize_stage2), because run_all dispatches
    // those through a base pointer.
 
    /// true if: task.prepare_owner_assignment(partition, nsubworld)
    template<typename Q>
    using has_prepare_owner_assignment_t =
        decltype(std::declval<Q&>().prepare_owner_assignment(
                     std::declval<const MacroTaskPartitioner::partitionT&>(), 0L));
    template<typename Q>
    static constexpr bool has_prepare_owner_assignment_v =
        madness::meta::is_detected_v<has_prepare_owner_assignment_t, Q>;
 
    /// true if: task.store_batches(world, subworld, cloud, argtuple, nsubworld)
    template<typename Q, typename ArgTuple>
    using has_store_batches_t =
        decltype(std::declval<Q&>().store_batches(
                     std::declval<World&>(), std::declval<World&>(), std::declval<Cloud&>(),
                     std::declval<const ArgTuple&>(), 0L));
    template<typename Q, typename ArgTuple>
    static constexpr bool has_store_batches_v =
        madness::meta::is_detected_v<has_store_batches_t, Q, ArgTuple>;
 
    /// true if: task.sym_pipeline_advance(subworld, ket, next_col, next_row, has_next)
    template<typename Q, typename VecT>
    using has_sym_pipeline_advance_t =
        decltype(std::declval<Q&>().sym_pipeline_advance(
                     std::declval<World&>(), std::declval<const VecT&>(),
                     std::declval<const Batch_1D&>(), std::declval<const Batch_1D&>(), true));
    template<typename Q, typename VecT>
    static constexpr bool has_sym_pipeline_advance_v =
        madness::meta::is_detected_v<has_sym_pipeline_advance_t, Q, VecT>;
 
    /// true if: task.accumulate_locally(subworld, result_subworld)
    template<typename Q, typename ResT>
    using has_accumulate_locally_t =
        decltype(std::declval<Q&>().accumulate_locally(
                     std::declval<World&>(), std::declval<const ResT&>()));
    template<typename Q, typename ResT>
    static constexpr bool has_accumulate_locally_v =
        madness::meta::is_detected_v<has_accumulate_locally_t, Q, ResT>;
 
    /// true if: task.wants_node_local_reduction()
    template<typename Q>
    using has_wants_node_local_reduction_t =
        decltype(std::declval<const Q&>().wants_node_local_reduction());
    template<typename Q>
    static constexpr bool has_wants_node_local_reduction_v =
        madness::meta::is_detected_v<has_wants_node_local_reduction_t, Q>;
 
    /// true if: task.finalize_stage1(subworld, nodeworld)
    template<typename Q>
    using has_finalize_stage1_t =
        decltype(std::declval<Q&>().finalize_stage1(std::declval<World&>(),
                                                    std::declval<World*>()));
    template<typename Q>
    static constexpr bool has_finalize_stage1_v =
        madness::meta::is_detected_v<has_finalize_stage1_t, Q>;
 
    /// true if: task.finalize_stage2(subworld, nodeworld, universe_result)
    template<typename Q, typename ResT>
    using has_finalize_stage2_t =
        decltype(std::declval<Q&>().finalize_stage2(std::declval<World&>(),
                                                    std::declval<World*>(),
                                                    std::declval<ResT&>()));
    template<typename Q, typename ResT>
    static constexpr bool has_finalize_stage2_v =
        madness::meta::is_detected_v<has_finalize_stage2_t, Q, ResT>;
 
    taskT task;
    bool debug=false;
    bool immediate_execution=false;
    World &world;
    std::shared_ptr<MacroTaskQ> taskq_ptr;
 
public:
 
    /// constructor takes the task, but no arguments to the task
    explicit MacroTask(World &world, taskT &task)
            : MacroTask(world,task, std::make_shared<MacroTaskQ>(MacroTaskQFactory(world))) {
        immediate_execution=true;
    }
 
    /// constructor takes task and a taskq factory for customization, immediate execution
    explicit MacroTask(World &world, taskT& task, const MacroTaskQFactory factory)
        : MacroTask(world,task, std::shared_ptr<MacroTaskQ>(new MacroTaskQ(factory))) {
        immediate_execution=true;
    }
 
    /// constructor takes the task, and a taskq, execution is not immediate
    explicit MacroTask(World &world, taskT &task, std::shared_ptr<MacroTaskQ> taskq_ptr)
            : task(task), world(world), taskq_ptr(taskq_ptr) {
 
        // someone might pass in a taskq nullptr
        immediate_execution=false;  // will be reset by the forwarding constructors
        if (this->taskq_ptr==0) {
            this->taskq_ptr=std::make_shared<MacroTaskQ>(MacroTaskQFactory(world));
            immediate_execution=true;
        }
 
        if (debug) this->taskq_ptr->set_printlevel(20);
        // set the cloud policies
        auto cloud_storage_policy = MacroTaskInfo::to_cloud_storage_policy(this->taskq_ptr->get_policy().storage_policy);
        this->taskq_ptr->cloud.set_storing_policy(cloud_storage_policy);
        this->taskq_ptr->cloud.set_replication_policy(this->taskq_ptr->get_policy().cloud_distribution_policy);
 
    }
 
    MacroTask& set_debug(const bool value) {
        debug=value;
        return *this;
    }
 
    std::shared_ptr<MacroTaskQ> get_taskq() const {
        return taskq_ptr;
    }
 
 
    /// this mimicks the original call to the task functor, called from the universe
 
    /// store all input to the cloud, create output Function<T,NDIM> in the universe,
    /// create the batched task and shove it into the taskq. Possibly execute the taskq.
    template<typename ... Ts>
    resultT operator()(const Ts &... args) {
 
        auto argtuple = std::tie(args...);
        static_assert(std::is_same<decltype(argtuple), argtupleT>::value, "type or number of arguments incorrect");
 
        // partition the argument vector into batches
        auto partitioner=task.partitioner;
        if (not partitioner) partitioner.reset(new MacroTaskPartitioner);
        partitioner->set_nsubworld(world.size());
        partitionT partition = partitioner->partition_tasks(argtuple);
 
        // if task has a compute_priority function implemented, recompute the priorities
        if constexpr (has_member_compute_priority<taskT, argtupleT>::value) {
            auto predicate = [&](const Batch& batch, const argtupleT& argtuple) {
                return task.compute_priority(batch,argtuple);
            };
            std::function<double(const Batch& batch, const argtupleT& argtuple)> predicate1 = predicate;
            MacroTaskPartitioner::recompute_priorities(partition,argtuple,predicate1);
        }
 
        // let the task assign batches to owners from the whole partition, then reorder it
        if constexpr (has_prepare_owner_assignment_v<taskT>)
            task.prepare_owner_assignment(partition, taskq_ptr->get_nsubworld());
 
        // store input and output: output being a pointer to a universe function (vector)
        if (debug and world.rank()==0) {
            print("MacroTask storage policy: ",taskq_ptr->get_policy().storage_policy);
            print("Cloud storage policy:     ",taskq_ptr->cloud.get_storing_policy());
        }
 
        recordlistT inputrecords = taskq_ptr->cloud.store(world, argtuple);
        // additionally store the inputs as owner-pinned cloud batches, if the task wants them
        if constexpr (has_store_batches_v<taskT, argtupleT>)
            task.store_batches(world, taskq_ptr->get_subworld(), taskq_ptr->cloud, argtuple,
                               taskq_ptr->get_nsubworld());
        resultT result = task.allocator(world, argtuple);
        auto outputrecords =prepare_output_records(taskq_ptr->cloud, result);
 
        // create tasks and add them to the taskq
        MacroTaskBase::taskqT vtask;
        for (const auto& batch_prio : partition) {
            const long owner_slot = task.owner_hint(batch_prio.first, taskq_ptr->get_nsubworld());
            vtask.push_back(
                    std::shared_ptr<MacroTaskBase>(new MacroTaskInternal(task, batch_prio, inputrecords, outputrecords, owner_slot)));
        }
        taskq_ptr->add_tasks(vtask);
        if (immediate_execution) taskq_ptr->run_all();
 
        return result;
    }
private:
 
 
    /// store *pointers* to the result WorldObject in the cloud and return the recordlist
    recordlistT prepare_output_records(Cloud &cloud, resultT& result) {
        if constexpr (is_tuple<resultT>::value) {
            static_assert(check_tuple_is_valid_task_result<resultT,0>(),
                            "tuple has invalid result type in prepare_output_records");
        } else {
            static_assert(is_valid_task_result_v<resultT>, "unknown result type in prepare_output_records");
        }
 
        if (debug) print("storing pointers to output in cloud");
        // store an element of the tuple only
        auto store_output_records = [&](const auto& result) {
            recordlistT outputrecords;
            typedef std::decay_t<decltype(result)> argT;
            if constexpr (is_madness_function<argT>::value) {
                outputrecords += cloud.store(world, result.get_impl().get()); // store pointer to FunctionImpl
            } else if constexpr (is_madness_function_vector<argT>::value) {
                outputrecords += cloud.store(world, get_impl(result));
            } else if constexpr (is_scalar_result<argT>::value) {
                outputrecords += cloud.store(world, result.get_impl());               // store pointer to ScalarResultImpl
            } else if constexpr (is_vector<argT>::value) {
                if (is_scalar_result<typename argT::value_type>::value) {
                    // argT = std::vector<ScalarResult<T>>
                    std::vector<std::shared_ptr<typename argT::value_type::implT>> v;
                    for (const auto& ptr : result) v.push_back(ptr.get_impl());
                    outputrecords+=cloud.store(world,v);
                } else {
                    MADNESS_EXCEPTION("\n\n  unknown vector result type in prepare_input ", 1);
                }
            } else {
                MADNESS_EXCEPTION("should not be here",1);
            }
            return outputrecords;
        };
 
        recordlistT outputrecords;
        if constexpr (is_tuple<resultT>::value) {
            // loop over tuple elements -- args is the individual tuple element
            std::apply([&](auto &&... args) {
                (( outputrecords+=store_output_records(args) ), ...);
            }, result);
        } else {
            outputrecords=store_output_records(result);
        }
        return outputrecords;
    }
 
 
    class MacroTaskInternal : public MacroTaskIntermediate<MacroTask> {
 
        typedef decay_tuple<typename taskT::argtupleT> argtupleT;   // removes const, &, etc
        typedef typename taskT::resultT resultT;
        recordlistT inputrecords;
        recordlistT outputrecords;
    public:
        taskT task;
        std::string get_name() const {
            if (task.name=="unknown_task") return typeid(task).name();
            return task.name;
        }
 
        MacroTaskInternal(const taskT &task, const std::pair<Batch,double> &batch_prio,
                          const recordlistT &inputrecords, const recordlistT &outputrecords,
                          const long owner_slot = -1)
      : inputrecords(inputrecords), outputrecords(outputrecords), task(task) {
            if constexpr (is_tuple<resultT>::value) {
                static_assert(check_tuple_is_valid_task_result<resultT,0>(),
                                "tuple has invalid result type in prepare_output_records");
            } else {
                static_assert(is_valid_task_result_v<resultT>, "unknown result type in prepare_output_records");
            }
            this->task.batch=batch_prio.first;
            this->priority=batch_prio.second;
            this->set_owner_slot(owner_slot);
        }
 
 
        void print_me(std::string s="") const override {
            print("this is task",get_name(),"with batch", task.batch,"priority",this->get_priority());
        }
 
        void print_me_as_table(std::string s="") const override {
            std::stringstream ss;
            std::string name=get_name();
            std::size_t namesize=std::min(std::size_t(28),name.size());
            name += std::string(28-namesize,' ');
 
            std::stringstream ssbatch;
            ssbatch << task.batch;
            std::string strbatch=ssbatch.str();
            int nspaces=std::max(int(0),35-int(ssbatch.str().size()));
            strbatch+=std::string(nspaces,' ');
 
            ss  << name
                << std::setw(10) << strbatch
                << this->print_priority_and_status_to_string();
            print(ss.str());
        }
 
        /// accumulate the result of the task into the final result living in the universe
        template<typename resultT1, std::size_t I=0>
        typename std::enable_if<is_tuple<resultT1>::value, void>::type
        accumulate_into_final_result(World &subworld, resultT1 &final_result, const resultT1 &tmp_result, const argtupleT& argtuple) {
            if constexpr(I < std::tuple_size_v<resultT1>) {
                using elementT = typename std::tuple_element<I, resultT>::type;// use decay types for determining a vector
                auto element_final=std::get<I>(final_result);
                auto element_tmp=std::get<I>(tmp_result);
                accumulate_into_final_result<elementT>(subworld, element_final, element_tmp, argtuple);
                accumulate_into_final_result<resultT1,I+1>(subworld, final_result, tmp_result, argtuple);
            }
        }
 
        /// accumulate the result of the task into the final result living in the universe
        template<typename resultT1>
        typename std::enable_if<not is_tuple<resultT1>::value, void>::type
        accumulate_into_final_result(World &subworld, resultT1 &result, const resultT1 &result_tmp, const argtupleT& argtuple) {
                if constexpr (is_madness_function<resultT1>::value) {
                    // gaxpy can be done in reconstructed or compressed mode
                    TreeState operating_state=result_tmp.get_impl()->get_tensor_type()==TT_FULL ? compressed : reconstructed;
                    result_tmp.change_tree_state(operating_state);
                    gaxpy(1.0,result,1.0, result_tmp);
                } else if constexpr(is_madness_function_vector<resultT1>::value) {
                    TreeState operating_state=result_tmp[0].get_impl()->get_tensor_type()==TT_FULL ? compressed : reconstructed;
                    change_tree_state(result_tmp,operating_state);
                    // compress(subworld, result_tmp);
                    // resultT1 tmp1=task.allocator(subworld,argtuple);
                    // tmp1=task.batch.template insert_result_batch(tmp1,result_tmp);
                    gaxpy(1.0,result,1.0,result_tmp,false);
                    // was using operator+=, but this requires a fence, which is not allowed here..
                    // result += tmp1;
                } else if constexpr (is_scalar_result<resultT1>::value) {
                    gaxpy(1.0, result, 1.0, result_tmp.get_local(), false);
                } else if constexpr (is_scalar_result_vector<resultT1>::value) {
                    // resultT1 tmp1=task.allocator(subworld,argtuple);
                    // tmp1=task.batch.template insert_result_batch(tmp1,result_tmp);
                    std::size_t sz=result.size();
                    for (size_t i=0; i<sz; ++i) {
                        gaxpy(1.0, result[i], 1.0, result_tmp[i].get_local(), false);
                    }
                }
 
        }
 
        /// Bridge the framework's virtual finalize to the task's hooks.
 
        /// A task that does not accumulate its own output has nothing to finalize, so the
        /// guard keeps every existing task on the queue's own accumulation path.
        bool wants_node_local_reduction() const override {
            if constexpr (has_wants_node_local_reduction_v<taskT>) return task.wants_node_local_reduction();
            else return false;
        }
 
        void finalize_stage1(World& subworld, World* nodeworld, Cloud& /*cloud*/) override {
            if (not task.accumulates_own_output()) return;
            if constexpr (has_finalize_stage1_v<taskT>) task.finalize_stage1(subworld, nodeworld);
        }
 
        void finalize_stage2(World& subworld, World* nodeworld, Cloud& cloud) override {
            if (not task.accumulates_own_output()) return;
            resultT result_universe = get_output(subworld, cloud);
            if constexpr (has_finalize_stage2_v<taskT, resultT>)
                task.finalize_stage2(subworld, nodeworld, result_universe);
        }
 
        /// the next task in taskq that this task's subworld also owns, or {-1, nullptr}
 
        /// A task that prefetches needs to know which batch it will be asked for next. With
        /// owner-pinned scheduling that is decidable locally: the queue order is fixed and
        /// ownership is already assigned, so the next owned entry is just the next match.
        std::pair<long, MacroTaskInternal*> find_next_owned_task(
                const MacroTaskBase::taskqT& taskq, const long element) const {
            const long owner_slot = this->get_owner_slot();
            if (owner_slot < 0) return {-1, nullptr};
            for (long next = element + 1; next < long(taskq.size()); ++next) {
                if (not taskq[next]->is_owned_by(owner_slot)) continue;
                auto next_task = std::dynamic_pointer_cast<MacroTaskInternal>(taskq[next]);
                if (next_task) return {next, next_task.get()};
            }
            return {-1, nullptr};
        }
 
        /// called by the MacroTaskQ when the task is scheduled
 
        /// Let this task start the transfer its successor will need.
 
        /// Called before the task body so the request overlaps this task's compute. The successor is
        /// the next task in the queue this rank owns; a task that declares no prefetch hook does
        /// nothing here.
        void prefetch_for_next_task(World& subworld, const argtupleT& argtuple,
                                    MacroTaskBase::taskqT& taskq, const long element) {
            // pre-compute: let the task prefetch what the next task it owns will need, so the
            // transfer overlaps this task's compute
            {
                auto [next_elem, next_ptr] = find_next_owned_task(taskq, element);
                const bool has_next = (next_ptr != nullptr) and (next_ptr->task.batch.input.size() > 1);
                const Batch_1D next_col = has_next ? next_ptr->task.batch.input[0] : Batch_1D();
                const Batch_1D next_row = has_next ? next_ptr->task.batch.input[1] : Batch_1D();
                // the ket is the third argument; a task with fewer cannot have the hook either
                if constexpr (std::tuple_size<argtupleT>::value >= 3) {
                    using ketT = std::decay_t<std::tuple_element_t<2, argtupleT>>;
                    if constexpr (has_sym_pipeline_advance_v<taskT, ketT>)
                        task.sym_pipeline_advance(subworld, std::get<2>(argtuple),
                                                  next_col, next_row, has_next);
                }
            }
        }
 
        /// Bring the operand coefficients into the subworld, unless the task does that itself.
        void copy_operands_into_subworld(World& subworld, Cloud& cloud, argtupleT& batched_argtuple,
                                         const MacroTaskInfo& policy, const bool debug) {
            // maybe move this block to the cloud?
            // A task that fetches its own operands must not have them copied in as well --
            // that would pay for the coefficients twice, once per task instead of once per
            // owned batch, which is the cost the owner-pinned path exists to avoid. Both pointer
            // policies store a pointer, so both would otherwise deep-copy here -- and which one is
            // in force comes from the queue, so a task sharing a queue with another operator (nemo
            // runs Coulomb and exchange through one) inherits that queue's policy, not its own.
            const bool need_auto_copy =
                (policy.storage_policy==MacroTaskInfo::StoreFunctionViaPointer or
                 policy.storage_policy==MacroTaskInfo::StorePointerToFunction)
                    and not task.handles_own_data_movement();
            if (need_auto_copy) {
                double cpu0=wall_time();
                Cloud::cloudtimer timer(subworld,cloud.copy_time);
                // the functions loaded from the cloud are pointers to the universe functions,
                // retrieve the function coefficients from the universe
                // aka: turn the current shallow copy into a deep copy
                if (debug) print("loading function coefficients from universe for task",get_name());
 
                // loop over the tuple -- copy the functions from the universe to the subworld
                auto copi = [&](auto& arg) {
                    typedef std::decay_t<decltype(arg)> argT;
                    if constexpr (is_madness_function<argT>::value) {
                        arg=copy(subworld, arg);
                    } else if constexpr (is_madness_function_vector<argT>::value) {
                        for (auto& f : arg) f=copy(subworld,f);
                    }
                };
 
                unary_tuple_loop(batched_argtuple,copi);
                double cpu1=wall_time();
                if (debug) {
                    io_redirect_cout io2;
                    print("copied coefficients for task",get_name(),"in",cpu1-cpu0,"seconds");
                }
            }
        }
 
        void run(World &subworld, Cloud &cloud, MacroTaskBase::taskqT &taskq, const long element, const bool debug,
            const MacroTaskInfo policy) override {
            // per-task files are a debugging aid and become thousands at fine granularity; a
            // failing task still reports on the terminal regardless, see the catch below
            const auto io_mode = debug ? io_redirect::Mode::File : io_redirect::Mode::Discard;
            io_redirect io(io_mode, element, get_name()+"_output", get_name()+"_task", debug);
            // The try covers the whole body, not just the compute: a throw while loading the
            // inputs or issuing a prefetch would otherwise escape uninstrumented and be
            // reported by the task backend with no task id and no what().
            try {
            const argtupleT argtuple = cloud.load<argtupleT>(subworld, inputrecords);
            argtupleT batched_argtuple = task.batch.copy_input_batch(argtuple);
 
            task.subworld_ptr=&subworld;
            // before the prefetch hook, which fetches through it
            task.cloud_ptr=&cloud;
            prefetch_for_next_task(subworld, argtuple, taskq, element);
 
            copy_operands_into_subworld(subworld, cloud, batched_argtuple, policy, debug);
 
                print("starting task no",element, ", '",get_name(),"', in subworld",subworld.id(),"with priority",this->get_priority(),"at time",wall_time());
                double cpu0=cpu_time();
                resultT result_batch = std::apply(task, batched_argtuple);      // lives in the subworld, is a batch of the full vector (if applicable)
                double cpu1=cpu_time();
                constexpr std::size_t bufsize=256;
                char buffer[bufsize];
                std::snprintf(buffer,bufsize,"completed task %3ld after %6.1fs at time %6.1fs\n",element,cpu1-cpu0,wall_time());
                print(std::string(buffer));
 
                // move the result from the batch to the final result, all still in subworld
                auto insert_batch = [&](auto& element1, auto& element2) {
                    typedef std::decay_t<decltype(element1)> decay_type;;
                    if constexpr (is_vector<decay_type>::value) {
                        element1=task.batch.insert_result_batch(element1,element2);
                    } else {
                        std::swap(element1,element2);
                    }
                };
                resultT result_subworld=task.allocator(subworld,argtuple);
                if constexpr (is_tuple<resultT>::value) {
                    binary_tuple_loop(result_subworld, result_batch, insert_batch);
                } else {
                    insert_batch(result_subworld,result_batch);
                }
 
                // Accumulate the batch result. A task that accumulates its own output keeps it in
                // a subworld-local buffer and drains it once in finalize_stage2, instead of one
                // subworld->universe gaxpy per batch.
                if (task.accumulates_own_output()) {
                    if constexpr (has_accumulate_locally_v<taskT, resultT>)
                        task.accumulate_locally(subworld, result_subworld);
                } else {
                    resultT result_universe=get_output(subworld, cloud);       // lives in the universe
                    accumulate_into_final_result<resultT>(subworld, result_universe, result_subworld, argtuple);
                }
 
            } catch (std::exception& e) {
                // RSS at the moment of failure separates a memory fault (near the node ceiling,
                // with a bad_alloc what()) from a logic or communication bug (well below it).
                const double rss_at_fail = get_rss_usage_in_GB();
                print("failing task no",element,"in subworld",subworld.id(),"at time",wall_time());
                print(e.what());
                print("RSS at failure (current resident, GB):", rss_at_fail);
                print_me_as_table();
                print("\n\n");
                {
                    // this task's stream may be a file or /dev/null, so repeat it on the terminal
                    io_redirect_cout io2;
                    print("failing task no",element,"in subworld",subworld.id(),"at time",wall_time());
                    print(e.what());
                    print("RSS at failure (current resident, GB):", rss_at_fail);
                }
                MADNESS_EXCEPTION("failing task",1);
            }
 
        };
 
        // this is called after all tasks have been executed and the taskq has ended
        void cleanup() override {
            task.cleanup();
        }
 
        template<typename T, std::size_t NDIM>
        static Function<T,NDIM> pointer2WorldObject(const std::shared_ptr<FunctionImpl<T,NDIM>> impl) {
            Function<T,NDIM> result;
            result.set_impl(impl);
            return result;
        }
 
        template<typename T, std::size_t NDIM>
        static std::vector<Function<T,NDIM>> pointer2WorldObject(const std::vector<std::shared_ptr<FunctionImpl<T,NDIM>>> v_impl) {
            std::vector<Function<T,NDIM>> vresult;
            vresult.resize(v_impl.size());
            set_impl(vresult,v_impl);
            return vresult;
        }
 
        template<typename T>
        static ScalarResult<T> pointer2WorldObject(const std::shared_ptr<ScalarResultImpl<T>> sr_impl) {
            return ScalarResult(sr_impl);
        }
 
        template<typename T>
        static std::vector<ScalarResult<T>> pointer2WorldObject(const std::vector<std::shared_ptr<ScalarResultImpl<T>>> v_sr_impl) {
            std::vector<ScalarResult<T>> vresult(v_sr_impl.size());
            for (size_t i=0; i<v_sr_impl.size(); ++i) {
                vresult[i].set_impl(v_sr_impl[i]);
            }
            return vresult;
        }
 
        /// return the WorldObjects or the result functions living in the universe
 
        /// read the pointers to the universe WorldObjects from the cloud,
        /// convert them to actual WorldObjects and return them
        resultT get_output(World &subworld, Cloud &cloud) const {
            resultT result;
 
            // save outputrecords, because they will be consumed by the cloud
            auto outputrecords1 = this->outputrecords;
 
            // turn an element of the tuple of pointers into an element of the tuple of WorldObjects
            auto doit = [&](auto& element) {
                typedef std::decay_t<decltype(element)> elementT;
 
                // load the elements from the cloud -- they contain pointers to WorldObjects
                if constexpr (is_madness_function_vector<elementT>::value) {
                    typedef typename elementT::value_type::implT implT;
                    auto ptr_element = cloud.consuming_load<std::vector<std::shared_ptr<implT>>>(
                        subworld, outputrecords1);
                    element = pointer2WorldObject(ptr_element);
                }
                else if constexpr (is_madness_function<elementT>::value) {
                    typedef typename elementT::implT implT;
                    auto ptr_element = cloud.consuming_load<std::shared_ptr<implT>>(subworld, outputrecords1);
                    element = pointer2WorldObject(ptr_element);
                }
                else if constexpr (is_scalar_result_vector<elementT>::value) {  // std::vector<ScalarResult<T>>
                    typedef typename elementT::value_type ScalarResultT;
                    typedef typename ScalarResultT::implT implT;
                    typedef std::vector<std::shared_ptr<implT>> vptrT;
                    auto ptr_element = cloud.consuming_load<vptrT>(subworld, outputrecords1);
                    element = pointer2WorldObject(ptr_element);
                }
                else if constexpr (is_scalar_result<elementT>::value) {
                    // elementT is a ScalarResultImpl<T>
                    // in cloud we store a std::shared_ptr<ScalarResultImpl<T>>
                    auto ptr_element = cloud.consuming_load<std::shared_ptr<typename elementT::implT>>(subworld, outputrecords1);
                    element = pointer2WorldObject(ptr_element);
                }
                else {
                    MADNESS_EXCEPTION("confused about the type of the result", 1);
                }
            };
            if constexpr (is_tuple<resultT>::value) {
                static_assert(check_tuple_is_valid_task_result<resultT, 0>(),
                              "invalid tuple task result -- must be vectors of functions");
                static_assert(is_tuple<resultT>::value, "is a tuple");
 
                // loop over all tuple elements
                //  1. load the pointers to the WorldObjects living in the universe
                //  2. create WorldObjects from the pointers and copy them into the tuple of type resultT
 
                // turn the tuple of pointers into a tuple of WorldObjects
                unary_tuple_loop(result,doit);
 
 
            } else {
                doit(result);
 
            }
            return result;
        }
 
    };
 
};
 
class MacroTaskOperationBase {
public:
    Batch batch;
    World* subworld_ptr=0;
    /// set by MacroTaskInternal::run before the task body runs, so a task that moves its
    /// own operands (see handles_own_data_movement) can fetch them from the cloud itself
    Cloud* cloud_ptr=0;
    std::string name="unknown_task";
    std::shared_ptr<MacroTaskPartitioner> partitioner=0;
    MacroTaskOperationBase() : batch(Batch(_, _, _)), partitioner(new MacroTaskPartitioner) {}
    MacroTaskOperationBase(const std::string name) : batch(Batch(_, _, _)),
                name(name), partitioner(new MacroTaskPartitioner) {}
    virtual ~MacroTaskOperationBase() {}
 
    /// which subworld should run this batch; -1 leaves the choice to the queue
 
    /// Reads whatever assignment prepare_owner_assignment computed, so the two go together:
    /// one decides the mapping for the whole partition, this one applies it per batch.
    virtual long owner_hint(const Batch&, const long /*nsubworld*/) const { return -1; }
 
    /// true if the task accumulates its own results and drains them in finalize_stage2,
    /// instead of the queue moving every batch result into the universe result itself
    virtual bool accumulates_own_output() const { return false; }
 
    /// true if the task moves its operand coefficients into the subworld itself
    virtual bool handles_own_data_movement() const { return false; }
 
    /// release whatever the task kept across its batches, called once the queue has run
 
    /// A task class that caches operands between its batches has to keep them in static
    /// state, since each batch is a separate task object. This is where that state must be
    /// dropped: it still refers to the subworld the batches were built in, and that
    /// subworld is destroyed after the queue finishes. Clearing it later means destroying
    /// function implementations whose world is already gone.
    virtual void cleanup() {}
};
 
 
} /* namespace madness */
 
#endif /* SRC_MADNESS_MRA_MACROTASKQ_H_ */