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 BOOST_SIMD_TOOLBOX_IEEE_FUNCTIONS_SCALAR_FAST_FREXP_HPP_INCLUDED
#define BOOST_SIMD_TOOLBOX_IEEE_FUNCTIONS_SCALAR_FAST_FREXP_HPP_INCLUDED
#include <boost/simd/toolbox/ieee/functions/fast_frexp.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/dispatch/meta/adapted_traits.hpp>
#include <boost/dispatch/meta/as_integer.hpp>
#include <boost/fusion/include/vector.hpp>
#include <boost/fusion/include/at.hpp>
#include <math.h>

namespace boost { namespace simd { namespace ext
{
  BOOST_SIMD_FUNCTOR_IMPLEMENTATION( boost::simd::tag::fast_frexp_, tag::cpu_, (A0)(A2)
                            , (scalar_< double_<A0> >)
                              (scalar_< double_<A0> >)
                              (scalar_< int64_<A2> >)
                            )
  {
    typedef int result_type;
    inline result_type operator()(A0 const& a0,A0 & a1,A2 & a2) const
    {
      int32_t tmp;
      a1 = ::frexp(a0, &tmp);
      a2 = tmp;
      return 0;
    }
  };

  BOOST_SIMD_FUNCTOR_IMPLEMENTATION( boost::simd::tag::fast_frexp_, tag::cpu_, (A0)(A2)
                            , (scalar_< single_<A0> >)
                              (scalar_< single_<A0> >)
                              (scalar_< int32_<A2> >)
                            )
  {
    typedef int result_type;
    inline result_type operator()(A0 const& a0,A0 & a1,A2 & a2) const
    {
      //      int32_t tmp;
      a1 = ::frexpf(a0, &a2);
      return 0;
    }
  };

  BOOST_SIMD_FUNCTOR_IMPLEMENTATION( boost::simd::tag::fast_frexp_, tag::cpu_, (A0)(A2)
                            , (scalar_< double_<A0> >)
                              (scalar_< int64_<A2> >)
                            )
  {
    typedef A0 result_type;
    inline result_type operator()(A0 const& a0,A2 & a2) const
    {
      typedef typename dispatch::meta::as_integer<A0, signed>::type      int_type;
      A0 a1;
      boost::simd::fast_frexp(a0, a1, a2);
      return a1;
    }
  };
  BOOST_SIMD_FUNCTOR_IMPLEMENTATION( boost::simd::tag::fast_frexp_, tag::cpu_, (A0)(A2)
                            , (scalar_< single_<A0> >)
                              (scalar_< int32_<A2> >)
                            )
  {
    typedef A0 result_type;
    inline result_type operator()(A0 const& a0,A2 & a2) const
    {
      typedef typename dispatch::meta::as_integer<A0, signed>::type      int_type;
      A0 a1;
      boost::simd::fast_frexp(a0, a1, a2);
      return a1;
    }
  };

  BOOST_SIMD_FUNCTOR_IMPLEMENTATION(boost::simd::tag::fast_frexp_, tag::cpu_,
                            (A0),
                            (scalar_ < double_<A0> > )
                           )
  {
    typedef typename boost::dispatch::meta::as_floating<A0>::type mantissa;
    typedef typename dispatch::meta::as_integer<A0,signed>::type                   exponent;
    typedef boost::fusion::vector<mantissa,exponent>                               result_type;
    BOOST_SIMD_FUNCTOR_CALL(1)
    {
      result_type res;
      int r1t;
      boost::fusion::at_c<0>(res) = ::frexp(a0, &r1t);
      boost::fusion::at_c<1>(res) = r1t;
      return res;
    }
  };

  BOOST_SIMD_FUNCTOR_IMPLEMENTATION(boost::simd::tag::fast_frexp_, tag::cpu_,
                            (A0),
                            (scalar_ < single_<A0> > )
                           )
  {
    typedef typename boost::dispatch::meta::as_floating<A0>::type mantissa;
    typedef typename dispatch::meta::as_integer<A0,signed>::type                   exponent;
    typedef boost::fusion::vector<mantissa,exponent>                            result_type;

    BOOST_SIMD_FUNCTOR_CALL(1)
    {
      result_type res;
      boost::fusion::at_c<0>(res) = ::frexpf(a0, &boost::fusion::at_c<1>(res));
      return res;
    }
  };
} } }
#endif


text:
/*
Copyright 2019 Association for the promotion of open - source insurance software and for the establishment of open interface standards in the insurance industry

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 <calculusprime/internal/parsing/Function.h>
#include <regex>
#include <sstream>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/throw_exception.hpp>
#include <calculusprime/internal/defs.h>
#include <calculusprime/FunctionIdBuilder.h>
#include <calculusprime/Logging.h>
#include <calculusprime/internal/ParsingContext.h>
#include <calculusprime/internal/parsing/FunctionArgument.h>
#include <calculusprime/internal/parsing/Parser.h>

namespace bal = boost::algorithm;

namespace {
const std::regex& getRegex()
{
    static const std::regex rgx("([a-zA-Z_]\\w*)\\s*\\(([^\\)]*)\\)");
    return rgx;
}
}

namespace CalculusPrime {

Function::Function(const std::string& header_, const std::string& implementation_)
    : m_implementation(implementation_)
{
    const std::regex& rgx(getRegex());
    std::smatch matches;
    if (std::regex_search(header_, matches, rgx)) {
        m_name = matches[1];
        std::string arguments(matches[2]);

        bal::split(m_arguments, arguments, [](char c) {
            return c == ',';
        });
        for (std::string& arg : m_arguments) {
            bal::trim(arg);
        }

        m_functionId = FunctionIdBuilder::createId(m_name, m_arguments.size());
    }
    else {
        BOOST_THROW_EXCEPTION(std::runtime_error("Invalid function header: " + header_));
    }
}

Function::~Function()
{
}

Value Function::execute(const std::vector<Value>& params_, const std::shared_ptr<IParsingContext>& parsingContext_)
{
    if (params_.size() != m_arguments.size()) {
        std::ostringstream msg;
        msg << "invalid number of arguments (" << params_.size() << ") for function " + m_name
            << ", expected: " << m_arguments.size();
        BOOST_THROW_EXCEPTION(std::runtime_error(msg.str()));
    }

    // create a copy of parsingContext_
    std::shared_ptr<ParsingContext> parsingContext(std::make_shared<ParsingContext>(std::static_pointer_cast<ParsingContext>(parsingContext_)));
    Parser parser(parsingContext);
    std::ostringstream oss;
    oss << m_name << '(';

    for (std::size_t i = 0; i < m_arguments.size(); ++i) {
        const Value& param = params_[i];
        parsingContext->assignParam(m_arguments[i], std::make_shared<FunctionArgument>(param));
        if (i != 0) {
            oss << ',';
        }
        oss << param;
    }
    oss << ')';
    log4cplus::NDCContextCreator ndc(oss.str());

    Value value = parser.parse(m_functionId, m_implementation);
    return value;
}

} // namespace CalculusPrime

ext: cpp

max_stars_repo_name: allianz/CalculusPrime-PnC

max_issues_repo_name: allianz/CalculusPrime-PnC

max_forks_repo_name: allianz/CalculusPrime-PnC


text:
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <boost/asio.hpp>

#include <stdlib.h>
#include <vector>

#include "NetworkDisplay.h"

NetworkDisplay::NetworkDisplay(NetworkDisplayConfig config) {
  mConfig = config;
  mTotalOutputPixels = 0;

#ifdef __linux__
  pthread_mutex_destroy(&mMutex);
  pthread_mutex_init(&mMutex, NULL);
#endif

  InitNetworkSegments();

  mFrameCount = 0;
  mFrameRate = config.frameRate;

  mInputScreenWidth = config.inputScreenWidth;
  mInputScreenHeight = config.inputScreenHeight;
  mTotalInputPixels = mInputScreenWidth * mInputScreenHeight;

  mInputBufferSize =  mTotalInputPixels * sizeof(uint16_t);

  mInputBuffer1 = (uint16_t *)malloc(mInputBufferSize);
  mInputBuffer2 = (uint16_t *)malloc(mInputBufferSize);
  mCurrInBuffer = mInputBuffer1;

  mOutputBufferSize = mTotalOutputPixels * sizeof(uint16_t);
  mOutputBuffer1 = (uint16_t *)malloc(mOutputBufferSize);
  mOutputBuffer2 = (uint16_t *)malloc(mOutputBufferSize);
  mCurrOutBuffer = mOutputBuffer1;

  mSinglePanelWidth = config.singlePanelWidth;
  mSinglePanelHeight = config.singlePanelHeight;

  mOutputScreenWidth  = config.totalPanelsWide * config.singlePanelWidth;
  mOutputScreenHeight = config.totalPanelsTall * config.singlePanelHeight;

  mSNow  = Milliseconds();
  mSNext = mSNow + 1000 / config.frameRate;

  StartThread();

#ifdef __USE_SDL2_VIDEO__
//  mSDL2Display = new SDL2Display(mInputScreenWidth, mInputScreenHeight);
  mSDL2Display = new SDL2Display(mOutputScreenWidth, mOutputScreenHeight);
#endif

}

void NetworkDisplay::InitNetworkSegments() {

  int ipFinalDigit = mConfig.destinationIpStartDigit;

  for (uint8_t i = 0; i < mConfig.totalSegments ; i++) {
    SegmentClientConfig segmentConfig;

    segmentConfig.segmentId = i;
    segmentConfig.singlePanelHeight = mConfig.singlePanelHeight;
    segmentConfig.singlePanelWidth = mConfig.singlePanelWidth;

    segmentConfig.numPanelsWide =  mConfig.segmentPanelsWide;
    segmentConfig.numPanelsTall = mConfig.segmentPanelsTall;

    segmentConfig.destinationPort = mConfig.destinationPort;

    char *destinationIp = (char *)malloc(strlen(mConfig.destinationIP));
    sprintf(destinationIp, mConfig.destinationIP, ipFinalDigit++);
    segmentConfig.destinationIP = destinationIp;

    mTotalOutputPixels += (segmentConfig.singlePanelWidth * segmentConfig.singlePanelHeight) * segmentConfig.numPanelsWide * segmentConfig.numPanelsTall;

    auto *segment = new SegmentClient(segmentConfig);
    mSegments.push_back(segment);

    segment->StartThread();
  }

  DescribeSegments();
}

void NetworkDisplay::ThreadFunction(NetworkDisplay *remoteDisplay) {
  uint16_t currentFrame = 0;
  uint16_t smallerSceen = mInputScreenWidth < mOutputScreenWidth ? mInputScreenWidth : mOutputScreenWidth;

  printf("Smaller screen is %s\n", (mInputScreenWidth < mOutputScreenWidth) ? "INPUT" : "OUTPUT");

  while (remoteDisplay->GetThreadRunnning()) {

    if (remoteDisplay->GetFrameCount() == currentFrame) {
      usleep(100);
      continue;
    }

    for (int segmentIdx = 0; segmentIdx < remoteDisplay->mSegments.size() - 1; segmentIdx++) {
      SegmentClient *segment = remoteDisplay->mSegments[segmentIdx];

      segment->LockMutex();
//      bzero(segment->GetInputBuffer(), segment->mTotalBytes);
      uint16_t startX = segmentIdx * segment->mSegmentWidth;

      for (uint16_t y = 0; y < segment->mSegmentHeight; y++) {
        uint16_t *screenBuffer = &mCurrOutBuffer[(y * smallerSceen) + (startX)];
        uint16_t *segmentBuffer = &segment->GetInputBuffer()[y * segment->mSegmentWidth];

        memcpy(segmentBuffer, screenBuffer, segment->mSegmentWidth * sizeof(uint16_t));
      }

      segment->UnlockMutex();
      segment->SwapBuffers();
      segment->IncrementFrameCount();
    }

    currentFrame = remoteDisplay->GetFrameCount();

  }

  printf("NetworkDisplay::ThreadFunction ended\n");
}
//uint3232_t  color = 0;
void NetworkDisplay::Update() {
//  SwapBuffers();
//  printf("frame  %i\n", mFrameCount);
  LockMutex();

  bzero(mCurrOutBuffer, mOutputBufferSize);
  size_t smallerBuffer = (mInputBufferSize < mOutputBufferSize) ? mInputBufferSize : mOutputBufferSize;
  memcpy(mCurrOutBuffer, mCurrInBuffer, smallerBuffer);

#ifdef __USE_SDL2_VIDEO__
  mSDL2Display->Update(mCurrInBuffer, mTotalInputPixels);
#endif

  UnlockMutex();

  mFrameCount++;
  NextFrameDelay();
}

void NetworkDisplay::DescribeSegments() {
  printf("I have %lu segments!\n", mSegments.size());
  for (int i = 0; i < mSegments.size(); i++) {
    mSegments[i]->Describe();
  }
}

uint16_t *NetworkDisplay::GetInputBuffer() {
  return mCurrInBuffer;
}

NetworkDisplay::~NetworkDisplay() {
  mThreadRunning = false;
  usleep(100);

  if (mThread.joinable()) {
    mThread.join();
  }

  for (int segmentIdx = 0; segmentIdx < mSegments.size(); segmentIdx++) {
    mSegments[segmentIdx]->StopThread();
  }

  delete mInputBuffer1;
  delete mInputBuffer2;
  delete mOutputBuffer1;
  delete mOutputBuffer2;
}

max_stars_repo_name: jaygarcia/network-rgb-matrix-display

max_issues_repo_name: jaygarcia/network-rgb-matrix-display

max_forks_repo_name: jaygarcia/network-rgb-matrix-display


text:
// Copyright 2019 Lawrence Livermore National Security, LLC and other Metall Project Developers.
// See the top-level COPYRIGHT file for details.
//
// SPDX-License-Identifier: (Apache-2.0 OR MIT)

#include <iostream>
#include <string>
#include <vector>
#include <cstddef>

#include <boost/interprocess/managed_external_buffer.hpp>
#include <boost/interprocess/segment_manager.hpp>
#include <boost/interprocess/managed_mapped_file.hpp>
#include <boost/interprocess/allocators/allocator.hpp>

#include <metall/detail/utility/file.hpp>
#include <metall/detail/utility/mmap.hpp>

#include "../data_structure/multithread_adjacency_list.hpp"
#include "bench_driver.hpp"

using namespace adjacency_list_bench;

using key_type = uint64_t;
using value_type = uint64_t;

namespace bip = boost::interprocess;
using manager_type = bip::managed_external_buffer;
using allocator_type = bip::allocator<void, typename manager_type::segment_manager>;
using adjacency_list_type =  data_structure::multithread_adjacency_list<key_type, value_type, allocator_type>;

void *map_file(const std::string &backing_file_name, const size_t file_size) {
  metall::detail::utility::create_file(backing_file_name);
  metall::detail::utility::extend_file_size(backing_file_name, file_size);

  auto ret = metall::detail::utility::map_file_write_mode(backing_file_name,
                                                          nullptr,
                                                          file_size, 0);
  if (ret.first == -1) {
    std::cerr << "Failed to map a file" << std::endl;
    std::abort();
  }

  if (!metall::detail::utility::os_close(ret.first)) {
    std::cerr << "Failed to close the file" << std::endl;
    std::abort();
  }

  return ret.second;
}

void *map_anonymous(const size_t file_size) {
  void *ret = metall::detail::utility::os_mmap(nullptr,
                                               file_size,
                                               PROT_READ | PROT_WRITE,
                                               MAP_PRIVATE | MAP_ANONYMOUS,
                                               -1, 0);
  if (!ret) {
    std::cerr << "Failed to map an anonymous region" << std::endl;
    std::abort();
  }
  return ret;
}

int main(int argc, char *argv[]) {
  bench_options option;
  if (!parse_options(argc, argv, &option)) {
    std::abort();
  }

  if (option.segment_size == 0) {
    std::cerr << "file size is required" << std::endl;
    std::abort();
  }

  void *addr = nullptr;
  if (option.datastore_path_list.empty()) {
    std::cout << "!!! Map ANONYMOUS region !!!" << std::endl;
    addr = map_anonymous(option.segment_size);
  } else {
    std::cout << "Map a file" << std::endl;
    bip::file_mapping::remove(option.datastore_path_list[0].c_str());
    addr = map_file(option.datastore_path_list[0], option.segment_size);
  }

  {
    manager_type manager(bip::create_only, addr, option.segment_size);
    auto adj_list = manager.construct<adjacency_list_type>(option.adj_list_key_name.c_str())
                                                          (manager.get_allocator<std::byte>());

    run_bench(option, adj_list);

    const auto start = metall::detail::utility::elapsed_time_sec();
    metall::detail::utility::os_msync(addr, option.segment_size, true);
    const auto elapsed_time = metall::detail::utility::elapsed_time_sec(start);
    std::cout << "sync_time (s)\t" << elapsed_time << std::endl;

    std::cout << "Segment usage (GB) "
              << static_cast<double>(manager.get_size() - manager.get_free_memory()) / (1ULL << 30) << std::endl;
    metall::detail::utility::munmap(addr, option.segment_size, !option.datastore_path_list.empty());
  }

  return 0;
}

max_stars_repo_name: correaa/metall

max_issues_repo_name: correaa/metall

max_forks_repo_name: correaa/metall


text:
/*

Copyright (c) 2005-2022, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifndef ABSTRACTCARDIACCELLFACTORY_HPP_
#define ABSTRACTCARDIACCELLFACTORY_HPP_

#include <boost/shared_ptr.hpp>

#include "AbstractCardiacCellInterface.hpp"
#include "AbstractTetrahedralMesh.hpp"
#include "AbstractIvpOdeSolver.hpp"
#include "EulerIvpOdeSolver.hpp"
#include "HeartGeometryInformation.hpp"
#include "HeartRegionCodes.hpp"
#include "ZeroStimulus.hpp"

/**
 * A factory to ease creating cardiac cell objects for use in a mono/bidomain simulation.
 *
 * The user should implement their own concrete class, in particular implementing
 * CreateCardiacCellForTissueNode(Node*), which should return the cell corresponding to a
 * given node. The user should also implement GetNumberOfCells() if this isn't equal
 * to the number of nodes. FinaliseCellCreation() can be used to (eg) add stimuli to
 * certain cells after they have been created.
 *
 * This class saves the user having to create cells in parallel, that work is done
 * by the pde instead.
 */
