From 12fdab54be153d1e63191396b9ab62b231eb3270 Mon Sep 17 00:00:00 2001 From: Oleg Bovykin Date: Fri, 13 Jun 2014 20:54:11 +0400 Subject: [PATCH 1/2] Total library rewrite --- .gitignore | 27 +- .rspec | 4 +- Gemfile | 4 + History.txt | 3 - README.txt => LICENSE.txt | 47 +-- Manifest.txt | 9 - README.md | 55 +++ Rakefile | 27 +- integration.gemspec | 27 ++ lib/integration.rb | 388 ++++-------------- lib/integration/gsl.rb | 49 +++ lib/integration/integrator.rb | 22 + .../integrator/adaptive_quadrature.rb | 26 ++ lib/integration/integrator/gauss.rb | 39 ++ lib/integration/integrator/monte_carlo.rb | 39 ++ lib/integration/integrator/rectangle.rb | 11 + lib/integration/integrator/romberg.rb | 33 ++ lib/integration/integrator/simpson.rb | 20 + lib/integration/integrator/trapezoid.rb | 15 + lib/integration/ruby.rb | 72 ++++ lib/integration/version.rb | 3 + spec/integration_spec.rb | 145 ++++--- spec/spec_helper.rb | 14 +- 23 files changed, 623 insertions(+), 456 deletions(-) create mode 100644 Gemfile delete mode 100644 History.txt rename README.txt => LICENSE.txt (53%) delete mode 100644 Manifest.txt create mode 100644 README.md create mode 100644 integration.gemspec create mode 100644 lib/integration/gsl.rb create mode 100644 lib/integration/integrator.rb create mode 100644 lib/integration/integrator/adaptive_quadrature.rb create mode 100644 lib/integration/integrator/gauss.rb create mode 100644 lib/integration/integrator/monte_carlo.rb create mode 100644 lib/integration/integrator/rectangle.rb create mode 100644 lib/integration/integrator/romberg.rb create mode 100644 lib/integration/integrator/simpson.rb create mode 100644 lib/integration/integrator/trapezoid.rb create mode 100644 lib/integration/ruby.rb create mode 100644 lib/integration/version.rb diff --git a/.gitignore b/.gitignore index 2500617..7178b3e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,28 @@ +*.gem +*.rbc +.bundle +.config +.yardoc +Gemfile.lock +InstalledFiles +_yardoc +coverage +doc/ +lib/bundler/man +pkg +rdoc +spec/reports +test/tmp +test/version_tmp +tmp +*.bundle +*.so +*.o +*.a +mkmf.log *~ coverage/* -pkg/* +*DS_Store +.rvmrc +.ruby-version +.ruby-gemset diff --git a/.rspec b/.rspec index d8ed776..0d786ba 100644 --- a/.rspec +++ b/.rspec @@ -1,3 +1,3 @@ --color --f s --b +--warnings +--require spec_helper diff --git a/Gemfile b/Gemfile new file mode 100644 index 0000000..0704604 --- /dev/null +++ b/Gemfile @@ -0,0 +1,4 @@ +source 'https://rubygems.org' + +# Specify your gem's dependencies in integration.gemspec +gemspec diff --git a/History.txt b/History.txt deleted file mode 100644 index 51ddfc2..0000000 --- a/History.txt +++ /dev/null @@ -1,3 +0,0 @@ -=== 0.1.0 / 2011-08-23 - -* Get code from http://rubyforge.org/projects/integral/ and converted to a gem. diff --git a/README.txt b/LICENSE.txt similarity index 53% rename from README.txt rename to LICENSE.txt index 884c4a1..32379bc 100644 --- a/README.txt +++ b/LICENSE.txt @@ -1,51 +1,17 @@ -= integration +Copyright (c) 2005 Beng + 2011 clbustos + 2014 arrowcircle -* http://github.com/clbustos/integration - -== DESCRIPTION: - -Numerical integration for Ruby, with a simple interface - -== FEATURES/PROBLEMS: - -* Use only one method: Integration.integrate -* Rectangular, Trapezoidal, Simpson, Adaptive quadrature, Monte Carlo and Romberg integration methods available on pure Ruby. -* If available, uses Ruby/GSL QNG non-adaptive Gauss-Kronrod integration and QAG adaptive integration, with support for Infinity+ and Infinity- - - -== SYNOPSIS: - - Integration.integrate(1,2,{:tolerance=>1e-10,:method=>:simpson}) {|x| x**2} => 2.333333 - - # Support for infinity bounds with GSL QAG adaptative integration - - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - Integration.integrate(Integration::MInfinity, 0 , {:tolerance=>1e-10}, &normal_pdf) => 0.5 - Integration.integrate(0, Integration::Infinity , {:tolerance=>1e-10}, &normal_pdf) => 0.5 - -== REQUIREMENTS: - -* Ruby/GSL, for better performance and support for infinite bounds - -== INSTALL: - - gem install integration - -== LICENSE: - - Copyright (c) 2005 Beng - 2011 clbustos - 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 @@ -53,8 +19,7 @@ Numerical integration for Ruby, with a simple interface 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. - + Except as contained in this notice, the name of the Beng shall not be used in advertising or otherwise to promote the sale, use or other dealings in this Software without prior written authorization from Beng. - diff --git a/Manifest.txt b/Manifest.txt deleted file mode 100644 index b23b9ea..0000000 --- a/Manifest.txt +++ /dev/null @@ -1,9 +0,0 @@ -.autotest -History.txt -Manifest.txt -README.txt -Rakefile -lib/integration.rb -spec/integration_spec.rb -spec/spec.opts -spec/spec_helper.rb diff --git a/README.md b/README.md new file mode 100644 index 0000000..d2daed8 --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ +# Integration + +Numerical integration for Ruby, with a simple interface + +## Installation + +Add gem `integration` to your Gemfile: + + gem 'integration' + +## Usage: + +Integrate have only one method: `integrate` + +Without GSL: + +``` +Integration.integrate(1, 2, {tolerance: 1e-10, method: :simpson}) { |x| x**2 } +=> 2.333333 +``` +With GSL (support for infinity bounds with GSL QAG adaptative integration): + +``` +func = ->(x) { (1/Math.sqrt(2 * Math::PI)) * Math.exp(-(x**2 / 2)) } + +Integration.integrate(Integration::MInfinity, 0, {tolerance: 1e-10 }, &func) +=> 0.5 + +Integration.integrate(0, Integration::Infinity, {tolerance: 1e-10}, &func) +=> 0.5 +``` + +## Available methods +Pure Ruby methods: + +* Simpson (`:simpson`, default method) +* Rectangular (`:rectangular`) +* Trapezoidal (`:trapezoidal`) +* Adaptive quadrature (`:adaptive_quadrature`) +* Romberg (`:romberg`) +* Monte Carlo (`:monte_carlo`, bad results) + +GSL methods: + +* QNG (`:qng`) +* QAG (`:qag`) + +## REQUIREMENTS: + +Integration works only with Ruby 1.9+ + +Integration depends on GSL ( GNU Scientific Library ) for infinity bounds and faster algoritms support. + +* For Mac OS X: `brew install gsl` +* For Ubuntu / Debian: `sudo apt-get install gsl` \ No newline at end of file diff --git a/Rakefile b/Rakefile index 1535c9f..809eb56 100644 --- a/Rakefile +++ b/Rakefile @@ -1,27 +1,2 @@ -# -*- ruby -*- -$:.unshift(File.expand_path(File.dirname(__FILE__)+"/lib/")) -require 'rubygems' -require 'hoe' -require 'integration' -require 'rubyforge' - -# Hoe.plugin :compiler -# Hoe.plugin :gem_prelude_sucks - Hoe.plugin :git -# Hoe.plugin :inline -# Hoe.plugin :racc -# Hoe.plugin :rubyforge - -Hoe.spec 'integration' do - - self.developer('Ben Gimpert', 'NO_EMAIL') - - self.developer('Claudio Bustos', 'clbustos_at_gmail.com') - - self.version=Integration::VERSION - self.extra_dev_deps << ["rspec",">=2.0"] << ["rubyforge",">=0"] - -end -# git log --pretty=format:"*%s[%cn]" v0.5.0..HEAD >> History.txt -# vim: syntax=ruby +require "bundler/gem_tasks" diff --git a/integration.gemspec b/integration.gemspec new file mode 100644 index 0000000..b802bee --- /dev/null +++ b/integration.gemspec @@ -0,0 +1,27 @@ +# coding: utf-8 +lib = File.expand_path('../lib', __FILE__) +$LOAD_PATH.unshift(lib) unless $LOAD_PATH.include?(lib) +require 'integration' + +Gem::Specification.new do |spec| + spec.name = 'integration' + spec.version = Integration::Version::VERSION + spec.authors = ['Ben Gimpert', 'Claudio Bustos', 'Oleg Bovykin'] + spec.email = ['clbustos_at_gmail.com', 'oleg.bovykin@gmail.com'] + spec.summary = %q{Integration methods, based on original work by Beng} + spec.description = %q{Numerical integration for Ruby, with a simple interface} + spec.homepage = '' + spec.license = 'MIT' + + spec.files = `git ls-files -z`.split("\x0") + spec.executables = spec.files.grep(%r{^bin/}) { |f| File.basename(f) } + spec.test_files = spec.files.grep(%r{^(test|spec|features)/}) + spec.require_paths = ['lib'] + + spec.add_development_dependency 'bundler', '~> 1.6' + spec.add_development_dependency 'rake' + spec.add_development_dependency 'rspec' + + spec.add_dependency 'rb-gsl' + spec.add_dependency 'activesupport' +end diff --git a/lib/integration.rb b/lib/integration.rb index 94d1ccc..1675478 100644 --- a/lib/integration.rb +++ b/lib/integration.rb @@ -1,324 +1,88 @@ -# Copyright (c) 2005 Beng (original code) -# 2011 clbustos -# -# 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 X CONSORTIUM 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. -# -# Except as contained in this notice, the name of the Beng shall not -# be used in advertising or otherwise to promote the sale, use or other dealings -# in this Software without prior written authorization from Beng. - -# Diverse integration methods -# Use Integration.integrate as wrapper to direct access to methods -# -# Method API -# +require 'active_support/core_ext/string' class Integration - VERSION = '0.1.0' - # Minus Infinity - MInfinity=:minfinity - # Infinity - Infinity=:infinity + require 'integration/version' + require 'integration/gsl' + require 'integration/ruby' + + M_INFINITY = :minfinity + INFINITY = :infinity + + RUBY_METHOD = [:rectangle, :trapezoid, :simpson, + :adaptive_quadrature, :gauss, :romberg, :monte_carlo] + GSL_METHOD = [:qng, :qag] + + DEFAULT_OPTIONS = { tolerance: 1e-10, initial_step: 16, step: 16 } + + attr_reader :lower_bound, :upper_bound, :options, :func + + def self.infinite?(bound) + bound == INFINITY || bound == M_INFINITY + end + + def initialize(lower, upper, options = {}, &func) + @lower_bound = lower + @upper_bound = upper + @options = default_options.merge(options) + @func = func + fails + end + + def fails + fail 'Non-numeric bounds' if invalid_bounds? + fail 'Ruby methods doesnt support infinity bounds' if inifinite_bounds? && ruby_method? + fail "Unknown integration method \"#{options[:method]}\"" unless (ruby_method? || gsl_method?) + end + + def result + return Integration::Ruby.new(lower_bound, upper_bound, options, &func).result if ruby_method? + Integration::Gsl.new(lower_bound, upper_bound, options, &func).result + end + + def default_method + self.class.has_gsl? ? :qag : :simpson + end + + def ruby_method? + RUBY_METHOD.include? options[:method] + end + + def gsl_method? + GSL_METHOD.include? options[:method] + end + + def default_options + DEFAULT_OPTIONS.merge(method: default_method) + end + + def inifinite_bounds? + (self.class.infinite?(lower_bound) || self.class.infinite?(upper_bound)) + end + + def invalid_bounds? + return false if inifinite_bounds? + !((lower_bound.is_a? Numeric) && (upper_bound.is_a? Numeric)) + end + class << self - - # Create a method 'has_' on Module - # which require a library and return true or false - # according to success of failure - def create_has_library(library) #:nodoc: + def create_has_library(library) define_singleton_method("has_#{library}?") do - cv="@@#{library}" - if !class_variable_defined? cv - begin - require library.to_s - class_variable_set(cv, true) - rescue LoadError - class_variable_set(cv, false) - end - end + cv = "@@#{library}" + assign_cv(library) unless class_variable_defined?(cv) class_variable_get(cv) end end - # Rectangle method - # +n+ implies number of subdivisions - # Source: - # * Ayres : Outline of calculus - def rectangle(t1, t2, n, &f) - d=(t2-t1) / n.to_f - n.times.inject(0) {|ac,i| - ac+f[t1+d*(i+0.5)] - }*d - end - alias_method :midpoint, :rectangle - # Trapezoid method - # +n+ implies number of subdivisions - # Source: - # * Ayres : Outline of calculus - def trapezoid(t1, t2, n, &f) - d=(t2-t1) / n.to_f - (d/2.0)*(f[t1]+ - 2*(1..(n-1)).inject(0){|ac,i| - ac+f[t1+d*i] - }+f[t2]) - end - # Simpson's rule - # +n+ implies number of subdivisions - # Source: - # * Ayres : Outline of calculus - def simpson(t1, t2, n, &f) - n += 1 unless n % 2 == 0 - d=(t2-t1) / n.to_f - out= (d / 3.0)*(f[t1.to_f].to_f+ - ((1..(n-1)).inject(0) {|ac,i| - ac+((i%2==0) ? 2 : 4)*f[t1+d*i] - })+f[t2.to_f].to_f) - out - end - - def adaptive_quadrature(a, b, tolerance) - h = (b.to_f - a) / 2 - fa = yield(a) - fc = yield(a + h) - fb = yield(b) - s = h * (fa + (4 * fc) + fb) / 3 - helper = Proc.new { |a, b, fa, fb, fc, h, s, level| - if level < 1/tolerance.to_f - fd = yield(a + (h / 2)) - fe = yield(a + (3 * (h / 2))) - s1 = h * (fa + (4.0 * fd) + fc) / 6 - s2 = h * (fc + (4.0 * fe) + fb) / 6 - if ((s1 + s2) - s).abs <= tolerance - s1 + s2 - else - helper.call(a, a + h, fa, fc, fd, h / 2, s1, level + 1) + - helper.call(a + h, b, fc, fb, fe, h / 2, s2, level + 1) - end - else - raise "Integral did not converge" - end - } - return helper.call(a, b, fa, fb, fc, h, s, 1) - end - - def gauss(t1, t2, n) - case n - when 1 - z = [0.0] - w = [2.0] - when 2 - z = [-0.57735026919, 0.57735026919] - w = [1.0, 1.0] - when 3 - z = [-0.774596669241, 0.0, 0.774596669241] - w = [0.555555555556, 0.888888888889, 0.555555555556] - when 4 - z = [-0.861136311594, -0.339981043585, 0.339981043585, 0.861136311594] - w = [0.347854845137, 0.652145154863, 0.652145154863, 0.347854845137] - when 5 - z = [-0.906179845939, -0.538469310106, 0.0, 0.538469310106, 0.906179845939] - w = [0.236926885056, 0.478628670499, 0.568888888889, 0.478628670499, 0.236926885056] - when 6 - z = [-0.932469514203, -0.661209386466, -0.238619186083, 0.238619186083, 0.661209386466, 0.932469514203] - w = [0.171324492379, 0.360761573048, 0.467913934573, 0.467913934573, 0.360761573048, 0.171324492379] - when 7 - z = [-0.949107912343, -0.741531185599, -0.405845151377, 0.0, 0.405845151377, 0.741531185599, 0.949107912343] - w = [0.129484966169, 0.279705391489, 0.381830050505, 0.417959183673, 0.381830050505, 0.279705391489, 0.129484966169] - when 8 - z = [-0.960289856498, -0.796666477414, -0.525532409916, -0.183434642496, 0.183434642496, 0.525532409916, 0.796666477414, 0.960289856498] - w = [0.10122853629, 0.222381034453, 0.313706645878, 0.362683783378, 0.362683783378, 0.313706645878, 0.222381034453, 0.10122853629] - when 9 - z = [-0.968160239508, -0.836031107327, -0.613371432701, -0.324253423404, 0.0, 0.324253423404, 0.613371432701, 0.836031107327, 0.968160239508] - w = [0.0812743883616, 0.180648160695, 0.260610696403, 0.31234707704, 0.330239355001, 0.31234707704, 0.260610696403, 0.180648160695, 0.0812743883616] - when 10 - z = [-0.973906528517, -0.865063366689, -0.679409568299, -0.433395394129, -0.148874338982, 0.148874338982, 0.433395394129, 0.679409568299, 0.865063366689, 0.973906528517] - w = [0.0666713443087, 0.149451349151, 0.219086362516, 0.26926671931, 0.295524224715, 0.295524224715, 0.26926671931, 0.219086362516, 0.149451349151, 0.0666713443087] - else - raise "Invalid number of spaced abscissas #{n}, should be 1-10" - end - sum = 0 - (0...n).each do |i| - t = ((t1.to_f + t2) / 2) + (((t2 - t1) / 2) * z[i]) - sum += w[i] * yield(t) - end - return ((t2 - t1) / 2.0) * sum - end - - def romberg(a, b, tolerance) - # NOTE one-based arrays are used for convenience - - h = b.to_f - a - m = 1 - close = 1 - r = [[], [], [], [], [], [], [], [], [], [], [], [], []]; - r[1][1] = (h / 2) * (yield(a) + yield(b)) - j = 1 - while j <= 11 && tolerance < close - j += 1 - r[j][0] = 0 - h /= 2 - sum = 0 - (1..m).each do |k| - sum += yield(a + (h * ((2 * k) - 1))) - end - m *= 2 - r[j][1] = r[j-1][1] / 2 + (h * sum) - (1..j-1).each do |k| - r[j][k+1] = r[j][k] + ((r[j][k] - r[j-1][k]) / ((4 ** k) - 1)) - end - close = (r[j][j] - r[j-1][j-1]) - end - return r[j][j] - end - - def monte_carlo(t1, t2, n) - width = (t2 - t1).to_f - height = nil - vals = [] - n.times do - t = t1 + (rand() * width) - ft = yield(t) - height = ft if height.nil? || ft > height - vals << ft - end - area_ratio = 0 - vals.each do |ft| - area_ratio += (ft / height.to_f) / n.to_f - end - return (width * height) * area_ratio - end - def is_infinite?(v) - v==Infinity or v==MInfinity - end - # Methods available on pure ruby - RUBY_METHOD=[:rectangle,:trapezoid,:simpson, :adaptive_quadrature , :gauss, :romberg, :monte_carlo] - # Methods available with Ruby/GSL library - GSL_METHOD=[:qng, :qag] - # Get the integral for a function +f+, with bounds +t1+ and - # +t2+ given a hash of +options+. - # If Ruby/GSL is available, you could use +Integration::Minfinity+ - # and +Integration::Infinity+ as bounds. Method - # Options are - # [:tolerance] Maximum difference between real and calculated integral. - # Default: 1e-10 - # [:initial_step] Initial number of subdivitions - # [:step] Subdivitions increment on each iteration - # [:method] Integration method. - # Methods are - # [:rectangle] for [:initial_step+:step*iteration] quadrilateral subdivisions - # [:trapezoid] for [:initial_step+:step*iteration] trapezoid-al subdivisions - # [:simpson] for [:initial_step+:step*iteration] parabolic subdivisions - # [:adaptive_quadrature] for recursive appoximations until error [tolerance] - # [:gauss] [:initial_step+:step*iteration] weighted subdivisons using translated -1 -> +1 endpoints - # [:romberg] extrapolation of recursion approximation until error < [tolerance] - # [:monte_carlo] make [:initial_step+:step*iteration] random samples, and check for above/below curve - # [:qng] GSL QNG non-adaptive Gauss-Kronrod integration - # [:qag] GSL QAG adaptive integration, with support for infinite bounds - def integrate(t1,t2,options=Hash.new, &f) - inf_bounds=(is_infinite?(t1) or is_infinite?(t2)) - raise "No function passed" unless block_given? - raise "Non-numeric bounds" unless ((t1.is_a? Numeric) and (t2.is_a? Numeric)) or inf_bounds - if(inf_bounds) - lower_bound=t1 - upper_bound=t2 - options[:method]=:qag if options[:method].nil? - else - lower_bound = [t1, t2].min - upper_bound = [t1, t2].max - end - def_method=(has_gsl?) ? :qag : :simpson - default_opts={:tolerance=>1e-10, :initial_step=>16, :step=>16, :method=>def_method} - options=default_opts.merge(options) - if RUBY_METHOD.include? options[:method] - raise "Ruby methods doesn't support infinity bounds" if inf_bounds - integrate_ruby(lower_bound,upper_bound,options,&f) - elsif GSL_METHOD.include? options[:method] - integrate_gsl(lower_bound,upper_bound,options,&f) - else - raise "Unknown integration method \"#{options[:method]}\"" - end - end - def integrate_gsl(lower_bound,upper_bound,options,&f) - - f = GSL::Function.alloc(&f) - method=options[:method] - tolerance=options[:tolerance] - - if(method==:qag) - w = GSL::Integration::Workspace.alloc() - if(is_infinite?(lower_bound) and is_infinite?(upper_bound)) - #puts "ambos" - val=f.qagi([tolerance,0.0], 1000, w) - elsif is_infinite?(lower_bound) - #puts "inferior #{upper_bound}" - val=f.qagil(upper_bound, [tolerance, 0], w) - elsif is_infinite?(upper_bound) - #puts "superior" - val=f.qagiu(lower_bound, [tolerance, 0], w) - else - - val=f.qag([lower_bound,upper_bound],[tolerance,0.0], GSL::Integration::GAUSS61, w) - end - elsif(method==:qng) - val=f.qng([lower_bound, upper_bound], [tolerance, 0.0]) - else - raise "Unknown integration method \"#{method}\"" - end - val[0] + + def assign_cv(library) + cv = "@@#{library}" + require library.to_s + class_variable_set(cv, true) + rescue LoadError + class_variable_set(cv, false) end - def integrate_ruby(lower_bound,upper_bound,options,&f) - method=options[:method] - tolerance=options[:tolerance] - initial_step=options[:initial_step] - step=options[:step] - - begin - method_obj = Integration.method(method.to_s.downcase) - rescue - raise "Unknown integration method \"#{method}\"" - end - current_step=initial_step - if(method==:adaptive_quadrature or method==:romberg or method==:gauss) - if(method==:gauss) - initial_step=10 if initial_step>10 - tolerance=initial_step - end - method_obj.call(lower_bound, upper_bound, tolerance, &f) - else - #puts "iniciando" - value=method_obj.call(lower_bound, upper_bound, current_step, &f) - previous=value+(tolerance*2) - diffs=[] - while((previous-value).abs > tolerance) do - #puts("Valor:#{value}, paso:#{current_step}") - #puts(current_step) - diffs.push((previous-value).abs) - #diffs.push value - current_step+=step - previous=value - #puts "Llamando al metodo" - - value=method_obj.call(lower_bound, upper_bound, current_step, &f) - end - #p diffs - - value - end + def integrate(t1, t2, options = {}, &f) + new(t1, t2, options, &f).result end end create_has_library :gsl diff --git a/lib/integration/gsl.rb b/lib/integration/gsl.rb new file mode 100644 index 0000000..3555d0f --- /dev/null +++ b/lib/integration/gsl.rb @@ -0,0 +1,49 @@ +class Integration::Gsl + M_INFINITY = :minfinity + INFINITY = :infinity + + def infinite?(bound) + bound == INFINITY || bound == M_INFINITY + end + + attr_reader :lower_bound, :upper_bound, :options, :func + def initialize(lower_bound, upper_bound, options, &func) + @lower_bound = lower_bound + @upper_bound = upper_bound + @options = options + @func = GSL::Function.alloc(&func) + end + + def result + return qag if qag? + qng + end + + def qng + func.qng([lower_bound, upper_bound], [tolerance, 0.0])[0] + end + + def qag + w = GSL::Integration::Workspace.alloc + return func.qagi([tolerance, 0.0], 1000, w)[0] if infinite?(lower_bound) && infinite?(upper_bound) + return func.qagil(upper_bound, [tolerance, 0], w)[0] if infinite?(lower_bound) + return func.qagiu(lower_bound, [tolerance, 0], w)[0] if infinite?(upper_bound) + func.qag([lower_bound, upper_bound], [tolerance, 0.0], GSL::Integration::GAUSS61, w)[0] + end + + def qag? + method == :qag + end + + def qng? + method == :qng + end + + def method + options[:method] + end + + def tolerance + options[:tolerance] + end +end diff --git a/lib/integration/integrator.rb b/lib/integration/integrator.rb new file mode 100644 index 0000000..7297b39 --- /dev/null +++ b/lib/integration/integrator.rb @@ -0,0 +1,22 @@ +class Integration::Integrator + require 'integration/integrator/simpson' + require 'integration/integrator/rectangle' + require 'integration/integrator/trapezoid' + require 'integration/integrator/adaptive_quadrature' + require 'integration/integrator/gauss' + require 'integration/integrator/romberg' + require 'integration/integrator/monte_carlo' + + attr_reader :lower_bound, :upper_bound, :n, :func + + def initialize(lower, upper, n, &func) + @lower_bound = lower + @upper_bound = upper + @n = n + @func = func + end + + def step + @step ||= (upper_bound - lower_bound) / n.to_f + end +end diff --git a/lib/integration/integrator/adaptive_quadrature.rb b/lib/integration/integrator/adaptive_quadrature.rb new file mode 100644 index 0000000..da9bad9 --- /dev/null +++ b/lib/integration/integrator/adaptive_quadrature.rb @@ -0,0 +1,26 @@ +class Integration::AdaptiveQuadrature < Integration::Integrator + def result + h = (upper_bound.to_f - lower_bound) / 2.0 + fa = func[lower_bound] + fc = func[lower_bound + h] + fb = func[upper_bound] + s = h * (fa + (4 * fc) + fb) / 3 + helper = Proc.new do |a, b, fa, fb, fc, h, s, level| + if level < 1 / n.to_f + fd = func[a + (h / 2)] + fe = func[a + (3 * (h / 2))] + s1 = h * (fa + (4.0 * fd) + fc) / 6 + s2 = h * (fc + (4.0 * fe) + fb) / 6 + if ((s1 + s2) - s).abs <= n + s1 + s2 + else + helper.call(a, a + h, fa, fc, fd, h / 2, s1, level + 1) + + helper.call(a + h, b, fc, fb, fe, h / 2, s2, level + 1) + end + else + fail 'Integral did not converge' + end + end + return helper.call(lower_bound, upper_bound, fa, fb, fc, h, s, 1) + end +end diff --git a/lib/integration/integrator/gauss.rb b/lib/integration/integrator/gauss.rb new file mode 100644 index 0000000..cc91eaa --- /dev/null +++ b/lib/integration/integrator/gauss.rb @@ -0,0 +1,39 @@ +class Integration::Gauss < Integration::Integrator + def result + fail "Invalid number of spaced abscissas #{n}, should be 1-10" if n > 10 || n < 1 + (0..n).to_a.reduce(0.0) do |memo, i| + t = ((lower_bound.to_f + upper_bound.to_f) / 2.0) + (step * z[i - 1]) + memo + w[i - 1] * func[t] + end * step + end + + def z + case n + when 1 then [0.0] + when 2 then [-0.57735026919, 0.57735026919] + when 3 then [-0.774596669241, 0.0, 0.774596669241] + when 4 then [-0.861136311594, -0.339981043585, 0.339981043585, 0.861136311594] + when 5 then [-0.906179845939, -0.538469310106, 0.0, 0.538469310106, 0.906179845939] + when 6 then [-0.932469514203, -0.661209386466, -0.238619186083, 0.238619186083, 0.661209386466, 0.932469514203] + when 7 then [-0.949107912343, -0.741531185599, -0.405845151377, 0.0, 0.405845151377, 0.741531185599, 0.949107912343] + when 8 then [-0.960289856498, -0.796666477414, -0.525532409916, -0.183434642496, 0.183434642496, 0.525532409916, 0.796666477414, 0.960289856498] + when 9 then [-0.968160239508, -0.836031107327, -0.613371432701, -0.324253423404, 0.0, 0.324253423404, 0.613371432701, 0.836031107327, 0.968160239508] + when 10 then [-0.973906528517, -0.865063366689, -0.679409568299, -0.433395394129, -0.148874338982, 0.148874338982, 0.433395394129, 0.679409568299, 0.865063366689, 0.973906528517] + end + end + + def w + case n + when 1 then [2.0] + when 2 then [1.0, 1.0] + when 3 then [0.555555555556, 0.888888888889, 0.555555555556] + when 4 then [0.347854845137, 0.652145154863, 0.652145154863, 0.347854845137] + when 5 then [0.236926885056, 0.478628670499, 0.568888888889, 0.478628670499, 0.236926885056] + when 6 then [0.171324492379, 0.360761573048, 0.467913934573, 0.467913934573, 0.360761573048, 0.171324492379] + when 7 then [0.129484966169, 0.279705391489, 0.381830050505, 0.417959183673, 0.381830050505, 0.279705391489, 0.129484966169] + when 8 then [0.10122853629, 0.222381034453, 0.313706645878, 0.362683783378, 0.362683783378, 0.313706645878, 0.222381034453, 0.10122853629] + when 9 then [0.0812743883616, 0.180648160695, 0.260610696403, 0.31234707704, 0.330239355001, 0.31234707704, 0.260610696403, 0.180648160695, 0.0812743883616] + when 10 then [0.0666713443087, 0.149451349151, 0.219086362516, 0.26926671931, 0.295524224715, 0.295524224715, 0.26926671931, 0.219086362516, 0.149451349151, 0.0666713443087] + end + end +end diff --git a/lib/integration/integrator/monte_carlo.rb b/lib/integration/integrator/monte_carlo.rb new file mode 100644 index 0000000..4bc67ae --- /dev/null +++ b/lib/integration/integrator/monte_carlo.rb @@ -0,0 +1,39 @@ +class Integration::MonteCarlo < Integration::Integrator + attr_accessor :height, :min, :max + + def result + area_ratio * width * height + end + + def get_vals + (0..n).to_a.reduce([]) do |memo, _i| + t = lower_bound + (rand * width) + ft = func[t] + @min ||= ft + @min = ft if ft < @min + @max ||= ft + @max = ft if ft > @max + @height ||= ft + @height = ft if ft > @height + memo + [ft] + end + end + + def vals + @vals ||= get_vals + end + + def area_ratio + vals.reduce(0) do |memo, ft| + memo + ft / (height.to_f * n.to_f) + end + end + + def width + (upper_bound - lower_bound).to_f + end + + def height + max - min + end +end diff --git a/lib/integration/integrator/rectangle.rb b/lib/integration/integrator/rectangle.rb new file mode 100644 index 0000000..e07bedf --- /dev/null +++ b/lib/integration/integrator/rectangle.rb @@ -0,0 +1,11 @@ +# Rectangle method +# +n+ implies number of subdivisions +# Source: +# * Ayres : Outline of calculus +class Integration::Rectangle < Integration::Integrator + def result + n.times.reduce(0) do |ac, i| + ac + func[lower_bound + step * (i + 0.5)] + end * step + end +end diff --git a/lib/integration/integrator/romberg.rb b/lib/integration/integrator/romberg.rb new file mode 100644 index 0000000..475a89b --- /dev/null +++ b/lib/integration/integrator/romberg.rb @@ -0,0 +1,33 @@ +class Integration::Romberg < Integration::Integrator + attr_accessor :h, :m, :close, :r, :j + + def initialize(lower_bound, upper_bound, n, &func) + super + @h = upper_bound.to_f - lower_bound + @close = 1 + @m = 1 + @r = [[], [], [], [], [], [], [], [], [], [], [], [], []] + @r[1][1] = (@h / 2) * (func[lower_bound] + func[upper_bound]) + @j = 1 + end + + def result + iterate while j <= 11 && n < close + r[j][j] + end + + def iterate + @j += 1 + r[j][0] = 0 + @h /= 2 + sum = (1..m).to_a.reduce(0) do |memo, k| + memo + func[lower_bound + (h * ((2 * k) - 1))] + end + @m *= 2 + r[j][1] = r[j - 1][1] / 2 + (h * sum) + (1..j - 1).each do |k| + @r[j][k + 1] = r[j][k] + ((r[j][k] - r[j - 1][k]) / ((4**k) - 1)) + end + @close = (r[j][j] - r[j - 1][j - 1]) + end +end diff --git a/lib/integration/integrator/simpson.rb b/lib/integration/integrator/simpson.rb new file mode 100644 index 0000000..309122f --- /dev/null +++ b/lib/integration/integrator/simpson.rb @@ -0,0 +1,20 @@ +# Simpson's rule +# +n+ implies number of subdivisions +# Source: +# * Ayres : Outline of calculus +class Integration::Simpson < Integration::Integrator + def initialize(lower_bound, upper_bound, n, &func) + super + @n += 1 unless @n.even? + end + + def result + (step / 3.0) * (func[lower_bound.to_f].to_f + func[upper_bound.to_f].to_f + iteration) + end + + def iteration + (1..(n - 1)).reduce(0) do |ac, i| + ac + (i.even? ? 2 : 4) * func[lower_bound + step * i] + end + end +end diff --git a/lib/integration/integrator/trapezoid.rb b/lib/integration/integrator/trapezoid.rb new file mode 100644 index 0000000..e1a7ac7 --- /dev/null +++ b/lib/integration/integrator/trapezoid.rb @@ -0,0 +1,15 @@ +# Trapezoid method +# +n+ implies number of subdivisions +# Source: +# * Ayres : Outline of calculus +class Integration::Trapezoid < Integration::Integrator + def result + (step / 2.0) * (func[lower_bound] + func[upper_bound] + iteration) + end + + def iteration + 2 * (1..(n - 1)).reduce(0) do |ac, i| + ac + func[lower_bound + step * i] + end + end +end diff --git a/lib/integration/ruby.rb b/lib/integration/ruby.rb new file mode 100644 index 0000000..d7b082a --- /dev/null +++ b/lib/integration/ruby.rb @@ -0,0 +1,72 @@ +require 'integration/integrator' +class Integration::Ruby + M_INFINITY = :minfinity + INFINITY = :infinity + + def infinite?(bound) + bound == INFINITY || bound == M_INFINITY + end + + attr_reader :lower_bound, :upper_bound, :options, :func, :current_step + def initialize(lower_bound, upper_bound, options, &func) + @lower_bound = lower_bound + @upper_bound = upper_bound + @options = options + @func = func + @current_step = initial_step + end + + def result + return non_iterative if non_iterative? + iterative + end + + def non_iterative + method_obj.new(lower_bound, upper_bound, tolerance, &func).result + end + + def iterative + value = method_obj.new(lower_bound, upper_bound, current_step, &func).result + previous = value + (tolerance * 2) + diffs = [] + while (previous - value).abs > tolerance + diffs.push((previous - value).abs) + @current_step += step + previous = value + value = method_obj.new(lower_bound, upper_bound, current_step, &func).result + end + value + end + + def non_iterative? + [:adaptive_quadrature, :romberg, :gauss].include? method + end + + def iterative? + !non_iterative? + end + + def method_obj + "Integration::#{method.to_s.camelize}".constantize + rescue + raise "Unknown integration method \"#{method}\"" + end + + def step + options[:step] + end + + def initial_step + return 10 if options[:initial_step] > 10 && method == :gauss + options[:initial_step] + end + + def method + options[:method] + end + + def tolerance + return initial_step if method == :gauss + options[:tolerance] + end +end diff --git a/lib/integration/version.rb b/lib/integration/version.rb new file mode 100644 index 0000000..78f7750 --- /dev/null +++ b/lib/integration/version.rb @@ -0,0 +1,3 @@ +module Integration::Version + VERSION = '0.1.1' +end diff --git a/spec/integration_spec.rb b/spec/integration_spec.rb index 824120d..3fb08ef 100644 --- a/spec/integration_spec.rb +++ b/spec/integration_spec.rb @@ -1,69 +1,112 @@ -require File.expand_path(File.dirname(__FILE__)+"/spec_helper.rb") +require 'spec_helper' + describe Integration do - a=lambda {|x| x**2} - # Integration over [1,2]=x^3/3=7/3 - methods=[:rectangle,:trapezoid, :simpson, :adaptive_quadrature, :romberg] - methods.each do |m| - it "should integrate correctly with ruby method #{m}" do - Integration.integrate(1,2,{:method=>m,:tolerance=>1e-8},&a).should be_within(1e-6).of(7.0 / 3 ) + let(:function) { ->(x) { x**2 } } + METHODS = [:rectangle, :trapezoid, :simpson, :adaptive_quadrature, :romberg] + let(:tolerance) { 1e-8 } + context 'over [1,2] = x^3/3 = 7/3' do + METHODS.each do |m| + it "calculates correct result with method #{m}" do + res = Integration.integrate(1, 2, { method: m, tolerance: tolerance }, &function) + expect(res).to be_within(1e-6).of(7.0 / 3) + end end end - - it "should return correct for trapezoid" do - a=rand() - b=rand()*10 - f=lambda {|x| x*a} - Integration.trapezoid(0,b,2,&f).should be_within(1e-14).of((a*b**2) / 2.0) - end - it "should return a correct value for a complex integration with ruby methods" do - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - Integration.integrate(0,1,{:tolerance=>1e-12,:method=>:simpson},&normal_pdf).should be_within(1e-11).of(0.341344746068) - Integration.integrate(0,1,{:tolerance=>1e-12,:method=>:adaptive_quadrature},&normal_pdf).should be_within(1e-11).of(0.341344746068) - end - it "should return a correct value for a complex integration with gsl methods" do - if Integration.has_gsl? - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - Integration.integrate(0,1,{:tolerance=>1e-12,:method=>:qng},&normal_pdf).should be_within(1e-11).of(0.341344746068) - Integration.integrate(0,1,{:tolerance=>1e-12,:method=>:qag},&normal_pdf).should be_within(1e-11).of(0.341344746068) - else - pending("GSL not available") + context 'Trapezoid' do + let(:a) { rand } + let(:b) { rand * 10 } + let(:function) { ->(x) { x * a } } + it 'calculates valid result' do + res = Integration.integrate(0, b, { method: :trapezoid }, &function) + expect(res).to be_within(1e-14).of((a * b**2) / 2.0) end end - - it "should return correct integration for infinity bounds" do - if Integration.has_gsl? - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} + context 'Ruby methods' do + let(:function) { ->(x) { (1.0 / Math.sqrt(2.0 * Math::PI)) * Math.exp(-(x**2.0 / 2.0)) } } + context 'Simpson' do + let(:method) { :simpson } + it 'calculates correct value' do + res = Integration.integrate(0, 1, { tolerance: 1e-12, method: method }, &function) + expect(res).to be_within(1e-11).of(0.341344746068) + end + end - Integration.integrate(Integration::MInfinity, Integration::Infinity,{:tolerance=>1e-10}, &normal_pdf).should be_within(1e-09).of(1) - else - pending("GSL not available") + context 'Adaptive quadrature' do + let(:method) { :adaptive_quadrature } + it 'calculates correct value' do + res = Integration.integrate(0, 1, { tolerance: 1e-12, method: method }, &function) + expect(res).to be_within(1e-11).of(0.341344746068) + end end end - it "should return correct integration for infinity lower bound" do - if Integration.has_gsl? - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - Integration.integrate(Integration::MInfinity, 0 , {:tolerance=>1e-10}, &normal_pdf).should be_within(1e-09).of(0.5) + context 'GSL methods' do + let(:function) { ->(x) { (1.0 / Math.sqrt(2.0 * Math::PI)) * Math.exp(-(x**2.0 / 2.0)) } } + let(:lower) { 0 } + let(:upper) { 1 } + let(:tolerance) { 1e-12 } + let(:imethod) { :qng } + let(:params) do + [ + ] + end + context 'QNG' do + let(:imethod) { :qng } + it 'calculates correct value' do + pending('GSL not available') unless Integration.has_gsl? + res = Integration.integrate(lower, upper, { tolerance: tolerance, method: imethod }, &function) + expect(res).to be_within(1e-11).of(0.341344746068) + end + end - else - pending("GSL not available") + context 'QAG' do + let(:imethod) { :qag } + it 'calculates correct value' do + pending('GSL not available') unless Integration.has_gsl? + res = Integration.integrate(lower, upper, { tolerance: tolerance, method: imethod }, &function) + expect(res).to be_within(1e-11).of(0.341344746068) + end end - end - it "should return correct integration for infinity upper bound" do - if Integration.has_gsl? - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - Integration.integrate(0,Integration::Infinity,{:tolerance=>1e-10}, &normal_pdf).should be_within(1e-09).of(0.5) + context 'Infinity' do + let(:imethod) { :qag } + let(:lower) { Integration::M_INFINITY } + let(:upper) { Integration::INFINITY } + context 'both' do + let(:tolerance) { 1e-10 } + it 'calculates correct value' do + pending('GSL not available') unless Integration.has_gsl? + res = Integration.integrate(lower, upper, { tolerance: tolerance, method: imethod }, &function) + expect(res).to be_within(1e-09).of(1) + end + end + + context 'lower' do + let(:upper) { 0 } + let(:tolerance) { 1e-10 } + it 'calculates correct value' do + pending('GSL not available') unless Integration.has_gsl? + res = Integration.integrate(lower, upper, { tolerance: tolerance, method: imethod }, &function) + expect(res).to be_within(1e-09).of(0.5) + end + end - else - pending("GSL not available") + context 'upper' do + let(:lower) { 0 } + let(:tolerance) { 1e-10 } + it 'calculates correct value' do + pending('GSL not available') unless Integration.has_gsl? + res = Integration.integrate(lower, upper, { tolerance: tolerance, method: imethod }, &function) + expect(res).to be_within(1e-09).of(0.5) + end + end + + it 'raises error if a ruby method is called with infinite bounds' do + res = -> { Integration.integrate(0, upper, { method: :simpson }, &function) } + expect(res).to raise_exception + end end end - it "should raise an error if a ruby methods is called with infinite bounds" do - normal_pdf=lambda {|x| (1/Math.sqrt(2*Math::PI))*Math.exp(-(x**2/2))} - lambda {Integration.integrate(0,Integration::Infinity,{:method=>:simpson}, &normal_pdf).should be_within(1e-09).of(0.5)}.should raise_exception() - end end - diff --git a/spec/spec_helper.rb b/spec/spec_helper.rb index 70f33ec..9b1f1ae 100644 --- a/spec/spec_helper.rb +++ b/spec/spec_helper.rb @@ -1,13 +1,5 @@ -$:.unshift(File.dirname(__FILE__)+"/../lib") -begin - require 'simplecov' - SimpleCov.start do - add_filter "/spec/" - add_group "Libraries", "lib" - end -rescue LoadError -end -require 'rspec' require 'integration' - +RSpec.configure do |config| + config.order = :random +end From cf7012e193ff658f069873fe39ae1fea222e77c2 Mon Sep 17 00:00:00 2001 From: Oleg Bovykin Date: Sun, 15 Jun 2014 19:15:19 +0400 Subject: [PATCH 2/2] Total library rewrite --- .autotest | 23 --------------- .travis.yml | 11 +++++++ README.md | 3 +- integration.gemspec | 8 ++--- lib/integration.rb | 4 +-- .../integrator/adaptive_quadrature.rb | 29 ++++++++++--------- lib/integration/integrator/monte_carlo.rb | 4 --- lib/integration/ruby.rb | 1 + lib/integration/version.rb | 2 +- 9 files changed, 35 insertions(+), 50 deletions(-) delete mode 100644 .autotest create mode 100644 .travis.yml diff --git a/.autotest b/.autotest deleted file mode 100644 index ef753ad..0000000 --- a/.autotest +++ /dev/null @@ -1,23 +0,0 @@ -# -*- ruby -*- - -require 'autotest/restart' - -# Autotest.add_hook :initialize do |at| -# at.extra_files << "../some/external/dependency.rb" -# -# at.libs << ":../some/external" -# -# at.add_exception 'vendor' -# -# at.add_mapping(/dependency.rb/) do |f, _| -# at.files_matching(/test_.*rb$/) -# end -# -# %w(TestA TestB).each do |klass| -# at.extra_class_map[klass] = "test/test_misc.rb" -# end -# end - -# Autotest.add_hook :run_command do |at| -# system "rake build" -# end diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..e6ad0f2 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,11 @@ +script: rspec +language: ruby +notifications: + email: false +rvm: + - 1.9 + - 2.0 + - 2.1 +before_install: + - sudo apt-get update -qq + - sudo apt-get install -y gsl-bin libgsl0-dev diff --git a/README.md b/README.md index d2daed8..a3c5dce 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,6 @@ # Integration +[![Build Status](https://travis-ci.org/arrowcircle/integration.svg?branch=master)](https://travis-ci.org/arrowcircle/integration) Numerical integration for Ruby, with a simple interface ## Installation @@ -33,7 +34,7 @@ Integration.integrate(0, Integration::Infinity, {tolerance: 1e-10}, &func) ## Available methods Pure Ruby methods: -* Simpson (`:simpson`, default method) +* Simpson (`:simpson`, default method) * Rectangular (`:rectangular`) * Trapezoidal (`:trapezoidal`) * Adaptive quadrature (`:adaptive_quadrature`) diff --git a/integration.gemspec b/integration.gemspec index b802bee..3333c6e 100644 --- a/integration.gemspec +++ b/integration.gemspec @@ -1,11 +1,11 @@ # coding: utf-8 lib = File.expand_path('../lib', __FILE__) $LOAD_PATH.unshift(lib) unless $LOAD_PATH.include?(lib) -require 'integration' +require 'integration/version' Gem::Specification.new do |spec| spec.name = 'integration' - spec.version = Integration::Version::VERSION + spec.version = IntegrationLib::VERSION spec.authors = ['Ben Gimpert', 'Claudio Bustos', 'Oleg Bovykin'] spec.email = ['clbustos_at_gmail.com', 'oleg.bovykin@gmail.com'] spec.summary = %q{Integration methods, based on original work by Beng} @@ -22,6 +22,6 @@ Gem::Specification.new do |spec| spec.add_development_dependency 'rake' spec.add_development_dependency 'rspec' - spec.add_dependency 'rb-gsl' - spec.add_dependency 'activesupport' + spec.add_runtime_dependency 'rb-gsl' + spec.add_runtime_dependency 'activesupport' end diff --git a/lib/integration.rb b/lib/integration.rb index 1675478..11e4d48 100644 --- a/lib/integration.rb +++ b/lib/integration.rb @@ -1,7 +1,5 @@ -require 'active_support/core_ext/string' - class Integration - require 'integration/version' + # require 'integration/version' require 'integration/gsl' require 'integration/ruby' diff --git a/lib/integration/integrator/adaptive_quadrature.rb b/lib/integration/integrator/adaptive_quadrature.rb index da9bad9..8784ca8 100644 --- a/lib/integration/integrator/adaptive_quadrature.rb +++ b/lib/integration/integrator/adaptive_quadrature.rb @@ -5,22 +5,23 @@ def result fc = func[lower_bound + h] fb = func[upper_bound] s = h * (fa + (4 * fc) + fb) / 3 - helper = Proc.new do |a, b, fa, fb, fc, h, s, level| - if level < 1 / n.to_f - fd = func[a + (h / 2)] - fe = func[a + (3 * (h / 2))] - s1 = h * (fa + (4.0 * fd) + fc) / 6 - s2 = h * (fc + (4.0 * fe) + fb) / 6 - if ((s1 + s2) - s).abs <= n - s1 + s2 - else - helper.call(a, a + h, fa, fc, fd, h / 2, s1, level + 1) + - helper.call(a + h, b, fc, fb, fe, h / 2, s2, level + 1) - end + calc(lower_bound, upper_bound, fa, fb, fc, h, s, 1) + end + + def calc(a, b, fa, fb, fc, h, s, level) + if level < 1 / n.to_f + fd = func[a + (h / 2)] + fe = func[a + (3 * (h / 2))] + s1 = h * (fa + (4.0 * fd) + fc) / 6 + s2 = h * (fc + (4.0 * fe) + fb) / 6 + if ((s1 + s2) - s).abs <= n + s1 + s2 else - fail 'Integral did not converge' + calc(a, a + h, fa, fc, fd, h / 2, s1, level + 1) + + calc(a + h, b, fc, fb, fe, h / 2, s2, level + 1) end + else + fail 'Integral did not converge' end - return helper.call(lower_bound, upper_bound, fa, fb, fc, h, s, 1) end end diff --git a/lib/integration/integrator/monte_carlo.rb b/lib/integration/integrator/monte_carlo.rb index 4bc67ae..b524222 100644 --- a/lib/integration/integrator/monte_carlo.rb +++ b/lib/integration/integrator/monte_carlo.rb @@ -32,8 +32,4 @@ def area_ratio def width (upper_bound - lower_bound).to_f end - - def height - max - min - end end diff --git a/lib/integration/ruby.rb b/lib/integration/ruby.rb index d7b082a..acaddfb 100644 --- a/lib/integration/ruby.rb +++ b/lib/integration/ruby.rb @@ -1,4 +1,5 @@ require 'integration/integrator' +require 'active_support/core_ext/string' class Integration::Ruby M_INFINITY = :minfinity INFINITY = :infinity diff --git a/lib/integration/version.rb b/lib/integration/version.rb index 78f7750..c19d5b5 100644 --- a/lib/integration/version.rb +++ b/lib/integration/version.rb @@ -1,3 +1,3 @@ -module Integration::Version +module IntegrationLib VERSION = '0.1.1' end