Line data    Source code

       1             : /*************************************************************************
       2             :  *
       3             :  *  Project
       4             :  *                         _____ _____  __  __ _____
       5             :  *                        / ____|  __ \|  \/  |  __ \
       6             :  *  ___  _ __   ___ _ __ | |  __| |__) | \  / | |__) |
       7             :  * / _ \| '_ \ / _ \ '_ \| | |_ |  ___/| |\/| |  ___/
       8             :  *| (_) | |_) |  __/ | | | |__| | |    | |  | | |
       9             :  * \___/| .__/ \___|_| |_|\_____|_|    |_|  |_|_|
      10             :  *      | |
      11             :  *      |_|
      12             :  *
      13             :  * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
      14             :  *
      15             :  * This software is licensed as described in the file LICENSE, which
      16             :  * you should have received as part of this distribution. The terms
      17             :  * among other details are referenced in the official documentation
      18             :  * seen here : https://akielaries.github.io/openGPMP/ along with
      19             :  * important files seen in this project.
      20             :  *
      21             :  * You may opt to use, copy, modify, merge, publish, distribute
      22             :  * and/or sell copies of the Software, and permit persons to whom
      23             :  * the Software is furnished to do so, under the terms of the
      24             :  * LICENSE file.
      25             :  *
      26             :  *
      27             :  *
      28             :  * This software is distributed on an AS IS basis, WITHOUT
      29             :  * WARRANTY OF ANY KIND, either express or implied.
      30             :  *
      31             :  ************************************************************************/
      32             : 
      33             : /*
      34             :  * Implementing various derivative related functions in C++
      35             :  */
      36             : #include <algorithm>
      37             : #include <cmath>
      38             : #include <cstdlib>
      39             : #include <iostream>
      40             : #include <openGPMP/calculus/differential.hpp>
      41             : #include <sstream>
      42             : #include <stdio.h>
      43             : #include <string>
      44             : 
      45           8 : void gpmp::Differential::add_term(double coefficient, int exponent) {
      46           8 :     terms.emplace_back(coefficient, exponent);
      47           8 : }
      48             : 
      49           0 : void gpmp::Differential::display() const {
      50           0 :     for (size_t i = 0; i < terms.size(); ++i) {
      51           0 :         const auto &term = terms[i];
      52           0 :         if (term.exponent > 1) {
      53           0 :             std::cout << term.coefficient << "*x^" << term.exponent;
      54           0 :         } else if (term.exponent == 1) {
      55           0 :             std::cout << term.coefficient << "*x";
      56             :         } else {
      57           0 :             std::cout << term.coefficient;
      58             :         }
      59             : 
      60           0 :         if (i < terms.size() - 1) {
      61           0 :             std::cout << " + ";
      62             :         }
      63             :     }
      64           0 :     std::cout << std::endl;
      65           0 : }
      66             : 
      67           1 : gpmp::Differential gpmp::Differential::power_rule() const {
      68           1 :     gpmp::Differential result;
      69           4 :     for (const auto &term : terms) {
      70           3 :         if (term.exponent > 0) {
      71           2 :             double new_coefficient = term.coefficient * term.exponent;
      72           2 :             int new_exponent = term.exponent - 1;
      73           2 :             result.add_term(new_coefficient, new_exponent);
      74             :         }
      75             :     }
      76           1 :     return result;
      77           0 : }
      78             : 
      79             : gpmp::Differential
      80           0 : gpmp::Differential::chain_rule(const gpmp::Differential &inner) const {
      81           0 :     gpmp::Differential result;
      82             : 
      83           0 :     for (const auto &outer_term : terms) {
      84             :         // Apply the chain rule to each term of the outer function
      85           0 :         gpmp::Differential inner_derivative = inner.power_rule();
      86             : 
      87             :         // Multiply each term of inner_derivative by the coefficient of the
      88             :         // outer term
      89           0 :         for (auto &inner_term : inner_derivative.terms) {
      90           0 :             inner_term.coefficient *= outer_term.coefficient;
      91             :         }
      92             : 
      93             :         // Multiply the inner derivative by the derivative of the outer term
      94           0 :         for (const auto &inner_term : inner_derivative.terms) {
      95           0 :             double new_coefficient = inner_term.coefficient;
      96           0 :             int new_exponent = outer_term.exponent + inner_term.exponent;
      97           0 :             result.add_term(new_coefficient, new_exponent);
      98             :         }
      99           0 :     }
     100             : 
     101           0 :     return result;
     102           0 : }
     103             : 
     104           0 : gpmp::Differential gpmp::Differential::nth_derivative(int n) const {
     105           0 :     gpmp::Differential result = *this;
     106           0 :     for (int i = 0; i < n; ++i) {
     107           0 :         result = result.power_rule();
     108             :     }
     109           0 :     return result;
     110           0 : }
     111             : 
     112           0 : double gpmp::Differential::eval(double x) const {
     113           0 :     double result = 0.0;
     114           0 :     for (const auto &term : terms) {
     115           0 :         result += term.coefficient * std::pow(x, term.exponent);
     116             :     }
     117           0 :     return result;
     118             : }
     119             : 
     120           0 : double gpmp::Differential::limit_at(double x) const {
     121           0 :     double result = 0.0;
     122             : 
     123           0 :     for (const auto &term : terms) {
     124           0 :         result += term.coefficient * std::pow(x, term.exponent);
     125             :     }
     126             : 
     127           0 :     return result;
     128             : }
     129             : 
     130           0 : double gpmp::Differential::limit_at_infinity() const {
     131             :     // Calculate the limit as x approaches infinity using the highest exponent
     132             :     // term
     133           0 :     if (!terms.empty()) {
     134             :         const auto &highest_term =
     135           0 :             *std::max_element(terms.begin(),
     136             :                               terms.end(),
     137           0 :                               [](const auto &a, const auto &b) {
     138           0 :                                   return a.exponent < b.exponent;
     139           0 :                               });
     140             : 
     141           0 :         if (highest_term.exponent > 0) {
     142             :             // If the highest term has a positive exponent, the limit is
     143             :             // infinity or negative infinity
     144           0 :             return (highest_term.coefficient > 0)
     145           0 :                        ? std::numeric_limits<double>::infinity()
     146           0 :                        : -std::numeric_limits<double>::infinity();
     147             :         } else {
     148             :             // If the highest term has a zero exponent, the limit is the
     149             :             // constant term
     150           0 :             return highest_term.coefficient;
     151             :         }
     152             :     }
     153             : 
     154             :     // If the differential is empty, the limit is undefined
     155           0 :     return std::numeric_limits<double>::quiet_NaN();
     156             : }