template<unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
class AbstractCardiacCellFactory
{
private:
    /** The mesh is automatically set in MonodomainProblem and BidomainProblem.
     *  This member variable should be accessed through GetMesh(), which will check if it has been set before
     *  and throw an exception otherwise.*/
    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;

    /**
     * A pointer to an HeartGeometryInformation information object
     * Can be accessed via get and set methods in this class.
     */
    HeartGeometryInformation<SPACE_DIM>* mpHeartGeometryInformation;

protected:
    /** For use at un-stimulated cells. */
    boost::shared_ptr<ZeroStimulus> mpZeroStimulus;
    /** The solver to give each of the cells */
    boost::shared_ptr<AbstractIvpOdeSolver> mpSolver;

public:
    /**
     * @return a newly created cell object for the given node.
     *
     * The default implementation checks whether the node is in the bath (in which
     * case a pointer to a (unique) fake cell is returned) and if not, calls
     * CreateCardiacCellForTissueNode (which must be defined by subclasses).
     *
     * @param pNode  Pointer to node object.
     */
    virtual AbstractCardiacCellInterface* CreateCardiacCellForNode(Node<SPACE_DIM>* pNode);

    /**
     * Must be overridden by subclasses to return a cell object for the given node.
     *
     * @param pNode  Pointer to node object.
     * @return a newly created cell object for the given tissue node.
     */
    virtual AbstractCardiacCellInterface* CreateCardiacCellForTissueNode(Node<SPACE_DIM>* pNode)=0;

    /**
     * May be overridden by subclasses to perform any necessary work after all cells
     * have been created.
     *
     * @param pCellsDistributed  Pointer to a vector of cardiac cell pointers.
     * @param lo  Lowest index owned by this process.
     * @param hi  Highest index owned by this process.
     */
    virtual void FinaliseCellCreation(std::vector< AbstractCardiacCellInterface* >* pCellsDistributed,
                                      unsigned lo, unsigned hi);

    /**
     * Method that fills in the vector of heterogeneity areas with the NodesLists
     * that correspond to a given layer (implemented in subclasses)
     */
    virtual void FillInCellularTransmuralAreas();

    /**
     * @return  The number of cells
     */
    virtual unsigned GetNumberOfCells();

    /**
     * Default constructor.
     *
     * @param pSolver  the ODE solver to use to simulate this cell.
     */
    AbstractCardiacCellFactory(boost::shared_ptr<AbstractIvpOdeSolver> pSolver = boost::shared_ptr<AbstractIvpOdeSolver>(new EulerIvpOdeSolver));
    /**
     * Destructor: free solver, zero stimulus and fake bath cell.
     */
    virtual ~AbstractCardiacCellFactory();

    /**
     * @param pMesh  the mesh for which to create cardiac cells.
     */
    virtual void SetMesh(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);

    /**
     * @return  the mesh used to create the cells.
     */
    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* GetMesh();

    /**
     * Set the HeartGeometryInformation object
     *
     * @param pHeartGeometryInformation the HeartGeometryInformation object that is to be set
     */
    void SetHeartGeometryInformation(HeartGeometryInformation<SPACE_DIM>* pHeartGeometryInformation);

    /**
     * @return the HeartGeometryInformation object
     */
    HeartGeometryInformation<SPACE_DIM>* GetHeartGeometryInformation();
};

#endif /*ABSTRACTCARDIACCELLFACTORY_HPP_*/

ext: hpp

max_stars_repo_name: stu-l/Chaste

max_issues_repo_name: stu-l/Chaste

max_forks_repo_name: stu-l/Chaste


text:
#ifndef BLUB_CORE_BYTEARRAY_HPP
#define BLUB_CORE_BYTEARRAY_HPP

#include "blub/core/globals.hpp"
#include "blub/serialization/access.hpp"
#include "blub/serialization/nameValuePair.hpp"

#include <boost/serialization/base_object.hpp>
#include <boost/serialization/vector.hpp>

#include <vector>

namespace blub
{

class byteArray : public std::vector<char>
{
public:
    typedef std::vector<char> t_base;

    byteArray()
    {
        ;
    }
    byteArray(const uint32& size)
        : t_base(size, 0)
    {;}

    byteArray(const byteArray &array);
    byteArray(const char *array, uint32 size);
    byteArray(const string &str);

    /**
     * @brief mid gets an copy from a specified index for a length
     * @param index
     * @param length if length is smaller zero, the method returns an array from index to the end
     * @return
     */
    byteArray mid(const uint32& index, const int32& length = -1) const;
    uint32 size() const
    {
        return static_cast<uint32>(t_base::size());
    }

    byteArray& operator +=(const byteArray& other);

private:
    BLUB_SERIALIZATION_ACCESS

    template <class formatType>
    void serialize(formatType & ar, const uint32& version)
    {
        (void)version;

        ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(t_base);
    }

};

}

#endif // BLUB_CORE_BYTEARRAY_HPP


text:
#include <ecto/ecto.hpp>
#include <boost/thread.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>

using ecto::tendrils;
using ecto::spore;

namespace ecto_test
{

  struct Accumulator
  {
    static void declare_io(const tendrils& p, tendrils& i, tendrils& o)
    {
      i.declare<double>(&Accumulator::left_,"left", "Left hand operand.");
      i.declare<double>(&Accumulator::right_,"right","Right hand operand.");
      o.declare<double>(&Accumulator::out_,"out","The current accumulation.", 0.0);
    }
    int process(const tendrils& inputs, const tendrils& /*outputs*/)
    {
      boost::mutex::scoped_lock lock(mutex);
      if ( inputs.find("left") != inputs.end() ) {
        std::cout << "  Left: " << *out_ << "+" << *left_ << "=" << *out_ + *left_ << std::endl;
        *out_ += *left_;
      }
      if ( inputs.find("right") != inputs.end() ) {
        std::cout << "  Right: " << *out_ << "+" << *right_ << "=" << *out_ + *right_ << std::endl;
        *out_ += *right_;
      }
      // sleep a bit to effectively test the sharing
      boost::this_thread::sleep(boost::posix_time::milliseconds(100));
      return ecto::OK;
    }
    spore<double> out_, left_, right_;
    mutable boost::mutex mutex;
  };
}

ECTO_CELL(ecto_test, ecto_test::Accumulator, "Accumulator", "Add inputs (potentially from different threads) to an incrementally accumulating sum.");

max_stars_repo_name: fujiehuang/ecto

max_issues_repo_name: fujiehuang/ecto

max_forks_repo_name: fujiehuang/ecto


text:
// Copyright 2019 Fred Hutchinson Cancer Research Center
// See the included LICENSE file for details on the licence that is granted to the user of this software.
#include <cytolib/GatingHierarchy.hpp>
#include <cytolib/global.hpp>
#include <boost/graph/graphviz.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/topological_sort.hpp>
#include <boost/graph/breadth_first_search.hpp>

