Skip to content

Instantly share code, notes, and snippets.

@MatteoRagni MatteoRagni/00_README.md
Last active Mar 13, 2017

Embed
What would you like to do?
Scripts for paper: M. Ragni - "Mr.CAS - A Minimalistic (pure) Ruby CAS for Fast Prototyping and Code Generation"

Mr.CAS - A Minimalistic (pure) Ruby CAS for Fast Prototyping and Code Generation

Abstract

There are complete Computer Algebra System (CAS) systems on the marketm with complete solutions for analysis of analytical models. But exporting a model, for optimization or any other research activity, requires a lot of work, even with a good software.

This work presents a Ruby library that exposes minimalistic CAS capabilities - i.e. simplifications, substitutions, evaluations, etc. Library aims at rapid prototyping of numerical interfaces and code generation for different target languages, separating mathematical expression from exportation rules - e.g. models from numerical conditioning best practices.

The library is implemented in pure Ruby language and compatible with all Ruby interpreter flavours.

Gist Repository

This Gist repository contains a working version of the scripts presented in the different listings of the paper. Due to restriction on number of words, some obvious part were omitted. All the scripts complete are presented here, numbered as in the paper.

Expected outputs in console are highlighted with =>

Tasks

There is a Rakefile that runs the scripts. For installing dependencies:

License for Gist Repository

Copyright (c) 2016 Matteo Ragni

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

This license covers all the file included in this repository

