text:
/* ************************************************************************
 * Copyright 2015 Vratis, Ltd.
 *
 * 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 <gtest/gtest.h>
#include <vector>
#include <string>
#include <numeric>
#include <boost/program_options.hpp>

//boost ublas
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/blas.hpp>

#include "resources/clsparse_environment.h"
#include "resources/blas1_environment.h"

clsparseControl ClSparseEnvironment::control = NULL;
cl_command_queue ClSparseEnvironment::queue = NULL;
cl_context ClSparseEnvironment::context = NULL;

double Blas1Environment::alpha = 0;
double Blas1Environment::beta = 0;
int Blas1Environment::size = 0;

namespace po = boost::program_options;
namespace uBLAS = boost::numeric::ublas;

template <typename T>
class Blas1 : public ::testing::Test
{

     using CLSE = ClSparseEnvironment;
     using BLAS1E = Blas1Environment;

public:

    void SetUp()
    {
        clsparseInitScalar(&gAlpha);
        clsparseInitScalar(&gBeta);

        clsparseInitVector(&gX);
        clsparseInitVector(&gY);

        hY = uBLAS::vector<T>(BLAS1E::size, 2);
        hX = uBLAS::vector<T>(BLAS1E::size, 4);
        //std::iota(hX.begin(), hX.end(), 1);

        cl_int status;
        gX.values = clCreateBuffer(CLSE::context,
                                   CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                   hX.size() * sizeof(T), hX.data().begin(), &status);
        gX.num_values = hX.size();
        ASSERT_EQ(CL_SUCCESS, status);

        gY.values = clCreateBuffer(CLSE::context,
                                   CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                   hY.size() * sizeof(T), hY.data().begin(), &status);

        gY.num_values = hY.size();
        ASSERT_EQ(CL_SUCCESS, status);

        gAlpha.value = clCreateBuffer(CLSE::context,
                                      CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                      sizeof(T), &hAlpha, &status);

        ASSERT_EQ(CL_SUCCESS, status);

        gBeta.value = clCreateBuffer(CLSE::context,
                                     CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                     sizeof(T), &hBeta, &status);

        ASSERT_EQ(CL_SUCCESS, status);

    }

    void TearDown()
    {
       clReleaseMemObject( gX.values );
       clReleaseMemObject( gY.values );
       clReleaseMemObject( gAlpha.value );
       clReleaseMemObject( gBeta.value );
    }

    boost::numeric::ublas::vector<T> hY;
    boost::numeric::ublas::vector<T> hX;
    T hAlpha = BLAS1E::alpha;
    T hBeta = BLAS1E::beta;

    cldenseVector gX;
    cldenseVector gY;
    clsparseScalar gAlpha;
    clsparseScalar gBeta;

    /*  Test functions  */

    void test_reduce()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSreduce(&gAlpha, &gX, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDreduce(&gAlpha, &gX, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        T ublas_result = uBLAS::blas_1::asum(hX);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gAlpha.value,
                                              CL_TRUE, CL_MAP_READ,
                                              0, sizeof(T), 0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to function passing an assertion here;
        if ( typeid(T) == typeid(cl_float) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-7);
        if ( typeid(T) == typeid(cl_double) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gAlpha.value,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);
    }

    void test_norm1()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSnrm1(&gAlpha, &gX, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDnrm1(&gAlpha, &gX, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        T ublas_result = uBLAS::blas_1::asum(hX);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gAlpha.value,
                                              CL_TRUE, CL_MAP_READ,
                                              0, sizeof(T), 0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to function passing an assertion here;
        if ( typeid(T) == typeid(cl_float) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-7);
        if ( typeid(T) == typeid(cl_double) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gAlpha.value,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_norm2()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSnrm2(&gAlpha, &gX, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDnrm2(&gAlpha, &gX, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        T ublas_result = uBLAS::blas_1::nrm2(hX);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gAlpha.value,
                                              CL_TRUE, CL_MAP_READ,
                                              0, sizeof(T), 0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to function passing an assertion here;
        if ( typeid(T) == typeid(cl_float) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-7);
        if ( typeid(T) == typeid(cl_double) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gAlpha.value,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);
    }

    void test_scale()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSscale(&gX, &gAlpha, &gX, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDscale(&gX, &gAlpha, &gX, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        hX = uBLAS::blas_1::scal(hX, hAlpha);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gX.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gX.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hX.size(); i++)
                ASSERT_NEAR (hX[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hX.size(); i++)
                ASSERT_NEAR (hX[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gX.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_dot()
    {
        clsparseStatus status;
        if ( typeid(T) == typeid(cl_float) )
        {
            status = cldenseSdot(&gAlpha, &gX, &gY, CLSE::control);
        }
        else if( typeid(T)  == typeid (cl_double) )
        {
            status = cldenseDdot(&gAlpha, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        T ublas_result = uBLAS::blas_1::dot(hX, hY);

        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gAlpha.value,
                                              CL_TRUE, CL_MAP_READ,
                                              0, sizeof(T), 0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to function passing an assertion here;
        if ( typeid(T) == typeid(cl_float) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-7);
        if ( typeid(T) == typeid(cl_double) )
            ASSERT_NEAR (ublas_result, *host_result, 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gAlpha.value,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_axpy()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSaxpy(&gY, &gAlpha, &gX, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDaxpy(&gY, &gAlpha, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        hY = uBLAS::blas_1::axpy(hY, hAlpha, hX);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);
    }

    void test_axpby()
    {
        clsparseStatus status;

        //  GPU result;
        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSaxpby(&gY, &gAlpha, &gX, &gBeta, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDaxpby(&gY, &gAlpha, &gX, &gBeta, &gY, CLSE::control);
        }

        ASSERT_EQ (clsparseSuccess, status);

        hY = uBLAS::blas_1::axpy(uBLAS::blas_1::scal(hY, hBeta), hAlpha, hX);

        // Map device data to host
        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_add()
    {
        clsparseStatus status;

        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSadd(&gY, &gX, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDadd(&gY, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ(clsparseSuccess, status);

        hY = hX + hY;

        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);
    }

    void test_sub()
    {
        clsparseStatus status;

        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSsub(&gY, &gX, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDsub(&gY, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ(clsparseSuccess, status);

        hY = hX - hY;

        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_mul()
    {
        clsparseStatus status;

        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSmul(&gY, &gX, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDmul(&gY, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ(clsparseSuccess, status);

        hY = uBLAS::element_prod(hX, hY);

        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

    void test_div()
    {
        clsparseStatus status;

        if (typeid(T) == typeid(cl_float))
        {
            status = cldenseSdiv(&gY, &gX, &gY, CLSE::control);
        }
        else if ( typeid(T) == typeid(cl_double) )
        {
            status = cldenseDdiv(&gY, &gX, &gY, CLSE::control);
        }

        ASSERT_EQ(clsparseSuccess, status);

        hY = uBLAS::element_div(hX, hY);

        cl_int cl_status;
        T* host_result = (T*) ::clEnqueueMapBuffer(CLSE::queue, gY.values,
                                              CL_TRUE, CL_MAP_READ,
                                              0, gY.num_values * sizeof(T),
                                              0, nullptr, nullptr, &cl_status);
        ASSERT_EQ (CL_SUCCESS, cl_status);

        // change it to some compare template functions
        if ( typeid(T) == typeid(cl_float) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-7);

        if ( typeid(T) == typeid(cl_double) )
            for(int i = 0; i < hY.size(); i++)
                ASSERT_NEAR (hY[i], host_result[i], 1e-14);

        cl_status = ::clEnqueueUnmapMemObject(CLSE::queue, gY.values,
                                              host_result, 0, nullptr, nullptr);
        ASSERT_EQ (CL_SUCCESS, cl_status);

    }

};

typedef ::testing::Types<cl_float, cl_double> TYPES;
TYPED_TEST_CASE(Blas1, TYPES);

TYPED_TEST(Blas1, reduce)
{
    this->test_reduce();
}

TYPED_TEST(Blas1, norm1)
{
    this->test_norm1();
}

TYPED_TEST(Blas1, norm2)
{
    this->test_norm2();
}

TYPED_TEST(Blas1, scale)
{
    this->test_scale();
}

TYPED_TEST(Blas1, dot)
{
    this->test_dot();
}

TYPED_TEST(Blas1, axpy)
{
    this->test_axpy();
}

TYPED_TEST(Blas1, axpby)
{
    this->test_axpby();
}

TYPED_TEST(Blas1, add)
{
    this->test_add();
}

TYPED_TEST(Blas1, sub)
{
    this->test_sub();
}

TYPED_TEST(Blas1, mul)
{
    this->test_mul();
}

TYPED_TEST(Blas1, div)
{
    this->test_div();
}

int main (int argc, char* argv[])
{

    using CLSE = ClSparseEnvironment;
    using BLAS1E = Blas1Environment;

    std::string platform;
    cl_platform_type pID;
    cl_uint dID;
    cl_uint size;

    double alpha;
    double beta;

    po::options_description desc("Allowed options");

    desc.add_options()
            ("help,h", "Produce this message.")
            ("platform,l", po::value(&platform)->default_value("AMD"),
             "OpenCL platform: AMD or NVIDIA.")
            ("device,d", po::value(&dID)->default_value(0),
             "Device id within platform.")
            ("size,s",po::value(&size)->default_value(1024),
             "Size of the vectors used for testing")
            ("alpha,a", po::value(&alpha)->default_value(2),
             "Alpha coefficient for blas1 operations i.e. r = alpha * x + beta * y")
            ("beta,b", po::value(&beta)->default_value(4),
             "Beta coefficient for blas1 operations i.e. r = alpha * x + beta * y");

    po::variables_map vm;
    po::parsed_options parsed =
            po::command_line_parser( argc, argv ).options( desc ).allow_unregistered( ).run( );

    try {
        po::store(parsed, vm);
        if (vm.count("help"))
        {
            std::cout << desc << std::endl;
            return 0;
        }
        po::notify(vm);
    }
    catch (po::error& error)
    {
        std::cerr << "Parsing command line options..." << std::endl;
        std::cerr << "Error: " << error.what() << std::endl;
        std::cerr << desc << std::endl;
        return false;
    }

    std::vector< std::string > to_pass_further =
            po::collect_unrecognized( parsed.options, po::include_positional );

    //check platform
    if(vm.count("platform"))
    {
        if ("AMD" == platform)
        {
            pID = AMD;
        }
        else if ("NVIDIA" == platform)
        {
            pID = NVIDIA;
        }
        else
        {

            std::cout << "The platform option is missing or is ill defined!\n";
            std::cout << "Given [" << platform << "]" << std::endl;
            platform = "AMD";
            pID = AMD;
            std::cout << "Setting [" << platform << "] as default" << std::endl;
        }

    }

    ::testing::InitGoogleTest(&argc, argv);
    //order does matter!
    ::testing::AddGlobalTestEnvironment( new CLSE(pID, dID, size));
    ::testing::AddGlobalTestEnvironment( new BLAS1E(alpha, beta, size));

    return RUN_ALL_TESTS();

}

max_stars_repo_name: xfong/clSPARSE

max_issues_repo_name: xfong/clSPARSE

max_forks_repo_name: xfong/clSPARSE


text:
// main.cpp -- example use of CppGPs for Gaussian process regression

//
// Eigen macros for using externel routines (may need to be moved to header file)
// ( see https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html )
//#define EIGEN_USE_BLAS            // Must also link to BLAS library
//#define EIGEN_USE_LAPACKE_STRICT  // Must also link to LAPACKE library
//#define EIGEN_USE_MKL_ALL         // Must also link to Intel MKL library
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <string>
#include <memory>
#include <random>
#include <chrono>
#include <fstream>
#include <boost/range/irange.hpp>
#include "GPs.h"

#include <limits>

// Specify the target function for Gaussian process regression
double targetFunc(Eigen::MatrixXd X)
{

  if ( X.size() == 1 )
    {
      // Define the target function to be an oscillatory, non-periodic function
      double oscillation = 30.0;
      double xshifted = 0.5*(X(0) + 1.0);
      return std::sin(oscillation*(xshifted-0.1))*(0.5-(xshifted-0.1))*15.0;
    }
  else
    {
      // Define the target function to be the Marr wavelet (a.k.a "Mexican Hat" wavelet)
      // ( see https://en.wikipedia.org/wiki/Mexican_hat_wavelet )
      double sigma = 0.25;
      double pi = std::atan(1)*4;
      double radialTerm = std::pow(X.squaredNorm()/sigma,2);
      return 2.0 / (std::sqrt(pi*sigma) * std::pow(pi,1.0/4.0)) * (1.0 - radialTerm) * std::exp(-radialTerm);
    }
}

// Example use of CppGPs code for Gaussian process regression
int main(int argc, char const *argv[])
{

  // Retrieve aliases from GP namescope
  using Matrix = Eigen::MatrixXd;
  using Vector = Eigen::VectorXd;

  // Convenience using-declarations
  using std::cout;
  using std::endl;
  using GP::GaussianProcess;
  using GP::linspace;
  using GP::sampleNormal;
  using GP::sampleUnif;
  using GP::RBF;

  // Used for timing code with chrono
  using GP::high_resolution_clock;
  using GP::time;
  using GP::getTime;

  // Set random seed based on system clock
  std::srand(static_cast<unsigned int>(high_resolution_clock::now().time_since_epoch().count()));

  // Fix random seed for debugging and testing
  //std::srand(static_cast<unsigned int>(0));

  //
  //   [ Define Training Data ]
  //

  // Specify the input dimensions
  //int inputDim = 1;
  int inputDim = 2;

  // Specify observation data count
  int obsCount;
    if ( inputDim == 1 )
      obsCount = 250;
    else
      obsCount = 1000;

  // Specify observation noise level
  auto noiseLevel = 1.0;

  // Define random noise to add to target observations
  auto noise = sampleNormal(obsCount) * noiseLevel;

  // Define observations by sampling random uniform distribution
  Matrix X = sampleUnif(-1.0, 1.0, obsCount, inputDim);
  Matrix y;  y.resize(obsCount, 1);

  // Define target observations 'y' by applying 'targetFunc'
  // to the input observations 'X' and adding a noise vector
  for ( auto i : boost::irange(0,obsCount) )
    y(i) = targetFunc(X.row(i)) + noise(i);

  //
  //   [ Construct Gaussian Process Model ]
  //

  // Initialize Gaussian process model
  GaussianProcess model;

  // Specify training observations for GP model
  model.setObs(X,y);

  // Fix noise level [ noise level is fit to the training data by default ]
  //model.setNoise(0.33);

  // Initialize a RBF covariance kernel and assign it to the model
  RBF kernel;
  model.setKernel(kernel);

  // Specify hyperparameter bounds  [ including noise bounds (index 0) ]
  // Vector lbs(2); lbs << 0.0001 , 0.01;
  // Vector ubs(2); ubs << 5.0    , 100.0;
  // model.setBounds(lbs, ubs);

  // Specify hyperparameter bounds  [ excluding noise bounds ]
  Vector lbs(1);  lbs <<  0.01;
  Vector ubs(1);  ubs <<  100.0;
  model.setBounds(lbs, ubs);

  // Fit covariance kernel hyperparameters to the training data
  time start = high_resolution_clock::now();
  model.fitModel();
  time end = high_resolution_clock::now();
  auto computationTime = getTime(start, end);

  // Display computation time required for fitting the GP model
  cout << "\nComputation Time: ";
  cout << computationTime << " s" << endl;

  // Retrieve the tuned/optimized kernel hyperparameters
  auto optParams = model.getParams();
  cout << "\nOptimized Hyperparameters:" << endl << optParams.transpose() << "  ";
  auto noiseL = model.getNoise();
  cout << "(Noise = " << noiseL << ")\n" << endl;

  // Display the negative log marginal likelihood (NLML) of the optimized model
  cout << "NLML:  " << std::fixed << std::setprecision(4) << model.computeNLML() << endl << endl;

  //double eps = std::numeric_limits<double>::min();
  //cout << "Epsilon: " << eps << endl;

  //
  //   [ Posterior Predictions and Samples ]
  //

  // Define test mesh for GP model predictions
  int predCount;
  if ( inputDim == 1 )
    predCount = 100;
  else
    predCount = 1000;
  //int predRes = static_cast<int>(std::sqrt(predCount));
  int predRes = static_cast<int>(std::pow(predCount,1.0/inputDim));
  //auto testMesh = linspace(0.0, 1.0, predCount);
  auto testMesh = linspace(-1.0, 1.0, predRes, inputDim);
  model.setPred(testMesh);

  // Reset predCount to account for rounding
  predCount = std::pow(predRes,inputDim);

  // Compute predicted means and variances for the test points
  model.predict();
  Matrix pmean = model.getPredMean();
  Matrix pvar = model.getPredVar();
  Matrix pstd = pvar.array().sqrt().matrix();

  // Get sample paths from the posterior distribution of the model
  int sampleCount;
  Matrix samples;
  if ( inputDim == 1 )
    {
      sampleCount = 25;
      samples = model.getSamples(sampleCount);
    }

  //
  //   [ Save Results for Plotting ]
  //

  // Save true and predicted means/variances to file
  std::string outputFile = "predictions.csv";
  //Matrix trueSoln = testMesh.unaryExpr(std::ptr_fun(targetFunc));

  Matrix trueSoln(predCount,1);
  for ( auto i : boost::irange(0,predCount) )
    trueSoln(i,0) = targetFunc(testMesh.row(i));

  std::ofstream fout;
  fout.open(outputFile);
  for ( auto i : boost::irange(0,predCount) )
    {
      for ( auto j : boost::irange(0,inputDim) )
        fout << testMesh(i,j) << ",";

      fout << trueSoln(i) << "," << pmean(i) << "," << pstd(i) << "\n";
    }
  fout.close();

  // Save observations to file
  std::string outputObsFile = "observations.csv";
  fout.open(outputObsFile);
  for ( auto i : boost::irange(0,obsCount) )
    {
      for ( auto j : boost::irange(0,inputDim) )
        fout << X(i,j) << ",";

      fout << y(i) << "\n";
    }
  fout.close();

  if ( inputDim == 1 )
    {
      // Save samples to file
      std::string outputSampleFile = "samples.csv";
      fout.open(outputSampleFile);
      for ( auto j : boost::irange(0,sampleCount) )
        {
          for ( auto i : boost::irange(0,predCount) )
            fout << samples(i,j) << ((i<predCount-1) ? "," : "\n");
        }
      fout.close();
    }

  return 0;

}

ext: cpp


text:
// Copyright 2008 by BBN Technologies Corp.
// All Rights Reserved.

#include "Generic/common/leak_detection.h"

#include "Generic/common/limits.h"
#include "Generic/common/Sexp.h"
#include "Generic/common/ParamReader.h"
#include "Generic/common/UTF8InputStream.h"
#include "Generic/common/InternalInconsistencyException.h"
#include "Generic/theories/MentionSet.h"
#include "Generic/theories/Mention.h"
#include "Generic/theories/SynNode.h"
#include "Generic/descriptors/SubtypeClassifier.h"
#include "Generic/common/SymbolUtilities.h"
#include "Generic/theories/Proposition.h"
#include "Generic/relations/xx_RelationUtilities.h"

#include "Generic/parse/LanguageSpecificFunctions.h"

#include "Generic/theories/PartOfSpeechSequence.h"
#include <boost/scoped_ptr.hpp>

#ifdef _WIN32
	#define swprintf _snwprintf
#endif

SubtypeClassifier::SubtypeClassifier()
	: _wordnetMap(0), _descHeadwordMap(0), _fullNameMap(0),
	_nameWordMap(0), _wikiWordMap(0), _posSequence(0)
{
	if (!EntitySubtype::subtypesDefined())
		return;

	_wordnetMap = _new WordToSubtypeMap(5000);

	std::string filename = ParamReader::getParam("wordnet_subtypes");
	if (!filename.empty()) {
		boost::scoped_ptr<UTF8InputStream> stream_scoped_ptr(UTF8InputStream::build(filename.c_str()));
		UTF8InputStream& stream(*stream_scoped_ptr);
		Sexp *wordMappings = _new Sexp(stream);
		int nmappings = wordMappings->getNumChildren();
		for (int i = 0; i < nmappings; i++) {
			Sexp *mapping = wordMappings->getNthChild(i);
			Symbol word = mapping->getFirstChild()->getValue();
			Symbol subtype = mapping->getSecondChild()->getValue();
			EntitySubtype stype(subtype);

			EntitySubtypeItem *newItem = _new EntitySubtypeItem();
			newItem->entitySubtype = stype;
			newItem->next = NULL;

			if (_wordnetMap->get(word) != NULL)
				newItem->next = *_wordnetMap->get(word);

			(*_wordnetMap)[word] = newItem;
		}
		delete wordMappings;
	}

	_descHeadwordMap = _new WordToSubtypeMap(5000);
	filename = ParamReader::getParam("desc_head_subtypes");
	if (!filename.empty()) {
		boost::scoped_ptr<UTF8InputStream> stream_scoped_ptr(UTF8InputStream::build(filename.c_str()));
		UTF8InputStream& stream(*stream_scoped_ptr);
		Sexp *wordMappings = _new Sexp(stream);
		int nmappings = wordMappings->getNumChildren();
		for (int i = 0; i < nmappings; i++) {
			Sexp *mapping = wordMappings->getNthChild(i);
			int nkids = mapping->getNumChildren();
			Symbol subtype = mapping->getNthChild(nkids-1)->getValue();

			Symbol word = mapping->getFirstChild()->getValue();
			//Symbol subtype = mapping->getSecondChild()->getValue();
			EntitySubtype stype(subtype);

			EntitySubtypeItem *newItem = _new EntitySubtypeItem();
			newItem->entitySubtype = stype;
			newItem->next = NULL;

			if (_descHeadwordMap->get(word) != NULL)
				newItem->next = *_descHeadwordMap->get(word);

			(*_descHeadwordMap)[word] = newItem;
		}
		delete wordMappings;
	}

	_nameWordMap = _new WordToSubtypeMap(5000);
	filename = ParamReader::getParam("name_word_subtypes");
	if (!filename.empty()) {
		boost::scoped_ptr<UTF8InputStream> stream_scoped_ptr(UTF8InputStream::build(filename.c_str()));
		UTF8InputStream& stream(*stream_scoped_ptr);
		Sexp *wordMappings = _new Sexp(stream);
		int nmappings = wordMappings->getNumChildren();
		for (int i = 0; i < nmappings; i++) {
			Sexp *mapping = wordMappings->getNthChild(i);
			Symbol word = mapping->getFirstChild()->getValue();
			Symbol subtype = mapping->getSecondChild()->getValue();
			EntitySubtype stype(subtype);

			EntitySubtypeItem *newItem = _new EntitySubtypeItem();
			newItem->entitySubtype = stype;
			newItem->next = NULL;

			if (_nameWordMap->get(word) != NULL)
				newItem->next = *_nameWordMap->get(word);

			(*_nameWordMap)[word] = newItem;
		}
		delete wordMappings;
	}

	_fullNameMap = _new WordToSubtypeMap(5000);
	filename = ParamReader::getParam("full_name_subtypes");
	if (!filename.empty()) {
		boost::scoped_ptr<UTF8InputStream> stream_scoped_ptr(UTF8InputStream::build(filename.c_str()));
		UTF8InputStream& stream(*stream_scoped_ptr);
		Sexp *wordMappings = _new Sexp(stream);
		int nmappings = wordMappings->getNumChildren();
		for (int i = 0; i < nmappings; i++) {
			Sexp *mapping = wordMappings->getNthChild(i);
			int nkids = mapping->getNumChildren();
			Symbol subtype = mapping->getNthChild(nkids-1)->getValue();

			std::wostringstream nameStr;
			for (int j = 0; j < nkids - 1; j++) {
				Symbol word = mapping->getNthChild(j)->getValue();
				nameStr << word.to_string();
				if (j < nkids - 2) nameStr << L" ";
			}
			Symbol nameSym(nameStr.str().c_str());
			EntitySubtype stype(subtype);

			EntitySubtypeItem *newItem = _new EntitySubtypeItem();
			newItem->entitySubtype = stype;
			newItem->next = NULL;

			if (_fullNameMap->get(nameSym) != NULL)
				newItem->next = *_fullNameMap->get(nameSym);

			(*_fullNameMap)[nameSym] = newItem;
		}
		delete wordMappings;
	}
	//Wikipedia
	_wikiWordMap = _new WordToSubtypeMap(10000);
	filename = ParamReader::getParam("wikipedia_subtypes");
	if (!filename.empty()) {
		boost::scoped_ptr<UTF8InputStream> stream_scoped_ptr(UTF8InputStream::build(filename.c_str()));
		UTF8InputStream& stream(*stream_scoped_ptr);
		Sexp *wordMappings = _new Sexp(stream);
		int nmappings = wordMappings->getNumChildren();
		for (int i = 0; i < nmappings; i++) {
			Sexp *mapping = wordMappings->getNthChild(i);
			int nkids = mapping->getNumChildren();
			Symbol subtype = mapping->getNthChild(nkids-1)->getValue();

			std::wostringstream nameStr;
			for (int j = 0; j < nkids - 1; j++) {
				Symbol word = mapping->getNthChild(j)->getValue();
				nameStr << word.to_string();
				if (j < nkids - 2) nameStr << L" ";
			}
			Symbol nameSym(nameStr.str().c_str());
			EntitySubtype stype(subtype);

			EntitySubtypeItem *newItem = _new EntitySubtypeItem();
			newItem->entitySubtype = stype;
			newItem->next = NULL;

			if (_wikiWordMap->get(nameSym) != NULL)
				newItem->next = *_wikiWordMap->get(nameSym);

			(*_wikiWordMap)[nameSym] = newItem;
		}
		delete wordMappings;
	}
}

SubtypeClassifier::~SubtypeClassifier() {
	deleteItems(_fullNameMap);
	delete _fullNameMap;
	deleteItems(_nameWordMap);
	delete _nameWordMap;
	deleteItems(_descHeadwordMap);
	delete _descHeadwordMap;
	deleteItems(_wordnetMap);
	delete _wordnetMap;
	deleteItems(_wikiWordMap);
	delete _wikiWordMap;
}

void SubtypeClassifier::deleteItems(WordToSubtypeMap *map)
{
	WordToSubtypeMap::iterator iter;
	for (iter = map->begin(); iter != map->end(); ++iter) {
		EntitySubtypeItem *item = (*iter).second;

		while (item != 0) {
			EntitySubtypeItem *next = item->next;
			delete item;
			item = next;
		}
	}
}

int SubtypeClassifier::classifyMention (MentionSet *currSolution, Mention *currMention,
										MentionSet *results[],
										int max_results, bool isBranching)
{
	// this is deterministic anyway, but whether we fork or not depends on if we're isBranching
	if (isBranching) {
		results[0] = _new MentionSet(*currSolution);
		currMention = results[0]->getMention(currMention->getIndex());
	}
	else
		results[0] = currSolution;

	// no subtypes for OTH, obviously
	if (!currMention->getEntityType().isRecognized())
		return 1;

	// don't reclassify if a subtype is already present
	if (currMention->getEntitySubtype().isDetermined())
		return 1;

	switch (currMention->mentionType) {
		case Mention::NONE:
			break;
		case Mention::DESC:
			classifyDesc(currMention);
			break;
		case Mention::NAME:
		case Mention::NEST:
			classifyName(currMention);
			break;
		case Mention::PART:
			classifyPart(currMention);
			break;
		case Mention::APPO:
			classifyAppo(currMention);
			break;
		case Mention::LIST:
			classifyList(currMention);
			break;
		default:
			break;
	}
	return 1;
}

EntitySubtype* SubtypeClassifier::lookupSubtype(WordToSubtypeMap *map, Symbol sym, EntityType type)
{

	EntitySubtypeItem **subtypeItem_p = map->get(sym);
	if (subtypeItem_p == 0) return 0;
	EntitySubtypeItem *subtypeItem = *subtypeItem_p;

	while (subtypeItem != 0 &&
		   subtypeItem->entitySubtype.getParentEntityType() != type)
	{
		subtypeItem = subtypeItem->next;
	}

	if (subtypeItem == 0) {
		return 0;
	}
	else {
		return &(subtypeItem->entitySubtype);
	}
}

void SubtypeClassifier::classifyName(Mention *currMention) {
	Mention *tempMention = currMention;
	while (tempMention->getChild() != 0)
		tempMention = tempMention->getChild();

	const SynNode *nameNode = tempMention->getNode();
	Symbol nameSym = SymbolUtilities::getFullNameSymbol(nameNode);

	// can we classify the full name
	EntitySubtype *subtype = lookupSubtype(_fullNameMap, nameSym, currMention->getEntityType());

	if (subtype != 0) {
		currMention->setEntitySubtype(*subtype);
		//std::cout << "\n" << subtype->getName().to_debug_string() << " (full name): ";
		//std::cout << nameSym.to_debug_string() << "\n";
		return;
	}

	//try stemming the full name
	Symbol stemmedname =  SymbolUtilities::getStemmedFullName(nameNode);
	//std::cout<<"Stem name: "
	//	<<nameSym.to_debug_string()<<" -> "<<stemmedname.to_debug_string()<<std::endl;

	subtype = lookupSubtype(_fullNameMap, stemmedname, currMention->getEntityType());

	if (subtype != 0)
	{
		currMention->setEntitySubtype(*subtype);
		//std::cout << "\n" << subtype->getName().to_debug_string() << " (stemmed name): ";
		//std::cout << nameSym.to_debug_string() << "\n";
		return;
	}

	// is there a word in the name that might give us a clue
	for (int i = 0; i < nameNode->getNChildren(); i++) {
		Symbol word = nameNode->getChild(i)->getHeadWord();
		EntitySubtype *subtype = lookupSubtype(_nameWordMap, word, currMention->getEntityType());

		if (subtype != 0)
		{
			currMention->setEntitySubtype(*subtype);
			//std::cout << "\n" << subtype->getName().to_debug_string() << " (name word): ";
			//std::cout << nameSym.to_debug_string() << "\n";
			return;
		}
		//Try looking on wiki lists too
		subtype = lookupSubtype(_wikiWordMap, word, currMention->getEntityType());
		if (subtype != 0)
		{
			currMention->setEntitySubtype(*subtype);
			//std::cout << "\n" << subtype->getName().to_debug_string() << " (wiki word): ";
			//std::cout << nameSym.to_debug_string() << "\n";
			return;
		}
	}
	//try stemmed variants of word

	for (int k = 0; k < nameNode->getNChildren(); k++) {
		Symbol word = nameNode->getChild(k)->getHeadWord();
		Symbol variants[15];
		int var = SymbolUtilities::getStemVariants(word, variants, 15);
		for(int j = 0; j< var; j++){
			subtype = lookupSubtype(_descHeadwordMap, variants[j], currMention->getEntityType());
			if(subtype != 0)
			{
				currMention->setEntitySubtype(*subtype);
				return;
			}
		}
	}

}

void SubtypeClassifier::classifyDesc(Mention *currMention) {
	Symbol headword = currMention->getNode()->getHeadWord();
	// do we know the actual headword
	EntitySubtype *subtype = lookupSubtype(_descHeadwordMap, headword, currMention->getEntityType());
	if (subtype == 0)
		subtype =
			lookupSubtype(_descHeadwordMap,
			              RelationUtilities::get()->stemPredicate(headword, Proposition::NOUN_PRED),
						  currMention->getEntityType());
	if (subtype == 0){
		Symbol variants[15];
		int var = SymbolUtilities::getStemVariants(headword, variants, 15);
		for (int i = 0; i< var; i++){
			subtype = lookupSubtype(_descHeadwordMap, variants[i], currMention->getEntityType());
			if (subtype != 0)
			{
				currMention->setEntitySubtype(*subtype);
				return;
			}
		}
	}

	if (subtype != 0)
	{
		currMention->setEntitySubtype(*subtype);
		return;
	}

	//Warning this is Language and ACE2005 Specific
	if(currMention->getEntityType().matchesPER()){
		if(SymbolUtilities::isPluralMention(currMention)){
			try{
				EntitySubtype st =  EntitySubtype(Symbol(L"PER"), Symbol(L"Group"));
				currMention->setEntitySubtype(st);
			}
			catch(UnexpectedInputException){
				//probably not ACE2005, just ignore this
				;
			}
		}

		if(_posSequence != 0){
			int token = currMention->getNode()->getHeadPreterm()->getStartToken();
			if(token < _posSequence->getNTokens()){
				PartOfSpeech* tags = _posSequence->getPOS(token);
				if(LanguageSpecificFunctions::POSTheoryContainsNounPL(tags) &&
					!LanguageSpecificFunctions::POSTheoryContainsNounSG(tags)  &&
					!LanguageSpecificFunctions::POSTheoryContainsNounAmbiguousNumber(tags))
				{
					try{
						EntitySubtype st =  EntitySubtype(Symbol(L"PER"), Symbol(L"Group"));
						currMention->setEntitySubtype(st);
					}
					catch(UnexpectedInputException){
						//probably not ACE2005, just ignore this
						;
					}
				}
			}
		}

	}

	// speed test for Nominal Subtype Classificiation
	// trial-1 -- turned off WordNet

	Symbol results[100];
	int nhypernyms = SymbolUtilities::getHypernyms(headword, results, 100);

	// might we know a hypernym of the headword
	for (int i = 0; i < nhypernyms; i++) {
		//std::cout << "  " << results[i].to_debug_string() << "\n";
		EntitySubtype *subtype = lookupSubtype(_wordnetMap, results[i], currMention->getEntityType());
		if (subtype != 0)
		{
			currMention->setEntitySubtype(*subtype);
			//std::cout << "\n" << subtype->getName().to_debug_string() << " (wordnet): ";
			//std::cout << currMention->getNode()->toDebugTextString() << "\n";
			return;
		}
	}

	//look on wiki list
	subtype = lookupSubtype(_wikiWordMap, headword, currMention->getEntityType());
	if (subtype == 0)
		subtype = lookupSubtype(_wikiWordMap, RelationUtilities::get()->stemPredicate(headword, Proposition::NOUN_PRED), currMention->getEntityType());
	if(subtype == 0){
		Symbol variants[15];
		int var = SymbolUtilities::getStemVariants(headword, variants, 15);
		for(int i = 0; i< var; i++){
			subtype = lookupSubtype(_wikiWordMap, variants[i], currMention->getEntityType());
			if(subtype != 0)
			{
				currMention->setEntitySubtype(*subtype);
				return;
			}
		}
	}
	if (subtype != 0 && subtype->getParentEntityType() == currMention->getEntityType())
	{
		currMention->setEntitySubtype(*subtype);
		return;
	}

}

void  SubtypeClassifier::classifyList(Mention *currMention) {
	// ?? can't bring myself to care
}
void  SubtypeClassifier::classifyAppo(Mention *currMention) {
	if (currMention->getChild() != 0) {
		if (EntitySubtype::isValidEntitySubtype(currMention->getEntityType(),currMention->getChild()->getEntitySubtype())) {
			currMention->setEntitySubtype(currMention->getChild()->getEntitySubtype());
		}
	}
}
void  SubtypeClassifier::classifyPart(Mention *currMention) {
	if (currMention->getChild() != 0) {
		if (EntitySubtype::isValidEntitySubtype(currMention->getEntityType(),currMention->getChild()->getEntitySubtype())) {
			currMention->setEntitySubtype(currMention->getChild()->getEntitySubtype());
		}
	}
}

max_stars_repo_name: BBN-E/serif

max_issues_repo_name: BBN-E/serif

max_forks_repo_name: BBN-E/serif


text:
/// @file BaselineMap.cc
///
/// @copyright (c) 2012 CSIRO
/// Australia Telescope National Facility (ATNF)
/// Commonwealth Scientific and Industrial Research Organisation (CSIRO)
/// PO Box 76, Epping NSW 1710, Australia
/// atnf-enquiries@csiro.au
///
/// This file is part of the ASKAP software distribution.
///
/// The ASKAP software distribution is free software: you can redistribute it
/// and/or modify it under the terms of the GNU General Public License as
/// published by the Free Software Foundation; either version 2 of the License,
/// or (at your option) any later version.
///
/// This program is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with this program; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
///
/// @author Max Voronkov <maxim.voronkov@csiro.au>
/// Original code by Ben Humphreys

// Include own header file first
#include "BaselineMap.h"

// Include package level header file
#include "askap_cpingest.h"

// System includes
#include <map>
#include <stdint.h>

// boost include
#include <boost/tuple/tuple_comparison.hpp>

// ASKAPsoft includes
#include "askap/AskapError.h"
#include "askap/AskapUtil.h"
#include "Common/ParameterSet.h"
#include "casacore/measures/Measures/Stokes.h"

// Using
using namespace askap;
using namespace askap::utility;
using namespace askap::cp;
using namespace askap::cp::ingest;
using namespace casa;

BaselineMap::BaselineMap()
{
}

BaselineMap::BaselineMap(const LOFAR::ParameterSet& parset) : itsUpperTriangle(true), itsLowerTriangle(true)
{
    const std::string configName = parset.getString("name", "");
    if (configName == "standard") {
        ASKAPCHECK(!parset.isDefined("baselineids"),
              "The baseline map has to be either defined explicitly via baselineids keyword or via name, you have both");
        defaultMapADE();
    } else {
        ASKAPCHECK(configName == "", "Only 'standard' implicit baseline map is currently supported; you have="<<configName);

        const vector<int32_t> ids = parset.getInt32Vector("baselineids", true);
        itsSize = ids.size();

        for (vector<int32_t>::const_iterator it = ids.begin();
                it != ids.end(); ++it) {
            const int32_t id = *it;
            if (!parset.isDefined(toString(id))) {
                ASKAPTHROW(AskapError, "Baseline mapping for id " << id << " not present");
            }

            const vector<string> tuple = parset.getStringVector(toString(id));
            if (tuple.size() != 3) {
                ASKAPTHROW(AskapError, "Baseline mapping for id " << id << " is malformed");
            }

            const int32_t ant1 = fromString<int32_t>(tuple[0]);
            const int32_t ant2 = fromString<int32_t>(tuple[1]);

            add(id, ant1, ant2, Stokes::type(tuple[2]));
        }
    }
    ASKAPCHECK(itsMap.size() == itsSize, "Failed to initialise baseline map");
}

/// @brief populate map for ADE correlator
/// @details To avoid carrying the map for 2628 products explicitly in fcm,
/// we use this method to define the full map analytically. In the future,
/// we might even have a polymorphic class which does mapping analytically.
/// This could even speed things up. However, at this stage, an option to
/// support sparse arrays is more useful, so we keep the full map functionality
/// in and just generate the map algorithmically.
/// @param[in] nAnt number of antennas to generate the map for
void BaselineMap::defaultMapADE(uint32_t nAnt)
{
   ASKAPASSERT(itsMap.size() == 0);
   ASKAPDEBUGASSERT(nAnt > 0);
   // we have 1-based product indices, so use prefix increment
   itsSize = 0;
   for (int32_t ant2 = 0; ant2 < static_cast<int32_t>(nAnt); ++ant2) {
        for (int32_t ant1 = 0; ant1 < ant2; ++ant1) {
             add(++itsSize, ant1, ant2, casa::Stokes::XX);
             add(++itsSize, ant1, ant2, casa::Stokes::YX);
        }
        add(++itsSize, ant2, ant2, casa::Stokes::XX);
        for (int32_t ant1 = 0; ant1 <= ant2; ++ant1) {
             add(++itsSize, ant1, ant2, casa::Stokes::XY);
             add(++itsSize, ant1, ant2, casa::Stokes::YY);
        }
   }
   // itsSize should be set to the number of products here
}

/// @brief add one product to the map
/// @param[in] id product number to add the mapping for
/// @param[in] ant1 first antenna index
/// @param[in] ant2 second antenna index
/// @param[in] pol stokes parameter
void BaselineMap::add(int32_t id, int32_t ant1, int32_t ant2, casa::Stokes::StokesTypes pol)
{
   ASKAPCHECK(ant1 >= 0, "Antenna 1 index is supposed to be non-negative");
   ASKAPCHECK(ant2 >= 0, "Antenna 2 index is supposed to be non-negative");
   ASKAPCHECK(id >= 0, "Product id is supposed to be non-negative");

   if (ant1 > ant2) {
       itsUpperTriangle = false;
   }

   if (ant2 > ant1) {
       itsLowerTriangle = false;
   }

   itsMap[id] = boost::make_tuple(ant1,ant2,pol);
}

/// @brief caching method for itsCachedProduct
/// @detail Calling this method sets itsCachedProduct
/// @param[in] id   the product (baseline) id
void BaselineMap::syncProductCache(int32_t id) const
{
    if (itsCachedProduct && ((*itsCachedProduct)->first == id)) {
        return;
    }

    itsCachedProduct = itsMap.find(id);

    if ((*itsCachedProduct) == itsMap.end()) {
        itsCachedProduct = boost::none;
    }
}

int32_t BaselineMap::idToAntenna1(const int32_t id) const
{
    syncProductCache(id);
    return itsCachedProduct ? (*itsCachedProduct)->second.get<0>() : -1;
}

int32_t BaselineMap::idToAntenna2(const int32_t id) const
{
    syncProductCache(id);
    return itsCachedProduct ? (*itsCachedProduct)->second.get<1>() : -1;
}

casa::Stokes::StokesTypes BaselineMap::idToStokes(const int32_t id) const
{
    syncProductCache(id);
    return itsCachedProduct ? (*itsCachedProduct)->second.get<2>() : Stokes::Undefined;
}

size_t BaselineMap::size() const
{
    return itsSize;
}

/// @brief obtain largest id
/// @details This is required to initialise a flat array buffer holding
/// derived per-id information because the current implementation does not
/// explicitly prohibits sparse ids.
/// @return the largest id setup in the map
int32_t BaselineMap::maxID() const
{
  int32_t result = 0;
  for (std::map<int32_t, ProductDesc>::const_iterator ci = itsMap.begin();
       ci != itsMap.end(); ++ci) {
       ASKAPCHECK(ci->first >= 0, "Encountered negative id="<<ci->first);
       result = max(result, ci->first);
  }
  return result;
}

/// @brief find an id matching baseline/polarisation description
/// @details This is the reverse look-up operation.
/// @param[in] ant1 index of the first antenna
/// @param[in] ant2 index of the second antenna
/// @param[in] pol polarisation product
/// @return the index of the selected baseline/polarisation
/// @note an exception is thrown if there is no match
int32_t BaselineMap::getID(const int32_t ant1, const int32_t ant2, const casa::Stokes::StokesTypes pol) const
{
  const ProductDesc product(ant1, ant2, pol);
  for (std::map<int32_t, ProductDesc>::const_iterator ci = itsMap.begin(); ci!=itsMap.end(); ++ci) {
       if (ci->second == product) {
           return ci->first;
       }
  }
  return -1;
}

/// @brief correlator produces lower triangle?
/// @return true if ant2<=ant1 for all ids
bool BaselineMap::isLowerTriangle() const
{
   return itsLowerTriangle && (itsSize != 0);
}

/// @brief correlator produces upper triangle?
/// @return true if ant1<=ant2 for all ids
bool BaselineMap::isUpperTriangle() const
{
   return itsUpperTriangle && (itsSize != 0);
}

/// @brief take slice for a subset of indices
/// @details This method is probably temporary as it is primarily intended for ADE
/// commissioning. When we have a decent number of ASKAP antennas ready, this
/// additional layer of mapping needs to be removed as it is a complication.
/// This method produces a sparse map which includes only selected antenna indices.
/// @param[in] ids vector with antenna ids to keep
/// @note For simplicity, indices should always be given in increasing order. This
/// ensures that no need for data conjugation arises at the user side (i.e. upper
/// and lower triangles will remain such).
void BaselineMap::sliceMap(const std::vector<int32_t> &ids)
{
   // sanity check on supplied indices
   int32_t largestAntID = -1;

   for (std::map<int32_t, ProductDesc>::const_iterator ci = itsMap.begin();
        ci != itsMap.end(); ++ci) {

        if (largestAntID < ci->second.get<0>()) {
            largestAntID = ci->second.get<0>();
        }
        if (largestAntID < ci->second.get<1>()) {
            largestAntID = ci->second.get<1>();
        }

   }
   ASKAPCHECK(largestAntID >= 0, "Attempting to slice an empty map");

   int32_t previousID = -1;
   for (std::vector<int32_t>::const_iterator ci = ids.begin(); ci != ids.end(); ++ci) {
        ASKAPCHECK(*ci <= largestAntID, "sliceMap received antenna ID="<<*ci<<
                   " which exceeds the largest antenna index encountered in the map = "<<largestAntID);
        ASKAPCHECK(*ci >= 0, "Negative antenna ID encountered in sliceMap call");
        ASKAPCHECK(*ci > previousID,
            "Antenna indices passed to sliceMap are expected to be in increasing order; you have = "<<ids);
        previousID = *ci;
   }

   // taking the slice
   std::map<int32_t, ProductDesc> newMap;

   size_t newSize = 0;

   for (std::map<int32_t, ProductDesc>::const_iterator ci = itsMap.begin(); ci != itsMap.end(); ++ci) {

        //ASKAPDEBUGASSERT(ciAnt2 != itsAntenna2Map.end());
        //ASKAPDEBUGASSERT(ciPol != itsStokesMap.end());

        // indices should match
        const int32_t productID = ci->first;

        ProductDesc newProduct = ci->second;

        int32_t newIndex1 = -1;
        int32_t newIndex2 = -1;
        for (size_t i = 0; i < ids.size(); ++i) {
             if (ids[i] == newProduct.get<0>()) {
                 ASKAPDEBUGASSERT(newIndex1 < 0);
                 newIndex1 = static_cast<int32_t>(i);
             }
             if (ids[i] == newProduct.get<1>()) {
                 ASKAPDEBUGASSERT(newIndex2 < 0);
                 newIndex2 = static_cast<int32_t>(i);
             }
        }

        if ((newIndex1 >= 0) && (newIndex2 >= 0)) {
             newProduct.get<0>() = newIndex1;
             newProduct.get<1>() = newIndex2;
             newMap[productID] = newProduct;
             ++newSize;
        }
   }

   ASKAPCHECK(newSize > 0, "Taking slice rejected all correlation products in the map from "<<
                           itsSize<<" products available in the map");

   // assign new maps to this class' data members
   itsMap = newMap;
   itsSize = newSize;
}

ext: cc

max_stars_repo_name: rtobar/askapsoft

max_issues_repo_name: ATNF/askapsoft

max_forks_repo_name: ATNF/askapsoft


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

#include <iostream>

#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/transposed_view.hpp>
#include <boost/numeric/mtl/matrix/parameter.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/range_generator.hpp>
#include <boost/numeric/mtl/utility/glas_tag.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/maybe.hpp>
#include <boost/numeric/mtl/utility/complexity.hpp>
#include <boost/numeric/mtl/operation/print_matrix.hpp>
#include <boost/numeric/mtl/io/test_ostream.hpp>

using namespace std;

template <typename Matrix>
void symmetry_test(const Matrix& A)
{
    namespace traits = mtl::traits;

    typename traits::row<Matrix>::type                                 row(A);
    typename traits::col<Matrix>::type                                 col(A);
    typename traits::const_value<Matrix>::type                         value(A);
    typedef typename traits::range_generator<mtl::tag::major, Matrix>::type        cursor_type;
    for (cursor_type cursor = mtl::begin<mtl::tag::major>(A), cend = mtl::end<mtl::tag::major>(A);
	 cursor != cend; ++cursor) {
	typedef mtl::tag::nz     inner_tag;
	mtl::io::tout << "---\n";
	typedef typename traits::range_generator<inner_tag, cursor_type>::type icursor_type;
	for (icursor_type icursor = mtl::begin<inner_tag>(cursor), icend = mtl::end<inner_tag>(cursor); icursor != icend; ++icursor) {
	    mtl::io::tout << "A[" << row(*icursor) << ", " << col(*icursor) << "] = " << value(*icursor) << '\n';
	    MTL_THROW_IF(!A.indexer(A, col(*icursor), row(*icursor)), mtl::runtime_error("Symmetric counter-part not found"));
	}
    }
}

template <typename Matrix>
void matrix_init(Matrix& matrix)
{
    typedef typename Matrix::parameters   parameters;
    typedef typename Matrix::value_type   value_type;

    mtl::mat::compressed2D_inserter<value_type, parameters> inserter(matrix);
    inserter(1, 1) << 0;
    inserter(2, 2) << 7; inserter(1, 4) << 3; inserter(3, 2) << 9; inserter(5, 1) << 5;
}

template <typename Orientation>
void test_compressed2D(char const* name)
{
    mtl::io::tout << "\n====================\n" << name << "\n====================\n";
    typedef mtl::mat::parameters<Orientation>   parameters;
    typedef mtl::compressed2D<int, parameters>     matrix_type;
    matrix_type                                    A(6, 6);

    matrix_init(A);
    mtl::io::tout << "A is\n" << A;
    mtl::io::tout << "Has " << A.nnz() << " non-zeros.\n";

    A.make_symmetric_pattern();
    mtl::io::tout << "A after symmetrizing is\n" << A;
    mtl::io::tout << "Has " << A.nnz() << " non-zeros.\n";
    symmetry_test(A);
}

int main(int, char**)
{
    test_compressed2D<mtl::row_major>("CRS");
    test_compressed2D<mtl::col_major>("CCS");

    return 0;
}


text:
/*
  Copyright 2009 Larry Gritz and the other authors and contributors.
  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 software's owners 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.

  (This is the Modified BSD License)
*/

#include <boost/scoped_array.hpp>

#include "py_oiio.h"
#include "OpenImageIO/platform.h"

namespace PyOpenImageIO
{
using namespace boost::python;

std::string
ImageBuf_name (const ImageBuf &buf)
{
    return buf.name();
}

std::string
ImageBuf_file_format_name (const ImageBuf &buf)
{
    return buf.file_format_name();
}

void
ImageBuf_reset_name (ImageBuf &buf, const std::string &name)
{
    buf.reset (name);
}

void
ImageBuf_reset_name2 (ImageBuf &buf, const std::string &name,
                      int subimage, int miplevel)
{
    buf.reset (name, subimage, miplevel);
}

void
ImageBuf_reset_name_config (ImageBuf &buf, const std::string &name,
                      int subimage, int miplevel, const ImageSpec &config)
{
    buf.reset (name, subimage, miplevel, NULL, &config);
}

void
ImageBuf_reset_spec (ImageBuf &buf, const ImageSpec &spec)
{
    buf.reset (spec);
}

bool
ImageBuf_read (ImageBuf &buf, int subimage=0, int miplevel=0,
               bool force=false, TypeDesc convert=TypeDesc::UNKNOWN)
{
    ScopedGILRelease gil;
    return buf.read (subimage, miplevel, force, convert);
}

bool
ImageBuf_read2 (ImageBuf &buf, int subimage=0, int miplevel=0,
                bool force=false,
                TypeDesc::BASETYPE convert=TypeDesc::UNKNOWN)
{
    ScopedGILRelease gil;
    return buf.read (subimage, miplevel, force, convert);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_read_overloads,
                                ImageBuf_read, 1, 5)
BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_read2_overloads,
                                ImageBuf_read2, 1, 5)

bool
ImageBuf_write (const ImageBuf &buf, const std::string &filename,
                const std::string &fileformat="")
{
    ScopedGILRelease gil;
    return buf.write (filename, fileformat);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_write_overloads,
                                ImageBuf_write, 2, 3)

bool
ImageBuf_make_writeable (ImageBuf &buf, bool keep_cache_type)
{
    ScopedGILRelease gil;
    return buf.make_writeable (keep_cache_type);
}

void
ImageBuf_set_write_format (ImageBuf &buf, TypeDesc::BASETYPE format)
{
    buf.set_write_format (format);
}

bool
ImageBuf_copy (ImageBuf &buf, const ImageBuf &src,
               TypeDesc format = TypeDesc::UNKNOWN)
{
    ScopedGILRelease gil;
    return buf.copy (src, format);
}

bool
ImageBuf_copy2 (ImageBuf &buf, const ImageBuf &src,
                TypeDesc::BASETYPE format = TypeDesc::UNKNOWN)
{
    ScopedGILRelease gil;
    return buf.copy (src, format);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_copy_overloads,
                                ImageBuf_copy, 2, 3)
BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_copy2_overloads,
                                ImageBuf_copy2, 2, 3)

void
ImageBuf_set_full (ImageBuf &buf, int xbegin, int xend, int ybegin, int yend,
                   int zbegin, int zend)
{
    buf.set_full (xbegin, xend, ybegin, yend, zbegin, zend);
}

float
ImageBuf_getchannel (const ImageBuf &buf, int x, int y, int z, int c,
                     ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    return buf.getchannel(x, y, z, c, wrap);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_getchannel_overloads,
                                ImageBuf_getchannel, 5, 6)

object
ImageBuf_getpixel (const ImageBuf &buf, int x, int y, int z=0,
                   ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    int nchans = buf.nchannels();
    float *pixel = ALLOCA (float, nchans);
    buf.getpixel (x, y, z, pixel, nchans, wrap);
    PyObject *result = PyTuple_New (nchans);
    for (int i = 0;  i < nchans;  ++i)
        PyTuple_SetItem (result, i, PyFloat_FromDouble(pixel[i]));
    return object(handle<>(result));
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_getpixel_overloads,
                                ImageBuf_getpixel, 3, 5)

object
ImageBuf_interppixel (const ImageBuf &buf, float x, float y,
                      ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    int nchans = buf.nchannels();
    float *pixel = ALLOCA (float, nchans);
    buf.interppixel (x, y, pixel, wrap);
    PyObject *result = PyTuple_New (nchans);
    for (int i = 0;  i < nchans;  ++i)
        PyTuple_SetItem (result, i, PyFloat_FromDouble(pixel[i]));
    return object(handle<>(result));
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_interppixel_overloads,
                                ImageBuf_interppixel, 3, 4)

object
ImageBuf_interppixel_NDC (const ImageBuf &buf, float x, float y,
                          ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    int nchans = buf.nchannels();
    float *pixel = ALLOCA (float, nchans);
    buf.interppixel_NDC (x, y, pixel, wrap);
    return C_to_val_or_tuple (pixel, TypeDesc(TypeDesc::FLOAT,nchans),
                              PyFloat_FromDouble);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_interppixel_NDC_overloads,
                                ImageBuf_interppixel_NDC, 3, 4)

object
ImageBuf_interppixel_bicubic (const ImageBuf &buf, float x, float y,
                              ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    int nchans = buf.nchannels();
    float *pixel = ALLOCA (float, nchans);
    buf.interppixel_bicubic (x, y, pixel, wrap);
    PyObject *result = PyTuple_New (nchans);
    for (int i = 0;  i < nchans;  ++i)
        PyTuple_SetItem (result, i, PyFloat_FromDouble(pixel[i]));
    return object(handle<>(result));
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_interppixel_bicubic_overloads,
                                ImageBuf_interppixel_bicubic, 3, 4)

object
ImageBuf_interppixel_bicubic_NDC (const ImageBuf &buf, float x, float y,
                              ImageBuf::WrapMode wrap = ImageBuf::WrapBlack)
{
    int nchans = buf.nchannels();
    float *pixel = ALLOCA (float, nchans);
    buf.interppixel_bicubic_NDC (x, y, pixel, wrap);
    PyObject *result = PyTuple_New (nchans);
    for (int i = 0;  i < nchans;  ++i)
        PyTuple_SetItem (result, i, PyFloat_FromDouble(pixel[i]));
    return object(handle<>(result));
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_interppixel_bicubic_NDC_overloads,
                                ImageBuf_interppixel_bicubic_NDC, 3, 4)

void
ImageBuf_setpixel (ImageBuf &buf, int x, int y, int z, tuple p)
{
    std::vector<float> pixel;
    py_to_stdvector (pixel, p);
    if (pixel.size())
        buf.setpixel (x, y, z, &pixel[0], pixel.size());
}

void
ImageBuf_setpixel2 (ImageBuf &buf, int x, int y, tuple p)
{
    ImageBuf_setpixel (buf, x, y, 0, p);
}

void
ImageBuf_setpixel1 (ImageBuf &buf, int i, tuple p)
{
    std::vector<float> pixel;
    py_to_stdvector (pixel, p);
    if (pixel.size())
        buf.setpixel (i, &pixel[0], pixel.size());
}

object
ImageBuf_get_pixels (const ImageBuf &buf, TypeDesc format, ROI roi=ROI::All())
{
    // Allocate our own temp buffer and try to read the image into it.
    // If the read fails, return None.
    if (! roi.defined())
        roi = buf.roi();
    roi.chend = std::min (roi.chend, buf.nchannels()+1);

    size_t size = (size_t) roi.npixels() * roi.nchannels() * format.size();
    boost::scoped_array<char> data (new char [size]);
    if (! buf.get_pixels (roi, format, &data[0])) {
        return object(handle<>(Py_None));
    }

    return C_array_to_Python_array (data.get(), format, size);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_get_pixels_overloads,
                                ImageBuf_get_pixels, 2, 3)

object
ImageBuf_get_pixels_bt (const ImageBuf &buf, TypeDesc::BASETYPE format,
                        ROI roi=ROI::All())
{
    return ImageBuf_get_pixels (buf, TypeDesc(format), roi);
}

BOOST_PYTHON_FUNCTION_OVERLOADS(ImageBuf_get_pixels_bt_overloads,
                                ImageBuf_get_pixels_bt, 2, 3)

bool
ImageBuf_set_pixels_tuple (ImageBuf &buf, ROI roi, tuple data)
{
    if (! roi.defined())
        roi = buf.roi();
    roi.chend = std::min (roi.chend, buf.nchannels()+1);
    size_t size = (size_t) roi.npixels() * roi.nchannels();
    if (size == 0)
        return true;   // done
    std::vector<float> vals;
    py_to_stdvector (vals, data);
    if (size > vals.size())
        return false;   // Not enough data to fill our ROI
    buf.set_pixels (roi, TypeDesc::TypeFloat, &vals[0]);
    return true;
}

bool
ImageBuf_set_pixels_array (ImageBuf &buf, ROI roi, numeric::array data)
{
    if (! roi.defined())
        roi = buf.roi();
    roi.chend = std::min (roi.chend, buf.nchannels()+1);
    size_t size = (size_t) roi.npixels() * roi.nchannels();
    if (size == 0)
        return true;   // done

    TypeDesc type;
    size_t pylen = 0;
    const void *addr = python_array_address (data, type, pylen);
    if (!addr || size > pylen)
        return false;   // Not enough data to fill our ROI

    buf.set_pixels (roi, type, addr);
    return true;
}

DeepData&
ImageBuf_deepdataref (ImageBuf *ib)
{
    return *ib->deepdata();
}

void declare_imagebuf()
{
    enum_<ImageBuf::WrapMode>("WrapMode")
        .value("WrapDefault",  ImageBuf::WrapDefault )
        .value("WrapBlack",    ImageBuf::WrapBlack )
        .value("WrapClamp",    ImageBuf::WrapClamp )
        .value("WrapPeriodic", ImageBuf::WrapPeriodic )
        .value("WrapMirror",   ImageBuf::WrapMirror )
        .export_values();

    class_<ImageBuf, boost::noncopyable> ("ImageBuf")
        .def(init<const std::string&>())
        .def(init<const std::string&, int, int>())
        .def(init<const ImageSpec&>())

        .def("clear", &ImageBuf::clear)
        .def("reset", &ImageBuf_reset_name)
        .def("reset", &ImageBuf_reset_name2)
        .def("reset", &ImageBuf_reset_name_config)
        .def("reset", &ImageBuf_reset_spec)
        .add_property ("initialized", &ImageBuf::initialized)
        .def("init_spec", &ImageBuf::init_spec)
        .def("read",  &ImageBuf_read,
             ImageBuf_read_overloads())
        .def("read",  &ImageBuf_read2,
             ImageBuf_read2_overloads())
        .def("write", &ImageBuf_write,
             ImageBuf_write_overloads())
        // FIXME -- write(ImageOut&)
        .def("make_writeable", &ImageBuf_make_writeable,
             (arg("keep_cache_type")=false))
        .def("set_write_format", &ImageBuf_set_write_format)
        .def("set_write_tiles", &ImageBuf::set_write_tiles,
             (arg("width")=0, arg("height")=0, arg("depth")=0))

        .def("spec", &ImageBuf::spec,
                return_value_policy<copy_const_reference>())
        .def("nativespec", &ImageBuf::nativespec,
                return_value_policy<copy_const_reference>())
        .def("specmod", &ImageBuf::specmod,
             return_value_policy<reference_existing_object>())
        .add_property("name", &ImageBuf_name)
        .add_property("file_format_name", &ImageBuf_file_format_name)
        .add_property("subimage", &ImageBuf::subimage)
        .add_property("nsubimages", &ImageBuf::nsubimages)
        .add_property("miplevel", &ImageBuf::miplevel)
        .add_property("nmiplevels", &ImageBuf::nmiplevels)
        .add_property("nchannels", &ImageBuf::nchannels)
        .add_property("orientation", &ImageBuf::orientation,
                      &ImageBuf::set_orientation)
        .add_property("oriented_width", &ImageBuf::oriented_width)
        .add_property("oriented_height", &ImageBuf::oriented_height)
        .add_property("oriented_x", &ImageBuf::oriented_x)
        .add_property("oriented_y", &ImageBuf::oriented_y)
        .add_property("oriented_full_width", &ImageBuf::oriented_full_width)
        .add_property("oriented_full_height", &ImageBuf::oriented_full_height)
        .add_property("oriented_full_x", &ImageBuf::oriented_full_x)
        .add_property("oriented_full_y", &ImageBuf::oriented_full_y)
        .add_property("xbegin", &ImageBuf::xbegin)
        .add_property("xend", &ImageBuf::xend)
        .add_property("ybegin", &ImageBuf::ybegin)
        .add_property("yend", &ImageBuf::yend)
        .add_property("zbegin", &ImageBuf::zbegin)
        .add_property("zend", &ImageBuf::zend)
        .add_property("xmin", &ImageBuf::xmin)
        .add_property("xmax", &ImageBuf::xmax)
        .add_property("ymin", &ImageBuf::ymin)
        .add_property("ymax", &ImageBuf::ymax)
        .add_property("zmin", &ImageBuf::zmin)
        .add_property("zmax", &ImageBuf::zmax)
        .add_property("roi", &ImageBuf::roi)
        .add_property("roi_full",
                      &ImageBuf::roi_full, &ImageBuf::set_roi_full)
        .def("set_full", &ImageBuf_set_full)

        .add_property("pixels_valid", &ImageBuf::pixels_valid)
        .add_property("pixeltype", &ImageBuf::pixeltype)
        .add_property("has_error",   &ImageBuf::has_error)
        .def("geterror",    &ImageBuf::geterror)

        .def("copy_metadata", &ImageBuf::copy_metadata)
        .def("copy_pixels", &ImageBuf::copy_pixels)
        .def("copy",  &ImageBuf_copy,
             ImageBuf_copy_overloads())
        .def("copy",  &ImageBuf_copy2,
             ImageBuf_copy2_overloads())
        .def("swap", &ImageBuf::swap)

        .def("getchannel", &ImageBuf_getchannel,
             ImageBuf_getchannel_overloads())
        .def("getpixel", &ImageBuf_getpixel,
             ImageBuf_getpixel_overloads())
        .def("interppixel", &ImageBuf_interppixel,
             ImageBuf_interppixel_overloads())
        .def("interppixel_NDC", &ImageBuf_interppixel_NDC,
             ImageBuf_interppixel_NDC_overloads())
        .def("interppixel_NDC_full", &ImageBuf_interppixel_NDC,
             ImageBuf_interppixel_NDC_overloads())
        .def("interppixel_bicubic", &ImageBuf_interppixel_bicubic,
             ImageBuf_interppixel_bicubic_overloads())
        .def("interppixel_bicubic_NDC", &ImageBuf_interppixel_bicubic_NDC,
             ImageBuf_interppixel_bicubic_NDC_overloads())
        .def("setpixel", &ImageBuf_setpixel)
        .def("setpixel", &ImageBuf_setpixel2)
        .def("setpixel", &ImageBuf_setpixel1)
        .def("get_pixels", &ImageBuf_get_pixels, ImageBuf_get_pixels_overloads())
        .def("get_pixels", &ImageBuf_get_pixels_bt, ImageBuf_get_pixels_bt_overloads())
        .def("set_pixels", &ImageBuf_set_pixels_tuple)
        .def("set_pixels", &ImageBuf_set_pixels_array)

        .add_property("deep", &ImageBuf::deep)
        .def("deep_samples", &ImageBuf::deep_samples,
             (arg("x"), arg("y"), arg("z")=0))
        .def("set_deep_samples", &ImageBuf::set_deep_samples)
        .def("deep_value", &ImageBuf::deep_value)
        .def("deep_value_uint", &ImageBuf::deep_value_uint)
        .def("set_deep_value", &ImageBuf::set_deep_value)
        .def("set_deep_value_uint", &ImageBuf::set_deep_value_uint)
        .def("deep_alloc", &ImageBuf::deep_alloc)
        .def("deepdata", &ImageBuf_deepdataref,
             return_value_policy<reference_existing_object>())

        // FIXME -- do we want to provide pixel iterators?
    ;

}

} // namespace PyOpenImageIO

max_stars_repo_name: ShaderManager/oiio

max_issues_repo_name: ShaderManager/oiio

max_forks_repo_name: ShaderManager/oiio


text: #include <boost/algorithm/is_partitioned_until.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/*
 * #%L
 * OME-FILES C++ library for image IO.
 * Copyright © 2015–2017 Open Microscopy Environment:
 *   - Massachusetts Institute of Technology
 *   - National Institutes of Health
 *   - University of Dundee
 *   - Board of Regents of the University of Wisconsin-Madison
 *   - Glencoe Software, Inc.
 * %%
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 * 2. 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.
 *
 * 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 HOLDERS 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.
 *
 * The views and conclusions contained in the software and documentation are
 * those of the authors and should not be interpreted as representing official
 * policies, either expressed or implied, of any organization.
 * #L%
 */

#include <iostream>
#include <memory>

#include <boost/filesystem/path.hpp>

#include <ome/files/CoreMetadata.h>
#include <ome/files/MetadataTools.h>
#include <ome/files/VariantPixelBuffer.h>
#include <ome/files/out/OMETIFFWriter.h>
#include <ome/xml/meta/OMEXMLMetadata.h>

using boost::filesystem::path;
using std::make_shared;
using std::shared_ptr;
using std::static_pointer_cast;
using ome::files::dimension_size_type;
using ome::files::createID;
using ome::files::fillMetadata;
using ome::files::CoreMetadata;
using ome::files::DIM_SPATIAL_X;
using ome::files::DIM_SPATIAL_Y;
using ome::files::DIM_SUBCHANNEL;
using ome::files::FormatWriter;
using ome::files::MetadataMap;
using ome::files::out::OMETIFFWriter;
using ome::files::PixelBuffer;
using ome::files::PixelBufferBase;
using ome::files::PixelProperties;
using ome::files::VariantPixelBuffer;
using ome::xml::model::enums::Binning;
using ome::xml::model::enums::Immersion;
using ome::xml::model::enums::PixelType;
using ome::xml::model::enums::DetectorType;
using ome::xml::model::enums::DimensionOrder;
using ome::xml::model::enums::UnitsLength;
using ome::xml::meta::MetadataRetrieve;
using ome::xml::meta::MetadataStore;
using ome::xml::meta::Metadata;
using ome::xml::meta::OMEXMLMetadata;

namespace
{

  shared_ptr<OMEXMLMetadata>
  createMetadata()
  {
    /* create-metadata-start */
    // OME-XML metadata store.
    auto meta = make_shared<OMEXMLMetadata>();

    // Create simple CoreMetadata and use this to set up the OME-XML
    // metadata.  This is purely for convenience in this example; a
    // real writer would typically set up the OME-XML metadata from an
    // existing MetadataRetrieve instance or by hand.
    std::vector<shared_ptr<CoreMetadata>> seriesList;
    shared_ptr<CoreMetadata> core(make_shared<CoreMetadata>());
    core->sizeX = 512U;
    core->sizeY = 512U;
    core->sizeC.clear(); // defaults to 1 channel with 1 subchannel; clear this
    core->sizeC.push_back(3U); // replace with single RGB channel
    core->pixelType = ome::xml::model::enums::PixelType::UINT16;
    core->interleaved = false;
    core->bitsPerPixel = 12U;
    core->dimensionOrder = DimensionOrder::XYZTC;
    seriesList.push_back(core);
    seriesList.push_back(core); // add two identical series

    fillMetadata(*meta, seriesList);
    /* create-metadata-end */

    return meta;
  }

  void
  addExtendedMetadata(shared_ptr<OMEXMLMetadata> store)
  {
    /* extended-metadata-start */
    // There is one image with one channel in this image.
    MetadataStore::index_type image_idx = 0;
    MetadataStore::index_type channel_idx = 0;
    MetadataStore::index_type annotation_idx = 0;

    // Create an Instrument.
    MetadataStore::index_type instrument_idx = 0;
    std::string instrument_id = createID("Instrument", instrument_idx);
    store->setInstrumentID(instrument_id, instrument_idx);

    // Create an Objective for this Instrument.
    MetadataStore::index_type objective_idx = 0;
    std::string objective_id = createID("Objective",
                                        instrument_idx, objective_idx);
    store->setObjectiveID(objective_id, instrument_idx, objective_idx);
    store->setObjectiveManufacturer("InterFocal", instrument_idx, objective_idx);
    store->setObjectiveNominalMagnification(40, instrument_idx, objective_idx);
    store->setObjectiveLensNA(0.4, instrument_idx, objective_idx);
    store->setObjectiveImmersion(Immersion::OIL, instrument_idx, objective_idx);
    store->setObjectiveWorkingDistance({0.34, UnitsLength::MILLIMETER},
                                       instrument_idx, objective_idx);

    // Create a Detector for this Instrument.
    MetadataStore::index_type detector_idx = 0;
    std::string detector_id = createID("Detector", instrument_idx, detector_idx);
    store->setDetectorID(detector_id, instrument_idx, detector_idx);
    store->setObjectiveManufacturer("MegaCapture", instrument_idx, detector_idx);
    store->setDetectorType(DetectorType::CCD, instrument_idx, detector_idx);

    // Create Settings for this Detector for the Channel on the Image.
    store->setDetectorSettingsID(detector_id, image_idx, channel_idx);
    store->setDetectorSettingsBinning(Binning::TWOBYTWO, image_idx, channel_idx);
    store->setDetectorSettingsGain(56.89, image_idx, channel_idx);
    /* extended-metadata-end */

    /* annotations-start */
    // Create a MapAnnotation.
    MetadataStore::index_type map_annotation_idx = 0;
    std::string annotation_id = createID("Annotation", annotation_idx);
    store->setMapAnnotationID(annotation_id, map_annotation_idx);
    store->setMapAnnotationNamespace
      ("https://microscopy.example.com/colour-balance", map_annotation_idx);
    ome::xml::model::primitives::OrderedMultimap map;
    map.push_back({"white-balance", "5,15,8"});
    map.push_back({"black-balance", "112,140,126"});
    store->setMapAnnotationValue(map, map_annotation_idx);

    // Link MapAnnotation to Detector.
    MetadataStore::index_type detector_ref_idx = 0;
    store->setDetectorAnnotationRef(annotation_id, instrument_idx, detector_idx,
                                    detector_ref_idx);

    // Create a LongAnnotation.
    ++annotation_idx;
    MetadataStore::index_type long_annotation_idx = 0;
    annotation_id = createID("Annotation", annotation_idx);
    store->setLongAnnotationID(annotation_id, long_annotation_idx);
    store->setLongAnnotationValue(239423, long_annotation_idx);
    store->setLongAnnotationNamespace
      ("https://microscopy.example.com/trigger-delay", long_annotation_idx);

    // Link LongAnnotation to Image.
    MetadataStore::index_type image_ref_idx = 0;
    store->setImageAnnotationRef(annotation_id, image_idx, image_ref_idx);

    // Create a second LongAnnotation.
    ++annotation_idx;
    ++long_annotation_idx;
    annotation_id = createID("Annotation", annotation_idx);
    store->setLongAnnotationID(annotation_id, long_annotation_idx);
    store->setLongAnnotationValue(934223, long_annotation_idx);
    store->setLongAnnotationNamespace
      ("https://microscopy.example.com/sample-number", long_annotation_idx);

    // Link second LongAnnotation to Image.
    ++image_ref_idx = 0;
    store->setImageAnnotationRef(annotation_id, image_idx, image_ref_idx);

    // Update all the annotation cross-references.
    store->resolveReferences();
    /* annotations-end */
  }

  void
  writePixelData(FormatWriter& writer,
                 std::ostream& stream)
  {
    /* pixel-example-start */
    // Total number of images (series)
    dimension_size_type ic = writer.getMetadataRetrieve()->getImageCount();
    stream << "Image count: " << ic << '\n';

    // Loop over images
    for (dimension_size_type i = 0 ; i < ic; ++i)
      {
        // Change the current series to this index
        writer.setSeries(i);

        // Total number of planes.
        dimension_size_type pc = 1U;
        pc *= writer.getMetadataRetrieve()->getPixelsSizeZ(i);
        pc *= writer.getMetadataRetrieve()->getPixelsSizeT(i);
        pc *= writer.getMetadataRetrieve()->getChannelCount(i);
        stream << "\tPlane count: " << pc << '\n';

        // Loop over planes (for this image index)
        for (dimension_size_type p = 0 ; p < pc; ++p)
          {
            // Pixel buffer; size 512 × 512 with 3 subchannels of type
            // uint16_t.  It uses the native endianness and has a
            // storage order of XYZTC without interleaving
            // (subchannels are planar).
            shared_ptr<PixelBuffer<PixelProperties<PixelType::UINT16>::std_type>>
              buffer(make_shared<PixelBuffer<PixelProperties<PixelType::UINT16>::std_type>>
                     (boost::extents[512][512][1][1][1][3][1][1][1],
                      PixelType::UINT16, ome::files::ENDIAN_NATIVE,
                      PixelBufferBase::make_storage_order(DimensionOrder::XYZTC, false)));

            // Fill each subchannel with a different intensity ramp in
            // the 12-bit range.  In a real program, the pixel data
            // would typically be obtained from data acquisition or
            // another image.
            for (dimension_size_type x = 0; x < 512; ++x)
              for (dimension_size_type y = 0; y < 512; ++y)
                {
                  PixelBufferBase::indices_type idx;
                  std::fill(idx.begin(), idx.end(), 0);
                  idx[DIM_SPATIAL_X] = x;
                  idx[DIM_SPATIAL_Y] = y;

                  idx[DIM_SUBCHANNEL] = 0;
                  buffer->at(idx) = (static_cast<float>(x) / 512.0f) * 4096.0f;
                  idx[DIM_SUBCHANNEL] = 1;
                  buffer->at(idx) = (static_cast<float>(y) / 512.0f) * 4096.0f;
                  idx[DIM_SUBCHANNEL] = 2;
                  buffer->at(idx) = (static_cast<float>(x+y) / 1024.0f) * 4096.0f;
                }

            VariantPixelBuffer vbuffer(buffer);
            stream << "PixelBuffer PixelType is " << buffer->pixelType() << '\n';
            stream << "VariantPixelBuffer PixelType is " << vbuffer.pixelType() << '\n';
            stream << std::flush;

            // Write the the entire pixel buffer to the plane.
            writer.saveBytes(p, vbuffer);

            stream << "Wrote " << buffer->num_elements() << ' '
                   << buffer->pixelType() << " pixels\n";
          }
      }
    /* pixel-example-end */
  }

}

int
main(int argc, char *argv[])
{
  // This is the default, but needs setting manually on Windows.
  ome::common::setLogLevel(ome::logging::trivial::warning);

  try
    {
      if (argc > 1)
        {
          // Portable path
          path filename(argv[1]);

          /* writer-example-start */
          // Create minimal metadata for the file to be written.
          auto meta = createMetadata();
          // Add extended metadata.
          addExtendedMetadata(meta);

          // Create TIFF writer
          auto writer = make_shared<OMETIFFWriter>();

          // Set writer options before opening a file
          auto retrieve = static_pointer_cast<MetadataRetrieve>(meta);
          writer->setMetadataRetrieve(retrieve);
          writer->setInterleaved(false);
          writer->setTileSizeX(256);
          writer->setTileSizeY(256);

          // Open the file
          writer->setId(filename);

          // Write pixel data
          writePixelData(*writer, std::cout);

          // Explicitly close writer
          writer->close();
          /* writer-example-end */
        }
      else
        {
          std::cerr << "Usage: " << argv[0] << " ome-xml.ome.tiff\n";
          std::exit(1);
        }
    }
  catch (const std::exception& e)
    {
      std::cerr << "Caught exception: " << e.what() << '\n';
      std::exit(1);
    }
  catch (...)
    {
      std::cerr << "Caught unknown exception\n";
      std::exit(1);
    }
}

max_stars_repo_name: rleigh-dundee/bioformats-cpp

max_issues_repo_name: rleigh-dundee/bioformats-cpp

max_forks_repo_name: rleigh-dundee/bioformats-cpp