namespace cytolib
{
	CytoFramePtr load_cytoframe(const string & uri, bool readonly, CytoCtx ctx)
	{
		 CytoFramePtr ptr;

		CytoVFS vfs(ctx);
		 auto fmt = uri_backend_type(uri, vfs);
		 bool is_exist = fmt == FileFormat::H5?vfs.is_file(uri):vfs.is_dir(uri);
		if(!is_exist)
		 throw(domain_error("cytoframe file missing for sample: " + uri));
		if(fmt == FileFormat::H5&&is_remote_path(uri))
		{

			 throw(domain_error("H5cytoframe doesn't support remote loading: " + uri));

		}
		else
		{

			if(fmt == FileFormat::H5)
				ptr.reset(new H5CytoFrame(uri, readonly));
			else
	#ifdef HAVE_TILEDB
				ptr.reset(new TileCytoFrame(uri, readonly, true, ctx));
	#else
			throw(domain_error("cytolib is not built with tiledb support!"));
	#endif

		}
		return ptr;
	}
	/**
	 * setter for channels (dedicated fro Rcpp API and  it won't throw on the unmatched old channel name)
	 * @param chnl_map
	 */
	void GatingHierarchy::set_channels(const CHANNEL_MAP & chnl_map)
	{
		//update comp
		comp.update_channels(chnl_map);

		//update gates

		VertexID_vec vertices=getVertices(0);

		for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
		{
			VertexID u=*it;
			nodeProperties & node=getNodeProperty(u);
			if(u!=0)
			{
				gatePtr g=node.getGate();
				if(g==NULL)
					throw(domain_error("no gate available for this node"));
				if(g_loglevel>=POPULATION_LEVEL)
					PRINT( "update channels for " +node.getName()+"\n");
				if(g->getType()!=BOOLGATE&&g->getType()!=LOGICALGATE&&g->getType()!=CLUSTERGATE)
					g->update_channels(chnl_map);
			}
		}

		//update local trans
		trans.update_channels(chnl_map);

		//update flag
		for(PARAM_VEC::iterator it = transFlag.begin(); it != transFlag.end(); it++)
			it->update_channels(chnl_map);

		//update fr
		frame_.set_channels(chnl_map);
	}

	/**
	 * compensate the data by the spillover provided by workspace
	 * or FCS TEXT keyword
	 * add prefix (e.g. Comp_ or <>) to channel name of the data
	 */
	void GatingHierarchy::compensate(CytoFrame & cytoframe)
	{
		if(comp.cid == "-2" || comp.cid == "")
		{
			if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("No compensation\n");
			return;
		}
		else if(comp.cid == "-1")
		{
			if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("Retrieve the Acquisition defined compensation matrix from FCS\n");
			/**
			 * compensate with spillover defined in keyword
			 */
			set_compensation(cytoframe.get_compensation(), false);
			//TODO:fix slash in compensation parameters for vX
			//currently I have not figured out the clean way to pass is_fix_slash flag from ws
			//this scenario may never occur so we won't bother the fix it until it bites us

		}

		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
			PRINT("Compensating...\n");

		cytoframe.compensate(comp);

		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
			PRINT("start prefixing data columns\n");

		for(const string & old : comp.marker)
		{
			cytoframe.set_channel(old, comp.prefix + old + comp.suffix);
		}

	}

	/**
	 * transform the data
	 * The reason we pass in MemCytoFrame is because the data member frame_ may not be finalized yet at this stage of parsing.
	 */
	void GatingHierarchy::transform_data(MemCytoFrame & cytoframe)
	{
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
			PRINT("start transforming data \n");
		if(cytoframe.n_rows()==0)
			throw(domain_error("data is not loaded yet!"));

	//	unsigned nEvents=fdata.nEvents;
	//	unsigned nChannls=fdata.nChannls;
		vector<string> channels=cytoframe.get_channels();
		int nEvents = cytoframe.n_rows();
		/*
		 * transforming each marker
		 */
		for(vector<string>::iterator it1=channels.begin();it1!=channels.end();it1++)
		{

			string curChannel=*it1;
			auto param_range = cytoframe.get_range(curChannel, ColType::channel, RangeType::instrument);
			TransPtr curTrans=trans.getTran(curChannel);

			if(curTrans)
			{
				if(curTrans->gateOnly())
					continue;

				EVENT_DATA_TYPE * x = cytoframe.get_data_memptr(curChannel, ColType::channel);
				if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				{
					string type;
					curTrans->getType(type);
					PRINT("transforming "+curChannel+" with func:"+type+"\n");
				}
				curTrans->transforming(x,nEvents);

				curTrans->transforming(&param_range.first, 1);

				curTrans->transforming(&param_range.second, 1);
			}

			cytoframe.set_range(curChannel, ColType::channel, param_range);

		}

	}

	void GatingHierarchy::calgate(MemCytoFrame & cytoframe, VertexID u, bool computeTerminalBool, INTINDICES &parentIndice)
	{
		nodeProperties & node=getNodeProperty(u);

		/*
		 * check if parent population is already gated
		 * because the boolgate's parent might be visited later than boolgate itself
		 */

	//	if(!parentNode.isGated())
	//	{
	//		if(g_loglevel>=POPULATION_LEVEL)
	//			PRINT("go to the ungated parent node:"+parentNode.getName()+"\n");
	//		calgate(cytoframe, pid, computeTerminalBool);
	//	}
	//

		if(g_loglevel>=POPULATION_LEVEL)
			PRINT("gating on:"+getNodePath(u)+"\n");

		gatePtr g=node.getGate();

		if(g==NULL)
			throw(domain_error("no gate available for this node"));

		/*
		 * calculate the indices for the current node
		 */

		switch(g->getType())
		{
		case BOOLGATE:
			{
				if(computeTerminalBool||getChildren(u).size()>0)
				{

					vector<bool> curIndices;
					curIndices=boolGating(cytoframe, u, computeTerminalBool);
					//combine with parent indices
					VertexID pid=getParent(u);
					nodeProperties & parentNode =getNodeProperty(pid);

					transform (curIndices.begin(), curIndices.end(), parentNode.getIndices().begin(), curIndices.begin(),logical_and<bool>());
					node.setIndices(curIndices);
				}
				else
				{
//					cout << computeTerminalBool << " : " << getChildren(u).size() <<  " : " << g->getType() << endl;
					return;
				}

				break;
			}
		case LOGICALGATE://skip any gating operation since the indice is already set once the gate is added
		case CLUSTERGATE:
			node.computeStats();
			return;
		default:
			{
				vector<unsigned> pind = parentIndice.getIndices_u();
				vector<unsigned> curIndices=g->gating(cytoframe, pind);
				node.setIndices(curIndices, parentIndice.getTotal());
			}

		}

		node.computeStats();
	}

	void GatingHierarchy::extendGate(MemCytoFrame & cytoframe, float extend_val){
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("\nstart extending Gates \n");

			if(cytoframe.n_rows()==0)
					throw(domain_error("data is not loaded yet!"));

			VertexID_vec vertices=getVertices(0);

			for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
			{
				VertexID u=*it;
				nodeProperties & node=getNodeProperty(u);
				if(u!=0)
				{
					gatePtr g=node.getGate();
					if(g==NULL)
						throw(domain_error("no gate available for this node"));
					if(g_loglevel>=POPULATION_LEVEL)
						PRINT(node.getName()+"\n");
					if(g->getType()!=BOOLGATE)
						g->extend(cytoframe,extend_val);
				}
			}
	}

	/**
	 *
	 * @param gh_pb
	 * @param cf_filename
	 * @param h5_opt
	 * @param is_skip_data whether to skip writing cytoframe data to pb. It is typically remain as default unless for debug purpose (e.g. re-writing gs that is loaded from legacy pb archive without actual data associated)
	 */
	void GatingHierarchy::convertToPb(pb::GatingHierarchy & gh_pb
			, string cf_filename
			, CytoFileOption h5_opt
			, bool is_skip_data
			, const CytoCtx & ctx){
		pb::populationTree * ptree = gh_pb.mutable_tree();
		/*
		 * cp tree
		 */
		VertexID_vec verIDs = getVertices();
		for(VertexID_vec::iterator it = verIDs.begin(); it != verIDs.end(); it++){
			VertexID thisVert = *it;
			nodeProperties & np = getNodeProperty(thisVert);

			pb::treeNodes * node = ptree->add_node();
			pb::nodeProperties * pb_node = node->mutable_node();
			bool isRoot = thisVert == 0;
			np.convertToPb(*pb_node, isRoot);
			if(!isRoot){
				node->set_parent(getParent(thisVert));
			}

		}
		//cp comp
		pb::COMP * comp_pb = gh_pb.mutable_comp();
		comp.convertToPb(*comp_pb);
		//cp trans
		pb::trans_local * trans_pb = gh_pb.mutable_trans();
		trans.convertToPb(*trans_pb);
		//cp trans flag
		BOOST_FOREACH(PARAM_VEC::value_type & it, transFlag){
			pb::PARAM * tflag_pb = gh_pb.add_transflag();
			it.convertToPb(*tflag_pb);
		}
		//fr
		if(!is_skip_data)
		{
			pb::CytoFrame * fr_pb = gh_pb.mutable_frame();
			bool flag = frame_.get_readonly();//get lock status
			frame_.set_readonly(false);//temporary unlock it
			frame_.flush_meta();
			frame_.set_readonly(flag);//restore the lock
			string ext;
			if(frame_.get_backend_type() == FileFormat::TILE)
				ext = ".tile";
			else
				ext = ".h5";

			frame_.convertToPb(*fr_pb, cf_filename + ext, h5_opt, ctx);
		}

	}

	//load legacy pb
	GatingHierarchy::GatingHierarchy(pb::GatingHierarchy & pb_gh, map<intptr_t, TransPtr> & trans_tbl){
		const pb::populationTree & tree_pb =  pb_gh.tree();
		int nNodes = tree_pb.node_size();

		tree = populationTree(nNodes);
		for(int i = 0; i < nNodes; i++){
			const pb::treeNodes & node_pb = tree_pb.node(i);
			const pb::nodeProperties & np_pb = node_pb.node();

			VertexID curChildID = i;
			tree[curChildID] = nodeProperties(np_pb);

			if(node_pb.has_parent()){
				VertexID parentID = node_pb.parent();
				boost::add_edge(parentID,curChildID,tree);
			}

		}
		//restore comp
		comp = compensation(pb_gh.comp());
		//restore trans flag
		for(int i = 0; i < pb_gh.transflag_size(); i++){
			transFlag.push_back(PARAM(pb_gh.transflag(i)));
		}
		//restore trans local
		trans = trans_local(pb_gh.trans(), trans_tbl);
	}
	/**
	 * add empty root node to the gating tree with the name set to 'root'
	 *
	 * \return the newly added root node Id
	 *
	 * For example:
	 * \code
	 * 	GatingHierarchy *curGh=new GatingHierarchy();
	 *	curGh->addRoot();
	 * \endcode
	 *
	 */
	VertexID GatingHierarchy::addRoot(){

		// Create  vertices in that graph
		VertexID u = boost::add_vertex(tree);
		nodeProperties & rootNode=tree[u];
		rootNode.setName("root");

		return(u);
	}

	/*
	 * this is for semi-automated pipeline to add population node associated with gate
	 * assuming gate split the parent population into two subpops, one of which is to keep
	 * depends on isNegate flag of the gate
	 */
	VertexID GatingHierarchy::addGate(gatePtr g,VertexID parentID,string popName)
	{

		//validity check
		int res = getChildren(parentID, popName);
		if( res >0 ){
			popName.append(" already exists!");
			throw(domain_error(popName));
		}else{
			VertexID curChildID = boost::add_vertex(tree);

			nodeProperties &curChild = tree[curChildID];
			curChild.setName(popName.c_str());
			curChild.setGate(g);
			if(g_loglevel>=POPULATION_LEVEL)
				PRINT("node created:"+curChild.getName()+"\n");
			//attach the populationNode to the boost node as property

			//add relation between current node and parent node
			boost::add_edge(parentID,curChildID,tree);
			return curChildID;
		}

	}
	/**
	 *
	 * It moves one node to the target parent.
	 *
	 * @param parent the target parent id
	 * @param child node id to be moved
	 */
	void GatingHierarchy::moveNode(string node, string parent){
		if(parent == node)
			throw(domain_error("Can't move the node to itself!"));

		VertexID cid = getNodeID(node), pid = getNodeID(parent);
		if(isDescendant(cid, pid))
			throw(domain_error("Can't move the node to its descendants!"));

		VertexID pid_old = getParent(cid);
		if(pid != pid_old)
		{
			boost::remove_edge(pid_old, cid, tree);
			boost::add_edge(pid, cid, tree);

		}

	}

	void GatingHierarchy::set_compensation(const compensation & _comp, bool is_update_prefix)
	{
		string prefix = comp.prefix;
		string suffix = comp.suffix;
		comp = _comp;
		//restore prefix
		if(!is_update_prefix)
		{
			comp.prefix = prefix;
			comp.suffix = suffix;
		}
	}
	void GatingHierarchy::set_compensation(compensation && _comp, bool is_update_prefix)
	{

		swap(comp, _comp);
		//restore prefix
		if(!is_update_prefix)
		{
			comp.prefix = _comp.prefix;
			comp.suffix = _comp.suffix;
		}
	}
	void GatingHierarchy::printLocalTrans(){
		PRINT("\nget trans from gating hierarchy\n");
		trans_map trans=this->trans.getTransMap();

		for (trans_map::iterator it=trans.begin();it!=trans.end();it++)
		{
			TransPtr curTrans=it->second;

			if(!curTrans->isInterpolated())
					throw(domain_error("non-interpolated calibration table:"+curTrans->getName()+curTrans->getChannel()+" from channel"+it->first));

			PRINT(it->first+curTrans->getName()+" "+curTrans->getChannel()+"\n");;

		}
	}