#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
z = CAS.vars 'z' # creates a variable
f = z ** 2 + 1 # define a symbolic expression
df = f.diff(z) # derivative w.r.t. z
puts df # => (((z)^((2 - 1)) * 2 * 1) + 0)
g = CAS.declare :g, f # creates implicit expression
dg = g.diff(z) # derivative w.r.t. z
puts dg # => ((((z)^((2 - 1)) * 2 * 1) + 0) * Dg[0](((z)^(2) + 1)))
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
x, y = CAS::vars 'x', 'y' # creates two variables
f = CAS.log( CAS.sin( y ) ) # symbolic expression
f.subs y => CAS.asin(CAS.exp(x)) # performs substitution
fs = f.simplify # simplifies expression
puts fs # => x
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
x = CAS.vars 'x' # creates a variable
f = x ** 2 + 1 # defines a symbolic expression
fc = f.call x => 2 # evaluates for x = 2
puts fc # => 5.0
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
x, y = CAS.vars 'x', 'y' # creates two variables
f = CAS.declare :f, x # declares implicit function :f
g = CAS.declare :g, x, y # declares implicit function :g
fg = f.greater_equal g # :f is greater than :g
puts fg # => (f(x) ≥ g(x, y))
max_fg = CAS::max f, g
puts max_fg # => ((f(x) ≥ g(x, y)) ? f(x) : g(x, y))
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
x = CAS::vars 'x' # creates a variable
f = CAS::log(CAS::sin(x)) # defines a symbolic function
proc = f.as_proc # compiles callable closure
r = proc.call('x' => Math::PI/2) # runs callable closure
puts r # => 0.0
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
# Rules definition for Fortran Language
module CAS
{
# . . .
CAS::Variable => Proc.new { "#{name}" }, # rule for Variable operation
CAS::Sin => Proc.new { "sin(#{x.to_fortran})" }, # rule for Sin operation
# . . .
}.each do |cls, prc|
cls.send(:define_method, :to_fortran, &prc) # Appends new to_fortran methods
end
end
# Usage
x = CAS.vars 'x' # creates a variable
code = (CAS.sin(x)).to_fortran # converts to target language
puts code # => sin(x)
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
%w|Mr.CAS
Mr.CAS/c-opt|.each do |gem|
raise LoadError, "Cannot load #{gem}" unless require gem
end
# Model
x, y = CAS.vars :x, :y # defines two variables
g = CAS.declare :g, x # declares implicit function g(x)
f = x ** y + g * CAS.log(CAS.sin(x ** y)) # f(x,y) = x^y + g(x) * log(sin(x^y))
# Code Generation
g.c_name = 'g_impl' # g(x) implementation token
clib = CAS::CLib.create "example" do
include_local "g_impl" # includes local g_impl header
implements_as "f_impl", f # token in interface of f(x,y)
end
puts clib.header
=begin
=>
// Header file for library: example.c
#ifndef example_H
#define example_H
// Standard Libraries
#include <math.h>
// Local Libraries
#include "g_impl"
// Definitions
// Functions
double f_impl(double x, double y);
#endif // example_H
=end
puts clib.source
=begin
=>
// Source file for library: example.c
#include "example.h"
double f_impl(double x, double y) {
double __t_0 = pow(x, y);
double __t_1 = g_impl(x);
double __t_2 = sin(__t_0);
double __t_3 = log(__t_2);
double __t_4 = (__t_1 + __t_3);
double __t_5 = (__t_0 + __t_4);
return __t_5;
}
// end of example.c
// Source file for library: example.c
#include "example.h"
double f_impl(double x, double y) {
double __t_0 = pow(x, y);
double __t_1 = g_impl(x);
double __t_2 = sin(__t_0);
double __t_3 = log(__t_2);
double __t_4 = (__t_1 + __t_3);
double __t_5 = (__t_0 + __t_4);
return __t_5;
}
// end of example.c
=end
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
%w|Mr.CAS
Mr.CAS/c-opt|.each do |gem|
raise LoadError, "Cannot load #{gem}" unless require gem
end
x = CAS::vars :x
# Parameters can be defined as implicit functions with no arguments
a = CAS.declare "PARAM_A"
g = (CAS.sqrt(x + a) - CAS.sqrt(x)) + CAS.sqrt(CAS::Pi + x)
# Particular Code Generation for difference between square roots.
# If the exported finds a difference between two square roots, it perform
# normalization. Please note that written in this way (for a sake of brevity)
# breaks optimized code generation flow for the argument of the square roots.
module CAS
class Diff
alias :__to_c_impl_old :__to_c_impl
def __to_c_impl(v)
if @x.is_a? CAS::Sqrt and @y.is_a? CAS::Sqrt
"(#{@x.x.__to_c(v)} + #{@y.x.__to_c(v)}) / " +
"( #{@x.__to_c(v)} + #{@y.__to_c(v)} )"
else
self.__to_c_impl_old(v)
end
end
end
end
clib = CAS::CLib.create "g_impl" do
define "PARAM_A()", 1.0 # defines the token PARAM_A()
# 1.0 is an arbitrary value for PARAM_A
define "M_PI", Math::PI # Defines the token M_PI, with value of
# Pi taken from interpreter
implements_as "g_impl", g # token in interface of f(x,y)
end
puts clib.header
=begin
=>
// Header file for library: g_impl.c
#ifndef g_impl_H
#define g_impl_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
#define PARAM_A() 1.0
#define M_PI 3.141592653589793
// Functions
double g_impl(double x);
#endif // g_impl_H
// Header file for library: g_impl.c
#ifndef g_impl_H
#define g_impl_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
#define PARAM_A() 1.0
#define M_PI 3.141592653589793
// Functions
double g_impl(double x);
#endif // g_impl_H
=end
puts clib.source
=begin
=>
// Source file for library: g_impl.c
#include "g_impl.h"
double g_impl(double x) {
double __t_0 = PARAM_A();
double __t_1 = (x + __t_0);
double __t_2 = sqrt(__t_1);
double __t_3 = sqrt(x);
double __t_4 = (__t_1 + x) / ( __t_2 + __t_3 );
double __t_5 = (M_PI + x);
double __t_6 = sqrt(__t_5);
double __t_7 = (__t_4 + __t_6);
return __t_7;
}
// end of g_impl.c
// Source file for library: g_impl.c
#include "g_impl.h"
double g_impl(double x) {
double __t_0 = PARAM_A();
double __t_1 = (x + __t_0);
double __t_2 = sqrt(__t_1);
double __t_3 = sqrt(x);
double __t_4 = (__t_1 + x) / ( __t_2 + __t_3 );
double __t_5 = (M_PI + x);
double __t_6 = sqrt(__t_5);
double __t_7 = (__t_4 + __t_6);
return __t_7;
}
// end of g_impl.c
=end
puts g
# => ((√((x + PARAM_A())) - √(x)) + √((π + x)))
exit 0
#!/usr/bin/env python3
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
import sympy
import math
##
# Evaluates integral of symbolic function f
# starting from a, ending in b, with step size
# h = (b - a)/n
# Returns a float
def integrate(f, a, b, n):
h = (b - a)/n
x = sympy.symbols('x')
func = sympy.lambdify((x), f)
sums = (func(a) +
func(b)) / 2.0
for i in range(1, n):
sums += func(a + i*h)
return sums * h
##
# Evaluates order using a function f,
# a starting point a and endpoint b with
# mesh size n
# Returns a float
def order(f, a, b, n):
x = sympy.symbols('x')
f_ab = sympy.Subs(f, (x), (b)).n()-\
sympy.Subs(f, (x), (a)).n()
df = f.diff(x)
f_1n = integrate(df, a, b, n)
f_2n = integrate(df, a, b, 2 * n)
return math.log(
(f_ab - f_1n) /
(f_ab - f_2n),
2)
x = sympy.symbols('x')
f = sympy.atan(x)
print(order(f, -1.0, 1.0, 100))
# => 1.9999999974244451
exit(0)
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
##
# Evaluates integral of symbolic function f
# starting from a, ending in b, with step size
# h = (b - a)/n
# Returns a Float
def integrate(f, a, b, n)
h = (b - a) / n
func = f.as_proc
sum = ((func.call 'x' => a) +
(func.call 'x' => b)) / 2.0
for i in (1...n)
sum += (func.call 'x' => (a + i*h))
end
return sum * h
end
##
# Evaluates order using a function f,
# a starting point a and endpoint b with
# mesh size n
# Returns a Float
def order(f, a, b, n)
x = CAS.vars 'x'
f_ab = (f.call x => b) -
(f.call x => a)
df = f.diff(x).simplify
f_1n = integrate(df, a, b, n)
f_2n = integrate(df, a, b, 2 * n)
return Math.log(
(f_ab - f_1n) /
(f_ab - f_2n),
2)
end
x = CAS.vars 'x' # Defines a variable
f = CAS.arctan x # Defines arctan(x)
puts(order f, -1.0, 1.0, 100) # => 1.9999999974244451
exit 0
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
%w|Mr.CAS
Mr.CAS/c-opt|.each do |gem|
raise LoadError, "Cannot load #{gem}" unless require gem
end
$x, $y, $h = CAS::vars :x, :y, :h
# Evaluates the factorial value of a
# CAS::Constant
def fact(n); (n < 2 ? 1 : n * fact(n - 1)); end
# Evaluates the vector of derivatives that will be used
# in the next function.
def coeff(f, n)
df = [f]
for _ in 2..n
df << df[-1].diff($x).simplify + (df[-1].diff($y).simplify * df[0])
end
return df
end
# Generates the symbolic graph of the Taylor step
# that can be used by an integrator.
def taylor(f, n)
df = coeff(f, n)
y = $y
for i in 0...df.size
y = y + (($h ** (i + 1))/(fact(i + 1)) * df[i])
end
return y.simplify
end
# In this case the generator is tested against the
# function f = x * y
f = $x * $y
# And it is exported in a C library
clib = CAS::CLib.create "taylor" do
implements_as "taylor_step", taylor(f, 4)
end
puts clib.header
=begin
=>
#ifndef taylor_H
#define taylor_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
// Functions
double taylor_step(double y, double h, double x);
#endif // taylor_H
=end
puts clib.source
=begin
=>
// Source file for library: taylor.c
#include "taylor.h"
double taylor_step(double y, double h, double x) {
double __t_0 = (x + y);
double __t_1 = (h + __t_0);
double __t_2 = pow(h, 2.0);
double __t_3 = (__t_2) / (2.0 + M_EPSILON);
double __t_4 = (x + __t_0);
double __t_5 = (y + __t_4);
double __t_6 = (__t_3 + __t_5);
double __t_7 = pow(h, 3);
double __t_8 = (__t_7) / (6 + M_EPSILON);
double __t_9 = (y + x);
double __t_10 = pow(x, 2.0);
double __t_11 = (__t_10 + 1.0);
double __t_12 = (__t_11 + __t_0);
double __t_13 = (__t_0 + __t_9 + __t_12);
double __t_14 = (__t_8 + __t_13);
double __t_15 = pow(h, 4);
double __t_16 = (__t_15) / (24 + M_EPSILON);
double __t_17 = (2.0 + y);
double __t_18 = (2.0 + x);
double __t_19 = (__t_18 + __t_0);
double __t_20 = (y + __t_11);
double __t_21 = (__t_19 + __t_20);
double __t_22 = (x + __t_11);
double __t_23 = (__t_18 + __t_22);
double __t_24 = (__t_23 + __t_0);
double __t_25 = (__t_17 + __t_21 + __t_24);
double __t_26 = (__t_16 + __t_25);
double __t_27 = (y + __t_1 + __t_6 + __t_14 + __t_26);
return __t_27;
}
// end of taylor.c
=end
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
# Taylor expansion example
# Factorial
def fact(n); n < 2 ? 1 : n * fact(n - 1); end
# Taylor function. Takes a function as input,
# extracts the argument, and use it to perform
# the Taylor expansion, where f is a function
# and n is the order of the series
def Taylor(f, n)
# We can use this method to intercept
# the argument of the function and create
# the central point for the Taylor series
x = f.args()[0]
x0 = CAS::vars (x.name.to_s + "0")
# We create the derivatives
df = [f]
n.times do; df << df[-1].diff(x).simplify; end
# and finally we accumulate all the argument
# of the series, using factorial function
tf = 0
for i in 0..n
tf += (df[i].subs(x => x0) / fact(i) * ((x - x0) ** i)).simplify
end
# The usage of simplifications allows us to
# reduce drastically the size of the summation
return tf.simplify
end
# Let's test it
z = CAS::vars "z"
f = CAS::sin(z)
tf = Taylor(f, 3)
puts tf
# => (sin(z0) + # order 0
# (cos(z0) * (z - z0)) + # order 1
# ((-sin(z0)) / (2) * ((z - z0))^(2)) + # order 2
# ((-cos(z0)) / (6) * ((z - z0))^(3))) # order 3
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
raise LoadError, 'Please install Mr.CAS first' unless require 'Mr.CAS'
# This module contains an exporter for Latex language.
# The exporter is not complete, but shows the capabilities
# and the semplicity on how a language exporter can be written
module CAS
{
# Terminal nodes
CAS::Constant => Proc.new { "#{x}" },
CAS::Variable => Proc.new { "#{name}" },
CAS::PI_CONSTANT => Proc.new { "\\pi" },
CAS::INFINITY_CONSTANT => Proc.new { "\\infty" },
CAS::NEG_INFINITY_CONSTANT => Proc.new { "-\\infty" },
# Base functions
CAS::Sum => Proc.new { "\\left( #{x.map(&:latex).join(" + ")} \\right)" },
CAS::Diff => Proc.new { "\\left( #{x.latex} - #{y.latex} \\right)" },
CAS::Prod => Proc.new { "\\left( #{x.map(&:latex).join(" ")} \\right)"},
CAS::Pow => Proc.new { "{#{x.latex}}^{#{y.latex}}" },
CAS::Div => Proc.new { "\\dfrac{#{x.latex}}{#{y.latex}}" },
CAS::Sqrt => Proc.new { "\\sqrt{#{x.latex}}" },
CAS::Invert => Proc.new { "-#{x.latex}" },
CAS::Abs => Proc.new { "\\left| #{x.latex} \\right|" },
# Trigonometric functions
CAS::Sin => Proc.new { "\\sin \\left( #{x.latex} \\right)" },
CAS::Asin => Proc.new { "\\arcsin \\left( #{x.latex} \\right)" },
CAS::Cos => Proc.new { "\\cos \\left( #{x.latex} \\right)" },
CAS::Acos => Proc.new { "\\arccos \\left( #{x.latex} \\right)" },
CAS::Tan => Proc.new { "\\tan \\left( #{x.latex} \\right)" },
CAS::Atan => Proc.new { "\\arctan \\left( #{x.latex} \\right)" },
# Trascendent functions
CAS::Exp => Proc.new { "e^#{x.latex}" },
CAS::Ln => Proc.new { "\\log \\left( #{x.latex} \\right)" },
# Box Conditions
CAS::BoxConditionOpen => Proc.new { "#{lower.latex} < #{x.latex} < #{upper.latex}" },
CAS::BoxConditionClosed => Proc.new { "#{lower.latex} \\leq #{x.latex} \\leq #{upper.latex}" },
CAS::BoxConditionUpperClosed => Proc.new { "#{lower.latex} < #{x.latex} \\leq #{upper.latex}" },
CAS::BoxConditionLowerClosed => Proc.new { "#{lower.latex} \\leq #{x.latex} < #{upper.latex}" },
# Conditions
CAS::Equal => Proc.new { "#{x.latex} = #{y.latex}" },
CAS::Smaller => Proc.new { "#{x.latex} < #{y.latex}" },
CAS::Greater => Proc.new { "#{x.latex} > #{y.latex}" },
CAS::SmallerEqual => Proc.new { "#{x.latex} \\leq #{y.latex}" },
CAS::GreaterEqual => Proc.new { "#{x.latex} \\geq #{y.latex}" },
# Piecewise
CAS::Piecewise => Proc.new {
"\\left\\{ " +
" \\begin{array}{lr} " +
" #{x.latex} & #{condition.latex} \\\\" +
" #{y.latex}" +
" \\end{array}" +
"\\right."
},
CAS::Max => Proc.new { "\\max \\left( #{x.latex}, \\, #{y.latex} \\right)" },
CAS::Min => Proc.new { "\\min \\left( #{x.latex}, \\, #{y.latex} \\right)" },
CAS::Function => Proc.new { "#{name}\\left( #{x.map(&:latex).join(", ")} \\right)" }
}.each do |cls, blk|
cls.send(:define_method, "latex", &blk)
end
end
x, y = CAS::vars "x", "y"
f = CAS::declare "f", x, y
g = CAS::declare "g", x, y
h = CAS::Piecewise.new(
f + g ** 2,
g - CAS::sqrt(f),
x.greater_equal(y))
puts h.latex
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
%w|Mr.CAS
Mr.CAS/c-opt|.each do |gem|
raise LoadError, "Cannot load #{gem}" unless require gem
end
# In this example a new Op (that implements the Sign function)
# will be implemented. A rule for the exportation in C is derived
# and also a version optimized with regularization that overloads
# the first one.
module CAS
class Sign < CAS::Op
# The fisrt function that is needed is the function
# that is able to evaluate it
def call(f)
(@x.call(f) >= 0 ? 1 : -1)
end
# Now we formulate a derivative for the sign function
# The derivative is quite questionable, and may be defined
# using delta of dirac, or defined by considering
# sign(x) = |x|/x for x != 0. We will go with this
# second alternative.
def diff(z)
dx = @x.diff(z)
(dx)/(abs(@x)) * (1 - (abs(dx) ** 2)/(@x ** 2))
end
# It is possible to give also a simplification
# routine. For the sign simplification. The only simplification
# we need is related to constant values. If we
# have a constant, we can decide to eliminate the sign function
def simplify
@x = @x.simplify
if @x.is_a? CAS::Constant
return @x if @x.call({}) >= 0
return -@x
end
return self
end
# Now we present a string representation for the
# new op
def to_s; "sign(#{@x})"; end
# And a code representation that can ge interpreted by Ruby
def to_code
"(#{@x.to_code} >= 0 ? 1 : -1)"
end
end # end
# It is possible to add a shortcut for calling rapidly the
# op we just created
def self.sign(x); CAS::Sign.new(x); end
end
# As additional step we can decide to append to our C exporter
# the sign function. Since it is only one Op to be added, we can go
# with a simpler solution with respect to the one presented in
# other examples. We also define the OPTIMIZED version, that requires
# three functions
module CAS
class Sign
# This function is the actual implementation
# for the C code of the Op
def __to_c_impl(v)
"(#{@x.__to_c(v)} ? 1 : -1)"
end
# This function is responsible of the
# creation of a new variable if does not exist.
# it also creates the identification
# The internal variable are always referred with __t_ID.
def __to_c(v)
sym = self.to_s.to_sym
v[sym] = {
def: "#{self.__to_c_impl(v)}",
var: "__t_#{v.keys.size}",
id: v.keys.size
} unless v[sym]
v[sym][:var]
end
# This is the user interface that is called by user.
def to_c
v = {}; self.__to_c(v); v
end
end
end
if ($0 == __FILE__)
x, y = CAS::vars "x", "y"
f = CAS.sign(x ** y)
clib = CAS::CLib.create "sign_example" do
implements_as "sign_example", f
end
puts clib.header
=begin
=>
// Header file for library: sign_example.c
#ifndef sign_example_H
#define sign_example_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
// Functions
double sign_example(double x, double y);
#endif // sign_example_H
=end
puts clib.source
=begin
=>
// Source file for library: sign_example.c
#include "sign_example.h"
double sign_example(double x, double y) {
double __t_0 = pow(x, y);
double __t_1 = (__t_0 ? 1 : -1);
return __t_1;
}
// end of sign_example.c
=end
end
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
if ($0 == __FILE__)
load "14_Example_07.rb"
# Let's try to overload the sign op that we have defined in the
# example 14_Example_07.rb for exporting a stabilized version of
# the sign(x) => sin(arctan(a * x))
x = CAS::vars "x"
f = CAS.sign(x)
# This is the version not stabilized
clib_a = CAS::CLib.create "sign_unstable" do
implements_as "sign", f
end
# Let's write an exportation rule stabilized
module CAS; class Sign
def __to_c_impl(v)
"sin(atan(ALPHA * #{x.__to_c(v)}))"
end
end; end
clib_b = CAS::CLib.create "sign_stable" do
define "ALPHA", 1e4
implements_as "sign", f
end
# Another way can be the substitution of
# the sign function with a new symbolic version
# but this CHANGES the graph of f
a = CAS.const 1e4
f.subs(CAS.sign(x) => CAS.sin(CAS.arctan(a * x)))
clib_c = CAS::CLib.create "sign_symbolic" do
implements_as "sign", f
end
puts clib_a.header, clib_a.source
=begin
=>
// Header file for library: sign_unstable.c
#ifndef sign_unstable_H
#define sign_unstable_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
// Functions
double sign(double x);
#endif // sign_unstable_H
// Source file for library: sign_unstable.c
#include "sign_unstable.h"
double sign(double x) {
double __t_0 = sin(atan(ALPHA * x));
return __t_0;
}
// end of sign_unstable.c
=end
puts clib_b.header, clib_b.source
=begin
=>
// Header file for library: sign_stable.c
#ifndef sign_stable_H
#define sign_stable_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
#define ALPHA 10000.0
// Functions
double sign(double x);
#endif // sign_stable_H
// Source file for library: sign_stable.c
#include "sign_stable.h"
double sign(double x) {
double __t_0 = sin(atan(ALPHA * x));
return __t_0;
}
// end of sign_stable.c
=end
puts clib_c.header, clib_c.source
=begin
=>
// Header file for library: sign_symbolic.c
#ifndef sign_symbolic_H
#define sign_symbolic_H
// Standard Libraries
#include <math.h>
// Local Libraries
// Definitions
// Functions
double sign(double x);
#endif // sign_symbolic_H
// Source file for library: sign_symbolic.c
#include "sign_symbolic.h"
double sign(double x) {
double __t_0 = sin(atan(ALPHA * x));
return __t_0;
}
// end of sign_symbolic.c
=end
end
#!/usr/bin/env ruby
# License MIT - 2016 Matteo Ragni :: refer to 0_README.md
%w|Mr.CAS
Mr.CAS/auto-diff|.each do |gem|
raise LoadError, "Cannot load #{gem}" unless require gem
end
def NewtonRaphson(f, x, x0, config = {})
# configurations for the algorithm
cfg = {maxsteps: 40, tolerance: 1e-12}
# Configuration overrides
config.keys.each { |k| cfg[k] = config[k] }
# Starting Point
xv = x0
0.upto(cfg[:maxsteps]) do
fv = f.auto_diff({x => xv})
# We return the result if we reached
# the zero of the function
return xv if fv.x.abs < cfg[:tolerance]
xv = xv - fv.x/fv.y
end
# If we cannot reach a minimum we
# raise an argument error
raise RuntimeError, "Not Found"
end
x, y = CAS::vars "x", "y"
a, b = CAS::const 1.0, 100.0
f = ((a - x) ** 2) + b * ((y - (x ** 2)) ** 2)
f.subs({y => (x ** 2)}).simplify
x0 = NewtonRaphson(f, x, 0.0)
puts "#{f} = #{f.call({x => x0})} =~ 0 (for x = #{x0})"
desc "Install ruby dependencies (Mr.CAS.) for current user"
task :dep_ruby do
`gem install Mr.CAS`
end
desc "Install ruby dependencies (Mr.CAS.) system wide"
task :dep_ruby_root do
`sudo gem install Mr.CAS`
end
desc "Install python dependencies (sympy) for current user"
task :dep_py do
`pip3 install sympy`
end
desc "Install python dependencies (sympy) system wide"
task :dep_py do
`pip3 install sympy`
end
require 'pry-byebug'
task_list = []
Dir.glob("*.rb").sort.each do |script|
name = script.gsub(/\d+_(List|Ex)\w+_(\d+)\.rb/, '\1\2')
desc "Run Ruby script #{script}"
task name.to_sym do
puts "== Running Ruby task :: #{name} :: #{script} =="
puts `ruby #{script}`
puts
end
task_list << name.to_sym
end
desc "Run Python3 script 10_Example_03.py"
task :Ex03py do
puts "== Running Python task :: Example03py :: 10_Example_03.py =="
puts `python3 10_Example_03.py`
puts
end
task_list << :Ex03py
desc "Runs all scripts in paper's order of appearance"
task :default => task_list
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.