diff --git a/stan/math/rev/mat.hpp b/stan/math/rev/mat.hpp index 82699a2173a..f421d51a51b 100644 --- a/stan/math/rev/mat.hpp +++ b/stan/math/rev/mat.hpp @@ -54,6 +54,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/mat/functor/de_integrator.hpp b/stan/math/rev/mat/functor/de_integrator.hpp new file mode 100644 index 00000000000..757b7fe859f --- /dev/null +++ b/stan/math/rev/mat/functor/de_integrator.hpp @@ -0,0 +1,147 @@ +/* + * The original code on which this file is based is from John Cook, + * which is licensed under the 2-clause BSD license, reproduced below: + * + * Copyright (c) 2015, John D. Cook + * 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. + * + * 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 STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_HPP +#define STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_HPP + +#include +#include +#include + +namespace stan { + +namespace math { + +/** + * Double Exponential Integrator. + * + * @tparam T Type of f. + * @param f a functor with signature double (double). + * @param a lower limit of integration, must be double type. + * @param b upper limit of integration, must be double type. + * @param tre target relative error. + * @param tae target absolute error. + * @return numeric integral of function f. + */ +template +inline double de_integrator(const F& f, double a, double b, double tre, + double tae) { + using std::fabs; + using std::log; + + // Apply the linear change of variables x = ct + d + // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$ + // c = (b-a)/2, d = (a+b)/2 + double c = 0.5 * (b - a); + double d = 0.5 * (a + b); + int num_function_evaluations; + + tae /= c; + + // Offsets to where each level's integration constants start. + // The last element is not a beginning but an end. + static int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193}; + int num_levels = sizeof(offsets) / sizeof(int) - 1; + + double new_contribution = 0.0; + double integral = 0.0; + double previous_integral = 0.0; + double error_estimate = DBL_MAX; + double h = 1.0; + double previous_delta, current_delta = DBL_MAX; + + integral = f(c * de_abcissas[0] + d) * de_weights[0]; + + int i; + for (i = offsets[0]; i != offsets[1]; ++i) + integral += de_weights[i] + * (f(c * de_abcissas[i] + d) + f(-c * de_abcissas[i] + d)); + + for (int level = 1; level != num_levels; ++level) { + h *= 0.5; + new_contribution = 0.0; + for (i = offsets[level]; i != offsets[level + 1]; ++i) + new_contribution + += de_weights[i] + * (f(c * de_abcissas[i] + d) + f(-c * de_abcissas[i] + d)); + new_contribution *= h; + + // difference in consecutive integral estimates + previous_delta = current_delta; + current_delta = fabs(0.5 * integral - new_contribution); + previous_integral = integral; + integral = 0.5 * integral + new_contribution; + + // Once convergence kicks in, error is approximately squared + // at each step. + // Determine whether we're in the convergent region by looking + // at the trend in the error. + if (level == 1) + // previous_delta meaningless, so cannot check + // convergence. + continue; + + // Exact comparison with zero is harmless here. Could possibly + // be replaced with a small positive upper limit on the size + // of current_delta, but determining that upper limit would be + // difficult. At worse, the loop is executed more times than + // necessary. But no infinite loop can result since there is + // an upper bound on the loop variable. + if (current_delta == 0.0) + break; + // previous_delta != 0 or would have been kicked out + // previously + double r = log(current_delta) / log(previous_delta); + + if (r > 1.9 && r < 2.1) { + // If convergence theory applied perfectly, r would be 2 in + // the convergence region. r close to 2 is good enough. We + // expect the difference between this integral estimate and + // the next one to be roughly delta^2. + error_estimate = current_delta * current_delta; + } else { + // Not in the convergence region. Assume only that error is + // decreasing. + error_estimate = current_delta; + } + + if (error_estimate < 0.1 * tae + && fabs(1 - integral / previous_integral) < tre) + break; + } + + num_function_evaluations = 2 * i - 1; + error_estimate *= c; + return c * integral; +} +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/mat/functor/de_integrator_constants.hpp b/stan/math/rev/mat/functor/de_integrator_constants.hpp new file mode 100644 index 00000000000..75680ffeb2e --- /dev/null +++ b/stan/math/rev/mat/functor/de_integrator_constants.hpp @@ -0,0 +1,202 @@ +/* + * The original code on which this file is based is from John Cook, + * which is licensed under the 2-clause BSD license, reproduced below: + * + * Copyright (c) 2015, John D. Cook + * 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. + * + * 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 STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_CONSTANTS_HPP +#define STAN_MATH_REV_MAT_FUNCTOR_DEINTEGRATOR_CONSTANTS_HPP + +namespace stan { + +namespace math { + +static const double de_abcissas[] = { + // 1st layer abcissas: transformed 0, 1, 2, 3 + 0.00000000000000000000, 0.95136796407274694573, 0.99997747719246159286, + 0.99999999999995705839, + // 2nd layer abcissas: transformed 1/2, 3/2, 5/2 + 0.67427149224843582608, 0.99751485645722438683, 0.99999998887566488198, + // 3rd layer abcissas: transformed 1/4, 3/4, ... + 0.37720973816403417379, 0.85956905868989663517, 0.98704056050737689169, + 0.99968826402835320905, 0.99999920473711471266, 0.99999999995285644818, + // 4th layer abcissas: transformed 1/8, 3/8, ... + 0.19435700332493543161, 0.53914670538796776905, 0.78060743898320029925, + 0.91487926326457461091, 0.97396686819567744856, 0.99405550663140214329, + 0.99906519645578584642, 0.99990938469514399984, 0.99999531604122052843, + 0.99999989278161241838, 0.99999999914270509218, 0.99999999999823216531, + // 5th layer abcissa: transformed 1/16, 3/16, ... + 0.097923885287832333262, 0.28787993274271591456, 0.46125354393958570440, + 0.61027365750063894488, 0.73101803479256151149, 0.82331700550640237006, + 0.88989140278426019808, 0.93516085752198468323, 0.96411216422354729193, + 0.98145482667733517003, 0.99112699244169880223, 0.99610866543750854254, + 0.99845420876769773751, 0.99945143443527460584, 0.99982882207287494166, + 0.99995387100562796075, 0.99998948201481850361, 0.99999801714059543208, + 0.99999969889415261122, 0.99999996423908091534, 0.99999999678719909830, + 0.99999999978973286224, 0.99999999999039393352, 0.99999999999970809734, + // 6th layer abcissas: transformed 1/32, 3/32, ... + 0.049055967305077886315, 0.14641798429058794053, 0.24156631953888365838, + 0.33314226457763809244, 0.41995211127844715849, 0.50101338937930910152, + 0.57558449063515165995, 0.64317675898520470128, 0.70355000514714201566, + 0.75669390863372994941, 0.80279874134324126576, 0.84221924635075686382, + 0.87543539763040867837, 0.90301328151357387064, 0.92556863406861266645, + 0.94373478605275715685, 0.95813602271021369012, 0.96936673289691733517, + 0.97797623518666497298, 0.98445883116743083087, 0.98924843109013389601, + 0.99271699719682728538, 0.99517602615532735426, 0.99688031812819187372, + 0.99803333631543375402, 0.99879353429880589929, 0.99928111192179195541, + 0.99958475035151758732, 0.99976797159956083506, 0.99987486504878034648, + 0.99993501992508242369, 0.99996759306794345976, 0.99998451990227082442, + 0.99999293787666288565, 0.99999693244919035751, 0.99999873547186590954, + 0.99999950700571943689, 0.99999981889371276701, 0.99999993755407837378, + 0.99999997987450320175, 0.99999999396413420165, 0.99999999832336194826, + 0.99999999957078777261, 0.99999999989927772326, 0.99999999997845533741, + 0.99999999999582460688, 0.99999999999927152627, 0.99999999999988636130, + // 7th layer abcissas: transformed 1/64, 3/64, ... + 0.024539763574649160379, 0.073525122985671294475, 0.12222912220155764235, + 0.17046797238201051811, 0.21806347346971200463, 0.26484507658344795046, + 0.31065178055284596083, 0.35533382516507453330, 0.39875415046723775644, + 0.44078959903390086627, 0.48133184611690504422, 0.52028805069123015958, + 0.55758122826077823080, 0.59315035359195315880, 0.62695020805104287950, + 0.65895099174335012438, 0.68913772506166767176, 0.71750946748732412721, + 0.74407838354734739913, 0.76886868676824658459, 0.79191549237614211447, + 0.81326360850297385168, 0.83296629391941087564, 0.85108400798784873261, + 0.86768317577564598669, 0.88283498824466895513, 0.89661425428007602579, + 0.90909831816302043511, 0.92036605303195280235, 0.93049693799715340631, + 0.93957022393327475539, 0.94766419061515309734, 0.95485549580502268541, + 0.96121861515111640753, 0.96682537031235585284, 0.97174454156548730892, + 0.97604156025657673933, 0.97977827580061576265, 0.98301279148110110558, + 0.98579936302528343597, 0.98818835380074264243, 0.99022624046752774694, + 0.99195566300267761562, 0.99341551316926403900, 0.99464105571251119672, + 0.99566407681695316965, 0.99651305464025377317, 0.99721334704346870224, + 0.99778739195890653083, 0.99825491617199629344, 0.99863314864067747762, + 0.99893703483351217373, 0.99917944893488591716, 0.99937140114093768690, + 0.99952223765121720422, 0.99963983134560036519, 0.99973076151980848263, + 0.99980048143113838630, 0.99985347277311141171, 0.99989338654759256426, + 0.99992317012928932869, 0.99994518061445869309, 0.99996128480785666613, + 0.99997294642523223656, 0.99998130127012072679, 0.99998722128200062811, + 0.99999136844834487344, 0.99999423962761663478, 0.99999620334716617675, + 0.99999752962380516793, 0.99999841381096473542, 0.99999899541068996962, + 0.99999937270733536947, 0.99999961398855024275, 0.99999976602333243312, + 0.99999986037121459941, 0.99999991800479471056, 0.99999995264266446185, + 0.99999997311323594362, 0.99999998500307631173, 0.99999999178645609907, + 0.99999999558563361584, 0.99999999767323673790, 0.99999999879798350040, + 0.99999999939177687583, 0.99999999969875436925, 0.99999999985405611550, + 0.99999999993088839501, 0.99999999996803321674, 0.99999999998556879008, + 0.99999999999364632387, 0.99999999999727404948, 0.99999999999886126543, + 0.99999999999953722654, 0.99999999999981720098, 0.99999999999992987953}; + +static const double de_weights[] = { + // First layer weights + 1.5707963267948966192, 0.230022394514788685, 0.00026620051375271690866, + 1.3581784274539090834e-12, + // 2nd layer weights + 0.96597657941230114801, 0.018343166989927842087, 2.1431204556943039358e-7, + // 3rd layer weights + 1.3896147592472563229, 0.53107827542805397476, 0.076385743570832304188, + 0.0029025177479013135936, 0.000011983701363170720047, + 1.1631165814255782766e-9, + // 4th layer weights + 1.5232837186347052132, 1.1934630258491569639, 0.73743784836154784136, + 0.36046141846934367417, 0.13742210773316772341, 0.039175005493600779072, + 0.0077426010260642407123, 0.00094994680428346871691, + 0.000062482559240744082891, 1.8263320593710659699e-6, + 1.8687282268736410132e-8, 4.9378538776631926964e-11, + // 5th layer weights + 1.5587733555333301451, 1.466014426716965781, 1.297475750424977998, + 1.0816349854900704074, 0.85017285645662006895, 0.63040513516474369106, + 0.44083323627385823707, 0.290240679312454185, 0.17932441211072829296, + 0.10343215422333290062, 0.055289683742240583845, 0.027133510013712003219, + 0.012083543599157953493, 0.0048162981439284630173, 0.0016908739981426396472, + 0.00051339382406790336017, 0.00013205234125609974879, + 0.000028110164327940134749, 4.8237182032615502124e-6, + 6.4777566035929719908e-7, 6.5835185127183396672e-8, + 4.8760060974240625869e-9, 2.5216347918530148572e-10, + 8.6759314149796046502e-12, + // 6th layer weights + 1.5677814313072218572, 1.5438811161769592204, 1.4972262225410362896, + 1.4300083548722996676, 1.3452788847662516615, 1.2467012074518577048, + 1.1382722433763053734, 1.0240449331118114483, 0.90787937915489531693, + 0.79324270082051671787, 0.68306851634426375464, 0.57967810308778764708, + 0.48475809121475539287, 0.39938474152571713515, 0.32408253961152890402, + 0.258904639514053516, 0.20352399885860174519, 0.15732620348436615027, + 0.11949741128869592428, 0.089103139240941462841, 0.065155533432536205042, + 0.046668208054846613644, 0.032698732726609031113, 0.022379471063648476483, + 0.014937835096050129696, 0.0097072237393916892692, 0.0061300376320830301252, + 0.0037542509774318343023, 0.0022250827064786427022, + 0.0012733279447082382027, 0.0007018595156842422708, + 0.00037166693621677760301, 0.00018856442976700318572, + 0.000091390817490710122732, 0.000042183183841757600604, + 0.000018481813599879217116, 7.6595758525203162562e-6, + 2.9916615878138787094e-6, 1.0968835125901264732e-6, + 3.7595411862360630091e-7, 1.1992442782902770219e-7, + 3.5434777171421953043e-8, 9.6498888961089633609e-9, + 2.4091773256475940779e-9, 5.482835779709497755e-10, + 1.1306055347494680536e-10, 2.0989335404511469109e-11, + 3.4841937670261059685e-12, + // 7th layer weights + 1.5700420292795931467, 1.5640214037732320999, 1.5520531698454121192, + 1.5342817381543034316, 1.5109197230741697127, 1.48224329788553807, + 1.4485862549613225916, 1.4103329714462590129, 1.3679105116808964881, + 1.3217801174437728579, 1.2724283455378627082, 1.2203581095793582207, + 1.1660798699324345766, 1.1101031939653403796, 1.0529288799552666556, + 0.99504180404613271514, 0.93690461274566793366, 0.87895234555278212039, + 0.82158803526696470334, 0.7651792989089561367, 0.71005590120546898385, + 0.65650824613162753076, 0.60478673057840362158, 0.55510187800363350959, + 0.5076251588319080997, 0.4624903980553677613, 0.41979566844501548066, + 0.37960556938665160999, 0.3419537959230168323, 0.30684590941791694932, + 0.27426222968906810637, 0.24416077786983990868, 0.21648020911729617038, + 0.19114268413342749532, 0.16805663794826916233, 0.14711941325785693248, + 0.12821973363120098675, 0.11123999898874453035, 0.096058391865189467849, + 0.082550788110701737654, 0.070592469906866999352, 0.060059642358636300319, + 0.05083075757257047107, 0.042787652157725676034, 0.035816505604196436523, + 0.029808628117310126969, 0.024661087314753282511, 0.020277183817500123926, + 0.016566786254247575375, 0.013446536605285730674, 0.010839937168255907211, + 0.0086773307495391815854, 0.0068957859690660035329, + 0.0054388997976239984331, 0.0042565295990178580165, + 0.0033044669940348302363, 0.0025440657675291729678, + 0.0019418357759843675814, 0.0014690143599429791058, + 0.0011011261134519383862, 0.00081754101332469493115, + 0.00060103987991147422573, 0.00043739495615911687786, + 0.00031497209186021200274, 0.00022435965205008549104, + 0.00015802788400701191949, 0.00011002112846666697224, + 0.000075683996586201477788, 0.000051421497447658802092, + 0.0000344921247593431977, 0.000022832118109036146591, + 0.000014908514031870608449, 9.5981941283784710776e-6, + 6.0899100320949039256e-6, 3.8061983264644899045e-6, + 2.3421667208528096843e-6, 1.4183067155493917523e-6, + 8.4473756384859863469e-7, 4.9458288702754198508e-7, + 2.8449923659159806339e-7, 1.6069394579076224911e-7, + 8.9071395140242387124e-8, 4.8420950198072369669e-8, 2.579956822953589238e-8, + 1.3464645522302038796e-8, 6.8784610955899001111e-9, + 3.4371856744650090511e-9, 1.6788897682161906807e-9, + 8.0099784479729665356e-10, 3.7299501843052790038e-10, + 1.6939457789411646876e-10, 7.4967397573818224522e-11, + 3.230446433325236576e-11, 1.3542512912336274432e-11, + 5.5182369468174885821e-12, 2.1835922099233609052e-12}; +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/mat/functor/integrate_1d_tsc.hpp b/stan/math/rev/mat/functor/integrate_1d_tsc.hpp new file mode 100644 index 00000000000..da1b19ce69d --- /dev/null +++ b/stan/math/rev/mat/functor/integrate_1d_tsc.hpp @@ -0,0 +1,191 @@ +#ifndef STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_tsc_HPP +#define STAN_MATH_PRIM_ARR_FUNCTOR_integrate_1d_tsc_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + +namespace math { + +/** + * Return the numerical integral of a function f given its gradient + * g. + * + * This function uses the algorithm of integrate_1d_tscg for + * numerical integration. + * + * @tparam T Type of f. + * @tparam G Type of g. + * @param f a functor with signature + * double (double, std::vector) or with signature + * double (double, T_param) where the first argument is one being + * integrated and the second one is either an extra scalar or vector + * being passed to f. + * @param g a functor with signature + * double (double, std::vector, int, std::ostream&) or with + * signature double (double, T_param, int, std::ostream&) where the + * first argument is being integrated and the second one is + * either an extra scalar or vector being passed to f and the + * third one selects which component of the gradient vector + * is to be returned. + * @param a lower limit of integration, must be a finite double. + * @param b upper limit of integration, must be a finite double. + * @param tre target relative error. + * @param tae target absolute error. + * @param param additional parameters to be passed to f. + * @param msgs stream. + * @return numeric integral of function f. + */ +template +inline typename scalar_type::type integrate_1d_tscg( + const F& f, const G& g, const double a, const double b, + const T_param& param, std::ostream& msgs, const double tre = 1e-6, + const double tae = 1e-6) { + check_finite("integrate_1d_tsc", "lower limit", a); + check_finite("integrate_1d_tsc", "upper limit", b); + + using std::placeholders::_1; + + // hard case, we want a normalizing factor + if (!is_constant_struct::value) { + size_t N = length(param); + std::vector grad(N); + + auto value_of_param = value_of(param); + + for (size_t n = 0; n < N; n++) + grad[n] = de_integrator( + std::bind(g, _1, value_of_param, static_cast(n + 1), + std::ref(msgs)), + a, b, tre, tae); + + double val_ = de_integrator( + std::bind(f, _1, value_of_param, std::ref(msgs)), a, b, tre, + tae); + + operands_and_partials ops_partials(param); + for (size_t n = 0; n < N; n++) + ops_partials.edge1_.partials_[n] += grad[n]; + + return ops_partials.build(val_); + // easy case, here we are calculating a normalizing constant, + // not a normalizing factor, so g doesn't matter at all + } + return de_integrator( + std::bind(f, _1, value_of(param), std::ref(msgs)), a, b, tre, + tae); +} + +/** + * Calculate gradient of f(x, param, std::ostream&) + * with respect to param_n (which must be an element of param) + */ +template +inline double gradient_of_f(const F& f, const double x, const T_param& param, + const var& param_n, std::ostream& msgs) { + set_zero_all_adjoints_nested(); + f(x, param, msgs).grad(); + return param_n.adj(); +} + +/** + * Return the numerical integral of a function f with its + * gradients being inferred automatically (but slowly). + * + * The numerical integration algorithm used is the Tanh-sinh + * quadrature method (also known as Double Exponential + * Transformation) which was proposed by Hidetosi Takahasi and + * Masatake Mori in 1974 + * + * The implementation of integration used is given John D. Cook with + * the slight modification of adding a relative tolerance error. See + * www.codeproject.com/kb/recipes/fastnumericalintegration.aspx + * www.johndcook.com/blog/double_exponential_integration/ + * + * Such implementation assumes f to be smooth with no discontinuity + * in the function nor in any of its derivatives. + * + * The integration is terminated when both relative and absolute + * tolerance are reached or when a predefined number of iterations + * is reached (such termination condition was designed to target + * floating point). + * + * @tparam T Type of f. + * @param f a functor with signature + * double (double, std::vector, std::ostream&) or with + * signature double (double, T_param, std::ostream&) where the first + * argument is one being integrated and the second one is either + * an extra scalar or vector being passed to f. + * @param a lower limit of integration, must be double type. + * @param b upper limit of integration, must be double type. + * @param tre target relative error. + * @param tae target absolute error. + * @param param additional parameters to be passed to f. + * @param msgs stream. + * @return numeric integral of function f. + */ +template +inline typename scalar_type::type integrate_1d_tsc( + const F& f, const double a, const double b, const T_param& param, + std::ostream& msgs, const double tre = 1e-6, const double tae = 1e-6) { + using std::placeholders::_1; + + stan::math::check_finite("integrate_1d_tsc", "lower limit", a); + stan::math::check_finite("integrate_1d_tsc", "upper limit", b); + + double val_ + = de_integrator(std::bind(f, _1, value_of(param), std::ref(msgs)), + a, b, tre, tae); + + if (!is_constant_struct::value) { + size_t N = stan::length(param); + std::vector results(N); + + try { + start_nested(); + + auto clean_param = to_var(value_of(param)); + typedef decltype(clean_param) clean_T_param; + + scalar_seq_view clean_param_vec(clean_param); + + for (size_t n = 0; n < N; n++) + results[n] = de_integrator( + std::bind(gradient_of_f, f, _1, + clean_param, clean_param_vec[n], std::ref(msgs)), + a, b, tre, tae); + } catch (const std::exception& e) { + recover_memory_nested(); + throw; + } + recover_memory_nested(); + + operands_and_partials ops_partials(param); + for (size_t n = 0; n < N; n++) + ops_partials.edge1_.partials_[n] += results[n]; + + return ops_partials.build(val_); + } else { + return val_; + } +} + +} // namespace math + +} // namespace stan + +#endif diff --git a/test/unit/math/prim/arr/functor/integrate_1d_tsc_test.cpp b/test/unit/math/prim/arr/functor/integrate_1d_tsc_test.cpp new file mode 100644 index 00000000000..d1da7da8d79 --- /dev/null +++ b/test/unit/math/prim/arr/functor/integrate_1d_tsc_test.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include +#include + +struct f1 { + template + inline typename stan::return_type::type operator()( + const T1& x, const T2& y, std::ostream& msgs) const { + return exp(x) + y; + } +}; + +struct f2 { + template + inline typename stan::return_type::type operator()( + const T1& x, const std::vector& y, std::ostream& msgs) const { + return exp(x) + pow(y[0], 2) + pow(y[1], 3); + } +}; + +struct f3 { + template + inline typename stan::return_type::type operator()( + const T1& x, const Eigen::Matrix& y, + std::ostream& msgs) const { + return exp(x) + pow(y(0), 2) + pow(y(1), 4) + 3 * y(2); + } +}; + +struct f4 { + template + inline typename stan::return_type::type operator()( + const T1& x, const Eigen::Matrix& y, + std::ostream& msgs) const { + return exp(x) + pow(y(0), 2) + pow(y(1), 5) + 3 * y(2); + } +}; + +std::ostringstream msgs; + +TEST(StanMath_integrate_1d_tsc, test1) { + using stan::math::integrate_1d_tsc; + + f1 if1; + + EXPECT_FLOAT_EQ(integrate_1d_tsc(if1, .2, .7, .5, msgs), 1.04235); + + f2 if2; + + EXPECT_FLOAT_EQ( + integrate_1d_tsc(if2, -.2, .7, std::vector(2, .4), msgs), + 1.396622); + f3 if3; + Eigen::VectorXd a(3); + a(0) = 4.0; + a(1) = 6.0; + a(2) = 5.1; + + EXPECT_FLOAT_EQ(integrate_1d_tsc(if3, -.2, 2.9, a, msgs), 4131.985); + + f4 if4; + Eigen::RowVectorXd b(3); + b(0) = 4.0; + b(1) = 6.0; + b(2) = 5.1; + + EXPECT_FLOAT_EQ(integrate_1d_tsc(if4, -.2, 2.9, b, msgs), 24219.99); +} diff --git a/test/unit/math/rev/arr/functor/integrate_1d_tsc_test.cpp b/test/unit/math/rev/arr/functor/integrate_1d_tsc_test.cpp new file mode 100644 index 00000000000..1573e6d721d --- /dev/null +++ b/test/unit/math/rev/arr/functor/integrate_1d_tsc_test.cpp @@ -0,0 +1,123 @@ +#include +#include +#include +#include +#include +#include +#include + +struct f1 { + template + inline typename stan::return_type::type operator()( + const T1& x, const T2& y, std::ostream& msgs) const { + return exp(x) + y; + } +}; + +struct f2 { + template + inline typename stan::return_type::type operator()( + const T1& x, const T2& y, std::ostream& msgs) const { + return exp(y * cos(2 * 3.141593 * x)) + y; + } +}; + +struct f3 { + template + inline typename stan::return_type::type operator()( + const T1& x, const std::vector& y, std::ostream& msgs) const { + return exp(x) + pow(y[0], 2.5) + 2 * pow(y[1], 3) + 2 * y[2]; + } +}; + +struct g3 { + template + inline typename stan::return_type::type operator()( + const T1& x, const std::vector& y, const int ii, + std::ostream& msgs) const { + if (ii == 1) + return 2.5 * pow(y[0], 1.5); + else if (ii == 2) + return 6 * pow(y[1], 2.); + else + return 2.; + } +}; + +std::ostringstream msgs; + +TEST(StanMath_integrate_1d_tsc, test1) { + using stan::math::integrate_1d_tsc; + + f1 if1; + + EXPECT_FLOAT_EQ( + integrate_1d_tsc(if1, .2, .7, stan::math::var(.5), msgs).val(), + 0.7923499 + .25); +} + +TEST(StanMath_integrate_1d_tsc, finite_diff) { + using stan::math::integrate_1d_tsc; + + { + f1 if1; + + AVAR a = .6; + AVAR f = integrate_1d_tsc(if1, .2, .7, a, msgs); + EXPECT_FLOAT_EQ(integrate_1d_tsc(if1, .2, .7, .6, msgs), f.val()); + + AVEC x = createAVEC(a); + VEC g; + f.grad(x, g); + + EXPECT_FLOAT_EQ((integrate_1d_tsc(if1, .2, .7, .6 + 1e-6, msgs) + - integrate_1d_tsc(if1, .2, .7, .6 - 1e-6, msgs)) + / 2e-6, + g[0]); + } + { + f2 if2; + + AVAR a = 0.68; + AVAR f = integrate_1d_tsc(if2, 0., 1.1, a, msgs); + EXPECT_FLOAT_EQ(integrate_1d_tsc(if2, 0., 1.1, .68, msgs), f.val()); + + AVEC x = createAVEC(a); + VEC g; + f.grad(x, g); + + EXPECT_FLOAT_EQ((integrate_1d_tsc(if2, 0., 1.1, .68 + 1e-6, msgs) + - integrate_1d_tsc(if2, 0., 1.1, .68 - 1e-6, msgs)) + / 2e-6, + g[0]); + } + { + f3 if3; + + AVAR a = 0.68; + AVAR b = 0.38; + AVAR c = 0.78; + AVEC vec = createAVEC(a, b, c); + AVAR f = integrate_1d_tsc(if3, 0., 1.1, vec, msgs); + + VEC g; + double p1; + double p2; + f.grad(vec, g); + + std::vector vecd = value_of(vec); + EXPECT_FLOAT_EQ(integrate_1d_tsc(if3, 0., 1.1, vecd, msgs), f.val()); + + vecd[0] += 1e-6; + p1 = integrate_1d_tsc(if3, 0., 1.1, vecd, msgs); + vecd[0] -= 2e-6; + p2 = integrate_1d_tsc(if3, 0., 1.1, vecd, msgs); + + EXPECT_FLOAT_EQ((p1 - p2) / 2e-6, g[0]); + + g3 ig3; + stan::math::set_zero_all_adjoints(); + integrate_1d_tscg(if3, ig3, 0., 1.1, vec, msgs).grad(); + EXPECT_FLOAT_EQ((p1 - p2) / 2e-6, a.adj()); + } +}