text:
// Boost.Geometry - gis-projections (based on PROJ4)

// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.

// This file was modified by Oracle on 2017, 2018.
// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.

// Use, modification and distribution is subject to 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)

// This file is converted from PROJ4, http://trac.osgeo.org/proj
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
// PROJ4 is maintained by Frank Warmerdam
// PROJ4 is converted to Boost.Geometry by Barend Gehrels

// Last updated version of proj: 5.0.0

// Original copyright notice:

// 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.

#ifndef BOOST_GEOMETRY_PROJECTIONS_OEA_HPP
#define BOOST_GEOMETRY_PROJECTIONS_OEA_HPP

#include <boost/math/special_functions/hypot.hpp>

#include <boost/geometry/srs/projections/impl/base_static.hpp>
#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
#include <boost/geometry/srs/projections/impl/projects.hpp>
#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
#include <boost/geometry/srs/projections/impl/aasincos.hpp>

namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace oea
    {
            template <typename T>
            struct par_oea
            {
                T    theta;
                T    m, n;
                T    two_r_m, two_r_n, rm, rn, hm, hn;
                T    cp0, sp0;
            };

            // template class, using CRTP to implement forward/inverse
            template <typename T, typename Parameters>
            struct base_oea_spheroid
                : public base_t_fi<base_oea_spheroid<T, Parameters>, T, Parameters>
            {
                par_oea<T> m_proj_parm;

                inline base_oea_spheroid(const Parameters& par)
                    : base_t_fi<base_oea_spheroid<T, Parameters>, T, Parameters>(*this, par)
                {}

                // FORWARD(s_forward)  sphere
                // Project coordinates from geographic (lon, lat) to cartesian (x, y)
                inline void fwd(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
                {
                    T Az, M, N, cp, sp, cl, shz;

                    cp = cos(lp_lat);
                    sp = sin(lp_lat);
                    cl = cos(lp_lon);
                    Az = aatan2(cp * sin(lp_lon), this->m_proj_parm.cp0 * sp - this->m_proj_parm.sp0 * cp * cl) + this->m_proj_parm.theta;
                    shz = sin(0.5 * aacos(this->m_proj_parm.sp0 * sp + this->m_proj_parm.cp0 * cp * cl));
                    M = aasin(shz * sin(Az));
                    N = aasin(shz * cos(Az) * cos(M) / cos(M * this->m_proj_parm.two_r_m));
                    xy_y = this->m_proj_parm.n * sin(N * this->m_proj_parm.two_r_n);
                    xy_x = this->m_proj_parm.m * sin(M * this->m_proj_parm.two_r_m) * cos(N) / cos(N * this->m_proj_parm.two_r_n);
                }

                // INVERSE(s_inverse)  sphere
                // Project coordinates from cartesian (x, y) to geographic (lon, lat)
                inline void inv(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
                {
                    T N, M, xp, yp, z, Az, cz, sz, cAz;

                    N = this->m_proj_parm.hn * aasin(xy_y * this->m_proj_parm.rn);
                    M = this->m_proj_parm.hm * aasin(xy_x * this->m_proj_parm.rm * cos(N * this->m_proj_parm.two_r_n) / cos(N));
                    xp = 2. * sin(M);
                    yp = 2. * sin(N) * cos(M * this->m_proj_parm.two_r_m) / cos(M);
                    cAz = cos(Az = aatan2(xp, yp) - this->m_proj_parm.theta);
                    z = 2. * aasin(0.5 * boost::math::hypot(xp, yp));
                    sz = sin(z);
                    cz = cos(z);
                    lp_lat = aasin(this->m_proj_parm.sp0 * cz + this->m_proj_parm.cp0 * sz * cAz);
                    lp_lon = aatan2(sz * sin(Az),
                        this->m_proj_parm.cp0 * cz - this->m_proj_parm.sp0 * sz * cAz);
                }

                static inline std::string get_name()
                {
                    return "oea_spheroid";
                }

            };

            // Oblated Equal Area
            template <typename Params, typename Parameters, typename T>
            inline void setup_oea(Params const& params, Parameters& par, par_oea<T>& proj_parm)
            {
                if (((proj_parm.n = pj_get_param_f<T, srs::spar::n>(params, "n", srs::dpar::n)) <= 0.) ||
                    ((proj_parm.m = pj_get_param_f<T, srs::spar::m>(params, "m", srs::dpar::m)) <= 0.)) {
                    BOOST_THROW_EXCEPTION( projection_exception(error_invalid_m_or_n) );
                } else {
                    proj_parm.theta = pj_get_param_r<T, srs::spar::theta>(params, "theta", srs::dpar::theta);
                    proj_parm.sp0 = sin(par.phi0);
                    proj_parm.cp0 = cos(par.phi0);
                    proj_parm.rn = 1./ proj_parm.n;
                    proj_parm.rm = 1./ proj_parm.m;
                    proj_parm.two_r_n = 2. * proj_parm.rn;
                    proj_parm.two_r_m = 2. * proj_parm.rm;
                    proj_parm.hm = 0.5 * proj_parm.m;
                    proj_parm.hn = 0.5 * proj_parm.n;
                    par.es = 0.;
                }
            }

    }} // namespace detail::oea
    #endif // doxygen

    /*!
        \brief Oblated Equal Area projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Miscellaneous
         - Spheroid
        \par Projection parameters
         - n (real)
         - m (real)
         - theta: Theta (degrees)
        \par Example
        \image html ex_oea.gif
    */
    template <typename T, typename Parameters>
    struct oea_spheroid : public detail::oea::base_oea_spheroid<T, Parameters>
    {
        template <typename Params>
        inline oea_spheroid(Params const& params, Parameters const& par)
            : detail::oea::base_oea_spheroid<T, Parameters>(par)
        {
            detail::oea::setup_oea(params, this->m_par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_oea, oea_spheroid, oea_spheroid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(oea_entry, oea_spheroid)

        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(oea_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(oea, oea_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_OEA_HPP


text:
#include <iostream>
#include <vector>
#include <fstream>
#include <boost/multiprecision/gmp.hpp>

#include "../../external_libs/json.hpp"
#include "../components/order.hpp"
#include "../components/token.hpp"
#include "../api/aux.hpp"
#include "./optimize.hpp"

/************************************** Global variables that are defined as static in various classes *************************************/

int Token::num_tokens{ 0 };                     // static integer variable of class Token, that counts number of tokens in our instance
std::map<std::string, int> Token::idx_tokens;   // static map of class Token, that maps the address of the token to a unique index

/*******************************************************************************************************************************************/

int main()
{

/********************************************************** VARIABLE DECLARATIONS **********************************************************/

    std::vector<Token> tokens;      // the vector containing all the tokens appearing in any order/AMM
    std::vector<Order> orders;      // the vector containing all the "regular" orders
    std::vector<CP_AMM> amms;       // the vector containing all the AMMs

/*******************************************************************************************************************************************/

    nlohmann::json json_file;
    std::ifstream input_file("./sample_instances/sample.json");
    input_file >> json_file;

    // Parsing the input
    parse_json_file(tokens, orders, amms, Token::num_tokens, Token::idx_tokens, json_file);

    // Just printing out the input
    std::cout.precision(std::numeric_limits<boost::multiprecision::mpf_float>::digits10);

    for (auto &i: tokens)
        print_token(i);

    for (auto &i: orders)
        print_order(i, tokens);

    for (auto &i: amms)
        print_cp_amm(i, tokens);

    // solve batch auction
    solve_auction(tokens, orders, amms);

    // return solution

    return 0;
}

ext: cpp

max_stars_repo_name: cowprotocol/simple-solver-template

max_issues_repo_name: cowprotocol/simple-solver-template

max_forks_repo_name: cowprotocol/simple-solver-template


text:
/*
	Copyright (C) 2003-2013 by Kristina Simpson <sweet.kristas@gmail.com>

	This software is provided 'as-is', without any express or implied
	warranty. In no event will the authors be held liable for any damages
	arising from the use of this software.

	Permission is granted to anyone to use this software for any purpose,
	including commercial applications, and to alter it and redistribute it
	freely, subject to the following restrictions:

	   1. The origin of this software must not be misrepresented; you must not
	   claim that you wrote the original software. If you use this software
	   in a product, an acknowledgment in the product documentation would be
	   appreciated but is not required.

	   2. Altered source versions must be plainly marked as such, and must not be
	   misrepresented as being the original software.

	   3. This notice may not be removed or altered from any source
	   distribution.
*/

#pragma once

#include <boost/property_tree/ptree.hpp>
#include "svg_attribs.hpp"
#include "svg_element.hpp"
#include "svg_length.hpp"
#include "svg_paint.hpp"
#include "svg_render.hpp"
#include "svg_transform.hpp"

namespace KRE
{
	namespace SVG
	{
		enum class GradientCoordSystem {
			USERSPACE_ON_USE,
			OBJECT_BOUNDING_BOX,
		};
		enum class GradientSpreadMethod {
			PAD,
			REFLECT,
			REPEAT,
		};

		class gradient_stop : public core_attribs// : public core_attribs, public presentation_attribs
		{
		public:
			gradient_stop(element* doc, const boost::property_tree::ptree& pt);
			virtual ~gradient_stop();
			void apply(render_context& ctx, cairo_pattern_t* pattern);
		private:
			double offset_;			// number or percent.
			paint_ptr color_;
			double opacity_;
			bool opacity_set_;
		};
		typedef std::shared_ptr<gradient_stop> gradient_stop_ptr;

		class gradient : public core_attribs
		{
		public:
			gradient(element* doc, const boost::property_tree::ptree& pt);
			virtual ~gradient();
			void set_source(render_context& ctx) const;
			void apply_stops(render_context& ctx, cairo_pattern_t* pattern) const;
			void apply_transforms(cairo_pattern_t* pattern) const;
		private:
			virtual void handle_set_source(render_context& ctx) const = 0;
			GradientCoordSystem coord_system_;
			GradientSpreadMethod spread_;
			std::string xlink_href_;
			std::vector<transform_ptr> transforms_;
			std::vector<gradient_stop_ptr> stops_;
		};
		typedef std::shared_ptr<gradient> gradient_ptr;

		class linear_gradient : public gradient
		{
		public:
			linear_gradient(element* doc, const boost::property_tree::ptree& pt);
			virtual ~linear_gradient();
		private:
			virtual void handle_set_source(render_context& ctx) const override;
			svg_length x1_;
			svg_length y1_;
			svg_length x2_;
			svg_length y2_;
			mutable std::shared_ptr<cairo_pattern_t> pattern_;
		};

		class radial_gradient : public gradient
		{
		public:
			radial_gradient(element* doc, const boost::property_tree::ptree& pt);
			virtual ~radial_gradient();
		private:
			virtual void handle_set_source(render_context& ctx) const override;
			svg_length cx_;
			svg_length cy_;
			svg_length r_;
			svg_length fx_;
			svg_length fy_;
			mutable std::shared_ptr<cairo_pattern_t> pattern_;
		};
	}
}

ext: hpp

max_stars_repo_name: sweetkristas/anura

max_issues_repo_name: sweetkristas/anura

max_forks_repo_name: sweetkristas/anura


text:
/*******************************************************************************
 *         Copyright 2003 & onward LASMEA UMR 6602 CNRS/Univ. Clermont II
 *         Copyright 2009 & onward LRI    UMR 8623 CNRS/Univ Paris Sud XI
 *
 *          Distributed under the Boost Software License, Version 1.0.
 *                 See accompanying file LICENSE.txt or copy at
 *                     http://www.boost.org/LICENSE_1_0.txt
 ******************************************************************************/
#ifndef NT2_SDK_SIMD_DETAILS_NATIVE_CONSTANTS_HPP_INCLUDED
#define NT2_SDK_SIMD_DETAILS_NATIVE_CONSTANTS_HPP_INCLUDED

#include <boost/tr1/functional.hpp>
#include <nt2/sdk/meta/as.hpp>
#include <nt2/sdk/simd/category.hpp>
#include <nt2/sdk/meta/from_bits.hpp>
#include <nt2/sdk/constant/boolean.hpp>
#include <nt2/sdk/simd/details/impl/splat.hpp>
#include <nt2/sdk/simd/details/impl/constants.hpp>

////////////////////////////////////////////////////////////////////////////////
// Forward all constant call to the simd version of themselves that splat
// the appropriate scalar constants into a proper SIMD vector.
////////////////////////////////////////////////////////////////////////////////
NT2_REGISTER_DISPATCH ( Tag, tag::cpu_
                      , (Tag)(A0)(X)
                      , ((target_< simd_< arithmetic_<A0>,X> >))
                      )

namespace nt2 { namespace ext
{
  template<class Tag, class X, class Dummy>
  struct  call< Tag(tag::target_< tag::simd_<tag::arithmetic_,X> >)
              , tag::cpu_, Dummy
              >
          : callable
  {
    template<class Sig> struct result;
    template<class This, class Target>
    struct result<This(Target)>
    {
      typedef typename meta::strip<Target>::type::type          target;
      typedef typename meta::scalar_of<target>::type            base;
      typedef typename
      std::tr1::result_of<functor<Tag>(nt2::meta::as_<base>)>::type  results;
      typedef typename target::template cast<results>::type      type;
    };

    NT2_FUNCTOR_CALL(1)
    {
      typedef typename meta::strip<A0>::type::type        target;
      typedef typename meta::scalar_of<target>::type      type;
      typedef typename NT2_RETURN_TYPE(1)::type           result_type;
      functor<Tag> callee;
      return splat<result_type>(callee( nt2::meta::as_<type>()));
    }
  };
} }

#endif


text:
#include "test.h"

#include "mjolnir/graphbuilder.h"

#include <algorithm>
#include <boost/optional.hpp>
#include <boost/property_tree/ptree.hpp>
#include <memory>
#include <sstream>
#include <string>
#include <thread>

using namespace std;
using namespace valhalla::mjolnir;

namespace {}

// TODO: sweet jesus add more tests of this class!

int main(int argc, char* argv[]) {
  testing::InitGoogleTest(&argc, argv);
  return RUN_ALL_TESTS();
}

ext: cc

max_stars_repo_name: mixvit/valhalla

max_issues_repo_name: Bolaxax/valhalla

max_forks_repo_name: Bolaxax/valhalla


text:
//  Boost common_factor_ct.hpp header file  ----------------------------------//

//  (C) Copyright John Maddock 2017.
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)

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

#ifndef BOOST_MATH_COMMON_FACTOR_CT_HPP
#define BOOST_MATH_COMMON_FACTOR_CT_HPP

#include <boost/integer/common_factor_ct.hpp>

namespace boost
{
namespace math
{

   using boost::integer::static_gcd;
   using boost::integer::static_lcm;
   using boost::integer::static_gcd_type;

}  // namespace math
}  // namespace boost

#endif  // BOOST_MATH_COMMON_FACTOR_CT_HPP


text:
#include <k52/optimization/hleborodov_rosenbrock_method.h>

#ifdef BUILD_WITH_MPI

#include <boost/mpi.hpp>
#include <k52/parallel/mpi/constants.h>

#endif

#include <cmath>
#include <stdexcept>
#include <algorithm>
#include <fstream>
#include <iostream>

#include <k52/optimization/params/i_continuous_parameters.h>

using ::std::vector;

namespace
{
const double kArgumentsIncreaseValue = 1.5;
const double kArgumentsDecreaseValue = -0.5;
}

namespace k52
{
namespace optimization
{

HleborodovRosenbrockMethod::HleborodovRosenbrockMethod(double precision, size_t max_iteration_number,
    double first_step, double max_step)
    : precision_(precision)
    , arguments_increase_(kArgumentsIncreaseValue)
    , arguments_decrease_(kArgumentsDecreaseValue)
    , first_step_(first_step)
    , max_step_(max_step)
    , max_iteration_number_(max_iteration_number)
{
}

HleborodovRosenbrockMethod* HleborodovRosenbrockMethod::Clone() const
{
    return new HleborodovRosenbrockMethod(precision_, max_iteration_number_, first_step_, max_step_);
}

std::string HleborodovRosenbrockMethod::get_name() const
{
    return "Hleborodov Rosenbrock Method";
}

#ifdef BUILD_WITH_MPI
void HleborodovRosenbrockMethod::Send(boost::mpi::communicator* communicator, int target) const
{
    communicator->send(target, k52::parallel::mpi::constants::kCommonTag, precision_);
    communicator->send(target, k52::parallel::mpi::constants::kCommonTag, first_step_);
    communicator->send(target, k52::parallel::mpi::constants::kCommonTag, max_step_);
    communicator->send(target, k52::parallel::mpi::constants::kCommonTag, max_iteration_number_);
}

void HleborodovRosenbrockMethod::Receive(boost::mpi::communicator* communicator, int source)
{
    communicator->recv(source, k52::parallel::mpi::constants::kCommonTag, precision_);
    communicator->recv(source, k52::parallel::mpi::constants::kCommonTag, first_step_);
    communicator->recv(source, k52::parallel::mpi::constants::kCommonTag, max_step_);
    communicator->recv(source, k52::parallel::mpi::constants::kCommonTag, max_iteration_number_);
}
#endif

vector<double> HleborodovRosenbrockMethod::FindOptimalParameters(const vector<double>& initial_parameters)
{
    vector<double> arguments = initial_parameters;
    dimension_ = arguments.size();
    vector<double> next_step_arguments(arguments);

    InitializeBasisAndSteps();

    for (size_t i = 1; i <= max_iteration_number_; i++)
    {
        // Research of next arguments, that minimize function
        MakeStep(next_step_arguments);

        if (IsExitCriteriaFulfilled(next_step_arguments, arguments))
        {
            break;
        }

        // Check number of iterations
        if (i == max_iteration_number_)
        {
            std::cout << " Minimun has not found";
        }

        // Make new basis by Gramm-Shmidt method
        CreateNewBasis();
        arguments = next_step_arguments;
    }

    return next_step_arguments;
}

void HleborodovRosenbrockMethod::InitializeBasisAndSteps()
{
    steps_array_ = vector<double> (dimension_,first_step_);

    basis_ = vector< vector<double> > (dimension_, vector<double> (dimension_, 0));

    // Initialize basis is Decart coordinate system
    for(size_t i = 0; i<dimension_; i++)
    {
        basis_[i][i] = 1;
    }
}

bool HleborodovRosenbrockMethod::IsExitCriteriaFulfilled(
    const vector<double>& arguments,
    const vector<double>& previous_step_arguments)
{
    bool escape = true;
    // check arguments increases to escape
    for(size_t j = 0; j<dimension_; j++)
    {
        if(std::abs(arguments[j] - previous_step_arguments[j]) > precision_)
        {
            escape = false;
        }
    }
    return escape;
}

void HleborodovRosenbrockMethod::MakeStep(vector<double> &arguments)
{
    double start_function_value = CountObjectiveFunctionValueToMinimize(arguments);

    vector<double> tested_arguments(arguments);
    vector<double> tested_arguments2(arguments);

    /// Sort steps_array_ by descending for beter method speed
    std::sort(steps_array_.begin(), steps_array_.end(), std::greater<double>());

    for(size_t i = 0; i<dimension_; i++)
    {
        // make test step to minimize function in current direction
        // make it by each of coordinates
        MakeStepPerCoordinate(i, &tested_arguments);

        double counted_value = CountObjectiveFunctionValueToMinimize(tested_arguments);

        if(start_function_value > counted_value)
        {
            steps_array_[i] *= arguments_increase_;
        }
        else
        {
            steps_array_[i] *= arguments_decrease_;
        }

        CorrectStep(i);

        MakeStepPerCoordinate(i, &tested_arguments2);

        tested_arguments = arguments;
    }

    arguments = tested_arguments2;
}

void HleborodovRosenbrockMethod::MakeStepPerCoordinate(int basis_index, vector<double>* arguments)
{
    for(size_t j = 0; j<dimension_; j++)
    {
        (*arguments)[j] += steps_array_[basis_index]*basis_[basis_index][j];
    }
}

void HleborodovRosenbrockMethod::CreateNewBasis()
{
    vector< vector<double> > matrix_of_coeff(dimension_, vector<double> (dimension_, 0));
    vector< vector<double> > unnormalized_basis(dimension_, vector<double> (dimension_, 0));

    for(size_t j = 0; j<dimension_; j++)
    {
        for(size_t i = j; i<dimension_; i++)
        {
            for(size_t k = 0; k<dimension_; k++)
            {
                matrix_of_coeff[j][k] += steps_array_[i]*basis_[i][k];
            }
        }
    }

    for(size_t i = 0; i<dimension_; i++)
    {
        for(size_t j = 0; j<dimension_; j++)
        {
            unnormalized_basis[i][j] = matrix_of_coeff[i][j];
        }
        basis_[0][i] = matrix_of_coeff[0][i];
    }

    NormalizeVector(basis_[0]);
    double temp;

    for(size_t j = 1; j<dimension_; j++)
    {
        for(size_t i = 0; i<j; i++)
        {
            temp = ScalarComposition(matrix_of_coeff[j],basis_[i]);
            for(size_t k = 0; k<dimension_; k++)
            {
                unnormalized_basis[j][k] -= temp*basis_[i][k];
            }
        }

        for(size_t k = 0; k<dimension_; k++)
        {
            basis_[j][k] = unnormalized_basis[j][k];
        }
        NormalizeVector(basis_[j]);
    }
}

void HleborodovRosenbrockMethod::NormalizeVector(vector<double> &target_vector)
{
    double temp = 0;
    for(size_t i = 0; i<dimension_; i++)
    {
        temp += target_vector[i]*target_vector[i];
    }
    temp = sqrt(temp);
    for(size_t i = 0; i<dimension_; i++)
    {
        if(temp != 0)
        {
            target_vector[i] /= temp;
        }
    }
}

double HleborodovRosenbrockMethod::ScalarComposition(
    const vector<double>& first_vector,
    const vector<double>& second_vector)
{
    double scalar_composition = 0;
    for(size_t i = 0; i<dimension_; i++)
    {
        scalar_composition += first_vector[i]*second_vector[i];
    }
    return scalar_composition;
}

void HleborodovRosenbrockMethod::CorrectStep(size_t step_index)
{
    if( std::abs(steps_array_[step_index]) > max_step_)
    {
        steps_array_[step_index] =
            max_step_*steps_array_[step_index]/std::abs(steps_array_[step_index]);
    }
}

}/* namespace optimization */
}/* namespace k52 */


text:
// Copyright (c) 2017 Graphcore Ltd. All rights reserved.
#include <boost/multi_array.hpp>
#include <cassert>
#include <poplibs_test/GeneralMatrixAdd.hpp>
#include <poplibs_test/GeneralMatrixMultiply.hpp>
#include <poplibs_test/Lstm.hpp>
#include <poplibs_test/NonLinearity.hpp>
#include <unordered_map>

// Fwd state array indices
#define LSTM_FWD_STATE_FORGET_GATE 2
#define LSTM_FWD_STATE_CAND_TANH 3
#define LSTM_FWD_STATE_INPUT_GATE 4
#define LSTM_FWD_STATE_OUTPUT_GATE 5
#define LSTM_FWD_STATE_OUTPUT_TANH 6

using IndexRange = boost::multi_array_types::index_range;
using Array1dRef = boost::multi_array_ref<double, 1>;
using Array2dRef = boost::multi_array_ref<double, 2>;
using Array2d = boost::multi_array<double, 2>;
using Array3dRef = boost::multi_array_ref<double, 3>;
using Array4dRef = boost::multi_array_ref<double, 4>;
using Array3d = boost::multi_array<double, 3>;

using namespace poplibs_test;

/**
 * Process a given unit type within an LSTM given its weights and biases.
 * The non-linearity is also specified although it may be derived from the unit
 */
static void processBasicLstmUnit(const Array2dRef prevOutput,
                                 const Array2dRef input,
                                 const Array3dRef weightsInput,
                                 const Array3dRef weightsOutput,
                                 const Array2dRef biases, Array2dRef output,
                                 unsigned lstmUnitOffset,
                                 popnn::NonLinearityType nonLinearityType) {
  const auto batchSize = prevOutput.shape()[0];
  const auto outputSize = prevOutput.shape()[1];

  /* split weight into two parts:
   * 1) part which weighs only the previous output
   * 2) part which weighs only the input
   */
  Array2d weightsOutputUnit = weightsOutput[lstmUnitOffset];
  Array2d weightsInputUnit = weightsInput[lstmUnitOffset];

  gemm::generalMatrixMultiply(prevOutput, weightsOutputUnit, output, output,
                              1.0, 0, false, false);
  gemm::generalMatrixMultiply(input, weightsInputUnit, output, output, 1.0, 1.0,
                              false, false);
  /* add bias */
  for (auto b = 0U; b != batchSize; ++b) {
    for (auto i = 0U; i != outputSize; ++i) {
      output[b][i] += biases[lstmUnitOffset][i];
    }
  }

  /* apply non-linearity */
  nonLinearity(nonLinearityType, output);
}

static std::unordered_map<BasicLstmCellUnit, unsigned>
getCellMapping(const std::vector<BasicLstmCellUnit> &cellOrder) {
  // build a mapping of the order that the gates are stored in.
  std::unordered_map<BasicLstmCellUnit, unsigned> cellMapping;
  for (unsigned i = 0; i < cellOrder.size(); ++i) {
    auto gate = cellOrder.at(i);
    cellMapping.insert(std::make_pair(gate, i));
  }

  return cellMapping;
}

void poplibs_test::lstm::basicLstmCellForwardPass(
    const Array3dRef input, const Array2dRef biases,
    const Array2dRef prevOutput, const Array3dRef weightsInput,
    const Array3dRef weightsOutput, Array2dRef prevCellState, Array4dRef state,
    const std::vector<BasicLstmCellUnit> &cellOrder) {
  const auto sequenceSize = state.shape()[1];
  const auto batchSize = state.shape()[2];
  const auto outputSize = state.shape()[3];
#ifndef NDEBUG
  const auto inputSize = input.shape()[2];
#endif
  assert(state.shape()[0] == LSTM_NUM_FWD_STATES);
  assert(weightsInput.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsInput.shape()[1] == inputSize);
  assert(weightsInput.shape()[2] == outputSize);
  assert(weightsOutput.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsOutput.shape()[1] == outputSize);
  assert(weightsOutput.shape()[2] == outputSize);
  assert(prevCellState.shape()[0] == batchSize);
  assert(prevCellState.shape()[1] == outputSize);
  assert(biases.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(biases.shape()[1] == outputSize);
  assert(prevOutput.shape()[0] == batchSize);
  assert(prevOutput.shape()[1] == outputSize);

  auto cellMapping = getCellMapping(cellOrder);

  for (auto s = 0U; s != sequenceSize; ++s) {
    Array2d ysm1 = s == 0 ? state[LSTM_FWD_STATE_ACTS_IDX][s]
                          : state[LSTM_FWD_STATE_ACTS_IDX][s - 1];
    Array2d csm1 = s == 0 ? state[LSTM_FWD_STATE_CELL_STATE_IDX][s]
                          : state[LSTM_FWD_STATE_CELL_STATE_IDX][s - 1];
    Array2d prevOutputThisStep = s == 0 ? prevOutput : ysm1;
    Array2d cellState = s == 0 ? prevCellState : csm1;
    Array2d inputThisStep = input[s];

    /* forget gate */
    Array2d forgetGate(boost::extents[batchSize][outputSize]);
    processBasicLstmUnit(prevOutputThisStep, inputThisStep, weightsInput,
                         weightsOutput, biases, forgetGate,
                         cellMapping.at(BASIC_LSTM_CELL_FORGET_GATE),
                         popnn::NonLinearityType::SIGMOID);
    state[LSTM_FWD_STATE_FORGET_GATE][s] = forgetGate;

    /* input gate */
    Array2d inputGate(boost::extents[batchSize][outputSize]);
    processBasicLstmUnit(prevOutputThisStep, inputThisStep, weightsInput,
                         weightsOutput, biases, inputGate,
                         cellMapping.at(BASIC_LSTM_CELL_INPUT_GATE),
                         popnn::NonLinearityType::SIGMOID);
    state[LSTM_FWD_STATE_INPUT_GATE][s] = inputGate;

    /* new candidate contribution to this cell */
    Array2d candidate(boost::extents[batchSize][outputSize]);
    processBasicLstmUnit(prevOutputThisStep, inputThisStep, weightsInput,
                         weightsOutput, biases, candidate,
                         cellMapping.at(BASIC_LSTM_CELL_CANDIDATE),
                         popnn::NonLinearityType::TANH);
    state[LSTM_FWD_STATE_CAND_TANH][s] = candidate;

    /* output gate */
    Array2d outputGate(boost::extents[batchSize][outputSize]);
    processBasicLstmUnit(prevOutputThisStep, inputThisStep, weightsInput,
                         weightsOutput, biases, outputGate,
                         cellMapping.at(BASIC_LSTM_CELL_OUTPUT_GATE),
                         popnn::NonLinearityType::SIGMOID);
    state[LSTM_FWD_STATE_OUTPUT_GATE][s] = outputGate;

    poplibs_test::gemm::hadamardProduct(forgetGate, cellState, cellState);
    poplibs_test::gemm::hadamardProduct(inputGate, candidate, candidate);
    poplibs_test::axpby::add(cellState, candidate, cellState);

    /* need to maintain the cell state for next step */
    Array2d outputThisStep = cellState;
    nonLinearity(popnn::NonLinearityType::TANH, outputThisStep);
    state[LSTM_FWD_STATE_OUTPUT_TANH][s] = outputThisStep;
    gemm::hadamardProduct(outputThisStep, outputGate, outputThisStep);

    state[LSTM_FWD_STATE_ACTS_IDX][s] = outputThisStep;
    state[LSTM_FWD_STATE_CELL_STATE_IDX][s] = cellState;
  }
}

static void computeGradients(const Array2dRef weightIn,
                             const Array2dRef weightPrev, const Array2dRef grad,
                             Array2dRef gradIn, Array2dRef gradPrev, bool acc) {
  double k = acc ? 1.0 : 0.0;
  gemm::generalMatrixMultiply(grad, weightIn, gradIn, gradIn, 1.0, k, false,
                              true);
  gemm::generalMatrixMultiply(grad, weightPrev, gradPrev, gradPrev, 1.0, k,
                              false, true);
}

void poplibs_test::lstm::basicLstmCellBackwardPass(
    const Array3dRef weightsInput, const Array3dRef weightsOutput,
    const Array3dRef gradsNextLayer, const Array2dRef prevCellState,
    const Array4dRef fwdState, Array4dRef bwdState, Array3dRef gradsPrevLayer,
    const std::vector<BasicLstmCellUnit> &cellOrder) {
  const auto sequenceSize = fwdState.shape()[1];
  const auto batchSize = fwdState.shape()[2];
  const auto outputSize = fwdState.shape()[3];
  const auto inputSize = gradsPrevLayer.shape()[2];

  assert(fwdState.shape()[0] == LSTM_NUM_FWD_STATES);
  assert(bwdState.shape()[0] == LSTM_NUM_BWD_STATES);
  assert(weightsInput.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsInput.shape()[1] == inputSize);
  assert(weightsInput.shape()[2] == outputSize);
  assert(weightsOutput.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsOutput.shape()[1] == outputSize);
  assert(weightsOutput.shape()[2] == outputSize);
  assert(prevCellState.shape()[0] == batchSize);
  assert(prevCellState.shape()[1] == outputSize);
  assert(fwdState.shape()[1] == sequenceSize);
  assert(fwdState.shape()[2] == batchSize);
  assert(fwdState.shape()[3] == outputSize);
  assert(bwdState.shape()[1] == sequenceSize);
  assert(bwdState.shape()[2] == batchSize);
  assert(bwdState.shape()[3] == outputSize);
  assert(gradsNextLayer.shape()[0] == sequenceSize);
  assert(gradsNextLayer.shape()[1] == batchSize);
  assert(gradsNextLayer.shape()[2] == outputSize);
  assert(gradsPrevLayer.shape()[0] == sequenceSize);
  assert(gradsPrevLayer.shape()[1] == batchSize);

  auto cellMapping = getCellMapping(cellOrder);

  // gradient of cell state for this step
  Array2d gradCellState(boost::extents[batchSize][outputSize]);
  for (auto it = gradCellState.data(),
            end = gradCellState.data() + gradCellState.num_elements();
       it != end; ++it) {
    *it = 0;
  }

  // gradient of output of this step
  Array2d gradOutput(boost::extents[batchSize][outputSize]);
  for (auto it = gradCellState.data(),
            end = gradCellState.data() + gradCellState.num_elements();
       it != end; ++it) {
    *it = 0;
  }

  for (auto i = sequenceSize; i != 0; --i) {
    const auto s = i - 1;

    Array2d sumGradOut(boost::extents[batchSize][outputSize]);
    Array2d gradOut = gradsNextLayer[s];
    axpby::add(gradOut, gradOutput, sumGradOut);

    Array2d actOutGate = fwdState[LSTM_FWD_STATE_OUTPUT_GATE][s];
    Array2d gradAtOTanhInp(boost::extents[batchSize][outputSize]);
    gemm::hadamardProduct(actOutGate, sumGradOut, gradAtOTanhInp);

    Array2d actTanhOutGate = fwdState[LSTM_FWD_STATE_OUTPUT_TANH][s];
    Array2d gradAtOutGate(boost::extents[batchSize][outputSize]);
    ;

    gemm::hadamardProduct(actTanhOutGate, sumGradOut, gradAtOutGate);

    bwdNonLinearity(popnn::NonLinearityType::TANH, actTanhOutGate,
                    gradAtOTanhInp);

    bwdNonLinearity(popnn::NonLinearityType::SIGMOID, actOutGate,
                    gradAtOutGate);

    Array2dRef gradAtCellStateSum = gradAtOTanhInp;
    axpby::add(gradAtOTanhInp, gradCellState, gradAtCellStateSum);

    Array2d actInpGate = fwdState[LSTM_FWD_STATE_INPUT_GATE][s];
    Array2d gradAtCand(boost::extents[batchSize][outputSize]);
    ;
    gemm::hadamardProduct(actInpGate, gradAtCellStateSum, gradAtCand);
    Array2d actCand = fwdState[LSTM_FWD_STATE_CAND_TANH][s];
    Array2d gradAtInpGate(boost::extents[batchSize][outputSize]);
    ;
    gemm::hadamardProduct(actCand, gradAtCellStateSum, gradAtInpGate);
    bwdNonLinearity(popnn::NonLinearityType::TANH, actCand, gradAtCand);
    bwdNonLinearity(popnn::NonLinearityType::SIGMOID, actInpGate,
                    gradAtInpGate);

    Array2d actForgetGate = fwdState[LSTM_FWD_STATE_FORGET_GATE][s];
    gemm::hadamardProduct(actForgetGate, gradAtCellStateSum, gradCellState);

    Array2d pCellAct(boost::extents[batchSize][outputSize]);

    if (s == 0) {
      pCellAct = prevCellState;
    } else {
      pCellAct = fwdState[LSTM_FWD_STATE_CELL_STATE_IDX][s - 1];
    }
    Array2d gradAtForgetGate(boost::extents[batchSize][outputSize]);
    ;

    gemm::hadamardProduct(pCellAct, gradAtCellStateSum, gradAtForgetGate);
    bwdNonLinearity(popnn::NonLinearityType::SIGMOID, actForgetGate,
                    gradAtForgetGate);

    Array2d gradIn(boost::extents[batchSize][inputSize]);
    ;
    Array2d weightsInUnit =
        weightsInput[cellMapping.at(BASIC_LSTM_CELL_FORGET_GATE)];
    Array2d weightsOutUnit =
        weightsOutput[cellMapping.at(BASIC_LSTM_CELL_FORGET_GATE)];
    computeGradients(weightsInUnit, weightsOutUnit, gradAtForgetGate, gradIn,
                     gradOutput, false);
    weightsInUnit = weightsInput[cellMapping.at(BASIC_LSTM_CELL_INPUT_GATE)];
    weightsOutUnit = weightsOutput[cellMapping.at(BASIC_LSTM_CELL_INPUT_GATE)];
    computeGradients(weightsInUnit, weightsOutUnit, gradAtInpGate, gradIn,
                     gradOutput, true);
    weightsInUnit = weightsInput[cellMapping.at(BASIC_LSTM_CELL_OUTPUT_GATE)];
    weightsOutUnit = weightsOutput[cellMapping.at(BASIC_LSTM_CELL_OUTPUT_GATE)];
    computeGradients(weightsInUnit, weightsOutUnit, gradAtOutGate, gradIn,
                     gradOutput, true);
    weightsInUnit = weightsInput[cellMapping.at(BASIC_LSTM_CELL_CANDIDATE)];
    weightsOutUnit = weightsOutput[cellMapping.at(BASIC_LSTM_CELL_CANDIDATE)];
    computeGradients(weightsInUnit, weightsOutUnit, gradAtCand, gradIn,
                     gradOutput, true);

    gradsPrevLayer[s] = gradIn;

    // save bwd state for weight update
    bwdState[BASIC_LSTM_CELL_FORGET_GATE][s] = gradAtForgetGate;
    bwdState[BASIC_LSTM_CELL_INPUT_GATE][s] = gradAtInpGate;
    bwdState[BASIC_LSTM_CELL_OUTPUT_GATE][s] = gradAtOutGate;
    bwdState[BASIC_LSTM_CELL_CANDIDATE][s] = gradAtCand;
  }
}

void poplibs_test::lstm::basicLstmCellParamUpdate(
    const Array3dRef prevLayerActs, const Array4dRef fwdState,
    const Array2dRef outputActsInit, const Array4dRef bwdState,
    Array3dRef weightsInputDeltas, Array3dRef weightsOutputDeltas,
    Array2dRef biasDeltas, const std::vector<BasicLstmCellUnit> &cellOrder) {
  const auto sequenceSize = prevLayerActs.shape()[0];
  const auto batchSize = prevLayerActs.shape()[1];
  const auto inputSize = prevLayerActs.shape()[2];
  const auto outputSize = fwdState.shape()[3];

  assert(fwdState.shape()[0] == LSTM_NUM_FWD_STATES);
  assert(fwdState.shape()[1] == sequenceSize);
  assert(fwdState.shape()[2] == batchSize);
  assert(outputActsInit.shape()[0] == batchSize);
  assert(outputActsInit.shape()[1] == outputSize);
  assert(bwdState.shape()[0] == LSTM_NUM_BWD_STATES);
  assert(bwdState.shape()[1] == sequenceSize);
  assert(bwdState.shape()[2] == batchSize);
  assert(bwdState.shape()[3] == outputSize);
  assert(weightsInputDeltas.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsInputDeltas.shape()[1] == inputSize);
  assert(weightsInputDeltas.shape()[2] == outputSize);
  assert(weightsOutputDeltas.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(weightsOutputDeltas.shape()[1] == outputSize);
  assert(weightsOutputDeltas.shape()[2] == outputSize);
  assert(biasDeltas.shape()[0] == BASIC_LSTM_CELL_NUM_UNITS);
  assert(biasDeltas.shape()[1] == outputSize);

  auto cellMapping = getCellMapping(cellOrder);

  for (auto it = weightsInputDeltas.data(),
            end = weightsInputDeltas.data() + weightsInputDeltas.num_elements();
       it != end; ++it) {
    *it = 0;
  }
  for (auto it = weightsOutputDeltas.data(),
            end =
                weightsOutputDeltas.data() + weightsOutputDeltas.num_elements();
       it != end; ++it) {
    *it = 0;
  }
  for (auto it = biasDeltas.data(),
            end = biasDeltas.data() + biasDeltas.num_elements();
       it != end; ++it) {
    *it = 0;
  }

  for (auto i = sequenceSize; i != 0; --i) {
    const auto s = i - 1;
    Array2d outActs(boost::extents[batchSize][outputSize]);
    if (s == 0) {
      outActs = outputActsInit;
    } else {
      outActs = fwdState[LSTM_FWD_STATE_ACTS_IDX][s - 1];
    }
    Array2d inActs = prevLayerActs[s];
    for (auto i = 0; i != BASIC_LSTM_CELL_NUM_UNITS; ++i) {
      const auto unit = static_cast<BasicLstmCellUnit>(i);

      Array2d grad = bwdState[i][s];
      Array2d wInputDeltasUnit(boost::extents[inputSize][outputSize]);

      gemm::generalMatrixMultiply(inActs, grad, wInputDeltasUnit,
                                  wInputDeltasUnit, 1.0, 0, true, false);
      for (auto ic = 0u; ic != inputSize; ++ic) {
        for (auto oc = 0u; oc != outputSize; ++oc) {
          weightsInputDeltas[cellMapping.at(unit)][ic][oc] +=
              wInputDeltasUnit[ic][oc];
        }
      }
      Array2d wOutputDeltasUnit(boost::extents[outputSize][outputSize]);

      gemm::generalMatrixMultiply(outActs, grad, wOutputDeltasUnit,
                                  wOutputDeltasUnit, 1.0, 0, true, false);
      for (auto oc1 = 0u; oc1 != outputSize; ++oc1) {
        for (auto oc2 = 0u; oc2 != outputSize; ++oc2) {
          weightsOutputDeltas[cellMapping.at(unit)][oc1][oc2] +=
              wOutputDeltasUnit[oc1][oc2];
        }
      }

      for (auto oc = 0u; oc != outputSize; ++oc) {
        for (auto b = 0u; b != batchSize; ++b) {
          biasDeltas[cellMapping.at(unit)][oc] += grad[b][oc];
        }
      }
    }
  }
}