	/*
	 * the version without the need of loading data
	 * by supplying the extend_to value
	 */
	void GatingHierarchy::extendGate(float extend_val, float extend_to){
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("\nstart extending Gates \n");

			VertexID_vec vertices=getVertices(0);

			for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
			{
				VertexID u=*it;
				nodeProperties & node=getNodeProperty(u);
				if(u!=0)
				{
					gatePtr g=node.getGate();
					if(g==NULL)
						throw(domain_error("no gate available for this node"));
					if(g_loglevel>=POPULATION_LEVEL)
						PRINT(node.getName()+"\n");
					if(g->getType()!=BOOLGATE)
						g->extend(extend_val,extend_to);
				}
			}
	}
	/*
	 * adjust gates by gains
	 */
	void GatingHierarchy::adjustGate(map<string,float> &gains){
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("\nstart rescale Gates by gains \n");

			VertexID_vec vertices=getVertices(0);

			for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
			{
				VertexID u=*it;
				nodeProperties & node=getNodeProperty(u);
				if(u!=0)
				{
					gatePtr g=node.getGate();
					if(g==NULL)
						throw(domain_error("no gate available for this node"));
					if(g_loglevel>=POPULATION_LEVEL)
						PRINT(node.getName()+"\n");
					if(g->getType()!=BOOLGATE)
						g->gain(gains);
				}
			}
	}

	/*
	 * transform gates
	 */
	void GatingHierarchy::transform_gate(){
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("\nstart transform Gates \n");

			VertexID_vec vertices=getVertices(0);

			for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
			{
				VertexID u=*it;
				nodeProperties & node=getNodeProperty(u);
				if(u!=0)
				{
					gatePtr g=node.getGate();
					if(g==NULL)
						throw(domain_error("no gate available for this node"));
					if(g_loglevel>=POPULATION_LEVEL)
						PRINT(node.getName()+"\n");
					unsigned short gateType= g->getType();
					if(gateType == CURLYQUADGATE)
					{
						CurlyQuadGate& curlyGate = dynamic_cast<CurlyQuadGate&>(*g);
						curlyGate.interpolate(trans);//the interpolated polygon is in raw scale
					}
					if(gateType!=BOOLGATE)
						g->transforming(trans);

				}
			}
	}

	/*
	 * Apply post-transformation shifts to gates (needed for magnetic gates)
	 */
	void GatingHierarchy::shift_gate(){
		if(g_loglevel>=GATING_HIERARCHY_LEVEL)
				PRINT("\nstart applying shifts to gates \n");

			VertexID_vec vertices=getVertices(0);

			for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
			{
				VertexID u=*it;
				nodeProperties & node=getNodeProperty(u);
				if(u!=0)
				{
					gatePtr g=node.getGate();
					if(g==NULL)
						throw(domain_error("no gate available for this node"));
					if(g_loglevel>=POPULATION_LEVEL)
						PRINT(node.getName()+"\n");
					unsigned short gateType= g->getType();
					if(gateType!=BOOLGATE && gateType!=CLUSTERGATE && gateType!=LOGICALGATE)
						g->shiftGate();

				}
			}
	}

	void GatingHierarchy::check_ungated_bool_node(VertexID u){
		nodeProperties & node = getNodeProperty(u);
		if(!node.isGated())
		{
			if(node.getGate()->getType()==BOOLGATE)
			{
				MemCytoFrame fr;
				gating(fr, u);//pass dummy frame since boolgating doesn't need it once the initial gating was completed thus all the ref nodes are guaranteed to be gated
			}
		}
	}
	/*
	 * traverse the tree to gate each pops
	 * assuming data have already been compensated and transformed
	 *
	 */
	void GatingHierarchy::gating(MemCytoFrame & cytoframe, VertexID u,bool recompute, bool computeTerminalBool, bool skip_faulty_node)
	{
		//get parent ind
		INTINDICES parentIndice;

		if(u>0)
		{
			VertexID pid = getParent(u);
			nodeProperties & node=getNodeProperty(pid);
			/*
			 * check if current population is already gated (by boolGate)
			 *
			 */
			if(!node.isGated())
				gating(cytoframe, u, recompute, computeTerminalBool, skip_faulty_node);

			parentIndice = INTINDICES(node.getIndices());

		}

		gating(cytoframe, u, recompute, computeTerminalBool, skip_faulty_node, parentIndice);
	}
	void GatingHierarchy::gating(MemCytoFrame & cytoframe, VertexID u,bool recompute, bool computeTerminalBool, bool skip_faulty_node, INTINDICES &parentIndice)
	{

	//	if(!isLoaded)
	//			throw(domain_error("data is not loaded yet!"));

		nodeProperties & node=getNodeProperty(u);
		if(u==0)
		{
			node.setIndices(cytoframe.n_rows());
			node.computeStats();
		}else
		{
			/*
			 * check if current population is already gated (by boolGate)
			 *
			 */
			if(recompute||!node.isGated())
			{
				try{
					calgate(cytoframe, u, computeTerminalBool, parentIndice);
				}
				catch(const std::exception & e)
				{
					if(skip_faulty_node)
					{
						PRINT(e.what());
						auto path = getNodePath(u, false);
						PRINT("\n Skipping the faulty node '" + path + "' and its descendants \n");
//						removeNode(path);//can't simply remove here since vertex ID will change and affect gating recursive call
						return;
					}
					else
						throw(domain_error(e.what()));
				}
			}
		}

		//recursively gate all the descendants of u
		if(node.isGated())
		{
			INTINDICES pind(node.getIndices_u(), node.getTotal());
			VertexID_vec children=getChildren(u);
			for(VertexID_vec::iterator it=children.begin();it!=children.end();it++)
			{
				//add boost node
				VertexID curChildID = *it;

				gating(cytoframe, curChildID,recompute, computeTerminalBool, skip_faulty_node, pind);
			}

		}

	}
	/*
	 * bool gating operates on the indices of reference nodes
	 * because they are global, thus needs to be combined with parent indices
	 * in cases of negated gate (i.e. !A & !B)
	 * @param u
	 * @return
	 */

	vector<bool> GatingHierarchy::boolGating(MemCytoFrame & cytoframe, VertexID u, bool computeTerminalBool){

		nodeProperties & node=getNodeProperty(u);
		gatePtr  g=node.getGate();

		//init the indices
	//	unsigned nEvents=fdata.getEventsCount();

	//	vector<bool> ind(nEvents,true);
		/*it is kinda of expensive to init a long bool vector
		 *
		 */
		vector<bool> ind;
		/*
		 * combine the indices of reference populations
		 */

		vector<BOOL_GATE_OP> boolOpSpec=g->getBoolSpec();
		for(vector<BOOL_GATE_OP>::iterator it=boolOpSpec.begin();it!=boolOpSpec.end();it++)
		{
			/*
			 * find id of reference node
			 */
			VertexID nodeID;
			/*
			 * assume the reference node has already added during the parsing stage
			 */

			nodeID=getRefNodeID(u,it->path);

			nodeProperties & curPop=getNodeProperty(nodeID);
			//prevent self-referencing
			if(nodeID == u){
				string strErr = "The boolean gate is referencing to itself: ";
				strErr.append(curPop.getName());
				throw(domain_error(strErr));
			}

			if(!curPop.isGated())
			{
				if(g_loglevel>=POPULATION_LEVEL)
					PRINT("go to the ungated reference node:"+curPop.getName()+"\n");
				gating(cytoframe, nodeID, true, computeTerminalBool);
			}

			vector<bool> curPopInd=curPop.getIndices();
			if(it->isNot)
				curPopInd.flip();

			/*
			 * for the first reference node
			 * assign the indices directly without logical operation
			 */
			if(it==boolOpSpec.begin())
				ind=curPopInd;
			else
			{
				switch(it->op)
				{
					case '&':
						transform (ind.begin(), ind.end(), curPopInd.begin(), ind.begin(),logical_and<bool>());
						break;
					case '|':
						transform (ind.begin(), ind.end(), curPopInd.begin(), ind.begin(),logical_or<bool>());
						break;
					default:
						throw(domain_error("not supported operator!"));
				}
			}

		}

		if(g->isNegate())
			ind.flip();

		return ind;

	}
	/*
	 * external boolOpSpec can be provided .
	 * It is mainly used by openCyto rectRef gate
	 * (needs to be combined with parent indices)
	 *
	 * @param u
	 * @param boolOpSpec
	 * @return
	 */
	vector<bool> GatingHierarchy::boolGating(MemCytoFrame & cytoframe, vector<BOOL_GATE_OP> boolOpSpec, bool computeTerminalBool){

		vector<bool> ind;
		/*
		 * combine the indices of reference populations
		 */

		for(vector<BOOL_GATE_OP>::iterator it=boolOpSpec.begin();it!=boolOpSpec.end();it++)
		{
			/*
			 * find id of reference node
			 */
			VertexID nodeID;
			/*
			 * assume the reference node has already added during the parsing stage
			 */

			nodeID=getNodeID(it->path);//search ID by path

			nodeProperties & curPop=getNodeProperty(nodeID);

			if(!curPop.isGated())
			{
				if(g_loglevel>=POPULATION_LEVEL)
					PRINT("go to the ungated reference node:"+curPop.getName()+"\n");
				gating(cytoframe, nodeID, true, computeTerminalBool);
			}

			vector<bool> curPopInd=curPop.getIndices();
			if(it->isNot)
				curPopInd.flip();

			/*
			 * for the first reference node
			 * assign the indices directly without logical operation
			 */
			if(it==boolOpSpec.begin())
				ind=curPopInd;
			else
			{
				switch(it->op)
				{
					case '&':
						transform (ind.begin(), ind.end(), curPopInd.begin(), ind.begin(),logical_and<bool>());
						break;
					case '|':
						transform (ind.begin(), ind.end(), curPopInd.begin(), ind.begin(),logical_or<bool>());
						break;
					default:
						throw(domain_error("not supported operator!"));
				}
			}

		}

		return ind;

	}

	/*
	 * current output the graph in dot format
	 * and further covert it to gxl in order for Rgraphviz to read since it does not support dot directly
	 * right now the data exchange is through file system,it would be nice to do it in memory
	 */
	void GatingHierarchy::drawGraph(string output)
	{
		ofstream outputFile(output.c_str());

		boost::write_graphviz(outputFile,tree,OurVertexPropertyWriterR(tree));
		outputFile.close();

	}

	class custom_bfs_visitor : public boost::default_bfs_visitor
		{

		public:
			custom_bfs_visitor(VertexID_vec& v) : vlist(v) { }
			VertexID_vec & vlist;
		  template < typename Vertex, typename Graph >
		  void discover_vertex(Vertex u, const Graph & g) const
		  {
			  vlist.push_back(u);
		//	  v=u;
		  }

		};

	/**
	 * retrieve all the node IDs
	 *
	 * @param order accept 3 values: REGULAR(0) is the same original order by which nodes were added;
	 * 								 TSORT(1) topological order;
	 * 								 BFS(2) breadth first searching order
	 * @return a vector of node IDs
	 */

	VertexID_vec GatingHierarchy::getVertices(unsigned short order){

		VertexID_vec res, vertices;
		switch (order)
		{

			case REGULAR:
			{
				VertexIt it_begin,it_end;
				boost::tie(it_begin,it_end)=boost::vertices(tree);
				for(VertexIt it=it_begin;it!=it_end;it++)
					res.push_back((unsigned long)*it);
			}
			break;

			case TSORT:
			{
				boost::topological_sort(tree,back_inserter(vertices));
				for(VertexID_vec::reverse_iterator it=vertices.rbegin();it!=vertices.rend();it++)
					res.push_back(*it);
			}
			break;

			case BFS:
			{
				custom_bfs_visitor vis(res);
	//			vector<VertexID> p(num_vertices(tree));
	//			populationTree tree_copy(num_vertices(tree));
				boost::breadth_first_search(tree, vertex(0, tree)
											, boost::visitor(
													vis
	//												boost::make_bfs_visitor(boost::record_predecessors(&p[0]
	//																									 ,boost::on_tree_edge()
	//																									)
	//																					)
															)
											);
	//			res=vis.vlist;

			}
			break;

			default:
				throw(domain_error("not valid sort type for tree traversal!"));
		}

		return(res);

	}

