text:
/**************************************************************************************************
 Software License Agreement (BSD License)

 Copyright (c) 2011-2013, LAR toolkit developers - University of Aveiro - http://lars.mec.ua.pt
 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 University of Aveiro 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.
***************************************************************************************************/
/*! \brief Code for feature extraction
 *
 *  extract features from targets received from MTT wrt the robot
 *  the features are printed in the command line and should be stored
 *  in a file (txt) for further processing with matlab
 *  the features will be used for training a classifier
 *
 *  it is also responsible of publishing these features to a matlab
 *  classifier that will then output a message with the classification
 *  of the target based on its features
 */

#include "ros/ros.h"
#include "std_msgs/String.h"
#include "tf/transform_listener.h"
#include "mtt/TargetList.h"
#include <boost/circular_buffer.hpp>
#include <boost/bind.hpp>
#include <boost/ref.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/weighted_sum.hpp>
#include <boost/accumulators/statistics/variance.hpp>
#include <numeric>

using namespace boost::accumulators;

// static ros::Time start_time;
ros::Time time_last_msg(0);
ros::Duration time_elapsed;
// ros::Duration duration_btw_msg;
// ros::Publisher nfeatures_pub;
// ros::Publisher marker_pub;

tf::TransformListener *p_listener;
tf::StampedTransform transform;

geometry_msgs::Twist twist;
geometry_msgs::Twist features;
// geometry_msgs::PoseWithCovariance nfeatures;

// std::list<geometry_msgs::PoseWithCovariance> matlab_list;

//statics
boost::circular_buffer<double> robot_posex_buffer(30);
// accumulator_set<double, stats<tag::mean> > acc;
// accumulator_set<double, features< tag::sum >, int > acc;

double target_x, target_y;
double target_vel, target_theta;
double robot_x, robot_y;
double robot_vel, robot_theta;
double position_diff, heading_diff;
double angle_to_robot;
double velocity_diff;
uint leader_tag = 0;
uint target_id = 0;
uint single_target = 0;
// uint counter = 0;

void targetsCallback(const mtt::TargetList& list)
{
  //will print information that should be stored in a file
  //file format: id, good/bad tag, time, pos x, pos y, vel, theta
  //             position_diff, heading_diff, angle_to_robot, velocity_diff
  static ros::Time start_time = ros::Time::now();
  time_elapsed = ros::Time::now() - start_time;

  /// /// ROBOT PART //////
  //use transformations to extract robot features
  try{
    p_listener->lookupTransform("/map", "/base_link", ros::Time(0), transform);
    p_listener->lookupTwist("/map", "/base_link", ros::Time(0), ros::Duration(0.5), twist);
  }
  catch (tf::TransformException ex){
    ROS_ERROR("%s",ex.what());
  }

  robot_x = transform.getOrigin().x();

  robot_posex_buffer.push_back(robot_x);

  robot_y = transform.getOrigin().y();
  robot_theta = tf::getYaw(transform.getRotation());
  robot_vel = sqrt(pow(twist.linear.x, 2) + pow(twist.linear.y, 2));

  //robot output line
  /// uncomment the following for training!
  printf("%d,%d,%.10f,%.10f,%.10f,%.10f,%.10f,0,0,0,0\n",
         -1, leader_tag, time_elapsed.toSec(),
         robot_x, robot_y, robot_vel, robot_theta);

  //testing new features extraction
//   accumulator_set<double, stats<tag::variance> > acc;
//   for_each(robot_posex_buffer.begin(), robot_posex_buffer.end(), boost::bind<void>(boost::ref(acc), _1));
//   printf("%f,%f,%f\n", robot_x, mean(acc), sqrt(variance(acc)));

  /// /// TARGETS PART //////
  //sweeps target list and extract features
  for(uint i = 0; i < list.Targets.size(); i++){
    target_id = list.Targets[i].id;
    target_x = list.Targets[i].pose.position.x;
    target_y = list.Targets[i].pose.position.y;
    target_theta = tf::getYaw(list.Targets[i].pose.orientation);
    target_vel = sqrt(pow(list.Targets[i].velocity.linear.x,2)+
                      pow(list.Targets[i].velocity.linear.y,2));
    position_diff = sqrt(pow(robot_x - target_x,2)+
                        pow(robot_y - target_y,2));
    heading_diff = robot_theta - target_theta;
    angle_to_robot = -robot_theta + atan2(target_y - robot_y,
                                        target_x - robot_x );
    velocity_diff = robot_vel - target_vel;

    //target output (to be used in adaboost training)

    // % output file format:
    // % 1: id
    // % 2: good/bad tag
    // % 3: time
    // % 4: pos x
    // % 5: pos y
    // % 6: vel
    // % 7: theta
    // % 8: pos diff
    // % 9: head diff
    // %10: angle 2 robot
    // %11: velocity diff

    /// uncomment the following to generate training file!
    printf("%d,%d,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f\n",
      target_id, leader_tag, time_elapsed.toSec(),
      target_x, target_y, target_vel, target_theta,
      position_diff, heading_diff, angle_to_robot, velocity_diff);

//     if(position_diff < 6.0 && target_vel > 0.5){
//       //if inside boundaries (in meters)
//       //store features in a covariance struct to send to matlab
//       nfeatures.pose.position.x = target_x;
//       nfeatures.pose.position.y = target_y;
//       nfeatures.pose.position.z = target_id;
//
//       nfeatures.covariance[0] = target_vel;
//       nfeatures.covariance[1] = velocity_diff;
//       nfeatures.covariance[2] = heading_diff;
//       nfeatures.covariance[3] = angle_to_robot;
//       nfeatures.covariance[4] = position_diff;
//
//       matlab_list.push_back(nfeatures);
// //       counter++;
//     }
  }
//   printf("targets within range: %d\n",counter);
//   counter = 0;

  //check if enough time has passed
  //and send batch of msgs to matlab
//   duration_btw_msg = ros::Time::now() - time_last_msg;
//
//   if(duration_btw_msg.toSec() > 0.01){
//     while(!matlab_list.empty()){
//       nfeatures_pub.publish(matlab_list.front());
//       matlab_list.pop_front();
//       usleep(0.01e6); //has to sleep, otherwise matlab do not get the msg
//     }
//     time_last_msg = ros::Time::now();
//   }
}

void tagCallback(const std_msgs::Header& tag)
{
  //ROS_INFO("tag received");
  //only works for a transition good/bad leader
  leader_tag = 1;
}

// void drawCallback(const geometry_msgs::Pose& input)
// {
//   visualization_msgs::Marker marker;
//
//   marker.header.frame_id = "/map";
//   marker.header.stamp = ros::Time::now();
//
//   marker.ns = "quality";
//   marker.id = 0;
//
//   marker.type = visualization_msgs::Marker::CYLINDER;
//   marker.action = visualization_msgs::Marker::ADD;
//
//   // Set the pose of the marker.  This is a full 6DOF pose relative to the frame/time specified in the header
//   marker.pose.position.x = input.position.x;
//   marker.pose.position.y = input.position.y;
//   marker.pose.position.z = 0;
//   marker.pose.orientation.x = 0.0;
//   marker.pose.orientation.y = 0.0;
//   marker.pose.orientation.z = 0.0;
//   marker.pose.orientation.w = 1.0;
//
//   // Set the scale of the marker
//   marker.scale.x = 0.5;
//   marker.scale.y = 0.5;
//   marker.scale.z = 0.5;
//
//   //good or bad leader
//   if(input.position.z == -1){
//     marker.color.r = 0.0f;
//     marker.color.g = 1.0f;
//   }
//   else if (input.position.z == 1){
//     marker.color.r = 1.0f;
//     marker.color.g = 0.0f;
//   }
//
//   marker.color.b = 0.0f;
//   marker.color.a = 1.0;
//
//   marker.lifetime = ros::Duration();
//
//   // Publish the marker
//   marker_pub.publish(marker);
// }

int main(int argc, char **argv)
{
  ros::init(argc, argv, "process_target");
  ros::NodeHandle n;

  tf::TransformListener listener;
  p_listener=&listener;

  //for testing purposes
//   single_target = atoi(argv[1]);
//   ROS_INFO("chosen target:%d",single_target);

  ros::Subscriber sub = n.subscribe("/targets", 1000, targetsCallback);
  ros::Subscriber sub_tag = n.subscribe("/timetag", 1000, tagCallback);
//   ros::Subscriber sub_draw = n.subscribe("/pose_to_draw", 1000, drawCallback);
//   nfeatures_pub = n.advertise<geometry_msgs::PoseWithCovariance>("/example_topic", 100);
//   marker_pub = n.advertise<visualization_msgs::Marker>("/leader_quality", 1);

  ros::spin();

  return 0;
}

max_stars_repo_name: procopiostein/leader

max_issues_repo_name: procopiostein/leader

max_forks_repo_name: procopiostein/leader


text:
/**************************************************************
 *
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you 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.
 *
 *************************************************************/

#ifndef INCLUDED_SLIDESHOW_EFFECT_REWINDER_HXX
#define INCLUDED_SLIDESHOW_EFFECT_REWINDER_HXX

#include "animationnode.hxx"
#include "eventhandler.hxx"
#include "animationeventhandler.hxx"
#include "event.hxx"
#include "screenupdater.hxx"

#include <com/sun/star/presentation/XSlideShow.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/function.hpp>
#include <vector>

namespace css = ::com::sun::star;

namespace slideshow { namespace internal {

class EventMultiplexer;
class EventQueue;
class UserEventQueue;

/** Rewind single effects of the main effect sequence.  A rewind is
    initiated by calling the Rewind() method.  Part of the processing is
    done asynchronously.  Multiple EventQueue::update() calls may be
    necessary to finish a rewind.

    Remember to call SetRootAnimationNode() when switching to a different
    slide so that the EffectRewinder can determine the number of main
    sequence effects.
*/
class EffectRewinder
{
public:
    EffectRewinder (
        EventMultiplexer& rEventMultiplexer,
        EventQueue& rEventQueue,
        UserEventQueue& rUserEventQueue);
    ~EffectRewinder (void);

    /** Call Dispose() before the ownder of an EffectRewinder object dies so
        that the EffectRewinder can release all references to the owner.

    */
    void dispose (void);

    /** Store the root node of the animation tree.  It is used in
        CountMainSequenceEffects() to count the number of main sequence
        effects (or effect groups.)
    */
    void setRootAnimationNode (
        const css::uno::Reference<css::animations::XAnimationNode>& xRootNode);

