diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 4a2f4287ac..2d8abfabfd 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -153,6 +153,8 @@ set(lisbslic3r_sources GCode/PrintExtents.hpp GCode/RetractWhenCrossingPerimeters.cpp GCode/RetractWhenCrossingPerimeters.hpp + GCode/SmallAreaInfillFlowCompensator.cpp + GCode/SmallAreaInfillFlowCompensator.hpp GCode/SpiralVase.cpp GCode/SpiralVase.hpp GCode/SeamPlacer.cpp diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index a02abefd94..415417ea91 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -1965,6 +1965,9 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato } else m_enable_extrusion_role_markers = false; + if (!print.config().small_area_infill_flow_compensation_model.empty()) + m_small_area_infill_flow_compensator = make_unique(print.config()); + // if thumbnail type of BTT_TFT, insert above header // if not, it is inserted under the header in its normal spot const GCodeThumbnailsFormat m_gcode_thumbnail_format = print.full_print_config().opt_enum("thumbnails_format"); @@ -5190,15 +5193,25 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, for (const Line& line : path.polyline.lines()) { const double line_length = line.length() * SCALING_FACTOR; path_length += line_length; + auto dE = e_per_mm * line_length; + if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) { + auto oldE = dE; + dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role()); + + if (m_config.gcode_comments && oldE > 0 && oldE != dE) { + description += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length); + } + } gcode += m_writer.extrude_to_xy( this->point_to_gcode(line.b), - e_per_mm * line_length, + dE, GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion()); } } else { // BBS: start to generate gcode from arc fitting data which includes line and arc const std::vector& fitting_result = path.polyline.fitting_result; for (size_t fitting_index = 0; fitting_index < fitting_result.size(); fitting_index++) { + std::string tempDescription = description; switch (fitting_result[fitting_index].path_type) { case EMovePathType::Linear_move: { size_t start_index = fitting_result[fitting_index].start_point_index; @@ -5207,10 +5220,19 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, const Line line = Line(path.polyline.points[point_index - 1], path.polyline.points[point_index]); const double line_length = line.length() * SCALING_FACTOR; path_length += line_length; + auto dE = e_per_mm * line_length; + if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) { + auto oldE = dE; + dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role()); + + if (m_config.gcode_comments && oldE > 0 && oldE != dE) { + tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length); + } + } gcode += m_writer.extrude_to_xy( this->point_to_gcode(line.b), - e_per_mm * line_length, - GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion()); + dE, + GCodeWriter::full_gcode_comment ? tempDescription : "", path.is_force_no_extrusion()); } break; } @@ -5220,12 +5242,21 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, const double arc_length = fitting_result[fitting_index].arc_data.length * SCALING_FACTOR; const Vec2d center_offset = this->point_to_gcode(arc.center) - this->point_to_gcode(arc.start_point); path_length += arc_length; + auto dE = e_per_mm * arc_length; + if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) { + auto oldE = dE; + dE = m_small_area_infill_flow_compensator->modify_flow(arc_length, dE, path.role()); + + if (m_config.gcode_comments && oldE > 0 && oldE != dE) { + tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, arc_length); + } + } gcode += m_writer.extrude_arc_to_xy( this->point_to_gcode(arc.end_point), center_offset, - e_per_mm * arc_length, + dE, arc.direction == ArcDirection::Arc_Dir_CCW, - GCodeWriter::full_gcode_comment ? description : "", path.is_force_no_extrusion()); + GCodeWriter::full_gcode_comment ? tempDescription : "", path.is_force_no_extrusion()); break; } default: @@ -5247,6 +5278,7 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, pre_fan_enabled = check_overhang_fan(new_points[0].overlap, path.role()); for (size_t i = 1; i < new_points.size(); i++) { + std::string tempDescription = description; const ProcessedPoint &processed_point = new_points[i]; const ProcessedPoint &pre_processed_point = new_points[i-1]; Vec2d p = this->point_to_gcode_quantized(processed_point.p); @@ -5285,8 +5317,17 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, gcode += m_writer.set_speed(new_speed, "", comment); last_set_speed = new_speed; } + auto dE = e_per_mm * line_length; + if (m_small_area_infill_flow_compensator && m_config.small_area_infill_flow_compensation.value) { + auto oldE = dE; + dE = m_small_area_infill_flow_compensator->modify_flow(line_length, dE, path.role()); + + if (m_config.gcode_comments && oldE > 0 && oldE != dE) { + tempDescription += Slic3r::format(" | Old Flow Value: %0.5f Length: %0.5f",oldE, line_length); + } + } gcode += - m_writer.extrude_to_xy(p, e_per_mm * line_length, GCodeWriter::full_gcode_comment ? description : ""); + m_writer.extrude_to_xy(p, dE, GCodeWriter::full_gcode_comment ? tempDescription : ""); prev = p; diff --git a/src/libslic3r/GCode.hpp b/src/libslic3r/GCode.hpp index e3b6b65a88..935cb95fc9 100644 --- a/src/libslic3r/GCode.hpp +++ b/src/libslic3r/GCode.hpp @@ -23,6 +23,7 @@ #include "GCode/ExtrusionProcessor.hpp" #include "GCode/PressureEqualizer.hpp" +#include "GCode/SmallAreaInfillFlowCompensator.hpp" #include #include @@ -531,6 +532,8 @@ private: std::unique_ptr m_wipe_tower; + std::unique_ptr m_small_area_infill_flow_compensator; + // Heights (print_z) at which the skirt has already been extruded. std::vector m_skirt_done; // Has the brim been extruded already? Brim is being extruded only for the first object of a multi-object print. @@ -593,6 +596,7 @@ private: friend class WipeTowerIntegration; friend class PressureEqualizer; friend class Print; + friend class SmallAreaInfillFlowCompensator; }; std::vector sort_object_instances_by_model_order(const Print& print, bool init_order = false); diff --git a/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.cpp b/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.cpp new file mode 100644 index 0000000000..573347d5ec --- /dev/null +++ b/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.cpp @@ -0,0 +1,88 @@ +// Modify the flow of extrusion lines inversely proportional to the length of +// the extrusion line. When infill lines get shorter the flow rate will auto- +// matically be reduced to mitigate the effect of small infill areas being +// over-extruded. + +// Based on original work by Alexander Þór licensed under the GPLv3: +// https://github.com/Alexander-T-Moss/Small-Area-Flow-Comp + +#include +#include +#include + +#include "../libslic3r.h" +#include "../PrintConfig.hpp" + +#include "SmallAreaInfillFlowCompensator.hpp" + +namespace Slic3r { + +bool nearly_equal(double a, double b) +{ + return std::nextafter(a, std::numeric_limits::lowest()) <= b && std::nextafter(a, std::numeric_limits::max()) >= b; +} + +SmallAreaInfillFlowCompensator::SmallAreaInfillFlowCompensator(const Slic3r::GCodeConfig& config) +{ + for (auto& line : config.small_area_infill_flow_compensation_model.values) { + std::istringstream iss(line); + std::string value_str; + double eLength = 0.0; + + if (std::getline(iss, value_str, ',')) { + try { + eLength = std::stod(value_str); + if (std::getline(iss, value_str, ',')) { + eLengths.push_back(eLength); + flowComps.push_back(std::stod(value_str)); + } + } catch (...) { + std::stringstream ss; + ss << "Error parsing data point in small area infill compensation model:" << line << std::endl; + + throw Slic3r::InvalidArgument(ss.str()); + } + } + } + + for (int i = 0; i < eLengths.size(); i++) { + if (i == 0) { + if (!nearly_equal(eLengths[i], 0.0)) { + throw Slic3r::InvalidArgument("First extrusion length for small area infill compensation model must be 0"); + } + } else { + if (nearly_equal(eLengths[i], 0.0)) { + throw Slic3r::InvalidArgument("Only the first extrusion length for small area infill compensation model can be 0"); + } + if (eLengths[i] <= eLengths[i - 1]) { + throw Slic3r::InvalidArgument("Extrusion lengths for subsequent points must be increasing"); + } + } + } + + if (!flowComps.empty() && !nearly_equal(flowComps.back(), 1.0)) { + throw Slic3r::InvalidArgument("Final compensation factor for small area infill flow compensation model must be 1.0"); + } + + flowModel.set_points(eLengths, flowComps); +} + +double SmallAreaInfillFlowCompensator::flow_comp_model(const double line_length) +{ + if (line_length == 0 || line_length > max_modified_length()) { + return 1.0; + } + + return flowModel(line_length); +} + +double SmallAreaInfillFlowCompensator::modify_flow(const double line_length, const double dE, const ExtrusionRole role) +{ + if (role == ExtrusionRole::erSolidInfill || role == ExtrusionRole::erTopSolidInfill || role == ExtrusionRole::erBottomSurface) { + return dE * flow_comp_model(line_length); + } + + return dE; +} + +} // namespace Slic3r diff --git a/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.hpp b/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.hpp new file mode 100644 index 0000000000..7d804c9b47 --- /dev/null +++ b/src/libslic3r/GCode/SmallAreaInfillFlowCompensator.hpp @@ -0,0 +1,35 @@ +#ifndef slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_ +#define slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_ + +#include "../libslic3r.h" +#include "../PrintConfig.hpp" +#include "../ExtrusionEntity.hpp" +#include "spline/spline.h" + +namespace Slic3r { + +class SmallAreaInfillFlowCompensator +{ +public: + SmallAreaInfillFlowCompensator() = delete; + explicit SmallAreaInfillFlowCompensator(const Slic3r::GCodeConfig& config); + ~SmallAreaInfillFlowCompensator() = default; + + double modify_flow(const double line_length, const double dE, const ExtrusionRole role); + +private: + // Model points + std::vector eLengths; + std::vector flowComps; + + // TODO: Cubic Spline + tk::spline flowModel; + + double flow_comp_model(const double line_length); + + double max_modified_length() { return eLengths.back(); } +}; + +} // namespace Slic3r + +#endif /* slic3r_GCode_SmallAreaInfillFlowCompensator_hpp_ */ diff --git a/src/libslic3r/Preset.cpp b/src/libslic3r/Preset.cpp index db675a5193..8fe5b2bf33 100644 --- a/src/libslic3r/Preset.cpp +++ b/src/libslic3r/Preset.cpp @@ -817,6 +817,7 @@ static std::vector s_Preset_print_options { "wipe_tower_cone_angle", "wipe_tower_extra_spacing", "wipe_tower_extruder", "wiping_volumes_extruders","wipe_tower_bridging", "single_extruder_multi_material_priming", "wipe_tower_rotation_angle", "tree_support_branch_distance_organic", "tree_support_branch_diameter_organic", "tree_support_branch_angle_organic", "hole_to_polyhole", "hole_to_polyhole_threshold", "hole_to_polyhole_twisted", "mmu_segmented_region_max_width", "mmu_segmented_region_interlocking_depth", + "small_area_infill_flow_compensation", "small_area_infill_flow_compensation_model", }; static std::vector s_Preset_filament_options { diff --git a/src/libslic3r/PrintConfig.cpp b/src/libslic3r/PrintConfig.cpp index 980e1b63db..f1120f88cb 100644 --- a/src/libslic3r/PrintConfig.cpp +++ b/src/libslic3r/PrintConfig.cpp @@ -2668,6 +2668,26 @@ def = this->add("filament_loading_speed", coFloats); def->mode = comAdvanced; def->set_default_value(new ConfigOptionString("")); + def = this->add("small_area_infill_flow_compensation", coBool); + def->label = L("Enable Flow Compensation"); + def->tooltip = L("Enable flow compensation for small infill areas"); + def->mode = comAdvanced; + def->set_default_value(new ConfigOptionBool(false)); + + def = this->add("small_area_infill_flow_compensation_model", coStrings); + def->label = L("Flow Compensation Model"); + def->tooltip = L( + "Flow Compensation Model, used to adjust the flow for small infill " + "areas. The model is expressed as a comma separated pair of values for " + "extrusion length and flow correction factors, one per line, in the " + "following format: \"1.234,5.678\""); + def->mode = comAdvanced; + def->gui_flags = "serialized"; + def->multiline = true; + def->full_width = true; + def->height = 15; + def->set_default_value(new ConfigOptionStrings{"0,0", "\n0.2,0.4444", "\n0.4,0.6145", "\n0.6,0.7059", "\n0.8,0.7619", "\n1.5,0.8571", "\n2,0.8889", "\n3,0.9231", "\n5,0.9520", "\n10,1"}); + { struct AxisDefault { std::string name; diff --git a/src/libslic3r/PrintConfig.hpp b/src/libslic3r/PrintConfig.hpp index c665cee6e9..1fef03f572 100644 --- a/src/libslic3r/PrintConfig.hpp +++ b/src/libslic3r/PrintConfig.hpp @@ -920,6 +920,7 @@ PRINT_CONFIG_CLASS_DEFINE( ((ConfigOptionEnum, wall_sequence)) ((ConfigOptionBool, is_infill_first)) + ((ConfigOptionBool, small_area_infill_flow_compensation)) ) PRINT_CONFIG_CLASS_DEFINE( @@ -1067,6 +1068,8 @@ PRINT_CONFIG_CLASS_DEFINE( ((ConfigOptionBool, enable_filament_ramming)) ((ConfigOptionBool, support_multi_bed_types)) + // Small Area Infill Flow Compensation + ((ConfigOptionStrings, small_area_infill_flow_compensation_model)) ) // This object is mapped to Perl as Slic3r::Config::Print. diff --git a/src/libslic3r/PrintObject.cpp b/src/libslic3r/PrintObject.cpp index ce64205aae..d3a206123f 100644 --- a/src/libslic3r/PrintObject.cpp +++ b/src/libslic3r/PrintObject.cpp @@ -928,6 +928,10 @@ bool PrintObject::invalidate_state_by_config_options( || opt_key == "wipe_on_loops" || opt_key == "wipe_speed") { steps.emplace_back(posPerimeters); + } else if ( + opt_key == "small_area_infill_flow_compensation" + || opt_key == "small_area_infill_flow_compensation_model") { + steps.emplace_back(posSlice); } else if (opt_key == "gap_infill_speed" || opt_key == "filter_out_gap_fill" ) { // Return true if gap-fill speed has changed from zero value to non-zero or from non-zero value to zero. diff --git a/src/slic3r/GUI/ConfigManipulation.cpp b/src/slic3r/GUI/ConfigManipulation.cpp index f6bbcecf09..dcf2c96196 100644 --- a/src/slic3r/GUI/ConfigManipulation.cpp +++ b/src/slic3r/GUI/ConfigManipulation.cpp @@ -742,6 +742,10 @@ void ConfigManipulation::toggle_print_fff_options(DynamicPrintConfig *config, co apply(config, &new_conf); } toggle_line("timelapse_type", is_BBL_Printer); + + + bool have_small_area_infill_flow_compensation = config->opt_bool("small_area_infill_flow_compensation"); + toggle_line("small_area_infill_flow_compensation_model", have_small_area_infill_flow_compensation); } void ConfigManipulation::update_print_sla_config(DynamicPrintConfig* config, const bool is_global_config/* = false*/) diff --git a/src/slic3r/GUI/Tab.cpp b/src/slic3r/GUI/Tab.cpp index 5a67afa514..fac5d6dfec 100644 --- a/src/slic3r/GUI/Tab.cpp +++ b/src/slic3r/GUI/Tab.cpp @@ -1970,6 +1970,14 @@ void TabPrint::build() optgroup->append_single_option_line("only_one_wall_first_layer"); optgroup->append_single_option_line("reduce_crossing_wall"); optgroup->append_single_option_line("max_travel_detour_distance"); + + optgroup = page->new_optgroup(L("Small Area Infill Flow Compensation (experimental)"), L"param_advanced"); + optgroup->append_single_option_line("small_area_infill_flow_compensation"); + Option option = optgroup->get_option("small_area_infill_flow_compensation_model"); + option.opt.full_width = true; + option.opt.is_code = true; + option.opt.height = 15; + optgroup->append_single_option_line(option); optgroup = page->new_optgroup(L("Bridging"), L"param_advanced"); optgroup->append_single_option_line("bridge_flow"); @@ -2200,7 +2208,7 @@ void TabPrint::build() optgroup->append_single_option_line("gcode_comments"); optgroup->append_single_option_line("gcode_label_objects"); optgroup->append_single_option_line("exclude_object"); - Option option = optgroup->get_option("filename_format"); + option = optgroup->get_option("filename_format"); // option.opt.full_width = true; option.opt.is_code = true; option.opt.multiline = true; diff --git a/src/spline/spline.h b/src/spline/spline.h new file mode 100644 index 0000000000..c8f08418fb --- /dev/null +++ b/src/spline/spline.h @@ -0,0 +1,951 @@ +/* + * spline.h + * + * simple cubic spline interpolation library without external + * dependencies + * + * --------------------------------------------------------------------- + * Copyright (C) 2011, 2014, 2016, 2021 Tino Kluge (ttk448 at gmail.com) + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * --------------------------------------------------------------------- + * + */ + + +#ifndef TK_SPLINE_H +#define TK_SPLINE_H + +#include +#include +#include +#include +#include +#ifdef HAVE_SSTREAM +#include +#include +#endif // HAVE_SSTREAM + +// not ideal but disable unused-function warnings +// (we get them because we have implementations in the header file, +// and this is because we want to be able to quickly separate them +// into a cpp file if necessary) +#if !defined(_MSC_VER) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-function" +#endif + +// unnamed namespace only because the implementation is in this +// header file and we don't want to export symbols to the obj files +namespace +{ + +namespace tk +{ + +// spline interpolation +class spline +{ +public: + // spline types + enum spline_type { + linear = 10, // linear interpolation + cspline = 30, // cubic splines (classical C^2) + cspline_hermite = 31 // cubic hermite splines (local, only C^1) + }; + + // boundary condition type for the spline end-points + enum bd_type { + first_deriv = 1, + second_deriv = 2, + not_a_knot = 3 + }; + +protected: + std::vector m_x,m_y; // x,y coordinates of points + // interpolation parameters + // f(x) = a_i + b_i*(x-x_i) + c_i*(x-x_i)^2 + d_i*(x-x_i)^3 + // where a_i = y_i, or else it won't go through grid points + std::vector m_b,m_c,m_d; // spline coefficients + double m_c0; // for left extrapolation + spline_type m_type; + bd_type m_left, m_right; + double m_left_value, m_right_value; + bool m_made_monotonic; + void set_coeffs_from_b(); // calculate c_i, d_i from b_i + size_t find_closest(double x) const; // closest idx so that m_x[idx]<=x + +public: + // default constructor: set boundary condition to be zero curvature + // at both ends, i.e. natural splines + spline(): m_type(cspline), + m_left(second_deriv), m_right(second_deriv), + m_left_value(0.0), m_right_value(0.0), m_made_monotonic(false) + { + ; + } + spline(const std::vector& X, const std::vector& Y, + spline_type type = cspline, + bool make_monotonic = false, + bd_type left = second_deriv, double left_value = 0.0, + bd_type right = second_deriv, double right_value = 0.0 + ): + m_type(type), + m_left(left), m_right(right), + m_left_value(left_value), m_right_value(right_value), + m_made_monotonic(false) // false correct here: make_monotonic() sets it + { + this->set_points(X,Y,m_type); + if(make_monotonic) { + this->make_monotonic(); + } + } + + + // modify boundary conditions: if called it must be before set_points() + void set_boundary(bd_type left, double left_value, + bd_type right, double right_value); + + // set all data points (cubic_spline=false means linear interpolation) + void set_points(const std::vector& x, + const std::vector& y, + spline_type type=cspline); + + // adjust coefficients so that the spline becomes piecewise monotonic + // where possible + // this is done by adjusting slopes at grid points by a non-negative + // factor and this will break C^2 + // this can also break boundary conditions if adjustments need to + // be made at the boundary points + // returns false if no adjustments have been made, true otherwise + bool make_monotonic(); + + // evaluates the spline at point x + double operator() (double x) const; + double deriv(int order, double x) const; + + // solves for all x so that: spline(x) = y + std::vector solve(double y, bool ignore_extrapolation=true) const; + + // returns the input data points + std::vector get_x() const { return m_x; } + std::vector get_y() const { return m_y; } + double get_x_min() const { assert(!m_x.empty()); return m_x.front(); } + double get_x_max() const { assert(!m_x.empty()); return m_x.back(); } + +#ifdef HAVE_SSTREAM + // spline info string, i.e. spline type, boundary conditions etc. + std::string info() const; +#endif // HAVE_SSTREAM + +}; + + + +namespace internal +{ + +// band matrix solver +class band_matrix +{ +private: + std::vector< std::vector > m_upper; // upper band + std::vector< std::vector > m_lower; // lower band +public: + band_matrix() {}; // constructor + band_matrix(int dim, int n_u, int n_l); // constructor + ~band_matrix() {}; // destructor + void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l + int dim() const; // matrix dimension + int num_upper() const + { + return (int)m_upper.size()-1; + } + int num_lower() const + { + return (int)m_lower.size()-1; + } + // access operator + double & operator () (int i, int j); // write + double operator () (int i, int j) const; // read + // we can store an additional diagonal (in m_lower) + double& saved_diag(int i); + double saved_diag(int i) const; + void lu_decompose(); + std::vector r_solve(const std::vector& b) const; + std::vector l_solve(const std::vector& b) const; + std::vector lu_solve(const std::vector& b, + bool is_lu_decomposed=false); + +}; + +double get_eps(); + +std::vector solve_cubic(double a, double b, double c, double d, + int newton_iter=0); + +} // namespace internal + + + + +// --------------------------------------------------------------------- +// implementation part, which could be separated into a cpp file +// --------------------------------------------------------------------- + +// spline implementation +// ----------------------- + +void spline::set_boundary(spline::bd_type left, double left_value, + spline::bd_type right, double right_value) +{ + assert(m_x.size()==0); // set_points() must not have happened yet + m_left=left; + m_right=right; + m_left_value=left_value; + m_right_value=right_value; +} + + +void spline::set_coeffs_from_b() +{ + assert(m_x.size()==m_y.size()); + assert(m_x.size()==m_b.size()); + assert(m_x.size()>2); + size_t n=m_b.size(); + if(m_c.size()!=n) + m_c.resize(n); + if(m_d.size()!=n) + m_d.resize(n); + + for(size_t i=0; i& x, + const std::vector& y, + spline_type type) +{ + assert(x.size()==y.size()); + assert(x.size()>=3); + // not-a-knot with 3 points has many solutions + if(m_left==not_a_knot || m_right==not_a_knot) + assert(x.size()>=4); + m_type=type; + m_made_monotonic=false; + m_x=x; + m_y=y; + int n = (int) x.size(); + // check strict monotonicity of input vector x + for(int i=0; i rhs(n); + for(int i=1; i2); + bool modified = false; + const int n=(int)m_x.size(); + // make sure: input data monotonic increasing --> b_i>=0 + // input data monotonic decreasing --> b_i<=0 + for(int i=0; i=m_y[i]) && (m_y[i]>=m_y[ip1]) && m_b[i]>0.0) ) { + modified=true; + m_b[i]=0.0; + } + } + // if input data is monotonic (b[i], b[i+1], avg have all the same sign) + // ensure a sufficient criteria for monotonicity is satisfied: + // sqrt(b[i]^2+b[i+1]^2) <= 3 |avg|, with avg=(y[i+1]-y[i])/h, + for(int i=0; i=0.0 && m_b[i+1]>=0.0 && avg>0.0) || + (m_b[i]<=0.0 && m_b[i+1]<=0.0 && avg<0.0) ) { + // input data is monotonic + double r = sqrt(m_b[i]*m_b[i]+m_b[i+1]*m_b[i+1])/std::fabs(avg); + if(r>3.0) { + // sufficient criteria for monotonicity: r<=3 + // adjust b[i] and b[i+1] + modified=true; + m_b[i] *= (3.0/r); + m_b[i+1] *= (3.0/r); + } + } + } + + if(modified==true) { + set_coeffs_from_b(); + m_made_monotonic=true; + } + + return modified; +} + +// return the closest idx so that m_x[idx] <= x (return 0 if x::const_iterator it; + it=std::upper_bound(m_x.begin(),m_x.end(),x); // *it > x + size_t idx = std::max( int(it-m_x.begin())-1, 0); // m_x[idx] <= x + return idx; +} + +double spline::operator() (double x) const +{ + // polynomial evaluation using Horner's scheme + // TODO: consider more numerically accurate algorithms, e.g.: + // - Clenshaw + // - Even-Odd method by A.C.R. Newbery + // - Compensated Horner Scheme + size_t n=m_x.size(); + size_t idx=find_closest(x); + + double h=x-m_x[idx]; + double interpol; + if(xm_x[n-1]) { + // extrapolation to the right + interpol=(m_c[n-1]*h + m_b[n-1])*h + m_y[n-1]; + } else { + // interpolation + interpol=((m_d[idx]*h + m_c[idx])*h + m_b[idx])*h + m_y[idx]; + } + return interpol; +} + +double spline::deriv(int order, double x) const +{ + assert(order>0); + size_t n=m_x.size(); + size_t idx = find_closest(x); + + double h=x-m_x[idx]; + double interpol; + if(xm_x[n-1]) { + // extrapolation to the right + switch(order) { + case 1: + interpol=2.0*m_c[n-1]*h + m_b[n-1]; + break; + case 2: + interpol=2.0*m_c[n-1]; + break; + default: + interpol=0.0; + break; + } + } else { + // interpolation + switch(order) { + case 1: + interpol=(3.0*m_d[idx]*h + 2.0*m_c[idx])*h + m_b[idx]; + break; + case 2: + interpol=6.0*m_d[idx]*h + 2.0*m_c[idx]; + break; + case 3: + interpol=6.0*m_d[idx]; + break; + default: + interpol=0.0; + break; + } + } + return interpol; +} + +std::vector spline::solve(double y, bool ignore_extrapolation) const +{ + std::vector x; // roots for the entire spline + std::vector root; // roots for each piecewise cubic + const size_t n=m_x.size(); + + // left extrapolation + if(ignore_extrapolation==false) { + root = internal::solve_cubic(m_y[0]-y,m_b[0],m_c0,0.0,1); + for(size_t j=0; j0) ? (m_x[i]-m_x[i-1]) : 0.0; + double eps = internal::get_eps()*512.0*std::min(h,1.0); + if( (-eps<=root[j]) && (root[j]0 && x.back()+eps > new_root) { + x.back()=new_root; // avoid spurious duplicate roots + } else { + x.push_back(new_root); + } + } + } + } + + // right extrapolation + if(ignore_extrapolation==false) { + root = internal::solve_cubic(m_y[n-1]-y,m_b[n-1],m_c[n-1],0.0,1); + for(size_t j=0; j0); + assert(n_u>=0); + assert(n_l>=0); + m_upper.resize(n_u+1); + m_lower.resize(n_l+1); + for(size_t i=0; i0) { + return (int)m_upper[0].size(); + } else { + return 0; + } +} + + +// defines the new operator (), so that we can access the elements +// by A(i,j), index going from i=0,...,dim()-1 +double & band_matrix::operator () (int i, int j) +{ + int k=j-i; // what band is the entry + assert( (i>=0) && (i=0) && (j diagonal, k<0 lower left part, k>0 upper right part + if(k>=0) return m_upper[k][i]; + else return m_lower[-k][i]; +} +double band_matrix::operator () (int i, int j) const +{ + int k=j-i; // what band is the entry + assert( (i>=0) && (i=0) && (j diagonal, k<0 lower left part, k>0 upper right part + if(k>=0) return m_upper[k][i]; + else return m_lower[-k][i]; +} +// second diag (used in LU decomposition), saved in m_lower +double band_matrix::saved_diag(int i) const +{ + assert( (i>=0) && (i=0) && (idim(); i++) { + assert(this->operator()(i,i)!=0.0); + this->saved_diag(i)=1.0/this->operator()(i,i); + j_min=std::max(0,i-this->num_lower()); + j_max=std::min(this->dim()-1,i+this->num_upper()); + for(int j=j_min; j<=j_max; j++) { + this->operator()(i,j) *= this->saved_diag(i); + } + this->operator()(i,i)=1.0; // prevents rounding errors + } + + // Gauss LR-Decomposition + for(int k=0; kdim(); k++) { + i_max=std::min(this->dim()-1,k+this->num_lower()); // num_lower not a mistake! + for(int i=k+1; i<=i_max; i++) { + assert(this->operator()(k,k)!=0.0); + x=-this->operator()(i,k)/this->operator()(k,k); + this->operator()(i,k)=-x; // assembly part of L + j_max=std::min(this->dim()-1,k+this->num_upper()); + for(int j=k+1; j<=j_max; j++) { + // assembly part of R + this->operator()(i,j)=this->operator()(i,j)+x*this->operator()(k,j); + } + } + } +} +// solves Ly=b +std::vector band_matrix::l_solve(const std::vector& b) const +{ + assert( this->dim()==(int)b.size() ); + std::vector x(this->dim()); + int j_start; + double sum; + for(int i=0; idim(); i++) { + sum=0; + j_start=std::max(0,i-this->num_lower()); + for(int j=j_start; joperator()(i,j)*x[j]; + x[i]=(b[i]*this->saved_diag(i)) - sum; + } + return x; +} +// solves Rx=y +std::vector band_matrix::r_solve(const std::vector& b) const +{ + assert( this->dim()==(int)b.size() ); + std::vector x(this->dim()); + int j_stop; + double sum; + for(int i=this->dim()-1; i>=0; i--) { + sum=0; + j_stop=std::min(this->dim()-1,i+this->num_upper()); + for(int j=i+1; j<=j_stop; j++) sum += this->operator()(i,j)*x[j]; + x[i]=( b[i] - sum ) / this->operator()(i,i); + } + return x; +} + +std::vector band_matrix::lu_solve(const std::vector& b, + bool is_lu_decomposed) +{ + assert( this->dim()==(int)b.size() ); + std::vector x,y; + if(is_lu_decomposed==false) { + this->lu_decompose(); + } + y=this->l_solve(b); + x=this->r_solve(y); + return x; +} + +// machine precision of a double, i.e. the successor of 1 is 1+eps +double get_eps() +{ + //return std::numeric_limits::epsilon(); // __DBL_EPSILON__ + return 2.2204460492503131e-16; // 2^-52 +} + +// solutions for a + b*x = 0 +std::vector solve_linear(double a, double b) +{ + std::vector x; // roots + if(b==0.0) { + if(a==0.0) { + // 0*x = 0 + x.resize(1); + x[0] = 0.0; // any x solves it but we need to pick one + return x; + } else { + // 0*x + ... = 0, no solution + return x; + } + } else { + x.resize(1); + x[0] = -a/b; + return x; + } +} + +// solutions for a + b*x + c*x^2 = 0 +std::vector solve_quadratic(double a, double b, double c, + int newton_iter=0) +{ + if(c==0.0) { + return solve_linear(a,b); + } + // rescale so that we solve x^2 + 2p x + q = (x+p)^2 + q - p^2 = 0 + double p=0.5*b/c; + double q=a/c; + double discr = p*p-q; + const double eps=0.5*internal::get_eps(); + double discr_err = (6.0*(p*p)+3.0*fabs(q)+fabs(discr))*eps; + + std::vector x; // roots + if(fabs(discr)<=discr_err) { + // discriminant is zero --> one root + x.resize(1); + x[0] = -p; + } else if(discr<0) { + // no root + } else { + // two roots + x.resize(2); + x[0] = -p - sqrt(discr); + x[1] = -p + sqrt(discr); + } + + // improve solution via newton steps + for(size_t i=0; i1e-8) { + x[i] -= f/f1; + } + } + } + + return x; +} + +// solutions for the cubic equation: a + b*x +c*x^2 + d*x^3 = 0 +// this is a naive implementation of the analytic solution without +// optimisation for speed or numerical accuracy +// newton_iter: number of newton iterations to improve analytical solution +// see also +// gsl: gsl_poly_solve_cubic() in solve_cubic.c +// octave: roots.m - via eigenvalues of the Frobenius companion matrix +std::vector solve_cubic(double a, double b, double c, double d, + int newton_iter) +{ + if(d==0.0) { + return solve_quadratic(a,b,c,newton_iter); + } + + // convert to normalised form: a + bx + cx^2 + x^3 = 0 + if(d!=1.0) { + a/=d; + b/=d; + c/=d; + } + + // convert to depressed cubic: z^3 - 3pz - 2q = 0 + // via substitution: z = x + c/3 + std::vector z; // roots of the depressed cubic + double p = -(1.0/3.0)*b + (1.0/9.0)*(c*c); + double r = 2.0*(c*c)-9.0*b; + double q = -0.5*a - (1.0/54.0)*(c*r); + double discr=p*p*p-q*q; // discriminant + // calculating numerical round-off errors with assumptions: + // - each operation is precise but each intermediate result x + // when stored has max error of x*eps + // - only multiplication with a power of 2 introduces no new error + // - a,b,c,d and some fractions (e.g. 1/3) have rounding errors eps + // - p_err << |p|, q_err << |q|, ... (this is violated in rare cases) + // would be more elegant to use boost::numeric::interval + const double eps = internal::get_eps(); + double p_err = eps*((3.0/3.0)*fabs(b)+(4.0/9.0)*(c*c)+fabs(p)); + double r_err = eps*(6.0*(c*c)+18.0*fabs(b)+fabs(r)); + double q_err = 0.5*fabs(a)*eps + (1.0/54.0)*fabs(c)*(r_err+fabs(r)*3.0*eps) + + fabs(q)*eps; + double discr_err = (p*p) * (3.0*p_err + fabs(p)*2.0*eps) + + fabs(q) * (2.0*q_err + fabs(q)*eps) + fabs(discr)*eps; + + // depending on the discriminant we get different solutions + if(fabs(discr)<=discr_err) { + // discriminant zero: one or two real roots + if(fabs(p)<=p_err) { + // p and q are zero: single root + z.resize(1); + z[0] = 0.0; // triple root + } else { + z.resize(2); + z[0] = 2.0*q/p; // single root + z[1] = -0.5*z[0]; // double root + } + } else if(discr>0) { + // three real roots: via trigonometric solution + z.resize(3); + double ac = (1.0/3.0) * acos( q/(p*sqrt(p)) ); + double sq = 2.0*sqrt(p); + z[0] = sq * cos(ac); + z[1] = sq * cos(ac-2.0*M_PI/3.0); + z[2] = sq * cos(ac-4.0*M_PI/3.0); + } else if (discr<0.0) { + // single real root: via Cardano's fromula + z.resize(1); + double sgnq = (q >= 0 ? 1 : -1); + double basis = fabs(q) + sqrt(-discr); + double C = sgnq * pow(basis, 1.0/3.0); // c++11 has std::cbrt() + z[0] = C + p/C; + } + for(size_t i=0; i1e-8) { + z[i] -= f/f1; + } + } + } + // ensure if a=0 we get exactly x=0 as root + // TODO: remove this fudge + if(a==0.0) { + assert(z.size()>0); // cubic should always have at least one root + double xmin=fabs(z[0]); + size_t imin=0; + for(size_t i=1; ifabs(z[i])) { + xmin=fabs(z[i]); + imin=i; + } + } + z[imin]=0.0; // replace the smallest absolute value with 0 + } + std::sort(z.begin(), z.end()); + return z; +} + + +} // namespace internal + + +} // namespace tk + + +} // namespace + +#if !defined(_MSC_VER) +#pragma GCC diagnostic pop +#endif + +#endif /* TK_SPLINE_H */