	/**
	 * retrieve the VertexID by the gating path
	 * @param gatePath single string containing full(or partial) gating path
	 *
	 * For example:
	 * \code
	 * gh->getNodeID("singlet");
	 * gh->getNodeID("CD3/CD4+");
	 * \endcode
	  */
	VertexID GatingHierarchy::getNodeID(string gatePath){

		deque<string> res;
		boost::split(res,gatePath,boost::is_any_of("/"));
		//remove the empty string
		res.erase(remove_if(res.begin(),res.end(), isEmpty), res.end());

		//prepend root if it is absolute path (starts with /) and root is not yet prepended yet
		if(gatePath[0] == '/' && res[0]!= "root")
			res.insert(res.begin(), "root");
		return (getNodeID(res));

	}
	/*
	 * retrieve the VertexID by the gating path.
	 * this serves as a parser to convert generic gating path into internal node ID
	 * and it doesn't allow ambiguity (throw the exception when multiple nodes match)
	 * @param gatePath a string vector of full(or partial) gating path
	 * @return node id
	 */
	VertexID GatingHierarchy::getNodeID(const deque<string> & gatePath){
		VertexID_vec nodeIDs = queryByPath(0,gatePath);
		unsigned nMatches = nodeIDs.size();
		if(nMatches == 1)
				return nodeIDs[0];
		else{
			string errMsg;
			for(unsigned i = 0; i < gatePath.size()-1; i ++)
				errMsg.append(gatePath[i] + "/");
			errMsg.append(gatePath[gatePath.size()-1]);
			if(nMatches == 0)
				throw(domain_error(errMsg + " not found!" ));
			else
				throw(domain_error(errMsg + " is ambiguous within the gating tree!"));
		}

	}

	/*
	 *  find the most immediate common ancestor
	 *
	 * @param nodeIDs input node IDs
	 * @param nDepths the depths of ancestor node. It is used to as the measurement to determine which reference node to win when multiple matches occur.
	 * 											   The deeper it is, the nearer the reference is to the boolean node.
	 * @return ancestor ID
	 */
	VertexID GatingHierarchy::getCommonAncestor(VertexID_vec nodeIDs, unsigned & nDepths){

		unsigned nSize = nodeIDs.size();
		vector<VertexID_vec> paths(nSize) ;
		VertexID CommonAncestor = 0;

		/*
		 * calculate the levels of going up in order to find their common ancestor
		 */

		/*
		 * trace each node back to the root
		 */
		unsigned minDepths = numeric_limits<unsigned>::max() ;
		for(unsigned i = 0; i < nSize; i++)
		{
			unsigned counter = 0;
			for(VertexID curNode = nodeIDs[i]; curNode != 0; curNode = getParent(curNode))
			{
				paths[i].push_back(curNode);
				counter++;
			}
			minDepths = min(counter, minDepths);
		}
		/*
		 * walking the paths (top-down) simultaneously and stop at the first diverging node
		 */

		for(nDepths = 0; nDepths < minDepths; nDepths++)
		{
			//check if all nodes are the same at this level
			unsigned j = 0;
			unsigned pos = paths[j].size()-nDepths -1;//root is at the end of vector
			VertexID u = paths[j][pos];
			//loop through the rest nodes to see if they equal to u
			for(j = 1; j < nSize; j++)
			{
				unsigned pos = paths[j].size()-nDepths -1;
				if(paths[j][pos] != u)
					break;
			}

			if(j == nSize)
				CommonAncestor = u;//update result if all are the same
			else
				break;//otherwise, stop the loop at this level
		}

		return CommonAncestor;

	}

	VertexID GatingHierarchy::getRefNodeID(VertexID u, const deque<string> & refPath){

		/*
		 * to save searching time, try the siblings first(go up one level)
		 * which represents most of scenarios (e.g. cytokine boolean gates in ICS assay
		 */
		VertexID boolParentID = getParent(u);
		VertexID_vec nodeIDs = queryByPath(boolParentID,refPath);
		unsigned nMatches = nodeIDs.size();
		if(nMatches == 1)
			return nodeIDs[0];
		else
		{
			/*
			 * if failed to find reference from siblings, then go to more general searching logic (which takes more time)
			 */

			 nodeIDs = queryByPath(0,refPath);
			 nMatches = nodeIDs.size();
			if(nMatches == 1)
					return nodeIDs[0];
			else{
				string errMsg = "Reference node: ";
				for(unsigned i = 0; i < refPath.size()-1; i ++)
					errMsg.append(refPath[i] + "/");
				errMsg.append(refPath[refPath.size()-1]);
				if(nMatches == 0)
					throw(domain_error(errMsg + " not found!" ));
				else
				{
					/*
					 * select the nearest one to u when multiple nodes matches
					 * using bfs search starting from u
					 */
					unordered_set<VertexID> candidates(nodeIDs.begin(), nodeIDs.end());
//					for(auto i : candidates)
//						cout << "candidate ref: " << getNodePath(i) << endl;
					queue<VertexID> toVisit;
					unordered_set<VertexID> visited;
					bool found = false;
					toVisit.push(u);
					VertexID_vec nearest;
					while(!toVisit.empty())
					{
						int nSize = toVisit.size();
//						cout << "# of nodes of the current level: " << nSize << endl;
						//visit the current level
						while(nSize-- > 0)
						{
							//pop from front and mark it
							VertexID v = toVisit.front();
							toVisit.pop();
							visited.insert(v);
//							cout << "pop " << getNodePath(v) << endl;
							if(candidates.find(v)!=candidates.end())
							{
								found = true;
								nearest.push_back(v);
//								cout << "found " << getNodePath(v) << endl;
							}

							//add all its adjacent nodes to the list
							if(v>0)//add parent
							{
								VertexID p = getParent(v);
								if(visited.find(p)==visited.end())
								{
//									cout << "push parent " << getNodePath(p) << endl;
									toVisit.push(p);
								}

							}
							if(v != u)//skip the descendants of the boolgate itself
							{

								//add children
								VertexID_vec children = getChildren(v);
								for(auto i : children)
								{
									if(visited.find(i)==visited.end())
									{
	//									cout << "push child " << getNodePath(i) << endl;
										toVisit.push(i);
									}

								}
							}
						}
						//if already match to refPath in the current level, then stop searching
						if(found)
							break;
					}

					if(nearest.size() > 1){
						errMsg += " can't be determined due to the multiple matches with the same distance to boolean node!\n";
						for(auto & v : nearest)
						{
							errMsg += getNodePath(v) + "\n";
						}
						throw(domain_error(errMsg));

					}
					else{
						if(g_loglevel>=POPULATION_LEVEL)
							PRINT("RefNode is matched to the nearest neighbor: " + getNodePath(nearest[0])+"\n");

						return nearest[0];
					}

				}

			}

		}

	}

	VertexID_vec GatingHierarchy::pathMatch(VertexID_vec leafIDs, const deque<string> & gatePath){
				VertexID_vec res;
				/*
				 * bottom-up searching each route from matched leaf nodes
				 */
				VertexID_vec::iterator it_leaf,it_matched;
				it_matched = leafIDs.end();

				for(it_leaf = leafIDs.begin(); it_leaf != leafIDs.end(); it_leaf++)
				{
					/*
					 * bottom up matching to the given gating path
					 *
					 */

					VertexID curLeafID = *it_leaf;

					// start from the parent of the leaf node
					VertexID curNodeID = curLeafID;
					deque<string>::const_reverse_iterator it;
					for(it = gatePath.rbegin()+1;it!=gatePath.rend();it++)
					{
						//get current parent from node path
						string parentNameFromPath = *it;

						/*
						 * retrieve the actual parent node from the tree
						 */
						VertexID parentID = getParent(curNodeID);
						string parentName = getNodeProperty(parentID).getName();
						//compare it to the parent node from the path
						if(parentName.compare(parentNameFromPath) != 0)
						{
							break; //not matched then exit current route
						}else{
							//move up to the next ancestor and continue the matching process
							curNodeID = parentID;
						}
					}

					//when it succeeds to the end of path
					if(it == gatePath.rend())
						res.push_back(curLeafID);

				}

				return res;

	}
	/*
	 * retrieve the VertexIDs by the gating path.(bottom-up searching)
	 * This routine allows multiple matches
	 * @param ancestorID when gatePath is partial path, this node ID narrow the searching range.
	 * @param gatePath input
	 * @return node IDs that matches to the query path
	 */
	VertexID_vec GatingHierarchy::queryByPath(VertexID ancestorID, const deque<string> & gatePath){
		/*
		 * search for the leaf node
		 */
		string leafName=gatePath[gatePath.size()-1];
		VertexID_vec leafIDs=getDescendants(ancestorID,leafName);
		return pathMatch(leafIDs, gatePath);

	}

	/**
	 * check if v is the descendant of u
	 * @param u
	 * @param v
	 * @return
	 */
	bool GatingHierarchy::isDescendant(VertexID u, VertexID v){
		VertexID_vec nodesTomatch;
		custom_bfs_visitor vis(nodesTomatch);
		boost::breadth_first_search(tree, u, boost::visitor(vis));

		for(auto & it : nodesTomatch)
		{
			if(it == v)
				return true;
		}
		return false;
	}

	/**
	 * search for all the nodes that matches the pop name given the ancestor node id
	 * @param u the ancestor node id to search from
	 * @param name the node name to search for
	 * @return the vector of node id that match
	 *
	 * For example:
	 * \code
	 * VertexID parentID = gh->getNodeID("CD3");
	 * //this may return two descendants: "CD3/CD4/CCR7+ 45RA+" and "CD3/CD8/CCR7+ 45RA+"
	 * gh->getDescendants(parentID, "CCR7+ 45RA+");
	 * \endcode

	 */
	VertexID_vec GatingHierarchy::getDescendants(VertexID u,string name){
		VertexID_vec nodesTomatch, res;
		custom_bfs_visitor vis(nodesTomatch);
		boost::breadth_first_search(tree, u, boost::visitor(vis));
		VertexID_vec::iterator it;
		for(it=nodesTomatch.begin();it!=nodesTomatch.end();it++)
		{
			u=*it;
			if(getNodeProperty(u).getName().compare(name)==0)
				res.push_back(u);
		}
	//	if(it==nodesTomatch.end())
	//	{
	//		if(g_loglevel>=POPULATION_LEVEL)
	//			PRINT(name+" not found under the node: "+boost::lexical_cast<string>(u)+". returning the root instead.\n");;
	//		u=0;
	//	}
		return res;
	}

	/*
	 * retrieve VertexID that matches population name given an ancestor node
	 * It is used to search for the first node in the gating path (full or partial).
	 * This is different from getRefNodeID in the way that pop name must be uniquely identifiable in the tree.
	 * @param u the ancestor node id
	 * @param popName the population name to match
	 * @return node ID
	 */
	VertexID GatingHierarchy::getDescendant(VertexID u,string popName){

		/*
		 * top-down searching from that ancestor
		 */
		VertexID_vec res = getDescendants(u,popName);
		unsigned nMatches = res.size();

		switch (nMatches){
		case 0:
				popName.append(" not found within the gating tree!");
				throw(domain_error(popName));
		case 1:
				return (res[0]);

		default:
				popName.append(" is ambiguous within the gating tree!");
				throw(domain_error(popName));
		}

	}

	/**
	 * retrieve all the population paths
	 *
	 * The assumption is each node only has one parent.
	 *
	 * @param order passed to getVertices function
	 * @param fullPath flag indicates whether to return full path or partial path
	 * @param showHidden whether to include the hidden nodes
	 * @return
	 */
	vector<string> GatingHierarchy::getNodePaths(unsigned short order,bool fullPath,bool showHidden){

		VertexID_vec vertices=getVertices(order);
		vector<string> res;
		for(VertexID_vec::iterator it=vertices.begin();it!=vertices.end();it++)
		{
			VertexID u=*it;
			if(!showHidden&&getNodeProperty(u).getHiddenFlag())
				continue;
			res.push_back(getNodePath(u, fullPath));

		}
		return res;
	}