    /** Rewind one effect of the main effect sequence.  When the current
        slide has not effects or no main sequence effect has yet been played
        then switch to the previous slide and replay all of its main
        sequence effects.
        The caller has to pass two functors that redisplay the current slide
        or switch to the previous slide so that it does not have to expose
        its internals to us.  Only one of the two functors is called.
        @param rpPaintLock
            This paint lock is released after the whole asynchronous
            procoess  of rewinding the current effect is completed.  It
            prevents intermediate repaints  that would show partial replay
            of effects.
        @param rSlideRewindFunctor
            This functor is called when the current slide is to be
            redisplayed.  When it is called then the other functor is not
            called.
        @param rPreviousSlideFunctor
            This functor is called to switch to the previous slide.  When it
            is called then the other functor is not called.
    */
    bool rewind (
        const ::boost::shared_ptr<ScreenUpdater::UpdateLock>& rpPaintLock,
        const ::boost::function<void(void)>& rSlideRewindFunctor,
        const ::boost::function<void(void)>& rPreviousSlideFunctor);

    /** Call this method after gotoPreviousEffect() triggered a slide change
        to the previous slide.
    */
    void skipAllMainSequenceEffects (void);

private:
    EventMultiplexer& mrEventMultiplexer;
    EventQueue& mrEventQueue;
    UserEventQueue& mrUserEventQueue;

    EventHandlerSharedPtr mpSlideStartHandler;
    EventHandlerSharedPtr mpSlideEndHandler;
    AnimationEventHandlerSharedPtr mpAnimationStartHandler;

    /** The number off main sequence effects so far.
    */
    sal_Int32 mnMainSequenceEffectCount;

    /** This is the currently scheduled event that executes the asynchronous
        part of the effect rewinding.  It is also used as flag that prevents
        nested rewinds.
    */
    EventSharedPtr mpAsynchronousRewindEvent;

    css::uno::Reference<css::animations::XAnimationNode> mxCurrentAnimationRootNode;
    ::boost::shared_ptr<ScreenUpdater::UpdateLock> mpPaintLock;

    bool mbNonUserTriggeredMainSequenceEffectSeen;

    void initialize (void);

    bool resetEffectCount (void);
    /** Called by listeners when an animation (not necessarily of a main
        sequence effect) starts.
    */
    bool notifyAnimationStart (const AnimationNodeSharedPtr& rpNode);

    /** Count the number of effects (or effect groups) in the main effect
        sequence.
    */
    sal_Int32 countMainSequenceEffects (void);

    /** Skip the next main sequence effect.
    */
    void skipSingleMainSequenceEffects (void);

    /** Skip the specified number of main sequence effects.
    */
    void skipSomeMainSequenceEffects (const sal_Int32 nSkipCount);

    /** Rewind the last effect of the main effect sequence by replaying all
        previous effects.
        @param nEffectCount
            The number of main sequence effects to replay.
        @param bRedisplayCurrentSlide
            When <TRUE/> then the current slide is redisplayed before the
            effects are replayed.
        @param rSlideRewindFunctor
            This functor is used to redisplay the current slide.
    */
    void asynchronousRewind (
        sal_Int32 nEffectCount,
        const bool bRedisplayCurrentSlide,
        const boost::function<void(void)>& rSlideRewindFunctor);

    /** Go to the previous slide and replay all of its main sequence effects
        (or effect groups).
        @param rPreviousSlideFunctor
            This functor is used to go to the previous slide.
    */
    void asynchronousRewindToPreviousSlide (
        const ::boost::function<void(void)>& rPreviousSlideFunctor);
};

} } // end of namespace ::slideshow::internal

#endif

ext: hxx

max_stars_repo_name: Grosskopf/openoffice

max_issues_repo_name: Grosskopf/openoffice

max_forks_repo_name: Grosskopf/openoffice


text:
#define NOHASH 1

#include <iostream>
#include <fstream>
#include <map>
#include <algorithm>

#include <CGAL/Surface_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
#include <CGAL/boost/graph/properties_Polyhedron_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Timer.h>
#include<boost/range/iterator_range.hpp>
#include <boost/unordered_map.hpp>
#include <unordered_map>

#include <boost/random/random_number_generator.hpp>
#include <boost/random/linear_congruential.hpp>

typedef CGAL::Simple_cartesian<double>  Kernel;
typedef Kernel::Point_3                 Point_3;
typedef Kernel::Vector_3                Vector_3;

typedef CGAL::Timer                     Timer;

