MatteoRagni/cas-rb

View on GitHub
README.md

Summary

Maintainability
Test Coverage
# Mr.CAS

[![Gem Version](https://badge.fury.io/rb/Mr.CAS.svg)](https://badge.fury.io/rb/Mr.CAS)
[![Code Climate](https://codeclimate.com/github/MatteoRagni/cas-rb/badges/gpa.svg)](https://codeclimate.com/github/MatteoRagni/cas-rb)
[![Build Status](https://travis-ci.org/MatteoRagni/cas-rb.svg?branch=master)](https://travis-ci.org/MatteoRagni/cas-rb)

## Introduction

An extremely simple graph, that handles only differentiation and simplification with one step ahead in the graph.

## Module Structure

![module structure image](https://gist.githubusercontent.com/MatteoRagni/639294457d97858bbb3cee2ffaf782dc/raw/2854e8877f2106307bfa190649e4253c70282139/graphviz.png)

## Example

Given the function of Rosenbrock, find the optimum of such a function.

## Execution

### Installation

First of all is necessary to install and load the `Mr.CAS` gem.

``` bash
gem install Mr.CAS
```

``` ruby
require 'Mr.CAS'
```

### Define a function

Now we can define the Rosenbrock function, that has two variable (`x` and `y`) and
two constants (`a` and `b`, with the default value on 1 and 100 in this case). The
variable `f` will contain our function:

``` ruby
x, y = CAS::vars :x, :y
a, b = CAS::const 1.0, 100.0

f = ((a - x) ** 2) + b * ((y - (x ** 2)) ** 2)
```

We can print this function as follows:

``` ruby
puts "f = #{f}"
# => ((1.0 - x)^2 + (100.0 * (y - x^2)^2))
```

### Derivatives

To find the minimum we need  to find the point in which the gradient of such a function is
equal to zero

The two derivatives are:

``` ruby
dfx = f.diff(x).simplify
dfy = f.diff(y).simplify

puts "df/dx = #{dfx}"
# => df/dx = ((((1 - x) * 2) * -1) + ((((y - x^2) * 2) * -(x * 2)) * 100.0))

puts "df/dy = #{dfy}"
# => df/dy = (((y - x^2) * 2) * 100.0)
```

### Substitutions

Now, from the second it is quite evident that

```
y = x^2 = g(x)
```
and thus we can perform on the first a substitution:

``` ruby
g = (x ** 2)
dfx = dfx.subs({y => g})

puts "df/dx = #{dfx}"
# => df/dx = (((1 - x) * 2) * -1)
```

### Meta-Programming (a Newton algorithm...)

That is something quite simple to solve (x = 1). Let's use this formulation for some meta-programming anyway! The stationary point is in the root of the function `dfx`, that depends
only on one variable:

``` ruby
puts "Arguments: #{dfx.args}"
```

We can write a Newton method. The solution is found iteratively solving the recursive equation:

```
x[k + 1] = x[k] - ( f(x[k]) / f'(x[k]) )
```

Let's write an algorithm that takes the symblic expression, and creates its own method
(e.g.: using `Proc` objects). Here an example with our function:

``` ruby
unused = (dfx/dfx.diff(x)).simplify
puts "unused(x) = #{unused}"
# => unused(x) = ((((1 - x) * 2) * -1)) / (((-1 * 2) * -1))

unused_proc = unused.as_proc
puts "unused(12.0) = #{unused.call({"x" => 12.0})}"
# => unused(12.0) = 11.0
```

First of all, let's write a function that contains the algorithm. We want to give as input
a symbolic function of one variable, it must reply with the value of the variable in which the function as value equal to zero:

``` ruby
def newton(f, init=3.0, tol=1E-8, kmax=100)
  k = 0

  x   = f.args[0]
  s   = {x.name => init}      # <-This is the solution dictionary

  fp  = f.as_proc             # <- This is a parsed object
  df  = f.diff(x).simplify    # <- This is still symbolic
  res = (f/df).as_proc        # <- This is a parsed object

  f0 = fp.call(s)
  k = 0

  loop do
    # Algorithm
    s[x.name] = s[x.name] - res.call(s)
    # Tolerance check
    f1        = fp.call(s)
    if (f1 - f0).abs < tol
      break
    else
      f0 = f1
    end

    # Iterations check
    k = k + 1
    break if k > kmax
  end
  puts "Solution after #{k} iterations: #{x} = #{s[x.name]}, f(#{x}) = #{f0}"
  return s[x.name]
end
```

**Transforming the symbolic function into a `Proc` makes the code more efficient.
It is not mandatory, but higlhy suggested in case of iterative algorithms**.

Let's call our new function, using as argument the derivative of the function
that we want to optimize:

``` ruby
x_opt = newton(dfx)
# => Solution after 1 iterations: x = 1.0, f(x) = -0.0
```

We can use the solution of `x`, to get the value of `y`:

``` ruby
puts "Optimum in #{x} = #{x_opt} and #{y} = #{g.call({x => x_opt})}"
# => Optimum in x = 1.0 and y = 1.0
```

## Disclaimer

This is a work in progress and mainly a proof of concept.
What really is missing is a ~~graph~~ a way to perform better simplifications.

## Full example code

``` ruby
require 'Mr.CAS'

# Define the function
x, y = CAS::vars "x", "y"
a, b = CAS::const 1.0, 100.0

f = ((a - x) ** 2) + b * ((y - (x ** 2)) ** 2)

puts "f = #{f}"

# Derivate wrt variables
dfx = f.diff(x).simplify
dfy = f.diff(y).simplify

puts "df/dx = #{dfx}"
puts "df/dy = #{dfy}"

# Perform substitutions
g = (x ** 2)
dfx = dfx.subs({y => g}).simplify
puts "df/dx = #{dfx}"

# Arguments of a function
puts "Arguments: #{dfx.args}"

# Metaprogramming: Create a Proc from symbolic math
unused = (dfx/dfx.diff(x)).simplify
puts "unused(x) = #{unused}"
unused_proc = unused.as_proc

# Testing the Proc. The feed must have string as key to identify
# the variable
puts "unused(12.0) = #{unused.call({"x" => 12.0})}"

# We will not use this function, instead we will let the
# algorithm to create its own


# Newton algorthm on the fly, that creates its own formula!
def newton(f, init=3.0, tol=1E-8, kmax=100)
  k = 0

  x   = f.args[0]
  s   = {x.name => init}      # <-This is the solution dictionary

  fp  = f.as_proc             # <- This is a parsed object
  df  = f.diff(x).simplify    # <- This is still symbolic
  res = (f/df).as_proc        # <- This is a parsed object

  f0 = fp.call(s)
  k = 0

  loop do
    # Algorithm
    s[x.name] = s[x.name] - res.call(s)
    # Tolerance check
    f1        = fp.call(s)
    if (f1 - f0).abs < tol
      break
    else
      f0 = f1
    end

    # Iterations check
    k = k + 1
    break if k > kmax
  end
  puts "Solution after #{k} iterations: #{x} = #{s[x.name]}, f(#{x}) = #{f0}"
  return s[x.name]
end

# Let's call our shining algorithm:
x_opt = newton(dfx)

# Let's see the final result for both:
puts "Optimum in #{x} = #{x_opt} and #{y} = #{g.call({x => x_opt})}"
```


## To do

#### Simple

 * [ ] substitute `CAS::Invert` with `CAS::Neg`
 * [ ] substitute `CAS::Condition` with `CAS::Expression`
 * [x] move code generation to external plugins that can be loaded on demand
 * [x] `CAS::Sum` and `CAS::Prod` inherits from `CAS::NaryOp` instead of `CAS::BinaryOp`

#### Complex

 * [ ] implement `CAS::VectorVariable` and `CAS::MatrixVariable`
 * [ ] Replace `CAS::Op` and `CAS::BinaryOp` with a single `CAS::Op` equal to `CAS::NaryOp`

##### Extremely complex

 * [ ] Implement the Risch algorithms for symbolic integrations