	/**
	 * Compute the depth of the given node
	 *
	 * @param u node ID
	 */
	unsigned GatingHierarchy::getNodeDepths(VertexID u){
		unsigned i = 0;
		while(u > 0){
			u = getParent(u);
			i++;
		}

		return i ;
	}
	/**
	 * Convert node Id to abs path
	 * @param u
	 * @return
	 */
	string GatingHierarchy::getNodePath(VertexID u,bool fullPath)
	{
		//init node path with the leaf
		string leafName=getNodeProperty(u).getName();
		string sNodePath = leafName;
		deque<string> nodePath;
		nodePath.push_front(sNodePath);

		//init searching routes
		VertexID_vec leafIDs=getDescendants(0,leafName);

		//start to trace back to ancestors
		while(u > 0)
		{
			//when full path is false, check if the current partial path is uniquely identifiable
			if(!fullPath)
			{
				//update searching routes
				leafIDs = pathMatch(leafIDs, nodePath);
				unsigned nMatches = leafIDs.size();
				if(nMatches == 1)
					break;//quit the path growing if unique
				else if(nMatches == 0)
					throw(domain_error(sNodePath + " not found!" ));

				// otherwise do nothing but continue to grow the path

			}

			sNodePath="/"+sNodePath;
			u=getParent(u);
			if(u>0)//don't append the root node
			{
				string pname = getNodeProperty(u).getName();
				nodePath.push_front(pname);
				sNodePath= pname + sNodePath;
			}

		}

		return sNodePath;

	}
	/**
	 * Get ancestor node for the given node
	 *
	 * Assume getParent only returns one parent node per GatingHierarchy
	 *
	 * @param u the given node ID
	 * @param level specify the distance from the given node
	 */
	VertexID GatingHierarchy::getAncestor(VertexID u,unsigned short level){

		for(unsigned short i=0;i<level;i++)
			u=getParent(u);
		return(u);
	}
	/*
	 * using boost in_edges out_edges to retrieve adjacent vertices
	 * assuming only one parent for each node
	 */
	EdgeID GatingHierarchy::getInEdges(VertexID target){
		vector<EdgeID> res;
		string err;
		err.append(boost::lexical_cast<string>(target));

		if(target<=boost::num_vertices(tree)-1)
		{

			boost::graph_traits<populationTree>::in_edge_iterator in_i, in_end;

			for (boost::tie(in_i, in_end) = in_edges(target,tree);
				         in_i != in_end; ++in_i)
			{
				EdgeID e = *in_i;
				res.push_back(e);
			}

		}
		else
			throw(domain_error(err+" :invalid vertexID!"));

		if(res.size()==0)
			throw(domain_error(err+" :parent not found!"));
		if(res.size()>1) //we only allow one parent per node
			throw(domain_error(err+" :multiple parent nodes found!"));

		return(res[0]);
	}

	/**
	 * Get parent node id for the given node
	 *
	 * @param target child ID
	 */
	VertexID GatingHierarchy::getParent(VertexID target){
		EdgeID e=getInEdges(target);
		return  boost::source(e, tree);
	}
	/**
	 * retrieve all children nodes
	 *
	 * @param source parent node ID
	 */
	VertexID_vec GatingHierarchy::getChildren(VertexID source){

		VertexID_vec res;
		if(source<=boost::num_vertices(tree)-1)
		{

			EdgeID e;
			boost::graph_traits<populationTree>::out_edge_iterator out_i, out_end;

			for (boost::tie(out_i, out_end) = out_edges(source,tree);
					 out_i != out_end; ++out_i)
				{
				  e = *out_i;
				  VertexID  targ = target(e, tree);
				  res.push_back(targ);
				}
		}
		else
		{
			PRINT("invalid vertexID:"+to_string(source)+"\n");
	//		res.push_back(0);
		}
		return(res);
	}

	/*
	 * retrieve single child node by parent id and child name.
	 * @param source id of the source node
	 * @param childName the child node name
	 * @return the child node id if succeeds; otherwise return -1.
	 */
	int GatingHierarchy::getChildren(VertexID source,string childName){

		int curNodeID;
		VertexID_vec children=getChildren(source);
		VertexID_vec::iterator it;
		for(it=children.begin();it!=children.end();it++)
		{
			curNodeID=*it;
			if(getNodeProperty(curNodeID).getName().compare(childName)==0)
				break;
		}
		if(it==children.end())
			curNodeID = -1;

		return(curNodeID);

	}

	/*
	 *
	 * make sure to use this API always since since it is safe way to access tree nodes due to the validity check
	 *
	 *since the vertex bundle should always exist as long as the  tree and node exist, thus it is safe
	 * to return the reference of it
	 */
	/**
	 * Retrieve the node properties
	 *
	 * It is the only way to access the gate, population indices and stats of the given node
	 * @param u node ID
	 * @return a reference to the nodeProperties object
	 */
	nodeProperties & GatingHierarchy::getNodeProperty(VertexID u){

		if(u<=boost::num_vertices(tree)-1)
			return(tree[u]);
		else
		{
			throw(out_of_range("returning empty node due to the invalid vertexID:" + boost::lexical_cast<std::string>(u)));

		}
	}

	GatingHierarchyPtr  GatingHierarchy::copy(bool is_copy_data, bool is_realize_data, const string & cf_filename) const{

		GatingHierarchyPtr res(new GatingHierarchy());

		res->comp=comp;
		res->tree=tree;
		res->transFlag = transFlag;
		res->trans = trans.copy();
		if(is_copy_data)
		{
			string ext;
			if(frame_.get_backend_type() == FileFormat::TILE)
				ext = ".tile";
			else
				ext = ".h5";

			if(is_realize_data)
				res->frame_ = frame_.copy_realized(cf_filename + ext);
			else
				res->frame_ = frame_.copy(cf_filename + ext);
		}
		else
			res->frame_ = frame_;
		return res;
	}

};

max_stars_repo_name: scignscape/PGVM

max_issues_repo_name: scignscape/PGVM

max_forks_repo_name: scignscape/PGVM


text:
/*********************************************************************
* Software License Agreement (BSD License)
*
*  Copyright (c) 2008, Willow Garage, Inc.
*  All rights reserved.
*
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following conditions
*  are met:
*
*   * Redistributions of source code must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*   * Redistributions in binary form must reproduce the above
*     copyright notice, this list of conditions and the following
*     disclaimer in the documentation and/or other materials provided
*     with the distribution.
*   * Neither the name of the Willow Garage nor the names of its
*     contributors may be used to endorse or promote products derived
*     from this software without specific prior written permission.
*
*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
*  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
*  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
*  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
*  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
*  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
*  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
*  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
*  POSSIBILITY OF SUCH DAMAGE.
*********************************************************************/
#include <boost/make_shared.hpp>
#include <boost/version.hpp>
#if ((BOOST_VERSION / 100) % 1000) >= 53
#include <boost/thread/lock_guard.hpp>
#endif

#include <ros/ros.h>
#include <nodelet/nodelet.h>
#include <image_transport/image_transport.h>
#include <sensor_msgs/image_encodings.h>
#include <dynamic_reconfigure/server.h>
#include <image_proc/DebayerConfig.h>

#include <opencv2/imgproc/imgproc.hpp>
// Until merged into OpenCV
#include "edge_aware.h"

#include <cv_bridge/cv_bridge.h>

namespace image_proc {

namespace enc = sensor_msgs::image_encodings;

class DebayerNodelet : public nodelet::Nodelet
{
  // ROS communication
  boost::shared_ptr<image_transport::ImageTransport> it_;
  image_transport::Subscriber sub_raw_;

  boost::mutex connect_mutex_;
  image_transport::Publisher pub_mono_;
  image_transport::Publisher pub_color_;

  // Dynamic reconfigure
  boost::recursive_mutex config_mutex_;
  typedef image_proc::DebayerConfig Config;
  typedef dynamic_reconfigure::Server<Config> ReconfigureServer;
  boost::shared_ptr<ReconfigureServer> reconfigure_server_;
  Config config_;

  virtual void onInit();

  void connectCb();

  void imageCb(const sensor_msgs::ImageConstPtr& raw_msg);

  void configCb(Config &config, uint32_t level);
};

void DebayerNodelet::onInit()
{
  ros::NodeHandle &nh         = getNodeHandle();
  ros::NodeHandle &private_nh = getPrivateNodeHandle();
  it_.reset(new image_transport::ImageTransport(nh));

  // Set up dynamic reconfigure
  reconfigure_server_.reset(new ReconfigureServer(config_mutex_, private_nh));
  ReconfigureServer::CallbackType f = boost::bind(&DebayerNodelet::configCb, this, _1, _2);
  reconfigure_server_->setCallback(f);

  // Monitor whether anyone is subscribed to the output
  typedef image_transport::SubscriberStatusCallback ConnectCB;
  ConnectCB connect_cb = boost::bind(&DebayerNodelet::connectCb, this);
  // Make sure we don't enter connectCb() between advertising and assigning to pub_XXX
  boost::lock_guard<boost::mutex> lock(connect_mutex_);
  pub_mono_  = it_->advertise("image_mono",  1, connect_cb, connect_cb);
  pub_color_ = it_->advertise("image_color", 1, connect_cb, connect_cb);
}

// Handles (un)subscribing when clients (un)subscribe
void DebayerNodelet::connectCb()
{
  boost::lock_guard<boost::mutex> lock(connect_mutex_);
  if (pub_mono_.getNumSubscribers() == 0 && pub_color_.getNumSubscribers() == 0)
    sub_raw_.shutdown();
  else if (!sub_raw_)
  {
    image_transport::TransportHints hints("raw", ros::TransportHints(), getPrivateNodeHandle());
    sub_raw_ = it_->subscribe("image_raw", 1, &DebayerNodelet::imageCb, this, hints);
  }
}

void DebayerNodelet::imageCb(const sensor_msgs::ImageConstPtr& raw_msg)
{
  int bit_depth = enc::bitDepth(raw_msg->encoding);
  //@todo Fix as soon as bitDepth fixes it
  if (raw_msg->encoding == enc::YUV422)
    bit_depth = 8;

  // First publish to mono if needed
  if (pub_mono_.getNumSubscribers())
  {
    if (enc::isMono(raw_msg->encoding))
      pub_mono_.publish(raw_msg);
    else
    {
      if ((bit_depth != 8) && (bit_depth != 16))
      {
        NODELET_WARN_THROTTLE(30,
                            "Raw image data from topic '%s' has unsupported depth: %d",
                            sub_raw_.getTopic().c_str(), bit_depth);
      } else {
        // Use cv_bridge to convert to Mono. If a type is not supported,
        // it will error out there
        sensor_msgs::ImagePtr gray_msg;
        try
        {
          if (bit_depth == 8)
            gray_msg = cv_bridge::toCvCopy(raw_msg, enc::MONO8)->toImageMsg();
          else
            gray_msg = cv_bridge::toCvCopy(raw_msg, enc::MONO16)->toImageMsg();
          pub_mono_.publish(gray_msg);
        }
        catch (cv_bridge::Exception &e)
        {
          NODELET_WARN_THROTTLE(30, "cv_bridge conversion error: '%s'", e.what());
        }
      }
    }
  }

  // Next, publish to color
  if (!pub_color_.getNumSubscribers())
    return;

  if (enc::isMono(raw_msg->encoding))
  {
    // For monochrome, no processing needed!
    pub_color_.publish(raw_msg);

    // Warn if the user asked for color
    NODELET_WARN_THROTTLE(30,
                            "Color topic '%s' requested, but raw image data from topic '%s' is grayscale",
                            pub_color_.getTopic().c_str(), sub_raw_.getTopic().c_str());
  }
  else if (enc::isColor(raw_msg->encoding))
  {
    pub_color_.publish(raw_msg);
  }
  else if (enc::isBayer(raw_msg->encoding)) {
    int type = bit_depth == 8 ? CV_8U : CV_16U;
    const cv::Mat bayer(raw_msg->height, raw_msg->width, CV_MAKETYPE(type, 1),
                        const_cast<uint8_t*>(&raw_msg->data[0]), raw_msg->step);

      sensor_msgs::ImagePtr color_msg = boost::make_shared<sensor_msgs::Image>();
      color_msg->header   = raw_msg->header;
      color_msg->height   = raw_msg->height;
      color_msg->width    = raw_msg->width;
      color_msg->encoding = bit_depth == 8? enc::BGR8 : enc::BGR16;
      color_msg->step     = color_msg->width * 3 * (bit_depth / 8);
      color_msg->data.resize(color_msg->height * color_msg->step);

      cv::Mat color(color_msg->height, color_msg->width, CV_MAKETYPE(type, 3),
                    &color_msg->data[0], color_msg->step);

      int algorithm;
      {
        boost::lock_guard<boost::recursive_mutex> lock(config_mutex_);
        algorithm = config_.debayer;
      }

      if (algorithm == Debayer_EdgeAware ||
          algorithm == Debayer_EdgeAwareWeighted)
      {
        // These algorithms are not in OpenCV yet
        if (raw_msg->encoding != enc::BAYER_GRBG8)
        {
          NODELET_WARN_THROTTLE(30, "Edge aware algorithms currently only support GRBG8 Bayer. "
                                "Falling back to bilinear interpolation.");
          algorithm = Debayer_Bilinear;
        }
        else
        {
          if (algorithm == Debayer_EdgeAware)
            debayerEdgeAware(bayer, color);
          else
            debayerEdgeAwareWeighted(bayer, color);
        }
      }
      if (algorithm == Debayer_Bilinear ||
          algorithm == Debayer_VNG)
      {
        int code = -1;
        if (raw_msg->encoding == enc::BAYER_RGGB8 ||
            raw_msg->encoding == enc::BAYER_RGGB16)
          code = cv::COLOR_BayerBG2BGR;
        else if (raw_msg->encoding == enc::BAYER_BGGR8 ||
                 raw_msg->encoding == enc::BAYER_BGGR16)
          code = cv::COLOR_BayerRG2BGR;
        else if (raw_msg->encoding == enc::BAYER_GBRG8 ||
                 raw_msg->encoding == enc::BAYER_GBRG16)
          code = cv::COLOR_BayerGR2BGR;
        else if (raw_msg->encoding == enc::BAYER_GRBG8 ||
                 raw_msg->encoding == enc::BAYER_GRBG16)
          code = cv::COLOR_BayerGB2BGR;

        if (algorithm == Debayer_VNG)
          code += cv::COLOR_BayerBG2BGR_VNG - cv::COLOR_BayerBG2BGR;

        cv::cvtColor(bayer, color, code);
      }

      pub_color_.publish(color_msg);
  }
  else if (raw_msg->encoding == enc::YUV422)
  {
    // Use cv_bridge to convert to BGR8
    sensor_msgs::ImagePtr color_msg;
    try
    {
      color_msg = cv_bridge::toCvCopy(raw_msg, enc::BGR8)->toImageMsg();
      pub_color_.publish(color_msg);
    }
    catch (cv_bridge::Exception &e)
    {
      NODELET_WARN_THROTTLE(30, "cv_bridge conversion error: '%s'", e.what());
    }
  }
  else if (raw_msg->encoding == enc::TYPE_8UC3)
  {
    // 8UC3 does not specify a color encoding. Is it BGR, RGB, HSV, XYZ, LUV...?
    NODELET_ERROR_THROTTLE(10,
                           "Raw image topic '%s' has ambiguous encoding '8UC3'. The "
                           "source should set the encoding to 'bgr8' or 'rgb8'.",
                           sub_raw_.getTopic().c_str());
  }
  else
  {
    NODELET_ERROR_THROTTLE(10, "Raw image topic '%s' has unsupported encoding '%s'",
                           sub_raw_.getTopic().c_str(), raw_msg->encoding.c_str());
  }
}

void DebayerNodelet::configCb(Config &config, uint32_t level)
{
  config_ = config;
}

} // namespace image_proc