template <typename G, typename Map>
void
run(const G& g)
{
  typedef typename boost::graph_traits<G>::vertex_descriptor vertex_descriptor;

  typedef typename boost::property_map<G,CGAL::vertex_point_t>::const_type  VPM;
  VPM vpm = get(CGAL::vertex_point,g);

  std::vector<vertex_descriptor> V, V2;
  std::vector<Point_3> P1, P2;

  for(vertex_descriptor vd : vertices(g)){
    V.push_back(vd);
    V2.push_back(vd);
  }

  boost::rand48 random;
  boost::random_number_generator<boost::rand48> rng(random);
  std::random_shuffle(V.begin(), V.end(), rng);

  Timer t;
#if 0
  t.start();
  Map vm;
  for(vertex_descriptor vd : V){
    vm[vd] = get(vpm,vd);
  }
  t.stop();  std::cerr << "Insertion:  " << t.time() << " sec.     " << std::endl;
#endif

  Vector_3 v(0,0,0);
  std::size_t st=0;

#if 0
  std::cerr << "foreach std::vector<vertex_descriptor)\n";
  t.reset(); t.start();
  for(int i=0; i<100; i++){
  for(vertex_descriptor vd : V2){
#ifdef NOHASH
    st += std::size_t(vd);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
  }

  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;
#endif

#if 1
  std::cerr << "foreach boost::iterator_range r = vertices(g))\n";

  t.reset(); t.start();
  for(int i=0; i<100; i++){
  boost::iterator_range<typename boost::graph_traits<G>::vertex_iterator> r = vertices(g);

  for(vertex_descriptor vd : r) {
#ifdef NOHASH
    st += std::size_t(vd);
    // v = v + (get(vpm,vd) - CGAL::ORIGIN);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
  }
  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;
#endif

#if 1
   std::cerr << "foreach CGAL::Iterator_range r = vertices(g))\n";
  t.reset(); t.start();
  for(int i=0; i<100; i++){
  CGAL::Iterator_range<typename boost::graph_traits<G>::vertex_iterator> ir = vertices(g);
  for(vertex_descriptor vd : ir) {
#ifdef NOHASH
    st += std::size_t(vd);
    //v = v + (get(vpm,vd) - CGAL::ORIGIN);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
  }
  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;

#endif

#if 1
  std::cerr << "foreach vertices(g))\n";
  t.reset(); t.start();
  for(int i=0; i<100; i++){
  for(vertex_descriptor vd : vertices(g)) {
#ifdef NOHASH
    st += std::size_t(vd);
    //v = v + (get(vpm,vd) - CGAL::ORIGIN);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
  }
  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;
#endif

#if 1
  std::cerr << "for vertices(g))\n";
  t.reset(); t.start();
  for(int i=0; i<100; i++){
    for(vertex_descriptor vd : vertices(g)) {
#ifdef NOHASH
    st += std::size_t(vd);
    //v = v + (get(vpm,vd) - CGAL::ORIGIN);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
  }
  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;

#endif

#if 0
  std::cerr << "boost::tie(vb,ve) = vertices(g);\n";
  t.reset(); t.start();
    for(int i=0; i<100; i++){
  typename boost::graph_traits<G>::vertex_iterator vb, ve;
  boost::tie(vb,ve) = vertices(g);
  for(; vb != ve; ++vb) {
    vertex_descriptor vd = *vb;
#ifdef NOHASH
    st += std::size_t(vd);
    //v = v + (get(vpm,vd) - CGAL::ORIGIN);
#else
    typename Map::iterator it = vm.find(vd);
    v = v + ((*it).second - CGAL::ORIGIN);
#endif
  }
    }
  t.stop();  std::cerr << "  " <<t.time() << " sec.     " << std::endl;
#endif

  std::cerr << "v = " << v << "  " << st << std::endl;

}

struct blob {
  int a, b, c, d, e, f, g;
};

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

  {
    typedef CGAL::Surface_mesh<Point_3>     Mesh;
    typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
    typedef std::map<vertex_descriptor,Point_3> SM;
    typedef std::unordered_map<vertex_descriptor,Point_3> SUM;
    typedef boost::unordered_map<vertex_descriptor,Point_3> BUM;

    Mesh m;
    std::ifstream input(argv[1]);
    input >> m;

    std::cerr << num_vertices(m) << " items\n";
    std::cerr << "\nSurface_mesh  std::map"<< std::endl;
    run<Mesh,SM>(m);
#if 0
    std::cerr << "\nSurface_mesh  std::unordered_map"<< std::endl;
    run<Mesh,SUM>(m);
    std::cerr << "\nSurface_mesh  boost::unordered_map"<< std::endl;
    run<Mesh,BUM>(m);
#endif
  }

#if 0
  {
    typedef CGAL::Polyhedron_3<Kernel>      Mesh;
    typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
    typedef std::map<vertex_descriptor,Point_3> SM;
    typedef std::unordered_map<vertex_descriptor,Point_3> SUM;
    typedef boost::unordered_map<vertex_descriptor,Point_3> BUM;

    Mesh m;
    std::ifstream input(argv[1]);
    input >> m;

    std::cerr << "\nPolyhedron_3  std::map" << std::endl;
    run<Mesh,SM>(m);
    std::cerr << "\nPolyhedron_3 std::unordered_map"<< std::endl;
    run<Mesh,SUM>(m);
    std::cerr << "\nPolyhedron_3 boost::unordered_map"<< std::endl;
    run<Mesh,BUM>(m);
  }

  {
    const int N = 3165798;
    std::cerr << "\nHashing "<< N << " pointers\n";
    std::vector<int> ints(N);

    std::vector<blob> data(N);
    for(int i =0; i <N ; i++){
      ints[i]=i;
    }
    std::random_shuffle(ints.begin(), ints.end());

    {
      boost::unordered_map<blob*,int> um;
      Timer t;
      t.start();
      for(int i= 0; i < N; i++){
        um[& (data[ints[i]])] = i;
      }
      t.stop();
      std::cerr << " boost::unordered_map: " <<  t.time() << " sec.\n";
    }

    {
      std::unordered_map<blob*,int> um;
      Timer t;
      t.start();
      for(int i= 0; i < N; i++){
        um[& (data[ints[i]])] = i;
      }
      t.stop();
      std::cerr  << " std::unordered_map: " << t.time() << " sec.\n";
    }
  }
#endif

  return 0;
}

ext: cpp

max_stars_repo_name: ffteja/cgal

max_issues_repo_name: ffteja/cgal

max_forks_repo_name: ffteja/cgal


text:
#include <thread>
#include <chrono>
#include <iostream>       // std::cout
#include <future>         // std::async, std::future
#include <boost/asio.hpp>
#include <csi-async/async.h>

void async_sleep(boost::asio::io_service& ios, int64_t ms, std::function<void(int ec)> cb) {
  auto h = std::make_shared<boost::asio::deadline_timer>(ios);
  h->expires_from_now(boost::posix_time::milliseconds(ms));
  std::cout << "sleeping for " << " (time : " << ms << " ms)" << std::endl;
  h->async_wait([h, ms, cb](const boost::system::error_code& ec) {
    std::cout << "sleep is over : " << " (time : " << ms << " ms)" << std::endl;
    cb(!ec ? 0 : -1);
  });
}

int main() {
  using namespace std::chrono_literals;
  boost::asio::io_service io_service;
  std::unique_ptr<boost::asio::io_service::work> keepalive_work = std::make_unique<boost::asio::io_service::work>(io_service);
  std::thread asio_thread([&] { io_service.run(); });

  auto outer(std::make_shared<csi::async::work<int>>(csi::async::SEQUENTIAL, csi::async::FIRST_FAIL));

  for (auto steps : {1,2,3}) {

    auto inner(std::make_shared<csi::async::work<int>>(csi::async::PARALLEL, csi::async::FIRST_SUCCESS)); // this will be each 100ms delay
    for (auto delay : {5,4,3,2,1}) {
      inner->push_back([&io_service, delay](csi::async::work<int>::callback cb) {
        async_sleep(io_service, delay * 100, cb);
      });
    }
    outer->push_back(*inner);
  }

  std::cout << "start ASYNC" << std::endl;
  auto start_ts = std::chrono::steady_clock::now();
  auto res = (*outer)();
  auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start_ts);
  std::cout << "done ASYNC duration: " << duration.count() << " ms" << std::endl;

  std::cout << "resetting dummy work" << std::endl;
  keepalive_work.reset();
  asio_thread.join();
  std::cout << "exiting" << std::endl;
  return 0;
}

max_stars_repo_name: bitbouncer/csi-async

max_issues_repo_name: bitbouncer/csi-async

max_forks_repo_name: bitbouncer/csi-async


text: #include <boost/spirit/home/x3/support/traits/move_to.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
// ***********************************************************************
//                      DO NOT EDIT THIS FILE !!!
// ***********************************************************************
//   This file is automatically generated from the grib_api templates. All
//   changes will be overridden. If you want to do permanent changes then
//   you should write them into the 'AreaProcessedCrossSectionProductImpl.*' files.
// ***********************************************************************

#include "AreaProcessedCrossSectionProduct.h"
#include "../../common/GeneralDefinitions.h"
#include "../../common/GeneralFunctions.h"
#include <boost/functional/hash.hpp>
#include <iostream>
#include <macgyver/Exception.h>

namespace SmartMet {
namespace GRIB2 {

/*! \brief The constructor of the class. */

AreaProcessedCrossSectionProduct::AreaProcessedCrossSectionProduct() {
  try {
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The destructor of the class. */

AreaProcessedCrossSectionProduct::~AreaProcessedCrossSectionProduct() {
}

/*! \brief The method reads and initializes all data related to the current object.

        \param memoryReader  This object controls the access to the memory mapped file.
*/

void AreaProcessedCrossSectionProduct::read(MemoryReader &memoryReader) {
  try {
    mParameter.read(memoryReader);
    mHorizontalDimensionProcessed = memoryReader.read_UInt8_opt();
    mTreatmentOfMissingData = memoryReader.read_UInt8_opt();
    mTypeOfStatisticalProcessing = memoryReader.read_UInt8_opt();
    mStartOfRange = memoryReader.read_UInt32_opt();
    mEndOfRange = memoryReader.read_UInt32_opt();
    mNumberOfDataValues = memoryReader.read_UInt16_opt();
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method writes all data related to the current object.

        \param dataWriter  This object is used for writing the object data.
*/

void AreaProcessedCrossSectionProduct::write(DataWriter &dataWriter) {
  try {
    mParameter.write(dataWriter);
    dataWriter << mHorizontalDimensionProcessed;
    dataWriter << mTreatmentOfMissingData;
    dataWriter << mTypeOfStatisticalProcessing;
    dataWriter << mStartOfRange;
    dataWriter << mEndOfRange;
    dataWriter << mNumberOfDataValues;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method is used for collecting the current class attributeList.

    \param prefix         The prefix of the each attributeList parameter.
    \param attributeList  The attributeList storage.
*/

void AreaProcessedCrossSectionProduct::getAttributeList(const std::string &prefix, T::AttributeList &attributeList) const {
  try {
    char name[300];
    sprintf(name, "%sAreaProcessedCrossSectionProduct.", prefix.c_str());
    mParameter.getAttributeList(name, attributeList);
    sprintf(name, "%sAreaProcessedCrossSectionProduct.HorizontalDimensionProcessed", prefix.c_str());
    attributeList.addAttribute(name, toString(mHorizontalDimensionProcessed));
    sprintf(name, "%sAreaProcessedCrossSectionProduct.TreatmentOfMissingData", prefix.c_str());
    attributeList.addAttribute(name, toString(mTreatmentOfMissingData));
    sprintf(name, "%sAreaProcessedCrossSectionProduct.TypeOfStatisticalProcessing", prefix.c_str());
    attributeList.addAttribute(name, toString(mTypeOfStatisticalProcessing));
    sprintf(name, "%sAreaProcessedCrossSectionProduct.StartOfRange", prefix.c_str());
    attributeList.addAttribute(name, toString(mStartOfRange));
    sprintf(name, "%sAreaProcessedCrossSectionProduct.EndOfRange", prefix.c_str());
    attributeList.addAttribute(name, toString(mEndOfRange));
    sprintf(name, "%sAreaProcessedCrossSectionProduct.NumberOfDataValues", prefix.c_str());
    attributeList.addAttribute(name, toString(mNumberOfDataValues));
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method is used for getting attribute values by their names.

    \param attributeName  The name of the attribute.
    \param attributeValue The value of the attribute (string).
*/

bool AreaProcessedCrossSectionProduct::getAttributeValue(const char *attributeName, std::string &attributeValue) const {
  try {
    if (attributeName == nullptr)
      return false;
    if (mParameter.getAttributeValue(attributeName, attributeValue))
      return true;
    if (strcasecmp(attributeName, "HorizontalDimensionProcessed") == 0) {
      attributeValue = toString(mHorizontalDimensionProcessed);
      return true;
    }
    if (strcasecmp(attributeName, "TreatmentOfMissingData") == 0) {
      attributeValue = toString(mTreatmentOfMissingData);
      return true;
    }
    if (strcasecmp(attributeName, "TypeOfStatisticalProcessing") == 0) {
      attributeValue = toString(mTypeOfStatisticalProcessing);
      return true;
    }
    if (strcasecmp(attributeName, "StartOfRange") == 0) {
      attributeValue = toString(mStartOfRange);
      return true;
    }
    if (strcasecmp(attributeName, "EndOfRange") == 0) {
      attributeValue = toString(mEndOfRange);
      return true;
    }
    if (strcasecmp(attributeName, "NumberOfDataValues") == 0) {
      attributeValue = toString(mNumberOfDataValues);
      return true;
    }
    return false;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method is used for checking if the attribute value matches to the given value.

    \param attributeName  The name of the attribute.
    \param attributeValue The value of the attribute (string).
*/

bool AreaProcessedCrossSectionProduct::hasAttributeValue(const char *attributeName, const char *attributeValue) const {
  try {
    if (attributeName == nullptr || attributeValue == nullptr)
      return false;
    if (mParameter.hasAttributeValue(attributeName, attributeValue))
      return true;
    if (strcasecmp(attributeName, "HorizontalDimensionProcessed") == 0 && strcasecmp(attributeValue, toString(mHorizontalDimensionProcessed).c_str()) == 0)
      return true;
    if (strcasecmp(attributeName, "TreatmentOfMissingData") == 0 && strcasecmp(attributeValue, toString(mTreatmentOfMissingData).c_str()) == 0)
      return true;
    if (strcasecmp(attributeName, "TypeOfStatisticalProcessing") == 0 && strcasecmp(attributeValue, toString(mTypeOfStatisticalProcessing).c_str()) == 0)
      return true;
    if (strcasecmp(attributeName, "StartOfRange") == 0 && strcasecmp(attributeValue, toString(mStartOfRange).c_str()) == 0)
      return true;
    if (strcasecmp(attributeName, "EndOfRange") == 0 && strcasecmp(attributeValue, toString(mEndOfRange).c_str()) == 0)
      return true;
    if (strcasecmp(attributeName, "NumberOfDataValues") == 0 && strcasecmp(attributeValue, toString(mNumberOfDataValues).c_str()) == 0)
      return true;
    return false;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method prints the content of the current object into the given stream.

    \param ostream      The output stream.
    \param level        The print level (used when printing multi-level structures).
    \param optionFlags  The printing options expressed in flag-bits.
*/

void AreaProcessedCrossSectionProduct::print(std::ostream &stream, uint level, uint optionFlags) const {
  try {
    stream << space(level) << "AreaProcessedCrossSectionProduct\n";
    mParameter.print(stream, level + 1, optionFlags);
    stream << space(level) << "- HorizontalDimensionProcessed = " << toString(mHorizontalDimensionProcessed) << "\n";
    stream << space(level) << "- TreatmentOfMissingData = " << toString(mTreatmentOfMissingData) << "\n";
    stream << space(level) << "- TypeOfStatisticalProcessing = " << toString(mTypeOfStatisticalProcessing) << "\n";
    stream << space(level) << "- StartOfRange = " << toString(mStartOfRange) << "\n";
    stream << space(level) << "- EndOfRange = " << toString(mEndOfRange) << "\n";
    stream << space(level) << "- NumberOfDataValues = " << toString(mNumberOfDataValues) << "\n";
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method counts the hash value of the current object. */

T::Hash AreaProcessedCrossSectionProduct::countHash() const {
  try {
    std::size_t seed = 0;
    if (mHorizontalDimensionProcessed)
      boost::hash_combine(seed, *mHorizontalDimensionProcessed);
    if (mTreatmentOfMissingData)
      boost::hash_combine(seed, *mTreatmentOfMissingData);
    if (mTypeOfStatisticalProcessing)
      boost::hash_combine(seed, *mTypeOfStatisticalProcessing);
    if (mStartOfRange)
      boost::hash_combine(seed, *mStartOfRange);
    if (mEndOfRange)
      boost::hash_combine(seed, *mEndOfRange);
    if (mNumberOfDataValues)
      boost::hash_combine(seed, *mNumberOfDataValues);
    return seed;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method return the template number of the current class. */

uint AreaProcessedCrossSectionProduct::getTemplateNumber() const {
  return 1002;
}

ProductDefinition *AreaProcessedCrossSectionProduct::createProductDefinition() const {
  try {
    return static_cast<ProductDefinition *>(new AreaProcessedCrossSectionProduct(*this));
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the pointer to the {@link mParameter} attribute. */

ParameterSettings *AreaProcessedCrossSectionProduct::getParameter() const {
  try {
    return static_cast<ParameterSettings *>(&mParameter);
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mHorizontalDimensionProcessed} attribute. */

const T::UInt8_opt &AreaProcessedCrossSectionProduct::getHorizontalDimensionProcessed() const {
  try {
    return mHorizontalDimensionProcessed;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mTreatmentOfMissingData} attribute. */

const T::UInt8_opt &AreaProcessedCrossSectionProduct::getTreatmentOfMissingData() const {
  try {
    return mTreatmentOfMissingData;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mTypeOfStatisticalProcessing} attribute. */

const T::UInt8_opt &AreaProcessedCrossSectionProduct::getTypeOfStatisticalProcessing() const {
  try {
    return mTypeOfStatisticalProcessing;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mStartOfRange} attribute. */

const T::UInt32_opt &AreaProcessedCrossSectionProduct::getStartOfRange() const {
  try {
    return mStartOfRange;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mEndOfRange} attribute. */

const T::UInt32_opt &AreaProcessedCrossSectionProduct::getEndOfRange() const {
  try {
    return mEndOfRange;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

/*! \brief The method returns the value of the {@link mNumberOfDataValues} attribute. */

const T::UInt16_opt &AreaProcessedCrossSectionProduct::getNumberOfDataValues() const {
  try {
    return mNumberOfDataValues;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setParameter(ParameterSettings &parameter) {
  try {
    mParameter = parameter;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setHorizontalDimensionProcessed(T::UInt8_opt horizontalDimensionProcessed) {
  try {
    mHorizontalDimensionProcessed = horizontalDimensionProcessed;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setTreatmentOfMissingData(T::UInt8_opt treatmentOfMissingData) {
  try {
    mTreatmentOfMissingData = treatmentOfMissingData;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setTypeOfStatisticalProcessing(T::UInt8_opt typeOfStatisticalProcessing) {
  try {
    mTypeOfStatisticalProcessing = typeOfStatisticalProcessing;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setStartOfRange(T::UInt32_opt startOfRange) {
  try {
    mStartOfRange = startOfRange;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setEndOfRange(T::UInt32_opt endOfRange) {
  try {
    mEndOfRange = endOfRange;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

void AreaProcessedCrossSectionProduct::setNumberOfDataValues(T::UInt16_opt numberOfDataValues) {
  try {
    mNumberOfDataValues = numberOfDataValues;
  } catch (...) {
    throw Fmi::Exception(BCP, "Operation failed", nullptr);
  }
}

} // namespace GRIB2
} // namespace SmartMet

max_stars_repo_name: fmidev/smartmet-library-grid-files

max_issues_repo_name: fmidev/smartmet-library-grid-files

max_forks_repo_name: fmidev/smartmet-library-grid-files


text:
/* Copyright 2003-2015 Joaquin M Lopez Munoz.
 * 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/libs/multi_index for library home page.
 */

#ifndef BOOST_MULTI_INDEX_DETAIL_INDEX_MATCHER_HPP
#define BOOST_MULTI_INDEX_DETAIL_INDEX_MATCHER_HPP

#if defined(_MSC_VER)
#pragma once
#endif

#include <algorithm>
#include <boost/config.hpp> /* keep it first to prevent nasty warns in MSVC */
#include <boost/multi_index/detail/auto_space.hpp>
#include <boost/multi_index/detail/raw_ptr.hpp>
#include <boost/noncopyable.hpp>
#include <cstddef>
#include <functional>

namespace boost
{

namespace multi_index
{

namespace detail
{

/* index_matcher compares a sequence of elements against a
 * base sequence, identifying those elements that belong to the
 * longest subsequence which is ordered with respect to the base.
 * For instance, if the base sequence is:
 *
 *   0 1 2 3 4 5 6 7 8 9
 *
 * and the compared sequence (not necesarilly the same length):
 *
 *   1 4 2 3 0 7 8 9
 *
 * the elements of the longest ordered subsequence are:
 *
 *   1 2 3 7 8 9
 *
 * The algorithm for obtaining such a subsequence is called
 * Patience Sorting, described in ch. 1 of:
 *   Aldous, D., Diaconis, P.: "Longest increasing subsequences: from
 *   patience sorting to the Baik-Deift-Johansson Theorem", Bulletin
 *   of the American Mathematical Society, vol. 36, no 4, pp. 413-432,
 *   July 1999.
 *   http://www.ams.org/bull/1999-36-04/S0273-0979-99-00796-X/
 *   S0273-0979-99-00796-X.pdf
 *
 * This implementation is not fully generic since it assumes that
 * the sequences given are pointed to by index iterators (having a
 * get_node() memfun.)
 */

namespace index_matcher
{

/* The algorithm stores the nodes of the base sequence and a number
 * of "piles" that are dynamically updated during the calculation
 * stage. From a logical point of view, nodes form an independent
 * sequence from piles. They are stored together so as to minimize
 * allocated memory.
 */

struct entry
{
    entry(void* node_, std::size_t pos_ = 0) : node(node_), pos(pos_)
    {
    }

    /* node stuff */

    void* node;
    std::size_t pos;
    entry* previous;
    bool ordered;

    struct less_by_node
    {
        bool operator()(const entry& x, const entry& y) const
        {
            return std::less<void*>()(x.node, y.node);
        }
    };

    /* pile stuff */

    std::size_t pile_top;
    entry* pile_top_entry;

    struct less_by_pile_top
    {
        bool operator()(const entry& x, const entry& y) const
        {
            return x.pile_top < y.pile_top;
        }
    };
};

/* common code operating on void *'s */

template <typename Allocator>
class algorithm_base : private noncopyable
{
  protected:
    algorithm_base(const Allocator& al, std::size_t size) :
        spc(al, size), size_(size), n_(0), sorted(false)
    {
    }

    void add(void* node)
    {
        entries()[n_] = entry(node, n_);
        ++n_;
    }

    void begin_algorithm() const
    {
        if (!sorted)
        {
            std::sort(entries(), entries() + size_, entry::less_by_node());
            sorted = true;
        }
        num_piles = 0;
    }

    void add_node_to_algorithm(void* node) const
    {
        entry* ent =
            std::lower_bound(entries(), entries() + size_, entry(node),
                             entry::less_by_node()); /* localize entry */
        ent->ordered = false;
        std::size_t n = ent->pos; /* get its position */

        entry dummy(0);
        dummy.pile_top = n;

        entry* pile_ent =    /* find the first available pile */
            std::lower_bound(/* to stack the entry            */
                             entries(), entries() + num_piles, dummy,
                             entry::less_by_pile_top());

        pile_ent->pile_top = n; /* stack the entry */
        pile_ent->pile_top_entry = ent;

        /* if not the first pile, link entry to top of the preceding pile */
        if (pile_ent > &entries()[0])
        {
            ent->previous = (pile_ent - 1)->pile_top_entry;
        }

        if (pile_ent == &entries()[num_piles])
        { /* new pile? */
            ++num_piles;
        }
    }

    void finish_algorithm() const
    {
        if (num_piles > 0)
        {
            /* Mark those elements which are in their correct position, i.e.
             * those belonging to the longest increasing subsequence. These are
             * those elements linked from the top of the last pile.
             */

            entry* ent = entries()[num_piles - 1].pile_top_entry;
            for (std::size_t n = num_piles; n--;)
            {
                ent->ordered = true;
                ent = ent->previous;
            }
        }
    }

    bool is_ordered(void* node) const
    {
        return std::lower_bound(entries(), entries() + size_, entry(node),
                                entry::less_by_node())
            ->ordered;
    }

  private:
    entry* entries() const
    {
        return raw_ptr<entry*>(spc.data());
    }

    auto_space<entry, Allocator> spc;
    std::size_t size_;
    std::size_t n_;
    mutable bool sorted;
    mutable std::size_t num_piles;
};

/* The algorithm has three phases:
 *   - Initialization, during which the nodes of the base sequence are added.
 *   - Execution.
 *   - Results querying, through the is_ordered memfun.
 */

template <typename Node, typename Allocator>
class algorithm : private algorithm_base<Allocator>
{
    typedef algorithm_base<Allocator> super;

  public:
    algorithm(const Allocator& al, std::size_t size) : super(al, size)
    {
    }

    void add(Node* node)
    {
        super::add(node);
    }

    template <typename IndexIterator>
    void execute(IndexIterator first, IndexIterator last) const
    {
        super::begin_algorithm();

        for (IndexIterator it = first; it != last; ++it)
        {
            add_node_to_algorithm(get_node(it));
        }

        super::finish_algorithm();
    }

    bool is_ordered(Node* node) const
    {
        return super::is_ordered(node);
    }

  private:
    void add_node_to_algorithm(Node* node) const
    {
        super::add_node_to_algorithm(node);
    }

    template <typename IndexIterator>
    static Node* get_node(IndexIterator it)
    {
        return static_cast<Node*>(it.get_node());
    }
};

} // namespace index_matcher

} // namespace detail

} /* namespace multi_index */

} /* namespace boost */

#endif

max_stars_repo_name: sotaoverride/backup

max_issues_repo_name: sotaoverride/backup

max_forks_repo_name: sotaoverride/backup


text:
#include <vector>
#include "boost/random.hpp"
#include "boost/generator_iterator.hpp"
#include <boost/random/normal_distribution.hpp>
#include <algorithm>

typedef boost::mt19937 RNGType; ///< mersenne twister generator

int main() {
    RNGType rng;
    boost::normal_distribution<> rdist(1.0,0.5); /**< normal distribution
                           with mean of 1.0 and standard deviation of 0.5 */

    boost::variate_generator< RNGType, boost::normal_distribution<> >
                    get_rand(rng, rdist);

    std::vector<double> v(1000);
    generate(v.begin(),v.end(),get_rand);
    return 0;
}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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 CRYPTSIMULATION2D_HPP_
#define CRYPTSIMULATION2D_HPP_

#include "ChasteSerialization.hpp"
#include <boost/serialization/base_object.hpp>

#include "WntConcentration.hpp"
#include "OffLatticeSimulation.hpp"
#include "VertexBasedCellPopulation.hpp"
#include "CryptSimulationBoundaryCondition.hpp"
#include "CryptCentreBasedDivisionRule.hpp"
#include "CryptVertexBasedDivisionRule.hpp"

/**
 * A 2D crypt simulation object. For more details on the crypt geometry, see the
 * papers by van Leeuwen et al (2009) [doi:10.1111/j.1365-2184.2009.00627.x] and
 * Osborne et al (2010) [doi:10.1098/rsta.2010.0173].
 */
class CryptSimulation2d : public OffLatticeSimulation<2>
{
    // Allow tests to access private members, in order to test computation of private functions e.g. DoCellBirth()
    friend class TestCryptSimulation2dWithMeshBasedCellPopulation;
    friend class TestCryptSimulation2dWithVertexBasedCellPopulation;

protected:

    /** Needed for serialization. */
    friend class boost::serialization::access;
    /**
     * Archive the simulation and member variable.
     *
     * @param archive the archive
     * @param version the current version of this class
     */
    template<class Archive>
    void serialize(Archive & archive, const unsigned int version)
    {
        archive & boost::serialization::base_object<OffLatticeSimulation<2> >(*this);

        SerializableSingleton<WntConcentration<2> >* p_wnt_wrapper = WntConcentration<2>::Instance()->GetSerializationWrapper();
        archive & p_wnt_wrapper;
    }

    /**
     * Helper member that stores whether we are using a MeshBasedCellPopulationWithGhostNodes
     * (if not, then we are assumed to be using a VertexBasedCellPopulation).
     */
    bool mUsingMeshBasedCellPopulation;

    /**
     * Overridden SetupSolve() method.
     *
     * Write initial beta catenin results to file if required.
     */
    void SetupSolve();

public:

    /**
     * Constructor.
     *
     * @param rCellPopulation A cell population object
     * @param deleteCellPopulationInDestructor Whether to delete the cell population on destruction to
     *     free up memory (defaults to false)
     * @param initialiseCells whether to initialise cells (defaults to true, set to false when loading from an archive)
     */
    CryptSimulation2d(AbstractCellPopulation<2>& rCellPopulation,
                      bool deleteCellPopulationInDestructor=false,
                      bool initialiseCells=true);

    /**
     * Destructor.
     *
     * This frees the CryptSimulationBoundaryCondition.
     */
    virtual ~CryptSimulation2d();

    /**
     * Set method for mUseJiggledBottomCells.
     */
    void UseJiggledBottomCells();

    /**
     * Sets the Ancestor index of all the cells at the bottom in order,
     * can be used to trace clonal populations.
     */
    void SetBottomCellAncestors();

    /**
     * Outputs simulation parameters to file
     *
     * As this method is pure virtual, it must be overridden
     * in subclasses.
     *
     * @param rParamsFile the file stream to which the parameters are output
     */
    void OutputSimulationParameters(out_stream& rParamsFile);
};

// Declare identifier for the serializer
#include "SerializationExportWrapper.hpp"
CHASTE_CLASS_EXPORT(CryptSimulation2d)

namespace boost
{
namespace serialization
{
/**
 * Serialize information required to construct a CryptSimulation2d.
 */
template<class Archive>
inline void save_construct_data(
    Archive & ar, const CryptSimulation2d * t, const unsigned int file_version)
{
    // Save data required to construct instance
    const AbstractCellPopulation<2>* p_cell_population = &(t->rGetCellPopulation());
    ar & p_cell_population;
}

/**
 * De-serialize constructor parameters and initialise a CryptSimulation2d.
 */
template<class Archive>
inline void load_construct_data(
    Archive & ar, CryptSimulation2d * t, const unsigned int file_version)
{
    // Retrieve data from archive required to construct new instance
    AbstractCellPopulation<2>* p_cell_population;
    ar & p_cell_population;

    // Invoke inplace constructor to initialise instance
    ::new(t)CryptSimulation2d(*p_cell_population, true, false);
}
}
} // namespace

#endif /*CRYPTSIMULATION2D_HPP_*/

max_stars_repo_name: stu-l/Chaste

max_issues_repo_name: stu-l/Chaste

max_forks_repo_name: stu-l/Chaste


text:
// ----------------------------------------------------------------------------
// Copyright (C) 2002-2006 Marcin Kalicinski
// Copyright (C) 2009 Sebastian Redl
//
// 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)
//
// For more information, see www.boost.org
// ----------------------------------------------------------------------------
#ifndef BOOST_PROPERTY_TREE_INI_PARSER_HPP_INCLUDED
#define BOOST_PROPERTY_TREE_INI_PARSER_HPP_INCLUDED

#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/detail/ptree_utils.hpp>
#include <boost/property_tree/detail/file_parser_error.hpp>
#include <fstream>
#include <string>
#include <sstream>
#include <stdexcept>
#include <locale>

namespace boost { namespace property_tree { namespace ini_parser
{

    /**
     * Determines whether the @c flags are valid for use with the ini_parser.
     * @param flags value to check for validity as flags to ini_parser.
     * @return true if the flags are valid, false otherwise.
     */
    inline bool validate_flags(int flags)
    {
        return flags == 0;
    }

    /** Indicates an error parsing INI formatted data. */
    class ini_parser_error: public file_parser_error
    {
    public:
        /**
         * Construct an @c ini_parser_error
         * @param message Message describing the parser error.
         * @param filename The name of the file being parsed containing the
         *                 error.
         * @param line The line in the given file where an error was
         *             encountered.
         */
        ini_parser_error(const std::string &message,
                         const std::string &filename,
                         unsigned long line)
            : file_parser_error(message, filename, line)
        {
        }
    };

    /**
     * Read INI from a the given stream and translate it to a property tree.
     * @note Clears existing contents of property tree. In case of error
     *       the property tree is not modified.
     * @throw ini_parser_error If a format violation is found.
     * @param stream Stream from which to read in the property tree.
     * @param[out] pt The property tree to populate.
     */
    template<class Ptree>
    void read_ini(std::basic_istream<
                    typename Ptree::key_type::value_type> &stream,
                  Ptree &pt)
    {
        typedef typename Ptree::key_type::value_type Ch;
        typedef std::basic_string<Ch> Str;
        const Ch semicolon = stream.widen(';');
        const Ch hash = stream.widen('#');
        const Ch lbracket = stream.widen('[');
        const Ch rbracket = stream.widen(']');

        Ptree local;
        unsigned long line_no = 0;
        Ptree *section = 0;
        Str line;

        // For all lines
        while (stream.good())
        {

            // Get line from stream
            ++line_no;
            std::getline(stream, line);
            if (!stream.good() && !stream.eof())
                BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                    "read error", "", line_no));

            // If line is non-empty
            line = property_tree::detail::trim(line, stream.getloc());
            if (!line.empty())
            {
                // Comment, section or key?
                if (line[0] == semicolon || line[0] == hash)
                {
                    // Ignore comments
                }
                else if (line[0] == lbracket)
                {
                    // If the previous section was empty, drop it again.
                    if (section && section->empty())
                        local.pop_back();
                    typename Str::size_type end = line.find(rbracket);
                    if (end == Str::npos)
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "unmatched '['", "", line_no));
                    Str key = property_tree::detail::trim(
                        line.substr(1, end - 1), stream.getloc());
                    if (local.find(key) != local.not_found())
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "duplicate section name", "", line_no));
                    section = &local.push_back(
                        std::make_pair(key, Ptree()))->second;
                }
                else
                {
                    Ptree &container = section ? *section : local;
                    typename Str::size_type eqpos = line.find(Ch('='));
                    if (eqpos == Str::npos)
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "'=' character not found in line", "", line_no));
                    if (eqpos == 0)
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "key expected", "", line_no));
                    Str key = property_tree::detail::trim(
                        line.substr(0, eqpos), stream.getloc());
                    Str data = property_tree::detail::trim(
                        line.substr(eqpos + 1, Str::npos), stream.getloc());
                    if (container.find(key) != container.not_found())
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "duplicate key name", "", line_no));
                    container.push_back(std::make_pair(key, Ptree(data)));
                }
            }
        }
        // If the last section was empty, drop it again.
        if (section && section->empty())
            local.pop_back();

        // Swap local ptree with result ptree
        pt.swap(local);

    }

    /**
     * Read INI from a the given file and translate it to a property tree.
     * @note Clears existing contents of property tree.  In case of error the
     *       property tree unmodified.
     * @throw ini_parser_error In case of error deserializing the property tree.
     * @param filename Name of file from which to read in the property tree.
     * @param[out] pt The property tree to populate.
     * @param loc The locale to use when reading in the file contents.
     */
    template<class Ptree>
    void read_ini(const std::string &filename,
                  Ptree &pt,
                  const std::locale &loc = std::locale())
    {
        std::basic_ifstream<typename Ptree::key_type::value_type>
            stream(filename.c_str());
        if (!stream)
            BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                "cannot open file", filename, 0));
        stream.imbue(loc);
        try {
            read_ini(stream, pt);
        }
        catch (ini_parser_error &e) {
            BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                e.message(), filename, e.line()));
        }
    }

    namespace detail
    {
        template<class Ptree>
        void check_dupes(const Ptree &pt)
        {
            if(pt.size() <= 1)
                return;
            const typename Ptree::key_type *lastkey = 0;
            typename Ptree::const_assoc_iterator it = pt.ordered_begin(),
                                                 end = pt.not_found();
            lastkey = &it->first;
            for(++it; it != end; ++it) {
                if(*lastkey == it->first)
                    BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                        "duplicate key", "", 0));
                lastkey = &it->first;
            }
        }

        template <typename Ptree>
        void write_keys(std::basic_ostream<
                                      typename Ptree::key_type::value_type
                                  > &stream,
                                  const Ptree& pt,
                                  bool throw_on_children)
        {
            typedef typename Ptree::key_type::value_type Ch;
            for (typename Ptree::const_iterator it = pt.begin(), end = pt.end();
                 it != end; ++it)
            {
                if (!it->second.empty()) {
                    if (throw_on_children) {
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "ptree is too deep", "", 0));
                    }
                    continue;
                }
                stream << it->first << Ch('=')
                    << it->second.template get_value<
                        std::basic_string<Ch> >()
                    << Ch('\n');
            }
        }

        template <typename Ptree>
        void write_top_level_keys(std::basic_ostream<
                                      typename Ptree::key_type::value_type
                                  > &stream,
                                  const Ptree& pt)
        {
            write_keys(stream, pt, false);
        }

        template <typename Ptree>
        void write_sections(std::basic_ostream<
                                typename Ptree::key_type::value_type
                            > &stream,
                            const Ptree& pt)
        {
            typedef typename Ptree::key_type::value_type Ch;
            for (typename Ptree::const_iterator it = pt.begin(), end = pt.end();
                 it != end; ++it)
            {
                if (!it->second.empty()) {
                    check_dupes(it->second);
                    if (!it->second.data().empty())
                        BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                            "mixed data and children", "", 0));
                    stream << Ch('[') << it->first << Ch(']') << Ch('\n');
                    write_keys(stream, it->second, true);
                }
            }
        }
    }

    /**
     * Translates the property tree to INI and writes it the given output
     * stream.
     * @pre @e pt cannot have data in its root.
     * @pre @e pt cannot have keys both data and children.
     * @pre @e pt cannot be deeper than two levels.
     * @pre There cannot be duplicate keys on any given level of @e pt.
     * @throw ini_parser_error In case of error translating the property tree to
     *                         INI or writing to the output stream.
     * @param stream The stream to which to write the INI representation of the
     *               property tree.
     * @param pt The property tree to tranlsate to INI and output.
     * @param flags The flags to use when writing the INI file.
     *              No flags are currently supported.
     */
    template<class Ptree>
    void write_ini(std::basic_ostream<
                       typename Ptree::key_type::value_type
                   > &stream,
                   const Ptree &pt,
                   int flags = 0)
    {
        BOOST_ASSERT(validate_flags(flags));
        (void)flags;

        if (!pt.data().empty())
            BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                "ptree has data on root", "", 0));
        detail::check_dupes(pt);

        detail::write_top_level_keys(stream, pt);
        detail::write_sections(stream, pt);
    }

    /**
     * Translates the property tree to INI and writes it the given file.
     * @pre @e pt cannot have data in its root.
     * @pre @e pt cannot have keys both data and children.
     * @pre @e pt cannot be deeper than two levels.
     * @pre There cannot be duplicate keys on any given level of @e pt.
     * @throw info_parser_error In case of error translating the property tree
     *                          to INI or writing to the file.
     * @param filename The name of the file to which to write the INI
     *                 representation of the property tree.
     * @param pt The property tree to tranlsate to INI and output.
     * @param flags The flags to use when writing the INI file.
     *              The following flags are supported:
     * @li @c skip_ini_validity_check -- Skip check if ptree is a valid ini. The
     *     validity check covers the preconditions but takes <tt>O(n log n)</tt>
     *     time.
     * @param loc The locale to use when writing the file.
     */
    template<class Ptree>
    void write_ini(const std::string &filename,
                   const Ptree &pt,
                   int flags = 0,
                   const std::locale &loc = std::locale())
    {
        std::basic_ofstream<typename Ptree::key_type::value_type>
            stream(filename.c_str());
        if (!stream)
            BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                "cannot open file", filename, 0));
        stream.imbue(loc);
        try {
            write_ini(stream, pt, flags);
        }
        catch (ini_parser_error &e) {
            BOOST_PROPERTY_TREE_THROW(ini_parser_error(
                e.message(), filename, e.line()));
        }
    }

} } }

namespace boost { namespace property_tree
{
    using ini_parser::ini_parser_error;
    using ini_parser::read_ini;
    using ini_parser::write_ini;
} }

#endif

max_stars_repo_name: rudylee/expo

max_issues_repo_name: rudylee/expo

max_forks_repo_name: rudylee/expo


text:
/*
For more information, please see: http://software.sci.utah.edu

The MIT License

Copyright (c) 2012 Scientific Computing and Imaging Institute,
University of Utah.

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/

//    File       : SolveInverseProblemWithTikhonov.cc
//    Author     : Moritz Dannhauer, Ayla Khan, Dan White
//    Date       : November 02th, 2012 (last update)

#include <boost/bind.hpp>
#include <boost/lexical_cast.hpp>

#include <Modules/Legacy/Inverse/SolveInverseProblemWithTikhonovImpl.h>

#include <Core/Datatypes/Matrix.h>
#include <Core/Datatypes/DenseMatrix.h>
#include <Core/Datatypes/DenseColumnMatrix.h>
#include <Core/Datatypes/SparseRowMatrix.h>
#include <Core/Datatypes/MatrixTypeConversions.h>

#include <Core/Logging/LoggerInterface.h>
#include <Core/Utils/Exception.h>

namespace BioPSE
{

using namespace SCIRun;
using namespace SCIRun::Core;
using namespace SCIRun::Core::Datatypes;
using namespace SCIRun::Core::Logging;

TikhonovAlgorithmImpl::TikhonovAlgorithmImpl(const DenseMatrixHandle& forwardMatrix,
                                             const DenseMatrixHandle& measuredData,
                                             AlgorithmChoice regularizationChoice,
                                             AlgorithmSolutionSubcase regularizationSolutionSubcase,
                                             AlgorithmResidualSubcase regularizationResidualSubcase,
                                             const DenseMatrixHandle sourceWeighting,
                                             const DenseMatrixHandle sensorWeighting,
                                             bool computeRegularizedInverse,
                                             LegacyLoggerInterface* pr)
: forwardMatrix_(forwardMatrix),
measuredData_(measuredData),
sourceWeighting_(sourceWeighting),
sensorWeighting_(sensorWeighting),
regularizationChoice_(regularizationChoice),
regularizationSolutionSubcase_(regularizationSolutionSubcase),
regularizationResidualSubcase_(regularizationResidualSubcase),
lambda_(0),
computeRegularizedInverse_(computeRegularizedInverse),
pr_(pr)
{
  //TODO: size checking here.
}

MatrixHandle TikhonovAlgorithmImpl::get_inverse_solution() const
{
  return inverseSolution_;
}

MatrixHandle TikhonovAlgorithmImpl::get_inverse_matrix() const
{
  return inverseMatrix_;
}

DenseColumnMatrixHandle TikhonovAlgorithmImpl::get_regularization_parameter() const
{
  return regularizationParameter_;
}

//! Find Corner, find the maximal curvature which corresponds to the L-curve corner
double
TikhonovAlgorithmImpl::FindCorner(const TikhonovAlgorithm::LCurveInput& input, int& lambda_index)
{
  const std::vector<double>& rho = input.rho_;
  const std::vector<double>& eta = input.eta_;
  const std::vector<double>& lambdaArray = input.lambdaArray_;
  int nLambda = input.nLambda_;

  std::vector<double> deta(nLambda);
  std::vector<double> ddeta(nLambda);
  std::vector<double> drho(nLambda);
  std::vector<double> ddrho(nLambda);
  std::vector<double> lrho(nLambda);
  std::vector<double> leta(nLambda);
  DenseColumnMatrix kapa(nLambda);

  double maxKapa = -1.0e10;
  for (int i = 0; i < nLambda; i++)
  {
    lrho[i] = std::log10(rho[i]);
    leta[i] = std::log10(eta[i]);
    if(i>0)
    {
      deta[i] = (leta[i]-leta[i-1]) / (lambdaArray[i]-lambdaArray[i-1]); // compute first derivative
      drho[i] = (lrho[i]-lrho[i-1]) / (lambdaArray[i]-lambdaArray[i-1]);
    }
    if(i>1)
    {
      ddeta[i] = (deta[i]-deta[i-1]) / (lambdaArray[i]-lambdaArray[i-1]); // compute second derivative from first
      ddrho[i] = (drho[i]-drho[i-1]) / (lambdaArray[i]-lambdaArray[i-1]);
    }
  }
  drho[0] = drho[1];
  deta[0] = deta[1];
  ddrho[0] = ddrho[2];
  ddrho[1] = ddrho[2];
  ddeta[0] = ddeta[2];
  ddeta[1] = ddeta[2];

  lambda_index = 0;
  for (int i = 0; i < nLambda; i++)
  {
    kapa[i] = std::abs((drho[i] * ddeta[i] - ddrho[i] * deta[i]) /  //compute curvature
                       std::sqrt(std::pow((deta[i]*deta[i]+drho[i]*drho[i]), 3.0)));
    if (kapa[i] > maxKapa) // find max curvature
    {
      maxKapa = kapa[i];
      lambda_index = i;
    }
  }

  return lambdaArray[lambda_index];
}

double
TikhonovAlgorithmImpl::LambdaLookup(const TikhonovAlgorithm::LCurveInput& input, double lambda, int& lambda_index, const double epsilon)
{
  const std::vector<double>& lambdaArray = input.lambdaArray_;
  int nLambda = input.nLambda_;

  for (int i = 0; i < nLambda-1; ++i)
  {
    if (i > 0 && (lambda < lambdaArray[i-1] || lambda > lambdaArray[i+1])) continue;

    double lambda_step_midpoint = std::abs(lambdaArray[i+1] - lambdaArray[i])/2;

    if (std::abs(lambda - lambdaArray[i]) <= epsilon)  // TODO: is this a reasonable comparison???
    {
      lambda_index = i;
      return lambdaArray[lambda_index];
    }

    if (std::abs(lambda - lambdaArray[i]) < lambda_step_midpoint)
    {
      lambda_index = i;
      return lambdaArray[lambda_index];
    }

    if (std::abs(lambda - lambdaArray[i+1]) < lambda_step_midpoint)
    {
      lambda_index = i+1;
      return lambdaArray[lambda_index];
    }
  }
  return -1;
}

namespace LinearAlgebra
{
  void solve_lapack(const DenseMatrix& A, const DenseColumnMatrix& b, DenseColumnMatrix& x)
  {
    x = A.lu().solve(b).eval();
  }

  class LapackError : public std::exception {};
}

void TikhonovAlgorithmImpl::run(const TikhonovAlgorithmImpl::Input& input)
{
  // TODO: use DimensionMismatch exception where appropriate
  // DIMENSION CHECK!!
  const int M = forwardMatrix_->nrows();
  const int N = forwardMatrix_->ncols();
  if (M != measuredData_->nrows())
  {
    BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Input matrix dimensions must agree."));
  }
  if (1 != measuredData_->ncols())
  {
    BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Measured data must be a vector"));
  }
  //decide used Tikhonov regularization formulation (underdetermined or overdetermined)
  //based purely on relationship of number of sensors compared to number of source reconstruction points
  //UNDERDETERMINED CASE
  if ( ((M < N) && (regularizationChoice_ == automatic)) || (regularizationChoice_ == underdetermined))
  {
    //.........................................................................
    // OPERATE ON DATA:
    // Compute X = R * R^T * A^T (A * R * R^T * A^T + LAMBDA * LAMBDA * C * C^T) * Y
    //.........................................................................
    DenseMatrix ARRtrAtr, RRtrAtr;
    double lambda=0, lambda_sq=0;
    DenseMatrix forward_transpose = forwardMatrix_->transpose();
    DenseMatrix regMat;
    #ifdef SCIRUN4_CODE_TO_BE_CONVERTED_LATER
    SparseRowMatrixHandle sourceWeighting_sparse;
    #endif
    DenseColumnMatrix solution(M);
    DenseColumnMatrix RRtrAtrsolution(N);

    const DenseMatrixHandle measuredDataRef = measuredData_;
    if (!sourceWeighting_) //check if a Source Space Weighting Matrix exist?
    {
      regMat = DenseMatrix::Identity(N, N);
      ARRtrAtr = *forwardMatrix_ * forward_transpose;
      RRtrAtr = forward_transpose;
    }
    else
    {
      if(((regularizationSolutionSubcase_==solution_constrained) && ((N != sourceWeighting_->nrows()) || (N != sourceWeighting_->ncols()))) ||
         ((regularizationSolutionSubcase_==solution_constrained_squared) && (N != sourceWeighting_->nrows())))
      {
        BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Solution Weighting Matrix (number of rows and columns) must fit matrix dimensions of Forward Matrix (number of columns) !"));
      }
      regMat = *sourceWeighting_;
      #ifdef SCIRUN4_CODE_TO_BE_CONVERTED_LATER
      sourceWeighting_sparse = sourceWeighting_->sparse();
      #endif

      if (regularizationSolutionSubcase_== solution_constrained)
      {
        RRtrAtr = *sourceWeighting_ * forward_transpose;
	      auto AR = forwardMatrix_;
	      ARRtrAtr = *AR * RRtrAtr;
      }
      else if (regularizationSolutionSubcase_== solution_constrained_squared)
      {
        auto AR = *forwardMatrix_ * *sourceWeighting_;
	      auto RtrAtr = AR.transpose().eval();
        RRtrAtr = *sourceWeighting_ * RtrAtr;
	      ARRtrAtr = AR * RtrAtr;
      }
    }
    DenseMatrix CCtr;
    if (!sensorWeighting_) //check if a Sensor Space (Noise Covariance) Weighting Matrix exist?
    {
      CCtr = DenseMatrix::Identity(M, M);
    }
    else
    {
      if (((regularizationResidualSubcase_ == residual_constrained) && ((M != sensorWeighting_->nrows()) || (M != sensorWeighting_->ncols()))) ||
          ((regularizationResidualSubcase_ == residual_constrained_squared) && (M != sensorWeighting_->nrows())))
      {
        BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Data Residual Weighting Matrix (number of rows and columns) must fit matrix dimensions of Forward Matrix (number of rows) !"));
      }
      else
      {
        if (regularizationResidualSubcase_ == residual_constrained)
          CCtr = *sensorWeighting_;
        else
          if ( regularizationResidualSubcase_ == residual_constrained_squared )
            CCtr = *sensorWeighting_ * sensorWeighting_->transpose();
      }
    }

    DenseMatrix regForMatrix(M, M);

    //Get Regularization parameter(s) : Lambda
    if ((input.regMethod_ == "single") || (input.regMethod_ == "slider"))
    {
      if (input.regMethod_ == "single")
      {
        // Use single fixed lambda value, entered in UI
        lambda = input.lambdaFromTextEntry_;
      }
      else if (input.regMethod_ == "slider")
      {
        // Use single fixed lambda value, select via slider
        lambda = input.lambdaSlider_;
      }
    }
    else if (input.regMethod_ == "lcurve")
    {
      const int nLambda = input.lambdaCount_;  //declare needed variables

      std::vector<double> lambdaArray(nLambda, 0.0);
      std::vector<double> rho(nLambda, 0.0);
      std::vector<double> eta(nLambda, 0.0);

      lambdaArray[0] = input.lambdaMin_;
      const double lam_step = pow(10.0, log10(input.lambdaMax_ / input.lambdaMin_) / (nLambda-1));

      DenseColumnMatrix Ax(M);
      DenseColumnMatrix Rx(N);

      int lambda_index = 0;
      const DenseMatrix &forwardMatrixRef = *forwardMatrix_;

      int nr_rows = CCtr.nrows();
      int nr_cols = CCtr.ncols();
      const int NR_BOUNDS_CHECK = nr_rows * nr_cols;

      double* AtrA = ARRtrAtr.data();
      double* LLtr = CCtr.data();

      for (int j = 0; j < nLambda; j++)
      {
        if (j)
        {
          lambdaArray[j] = lambdaArray[j-1] * lam_step;
        }

        lambda_sq = lambdaArray[j] * lambdaArray[j]; //generate current lambda

        double* rm = regForMatrix.data();

        for (int i = 0; i < NR_BOUNDS_CHECK; i++)
        {
          rm[i] = AtrA[i] + lambda_sq * LLtr[i];
        }

        try
        {
          LinearAlgebra::solve_lapack(regForMatrix, *measuredDataRef, solution);
        }
        catch (LinearAlgebra::LapackError&)
        {
          const std::string errorMessage("The Tikhonov linear system could not be solved for a regularization parameter in the Lambda Range of the L-curve. Use a higher Lambda Range ''From'' value for the L-Curve calculation.");
          if (pr_)
          {
            pr_->error(errorMessage);
          }
          else
          {
            std::cerr << errorMessage << std::endl;
          }
          throw;
        }
        catch(DimensionMismatch&)
        {
          const std::string errorMessage("Invalid matrix sizes are being used in the Tikhonov linear system.");
          if (pr_)
          {
            pr_->error(errorMessage);
          }
          else
          {
            std::cerr << errorMessage << std::endl;
          }
          throw;
        }

        RRtrAtrsolution = RRtrAtr * solution;

        if (sourceWeighting_)
        {
          if (RRtrAtrsolution.nrows() == sourceWeighting_->ncols())
            Rx = *sourceWeighting_ * RRtrAtrsolution;
          else if  (RRtrAtrsolution.nrows() == sourceWeighting_->nrows())
            Rx = sourceWeighting_->transpose() * RRtrAtrsolution;
          else
          {
            const std::string errorMessage(" Solution weighting matrix unexpectedly does not fit to compute the weighted solution norm. ");
            if (pr_)
            {
              pr_->error(errorMessage);
            }
            else
            {
              std::cerr << errorMessage << std::endl;
            }
          }
        }
        else
          Rx = RRtrAtrsolution;

        // Calculate the norm of Ax-b and Rx for L curve
        Ax = forwardMatrixRef * RRtrAtrsolution;

        rho[j]=0; eta[j]=0;
        for (int k = 0; k < Ax.nrows(); k++)
        {
          double T = Ax(k) - (*measuredDataRef)(k);
          rho[j] += T*T; //norm of the data fit term
        }

        for (int k = 0; k < Rx.nrows(); k++)
        {
          double T = Rx[k];
          eta[j] += T*T; //norm of the model term
        }
        // eta and rho needed to plot the Lcurve and determine the L corner
        rho[j] = sqrt(rho[j]);
        eta[j] = sqrt(eta[j]);
      }
      boost::shared_ptr<TikhonovAlgorithm::LCurveInput> lcurveInput(new TikhonovAlgorithm::LCurveInput(rho, eta, lambdaArray, nLambda));
      lcurveInput_handle_ = lcurveInput;
      lambda = FindCorner(*lcurveInput_handle_, lambda_index);

      if (input.updateLCurveGui_)
        input.updateLCurveGui_(lambda, *lcurveInput_handle_, lambda_index);
    }
    lambda_sq = lambda * lambda;
    //compute the solution with the selected regularization parameter: lambda
    regularizationParameter_.reset(new DenseColumnMatrix(1));
    (*regularizationParameter_)(0) = lambda;
    int nr_rows = regForMatrix.nrows();
    int nr_cols = regForMatrix.ncols();

    double* AtrA = ARRtrAtr.data();
    double* RtrR = CCtr.data();
    double* rm   = regForMatrix.data();
    for (int i=0; i<(int)(nr_rows*nr_cols); i++)
    {
      rm[i] = AtrA[i] + lambda_sq * RtrR[i];
    }

    if (computeRegularizedInverse_)
    {
      if (sourceWeighting_)
      {
        inverseMatrix_.reset(new DenseMatrix(RRtrAtr * regForMatrix.inverse()));
      }
      else
      {
        inverseMatrix_.reset(new DenseMatrix(forward_transpose * regForMatrix.inverse()));
      }
      inverseSolution_.reset(new DenseMatrix(*inverseMatrix_ * *measuredData_));
    }
    else
    {
      try
      {
        LinearAlgebra::solve_lapack(regForMatrix, *measuredDataRef, solution);
      }
      catch (LinearAlgebra::LapackError&)
      {
        const std::string errorMessage("The Tikhonov linear system could not be solved for a regularization parameter in the Lambda Range of the L-curve. Use a higher Lambda Range ''From'' value for the L-Curve calculation.");
        if (pr_)
        {
          pr_->error(errorMessage);
        }
        else
        {
          std::cerr << errorMessage << std::endl;
        }
        throw;
      }
      catch(DimensionMismatch&)
      {
        const std::string errorMessage("Invalid matrix sizes are being used in the Tikhonov linear system.");
        if (pr_)
        {
          pr_->error(errorMessage);
        }
        else
        {
          std::cerr << errorMessage << std::endl;
        }
        throw;
      }

      RRtrAtrsolution = RRtrAtr * solution;
      inverseSolution_.reset(new DenseMatrix(RRtrAtrsolution));
    }

  }
  else
    //OVERDETERMINED CASE,
    //similar procedure as underdetermined case (documentation comments similar, see above)
    if ( ((regularizationChoice_ == automatic) && (M>=N)) || (regularizationChoice_==overdetermined) )
    {
      //.........................................................................
      // OPERATE ON DATA:
      // Computes X = (A^T * C^T * C * A + LAMBDA * LAMBDA * R^T * R) * A^T * C^T * C *Y
      //.........................................................................
      // calculate A^T * A
      DenseMatrix forward_transpose = forwardMatrix_->transpose();
      DenseMatrix AtrA;
      DenseMatrix CCtr;
      DenseMatrix matrixNoiseCov_transpose;

      if (!sensorWeighting_)
      //check if a Sensor Space (Noise Covariance) Weighting Matrix exist?
      {
        // For large problems, allocating this matrix causes the module to hang.
        // The CCtr matrix is only used if sensorWeighting is available.
        //
        // In general, this code really should be rewritten to use sparse matrices (SCIRun 5)
        // or have replaced with more efficient calculations
        // (Eigen may have algorithms we can use).
        //CCtr = DenseMatrix::identity(M);
        AtrA = forward_transpose * *forwardMatrix_;
      }
      else
      {
        if (((regularizationResidualSubcase_ == residual_constrained) && ((M != sensorWeighting_->nrows()) || (M != sensorWeighting_->ncols()))) ||
            ((regularizationResidualSubcase_ == residual_constrained_squared) && (M != sensorWeighting_->ncols())))
        {
          BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Data Residual Weighting Matrix (number of rows and columns) must fit matrix dimensions of Forward Matrix (number of rows)!"));
        }

        matrixNoiseCov_transpose = sensorWeighting_->transpose();

        // regularizationResidualSubcase_ is either residual_constrained (default) or residual_constrained_squared
        if (regularizationResidualSubcase_ == residual_constrained)
          CCtr = *sensorWeighting_;
        else
          if (regularizationResidualSubcase_ == residual_constrained_squared)
            CCtr = matrixNoiseCov_transpose * *sensorWeighting_;

        AtrA = forward_transpose * CCtr * *forwardMatrix_;
      }
      // calculate R^T * R
      DenseMatrix RtrR;
      DenseMatrix regMat;

      if (!sourceWeighting_)
      {
        regMat = DenseMatrix::Identity(N, N);
        RtrR = DenseMatrix::Identity(N, N);
      }
      else
      {
        if(((regularizationSolutionSubcase_==solution_constrained) && ((N != sourceWeighting_->nrows()) || (N != sourceWeighting_->ncols()))) ||
           ((regularizationSolutionSubcase_==solution_constrained_squared) && (N != sourceWeighting_->ncols())))
        {
          BOOST_THROW_EXCEPTION(DimensionMismatch() << DimensionMismatchInfo("Solution Weighting Matrix (number of rows and columns) must fit matrix dimensions of Forward Matrix (number of columns) ! "));
        }
        regMat = *sourceWeighting_;

        if (regularizationSolutionSubcase_ == solution_constrained_squared)
          RtrR = sourceWeighting_->transpose() * regMat;
        else if (regularizationSolutionSubcase_ == solution_constrained)
          RtrR = regMat;
      }

      double lambda = 0, lambda_sq = 0;
      int lambda_index = 0;
      // calculate A^T * Y
      DenseColumnMatrix AtrY(N);
      DenseColumnMatrix measuredDataRef;
      if (sensorWeighting_)
      {
        measuredDataRef = CCtr * *measuredData_;
      }
      else
      {
        measuredDataRef = *measuredData_;
      }
      AtrY = forwardMatrix_->transpose() * measuredDataRef;
      DenseMatrix regForMatrix(N, N);
      DenseColumnMatrix solution(N);
      DenseColumnMatrix Ax(M);
      DenseColumnMatrix Rx(N);
      if ((input.regMethod_ == "single") || (input.regMethod_ == "slider"))
      {
        if (input.regMethod_ == "single")
        {
          // Use single fixed lambda value, entered in UI
          lambda = input.lambdaFromTextEntry_;
        }
        else if (input.regMethod_ == "slider")
        {
          // Use single fixed lambda value, select via slider
          lambda = input.lambdaSlider_;
        }
      }
      else if (input.regMethod_ == "lcurve")
      {
        // Use L-curve, lambda from corner of the L-curve
        const int nLambda = input.lambdaCount_;

        std::vector<double> lambdaArray(nLambda);
        std::vector<double> rho(nLambda);
        std::vector<double> eta(nLambda);

        lambdaArray[0] = input.lambdaMin_;
        const double lam_step = pow(10.0, log10(input.lambdaMax_ / input.lambdaMin_) / (nLambda-1));

        double* AtrAptr = AtrA.data();
        double* RtrRptr = RtrR.data();
        double* rm   = regForMatrix.data();

        int s = N*N;
        for (int j = 0; j < nLambda; j++)
        {
          if (j)
          {
            lambdaArray[j] = lambdaArray[j-1] * lam_step;
          }

          lambda_sq = lambdaArray[j] * lambdaArray[j];

          ///////////////////////////////////////
          ////Calculating the solution directly
          ///////////////////////////////////////
          for (int i = 0; i < s; i++)
          {
            rm[i] = AtrAptr[i] + lambda_sq * RtrRptr[i];
          }

          // Before, solution will be equal to (A^T * y)
          // After, solution will be equal to x_reg
          try
          {
            LinearAlgebra::solve_lapack(regForMatrix, AtrY, solution);
          }
          catch (LinearAlgebra::LapackError&)
          {
            const std::string errorMessage("The Tikhonov linear system could not be solved for a regularization parameter in the Lambda Range of the L-curve. Use a higher Lambda value for the L-Curve calculation.");
            if (pr_)
            {
              pr_->error(errorMessage);
            }
            else
            {
              std::cerr << errorMessage << std::endl;
            }
            throw;
          }
          catch(DimensionMismatch&)
          {
            const std::string errorMessage("Invalid matrix sizes are being used in the Tikhonov linear system.");
            if (pr_)
            {
              pr_->error(errorMessage);
            }
            else
            {
              std::cerr << errorMessage << std::endl;
            }
            throw;
          }

          ////////////////////////////////
          const DenseMatrix &forwardMatrixRef = *forwardMatrix_;

          Ax = forwardMatrixRef * solution;
          Rx = regMat * solution;
          // Calculate the norm of Ax-b and Rx
          rho[j]=0; eta[j]=0;
          for (int k = 0; k < M; k++)
          {
            double T = Ax(k) - (*measuredData_)(k);
            rho[j] += T*T;
          }

          for (int k = 0; k < N; k++)
          {
            double T = Rx(k);
            eta[j] += T*T;
          }

          rho[j] = sqrt(rho[j]);
          eta[j] = sqrt(eta[j]);
        }
        boost::shared_ptr<TikhonovAlgorithm::LCurveInput> lcurveInput(new TikhonovAlgorithm::LCurveInput(rho, eta, lambdaArray, nLambda));
        lcurveInput_handle_ = lcurveInput;
        lambda = FindCorner(*lcurveInput_handle_, lambda_index);

        if (input.updateLCurveGui_)
          input.updateLCurveGui_(lambda, *lcurveInput_handle_, lambda_index);

      } // END  else if (reg_method_.get() == "lcurve")
      lambda_sq = lambda * lambda;

      regularizationParameter_.reset(new DenseColumnMatrix(1));
      (*regularizationParameter_)[0] = lambda;
      int nr_rows = regForMatrix.nrows();
      int nr_cols = regForMatrix.ncols();
      const int NR_BOUNDS_CHECK = nr_rows * nr_cols;
      const double* AtrAptr = AtrA.data();
      const double* RtrRptr = RtrR.data();
      double* rm   = regForMatrix.data();

      for (int i = 0; i < NR_BOUNDS_CHECK; i++)
      {
        rm[i] = AtrAptr[i] + lambda_sq * RtrRptr[i];
      }
      if (computeRegularizedInverse_)
      {
        if (CCtr.empty())
          inverseMatrix_.reset(new DenseMatrix((regForMatrix.inverse() * forward_transpose)));
        else
          inverseMatrix_.reset(new DenseMatrix((regForMatrix.inverse() * forward_transpose) * CCtr));

        inverseSolution_.reset(new DenseMatrix(*inverseMatrix_ * *measuredData_));
      }
      else
      {
        try
        {
          LinearAlgebra::solve_lapack(regForMatrix, AtrY, solution);
        }
        catch (LinearAlgebra::LapackError&)
        {
          const std::string errorMessage("The Tikhonov linear system could not be solved for a regularization parameter in the Lambda Range of the L-curve. Use a higher Lambda value for the L-Curve calculation.");
          if (pr_)
          {
            pr_->error(errorMessage);
          }
          else
          {
            std::cerr << errorMessage << std::endl;
          }
          throw;
        }
        catch(DimensionMismatch&)
        {
          const std::string errorMessage("Invalid matrix sizes are being used in the Tikhonov linear system.");
          if (pr_)
          {
            pr_->error(errorMessage);
          }
          else
          {
            std::cerr << errorMessage << std::endl;
          }
          throw;
        }
        inverseSolution_.reset(new DenseMatrix(solution));
      }

    }
}

void TikhonovAlgorithmImpl::update_graph(const TikhonovAlgorithmImpl::Input& input, double lambda, int lambda_index, const double epsilon)
{
  if (lcurveInput_handle_ && input.updateLCurveGui_)
  {
    lambda = LambdaLookup(*lcurveInput_handle_, lambda, lambda_index, epsilon);
    if (lambda >= 0)
    {
      input.updateLCurveGui_(lambda, *lcurveInput_handle_, lambda_index);
    }
  }
}

#if 0

void SolveInverseProblemWithTikhonov::tcl_command(GuiArgs& args, void* userdata)
{
  if (args[1] == "updategraph" && args.count() == 4)
  {
    double lambda = boost::lexical_cast<double>(args[2]);
    int lambda_index = boost::lexical_cast<double>(args[3]);

    if (input_handle_.get() != 0 && algo_handle_.get() != 0)
    {
      algo_handle_->update_graph(*input_handle_, lambda, lambda_index, lambda_resolution_.get());
    }
  }
  else
  {
    // Relay data to the Module class
    Module::tcl_command(args, userdata);
  }
}

#endif

TikhonovAlgorithm::LCurveInput::LCurveInput(const std::vector<double>& rho, const std::vector<double>& eta, const std::vector<double>& lambdaArray, int nLambda)
: rho_(rho), eta_(eta), lambdaArray_(lambdaArray), nLambda_(nLambda)
{}

TikhonovAlgorithmImpl::Input::Input(const std::string& regMethod, double lambdaFromTextEntry, double lambdaSlider, int lambdaCount, double lambdaMin, double lambdaMax,
                                    lcurveGuiUpdate updateLCurveGui)
: regMethod_(regMethod), lambdaFromTextEntry_(lambdaFromTextEntry), lambdaSlider_(lambdaSlider), lambdaCount_(lambdaCount), lambdaMin_(lambdaMin), lambdaMax_(lambdaMax),
updateLCurveGui_(updateLCurveGui)
{}

} // End namespace BioPSE

ext: cc

max_stars_repo_name: benjaminlarson/SCIRunGUIPrototype

max_issues_repo_name: benjaminlarson/SCIRunGUIPrototype

max_forks_repo_name: benjaminlarson/SCIRunGUIPrototype


text:
#pragma once

#include "util.hpp"

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/ring.hpp>
#include <boost/geometry/algorithms/covered_by.hpp>

#ifndef HABITAT_H
#define HABITAT_H

typedef boost::geometry::model::point<double,2,boost::geometry::cs::cartesian> Point;
typedef boost::geometry::model::ring<Point> Ring;

/*
 The habitat is represented as a ring (a polygon without holes).
 The vertices of the ring are specified by the user in the file `datapath.outer`.
 For example, suppose that `datapath.outer` contains the following six lines:
 0 0
 3 3
 0 6
 12 6
 12 0
 0 0
 Then the habitat is given by (0,0) - (3,3) - (0,6) - (12,6) - (12,0) - (0,0).
 The vertices of the population graph (the demes) will fall inside the habitat by construction.

 However, EEMS will not check that the sampling locations fall inside the habitat -- instead,
 each sample will be assigned to the closest deme. Therefore, the user should specify a habitat
 that is sufficiently large to cover all sampling locations.
 */

class Habitat {
public:

    Habitat( );
    ~Habitat( );

    void generate_outer(const string &datapath);
    bool dlmwrite_outer(const string &mcmcpath) const;
    bool in_point(const double x, const double y) const;

    double get_area( ) const;
    double get_xmin( ) const;
    double get_xmax( ) const;
    double get_ymin( ) const;
    double get_ymax( ) const;
    double get_xspan( ) const;
    double get_yspan( ) const;

private:

    Ring domain;
    double xmin, xmax, xspan;
    double ymin, ymax, yspan;

};

#endif