// Register nodelet
#include <pluginlib/class_list_macros.h>
PLUGINLIB_EXPORT_CLASS( image_proc::DebayerNodelet, nodelet::Nodelet)


text:
// Copyright András Vukics 2006–2020. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
#include "ExampleMode.h"

#include "ElementLiouvillean.h"
#include "Tridiagonal.tcc"
#include "TridiagonalHamiltonian.h"

#include <boost/bind.hpp>

using boost::bind;

void aJump   (StateVectorLow&, double kappa_nPlus1);
void aDagJump(StateVectorLow&, double kappa_n     );

double aJumpRate   (const LazyDensityOperator&, double kappa_nPlus1);
double aDagJumpRate(const LazyDensityOperator&, double kappa_n     );

basic::PumpedLossyMode::PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
  : Free(cutoff,{},
         {
           CF{"(kappa*(2*nTh+1),delta)",dcomp(kappa*(2*nTh+1),delta),     cutoff },
           CF{"eta"                    ,eta                         ,sqrt(cutoff)}
         }),
    TridiagonalHamiltonian<1,false>(dcomp(-kappa*(2*nTh+1),delta)*nop(cutoff)
                                    +
                                    tridiagPlusHC_overI(conj(eta)*aop(cutoff))),
    ElementLiouvilleanStrategies<1,2>(JumpStrategies(bind(aJump   ,_1,kappa*(nTh+1)),
                                                     bind(aDagJump,_1,kappa* nTh   )),
                                      JumpRateStrategies(bind(aJumpRate   ,_1,kappa*(nTh+1)),
                                                         bind(aDagJumpRate,_1,kappa* nTh   )),
                                      "Mode",{"photon loss","photon absorption"}),
    ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"})
{
  getParsStream()<<"Pumped lossy mode";
}

auto basic::PumpedLossyMode::average_v(NoTime, const LazyDensityOperator& matrix) const -> const Averages
{
  auto averages(initializedAverages());

  averages(0)=aJumpRate(matrix,1);

  for (int n=1; n<int(matrix.getDimension()); n++) {
    double sqrtn=sqrt(n);
    dcomp offdiag(matrix(n)(n-1));
    averages(1)+=sqrtn*real(offdiag);
    averages(2)+=sqrtn*imag(offdiag);
  }

  return averages;

}

const Tridiagonal::Diagonal mainDiagonal(dcomp z, size_t cutoff);

basic::PumpedLossyModeIP::PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
  : Free(cutoff,
         {
           RF{"kappa*(2*nTh+1)",kappa*(2*nTh+1),cutoff},
           RF{"delta",delta,1}
         },
         {CF{"eta",eta,sqrt(cutoff)}}),
    FreeExact<false>(cutoff),
    TridiagonalHamiltonian<1,true>(furnishWithFreqs(tridiagPlusHC_overI(conj(eta)*aop(cutoff)),
                                                    mainDiagonal(dcomp(kappa*(2*nTh+1),-delta),cutoff))),
    ElementLiouvilleanStrategies<1,2>(JumpStrategies(bind(aJump   ,_1,kappa*(nTh+1)),
                                                     bind(aDagJump,_1,kappa* nTh   )),
                                      JumpRateStrategies(bind(aJumpRate   ,_1,kappa*(nTh+1)),
                                                         bind(aDagJumpRate,_1,kappa* nTh   )),
                                      "Mode",{"photon loss","photon absorption"}),
    ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"}),
    z_(kappa*(2*nTh+1),-delta)
{}

void basic::PumpedLossyModeIP::updateU(OneTime dtDid) const
{
  getDiagonal()=exp(-z_*(dtDid*blitz::tensor::i));
}

// PumpedLossyModeIP::average_v exactly the same as PumpedLossyMode::average_v above
auto basic::PumpedLossyModeIP::average_v(NoTime, const LazyDensityOperator& matrix) const -> const Averages
{
  auto averages(initializedAverages());

  averages(0)=aJumpRate(matrix,1);

  for (int n=1; n<int(matrix.getDimension()); n++) {
    double sqrtn=sqrt(n);
    dcomp offdiag(matrix(n)(n-1));
    averages(1)+=sqrtn*real(offdiag);
    averages(2)+=sqrtn*imag(offdiag);
  }

  return averages;

}

hierarchical::ModeBase::ModeBase(double kappa, double nTh, size_t cutoff)
  : Free(cutoff,{RF{"kappa*(2*nTh+1)",kappa*(2*nTh+1),cutoff}}),
    ElementLiouvillean<1,2>("Mode",{"photon loss","photon absorption"}),
    ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"}),
    kappa_(kappa), nTh_(nTh)
{}

void hierarchical::ModeBase::doActWithJ(NoTime, StateVectorLow& psi, LindbladNo<0>) const
{
  aJump(psi,kappa_*(nTh_+1));
}

void hierarchical::ModeBase::doActWithJ(NoTime, StateVectorLow& psi, LindbladNo<1>) const
{
  aDagJump(psi,kappa_*nTh_);
}

double hierarchical::ModeBase::rate(NoTime, const LazyDensityOperator& matrix, LindbladNo<0>) const
{
  return aJumpRate(matrix,kappa_*(nTh_+1));
}

double hierarchical::ModeBase::rate(NoTime, const LazyDensityOperator& matrix, LindbladNo<1>) const
{
  return aDagJumpRate(matrix,kappa_*nTh_);
}

// ModeBase::average_v exactly the same as PumpedLossyMode::average_v above
auto hierarchical::ModeBase::average_v (NoTime, const hierarchical::ModeBase::LazyDensityOperator& matrix) const -> const Averages
{
  auto averages(initializedAverages());

  averages(0)=aJumpRate(matrix,1);

  for (int n=1; n<int(matrix.getDimension()); n++) {
    double sqrtn=sqrt(n);
    dcomp offdiag(matrix(n)(n-1));
    averages(1)+=sqrtn*real(offdiag);
    averages(2)+=sqrtn*imag(offdiag);
  }

  return averages;

}

hierarchical::PumpedLossyMode::PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
  : ModeBase(kappa,nTh,cutoff),
    TridiagonalHamiltonian<1,false>(dcomp(-kappa*(2*nTh+1),delta)*nop(cutoff)
                                    +
                                    tridiagPlusHC_overI(conj(eta)*aop(cutoff)))
{}

hierarchical::PumpedLossyModeIP::PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
  : ModeBase(kappa,nTh,cutoff),
    FreeExact<false>(cutoff),
    TridiagonalHamiltonian<1,true>(furnishWithFreqs(tridiagPlusHC_overI(conj(eta)*aop(cutoff)),
                                                    mainDiagonal(dcomp(kappa*(2*nTh+1),-delta),cutoff))),
    z_(kappa*(2*nTh+1),-delta)
{}

// PumpedLossyModeIP::updateU exactly the same as above
void hierarchical::PumpedLossyModeIP::updateU(OneTime dtDid) const
{
  getDiagonal()=exp(-z_*(dtDid*blitz::tensor::i));
}

ext: cc

max_stars_repo_name: bartoszek/cppqed

max_issues_repo_name: bartoszek/cppqed

max_forks_repo_name: bartoszek/cppqed


text:
/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
// text_oarchive_impl.ipp:

// (C) Copyright 2002 Robert Ramey - http://www.rrsd.com .
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

//  See http://www.boost.org for updates, documentation, and revision history.

#include <string>
#include <boost/config.hpp>
#include <locale>
#include <cstddef> // size_t

#include <boost/config.hpp>
#if defined(BOOST_NO_STDC_NAMESPACE)
namespace std{
    using ::size_t;
} // namespace std
#endif

#ifndef BOOST_NO_CWCHAR
#include <cwchar>
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std{ using ::wcslen; }
#endif
#endif

#include <boost/archive/add_facet.hpp>
#include <boost/archive/text_oarchive.hpp>

namespace boost {
namespace archive {

//////////////////////////////////////////////////////////////////////
// implementation of basic_text_oprimitive overrides for the combination
// of template parameters used to create a text_oprimitive

template<class Archive>
BOOST_ARCHIVE_DECL void
text_oarchive_impl<Archive>::save(const char * s)
{
    const std::size_t len = std::ostream::traits_type::length(s);
    *this->This() << len;
    this->This()->newtoken();
    os << s;
}

template<class Archive>
BOOST_ARCHIVE_DECL void
text_oarchive_impl<Archive>::save(const std::string &s)
{
    const std::size_t size = s.size();
    *this->This() << size;
    this->This()->newtoken();
    os << s;
}

#ifndef BOOST_NO_CWCHAR
#ifndef BOOST_NO_INTRINSIC_WCHAR_T
template<class Archive>
BOOST_ARCHIVE_DECL void
text_oarchive_impl<Archive>::save(const wchar_t * ws)
{
    const std::size_t l = std::wcslen(ws);
    * this->This() << l;
    this->This()->newtoken();
    os.write((const char *)ws, l * sizeof(wchar_t)/sizeof(char));
}
#endif

#ifndef BOOST_NO_STD_WSTRING
template<class Archive>
BOOST_ARCHIVE_DECL void
text_oarchive_impl<Archive>::save(const std::wstring &ws)
{
    const std::size_t l = ws.size();
    * this->This() << l;
    this->This()->newtoken();
    os.write((const char *)(ws.data()), l * sizeof(wchar_t)/sizeof(char));
}
#endif
#endif // BOOST_NO_CWCHAR

template<class Archive>
BOOST_ARCHIVE_DECL
text_oarchive_impl<Archive>::text_oarchive_impl(
    std::ostream & os,
    unsigned int flags
) :
    basic_text_oprimitive<std::ostream>(
        os,
        0 != (flags & no_codecvt)
    ),
    basic_text_oarchive<Archive>(flags)
{
    if(0 == (flags & no_header))
        #if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3205))
        this->init();
        #else
        this->basic_text_oarchive<Archive>::init();
        #endif
}

template<class Archive>
BOOST_ARCHIVE_DECL void
text_oarchive_impl<Archive>::save_binary(const void *address, std::size_t count){
    put('\n');
    this->end_preamble();
    #if ! defined(__MWERKS__)
    this->basic_text_oprimitive<std::ostream>::save_binary(
    #else
    this->basic_text_oprimitive::save_binary(
    #endif
        address,
        count
    );
    this->delimiter = this->eol;
}

} // namespace archive
} // namespace boost

ext: ipp


text:
#include <tiny_functions.h>

#define BOOST_AUTO_TEST_MAIN
#include <boost/test/auto_unit_test.hpp>
#include <boost/test/unit_test_suite.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/test/test_tools.hpp>

BOOST_AUTO_TEST_SUITE(tiny_clamp);

  BOOST_AUTO_TEST_CASE(double_testing)
  {
    double const min_val = 0.1;
    double const max_val = 2.9;
    double clamped_val;

    // normal clamp tests
    clamped_val = tiny::clamp(1.5, min_val, max_val);
    BOOST_CHECK( clamped_val == 1.5 );

    clamped_val = tiny::clamp(-1.5, min_val, max_val);
    BOOST_CHECK( clamped_val == min_val );

    clamped_val = tiny::clamp(3.0, min_val, max_val);
    BOOST_CHECK( clamped_val == max_val );

    // clamp minimum tests
    clamped_val = tiny::clamp_min(1.5, min_val);
    BOOST_CHECK( clamped_val == 1.5 );

    clamped_val = tiny::clamp_min(-1.5, min_val);
    BOOST_CHECK( clamped_val == min_val );

    clamped_val = tiny::clamp_min(3.0, min_val);
    BOOST_CHECK( clamped_val == 3.0 );

    // clamp maximum tests
    clamped_val = tiny::clamp_max(1.5, max_val);
    BOOST_CHECK( clamped_val == 1.5 );

    clamped_val = tiny::clamp_max(-1.5, max_val);
    BOOST_CHECK( clamped_val == -1.5 );

    clamped_val = tiny::clamp_max(3.0, max_val);
    BOOST_CHECK( clamped_val == max_val );

    // clamp zero..one tests
    clamped_val = tiny::clamp_zero_one(0.5);
    BOOST_CHECK( clamped_val == 0.5 );

    clamped_val = tiny::clamp_zero_one(-1.5);
    BOOST_CHECK( clamped_val == 0.0 );

    clamped_val = tiny::clamp_zero_one(3.0);
    BOOST_CHECK( clamped_val == 1.0 );
  }

  BOOST_AUTO_TEST_CASE(float_testing)
  {
    float const min_val = 0.1f;
    float const max_val = 2.9f;
    float clamped_val;

    // normal clamp tests
    clamped_val = tiny::clamp(1.5f, min_val, max_val);
    BOOST_CHECK( clamped_val == 1.5f );

    clamped_val = tiny::clamp(-1.5f, min_val, max_val);
    BOOST_CHECK( clamped_val == min_val );

    clamped_val = tiny::clamp(3.0f, min_val, max_val);
    BOOST_CHECK( clamped_val == max_val );

    // clamp minimum tests
    clamped_val = tiny::clamp_min(1.5f, min_val);
    BOOST_CHECK( clamped_val == 1.5f );

    clamped_val = tiny::clamp_min(-1.5f, min_val);
    BOOST_CHECK( clamped_val == min_val );

    clamped_val = tiny::clamp_min(3.0f, min_val);
    BOOST_CHECK( clamped_val == 3.0f );

    // clamp maximum tests
    clamped_val = tiny::clamp_max(1.5f, max_val);
    BOOST_CHECK( clamped_val == 1.5f );

    clamped_val = tiny::clamp_max(-1.5f, max_val);
    BOOST_CHECK( clamped_val == -1.5f );

    clamped_val = tiny::clamp_max(3.0f, max_val);
    BOOST_CHECK( clamped_val == max_val );

    // clamp zero..one tests
    clamped_val = tiny::clamp_zero_one(0.5f);
    BOOST_CHECK( clamped_val == 0.5f );

    clamped_val = tiny::clamp_zero_one(-1.5f);
    BOOST_CHECK( clamped_val == 0.0f );

    clamped_val = tiny::clamp_zero_one(3.0f);
    BOOST_CHECK( clamped_val == 1.0f );
  }

BOOST_AUTO_TEST_SUITE_END();

max_stars_repo_name: diku-dk/PROX

max_issues_repo_name: erleben/matchstick

max_forks_repo_name: erleben/matchstick


text:
// Copyright 2020 Hans Dembinski
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_HISTOGRAM_DETAIL_REDUCE_COMMAND_HPP
#define BOOST_HISTOGRAM_DETAIL_REDUCE_COMMAND_HPP

#include <boost/histogram/detail/span.hpp>
#include <boost/histogram/fwd.hpp>
#include <boost/throw_exception.hpp>
#include <cassert>
#include <stdexcept>
#include <string>

namespace boost {
namespace histogram {
namespace detail {

struct reduce_command {
  static constexpr unsigned unset = static_cast<unsigned>(-1);
  unsigned iaxis = unset;
  enum class range_t : char {
    none,
    indices,
    values,
  } range = range_t::none;
  union {
    axis::index_type index;
    double value;
  } begin{0}, end{0};
  unsigned merge = 0; // default value indicates unset option
  bool crop = false;
  // for internal use by the reduce algorithm
  bool is_ordered = true;
  bool use_underflow_bin = true;
  bool use_overflow_bin = true;
};

// - puts commands in correct axis order
// - sets iaxis for positional commands
// - detects and fails on invalid settings
// - fuses merge commands with non-merge commands
inline void normalize_reduce_commands(span<reduce_command> out,
                                      span<const reduce_command> in) {
  unsigned iaxis = 0;
  for (const auto& o_in : in) {
    assert(o_in.merge > 0);
    if (o_in.iaxis != reduce_command::unset && o_in.iaxis >= out.size())
      BOOST_THROW_EXCEPTION(std::invalid_argument("invalid axis index"));
    auto& o_out = out.begin()[o_in.iaxis == reduce_command::unset ? iaxis : o_in.iaxis];
    if (o_out.merge == 0) {
      o_out = o_in;
    } else {
      // Some command was already set for this axis, try to fuse commands.
      if (!((o_in.range == reduce_command::range_t::none) ^
            (o_out.range == reduce_command::range_t::none)) ||
          (o_out.merge > 1 && o_in.merge > 1))
        BOOST_THROW_EXCEPTION(std::invalid_argument(
            "multiple conflicting reduce commands for axis " +
            std::to_string(o_in.iaxis == reduce_command::unset ? iaxis : o_in.iaxis)));
      if (o_in.range != reduce_command::range_t::none) {
        o_out.range = o_in.range;
        o_out.begin = o_in.begin;
        o_out.end = o_in.end;
      } else {
        o_out.merge = o_in.merge;
      }
    }
    ++iaxis;
  }

  iaxis = 0;
  for (auto& o : out) o.iaxis = iaxis++;
}

} // namespace detail
} // namespace histogram
} // namespace boost

#endif


text:
//
// Copyright (c) 2016-2017 Vinnie Falco (vinnie dot falco at gmail dot com)
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//
// Official repository: https://github.com/boostorg/beast
//

#ifndef BOOST_BEAST_HTTP_IMPL_VERB_IPP
#define BOOST_BEAST_HTTP_IMPL_VERB_IPP

#include <boost/beast/core/detail/config.hpp>
#include <boost/throw_exception.hpp>
#include <stdexcept>

namespace boost {
namespace beast {
namespace http {

namespace detail {

template<class = void>
inline
string_view
verb_to_string(verb v)
{
    switch(v)
    {
    case verb::delete_:       return "DELETE";
    case verb::get:           return "GET";
    case verb::head:          return "HEAD";
    case verb::post:          return "POST";
    case verb::put:           return "PUT";
    case verb::connect:       return "CONNECT";
    case verb::options:       return "OPTIONS";
    case verb::trace:         return "TRACE";

    case verb::copy:          return "COPY";
    case verb::lock:          return "LOCK";
    case verb::mkcol:         return "MKCOL";
    case verb::move:          return "MOVE";
    case verb::propfind:      return "PROPFIND";
    case verb::proppatch:     return "PROPPATCH";
    case verb::search:        return "SEARCH";
    case verb::unlock:        return "UNLOCK";
    case verb::bind:          return "BIND";
    case verb::rebind:        return "REBIND";
    case verb::unbind:        return "UNBIND";
    case verb::acl:           return "ACL";

    case verb::report:        return "REPORT";
    case verb::mkactivity:    return "MKACTIVITY";
    case verb::checkout:      return "CHECKOUT";
    case verb::merge:         return "MERGE";

    case verb::msearch:       return "M-SEARCH";
    case verb::notify:        return "NOTIFY";
    case verb::subscribe:     return "SUBSCRIBE";
    case verb::unsubscribe:   return "UNSUBSCRIBE";

    case verb::patch:         return "PATCH";
    case verb::purge:         return "PURGE";

    case verb::mkcalendar:    return "MKCALENDAR";

    case verb::link:          return "LINK";
    case verb::unlink:        return "UNLINK";

    case verb::unknown:
        return "<unknown>";
    }

    BOOST_THROW_EXCEPTION(std::invalid_argument{"unknown verb"});
}

template<class = void>
verb
string_to_verb(string_view v)
{
/*
    ACL
    BIND
    CHECKOUT
    CONNECT
    COPY
    DELETE
    GET
    HEAD
    LINK
    LOCK
    M-SEARCH
    MERGE
    MKACTIVITY
    MKCALENDAR
    MKCOL
    MOVE
    NOTIFY
    OPTIONS
    PATCH
    POST
    PROPFIND
    PROPPATCH
    PURGE
    PUT
    REBIND
    REPORT
    SEARCH
    SUBSCRIBE
    TRACE
    UNBIND
    UNLINK
    UNLOCK
    UNSUBSCRIBE
*/
    if(v.size() < 3)
        return verb::unknown;
    // s must be null terminated
    auto const eq =
        [](string_view sv, char const* s)
        {
            auto p = sv.data();
            for(;;)
            {
                if(*s != *p)
                    return false;
                ++s;
                ++p;
                if(! *s)
                    return p == sv.end();
            }
        };
    auto c = v[0];
    v.remove_prefix(1);
    switch(c)
    {
    case 'A':
        if(v == "CL")
            return verb::acl;
        break;

    case 'B':
        if(v == "IND")
            return verb::bind;
        break;

    case 'C':
        c = v[0];
        v.remove_prefix(1);
        switch(c)
        {
        case 'H':
            if(eq(v, "ECKOUT"))
                return verb::checkout;
            break;

        case 'O':
            if(eq(v, "NNECT"))
                return verb::connect;
            if(eq(v, "PY"))
                return verb::copy;
            BOOST_BEAST_FALLTHROUGH;

        default:
            break;
        }
        break;

    case 'D':
        if(eq(v, "ELETE"))
            return verb::delete_;
        break;

    case 'G':
        if(eq(v, "ET"))
            return verb::get;
        break;

    case 'H':
        if(eq(v, "EAD"))
            return verb::head;
        break;

    case 'L':
        if(eq(v, "INK"))
            return verb::link;
        if(eq(v, "OCK"))
            return verb::lock;
        break;

    case 'M':
        c = v[0];
        v.remove_prefix(1);
        switch(c)
        {
        case '-':
            if(eq(v, "SEARCH"))
                return verb::msearch;
            break;

        case 'E':
            if(eq(v, "RGE"))
                return verb::merge;
            break;

        case 'K':
            if(eq(v, "ACTIVITY"))
                return verb::mkactivity;
            if(v[0] == 'C')
            {
                v.remove_prefix(1);
                if(eq(v, "ALENDAR"))
                    return verb::mkcalendar;
                if(eq(v, "OL"))
                    return verb::mkcol;
                break;
            }
            break;

        case 'O':
            if(eq(v, "VE"))
                return verb::move;
            BOOST_BEAST_FALLTHROUGH;

        default:
            break;
        }
        break;

    case 'N':
        if(eq(v, "OTIFY"))
            return verb::notify;
        break;

    case 'O':
        if(eq(v, "PTIONS"))
            return verb::options;
        break;

    case 'P':
        c = v[0];
        v.remove_prefix(1);
        switch(c)
        {
        case 'A':
            if(eq(v, "TCH"))
                return verb::patch;
            break;

        case 'O':
            if(eq(v, "ST"))
                return verb::post;
            break;

        case 'R':
            if(eq(v, "OPFIND"))
                return verb::propfind;
            if(eq(v, "OPPATCH"))
                return verb::proppatch;
            break;

        case 'U':
            if(eq(v, "RGE"))
                return verb::purge;
            if(eq(v, "T"))
                return verb::put;
            BOOST_BEAST_FALLTHROUGH;

        default:
            break;
        }
        break;

    case 'R':
        if(v[0] != 'E')
            break;
        v.remove_prefix(1);
        if(eq(v, "BIND"))
            return verb::rebind;
        if(eq(v, "PORT"))
            return verb::report;
        break;

    case 'S':
        if(eq(v, "EARCH"))
            return verb::search;
        if(eq(v, "UBSCRIBE"))
            return verb::subscribe;
        break;

    case 'T':
        if(eq(v, "RACE"))
            return verb::trace;
        break;

    case 'U':
        if(v[0] != 'N')
            break;
        v.remove_prefix(1);
        if(eq(v, "BIND"))
            return verb::unbind;
        if(eq(v, "LINK"))
            return verb::unlink;
        if(eq(v, "LOCK"))
            return verb::unlock;
        if(eq(v, "SUBSCRIBE"))
            return verb::unsubscribe;
        break;

    default:
        break;
    }

    return verb::unknown;
}

} // detail

inline
string_view
to_string(verb v)
{
    return detail::verb_to_string(v);
}

inline
verb
string_to_verb(string_view s)
{
    return detail::string_to_verb(s);
}

} // http
} // beast
} // boost

#endif

max_stars_repo_name: fictheader/fcolorwheel

max_issues_repo_name: fictheader/fcolorwheel

max_forks_repo_name: fictheader/fcolorwheel