fixed errors in rosenbrock

Author: MatteoRagni

.gitignore

@@ -1,4 +1,4 @@
 *.o
 *.so
-fuse_problem
+fuse.problem
 .fuse_hidden*

Makefile

@@ -36,7 +36,7 @@ all: fuse.problem
 	$(GXX) $(LDFLAGS) -o fuse.problem problem.o fuse.problem.o
 
 fuse.problem: problem
-	$(GCC) -c main.c $(CFLAGS) $(CFUSE) -o fuse.problem.o
+	$(GCC) -c src/main.c $(CFLAGS) $(CFUSE) -o fuse.problem.o
 
 problem:
 	$(GXX) -c src/problem.cpp $(CFLAGS) -o problem.o

README.md

@@ -85,5 +85,39 @@ File system **must** always be unmounted, or it will continue to run in backgrou
 
 ## Example
 
-Let's run `librosenbrock` to test usage, with a ruby script that reads and writes files.
-TO BE FINISHED
+Let's run `librosenbrock` to test usage, with a ruby script that reads and writes files. The script implements a Newton method on the Rosenbrock function.
+
+The necessary file system is loaded as:
+```sh
+make all
+make rosenbrock
+
+./fuse.problem -s /tmp/mount shared/librosenbrock.so
+```
+
+Then `rosenbrock.rb` can be run as test problem. The core elements to interface file system with Ruby language are the three functions:
+```ruby
+def write(x)
+  File.open($x, "w") { |xf| xf.puts x.map(&:to_s).join("\n") }
+end
+```
+for writing a new value (it also triggers a re-evaluation),
+```ruby
+def write_p(p)
+  File.open($p, "w") { |pf|  pf.puts p.map(&:to_s).join("\n") }
+end
+```
+for writing a parameter, and
+```ruby
+def read()
+  return [
+    File.open($df, "r").read.chomp.split("\n").map(&:to_f),
+    File.open($ddf, "r").read.chomp.split("\n").map(&:to_f)
+  ]
+end
+
+# [...]
+
+File.open($f, "r").read.chomp
+```
+for reading gradient, hessian (function `read`) and reading the function value (the last line).

include/problem.h

@@ -120,7 +120,7 @@ template <class FP> class Problem {
   map<string, int> paths; /**< path to internal representation (enums) mapping */
 
   size_t x_s; /**< declared dimension for input x */
-  size_t p_s  /**< declared dimension for parameter p */
+  size_t p_s; /**< declared dimension for parameter p */
   size_t f_s; /**< number of function declared in shared object */
 
   vector<Function> f; /**< vector of function pointer */

shared/librosenbrock.cpp

@@ -42,8 +42,8 @@
  * $$
  * \\nabla f = \\left[
  * \\begin{array}{c}
- * 4 a x^3 - 4 b x y - 2 a + 2 x \\\\
- * 3 b (y - x^2)
+ * 2 x - 2 a - 4 b x (-x^2 + y) \\\\
+ * 2 b (y - x^2)
  * \\end{array}
  * \\right]
  * $$
@@ -53,8 +53,10 @@
  * $$
  * \\nabla ^2 f = \\left[
  * \\begin{array}{cc}
- * 2 + (12 x - 4 y) b & -4 b x \\\\
- * -4 b x & 2 b
+ * 2 + 8 b x^2 - 4 b (-x^2 + y) - 4 b x &
+ * -4 b x \\\\
+ * -4 b x &
+ * 2 b
  * \\end{array}
  * \\right]
  * $$
@@ -97,7 +99,7 @@ FUNCTION_DEFINE(f_rosenbrock, {
 FUNCTION_DEFINE(df_rosenbrock, {
   Y().clear();
   Y().push_back(
-    4.0 * P()[1] * pow(X()[0], 3.0) - 4.0 * P()[1] * X()[0] * X()[1] - 2.0 * P()[0] + 2.0 * X()[0]
+    2 * X()[0] - 2 * P()[0] - 4 * P()[1] * X()[0] * (-pow(X()[0], 2.0) + X()[1])
   );
   Y().push_back(
     2.0 * P()[1] * (X()[1] - pow(X()[0], 2.0))
@@ -110,7 +112,7 @@ FUNCTION_DEFINE(ddf_rosenbrock, {
   PROBLEM_PRECISION ddf01 = -4.0 * P()[1] * X()[0];
 
   Y().push_back(
-    2.0 + (12.0 * pow(X()[0], 2.0) - 4.0 * X()[1]) * P()[1]
+    2.0 + 8.0 * P()[1] * pow(X()[0], 2) - 4 * P()[1] * (X()[1] - pow(X()[0], 2))
   );
   Y().push_back(ddf01);
   Y().push_back(ddf01);

shared/rosenbrock.rb

@@ -5,39 +5,75 @@
 # dimensions and we already know that f returns only a
 # value, while gradient returns a vector.
 
-$mount = "/tmp/mount/"
+$mount = "/tmp/mount/"  # base mountpoint
 $x     = $mount + "x"   # path for input
 $p     = $mount + "p"   # path for parameters
+$f     = $mount + "y0"  # path for function
 $df    = $mount + "y1"  # path for gradient
 $ddf   = $mount + "y2"  # path for hessian
 
 # The next function only writes new parameters
 # in parameter file
 def write_p(p)
-  File.open($p, "w") do |pf|
-    pf.puts p.map(&:to_s).join("\n")
-  end
+  File.open($p, "w") { |pf|  pf.puts p.map(&:to_s).join("\n") }
 end
 
 # The next function only writes new input to the input
 # file
 def write(x)
-  File.open($x, "w") do |xf|
-    xf.puts x.map(&:to_s).join("\n")
-  end
+  File.open($x, "w") { |xf| xf.puts x.map(&:to_s).join("\n") }
 end
 
 # The next function reads the gradient and hessian, thus returns
-# an array of array of values
+# an array of array of values. This is what we need to evaluate a
+# Newton algorithm
 def read()
-  df = Files($df, "r").read.chomp.split("\n").map(&:to_f)
-  ddf = Files($ddf, "r").read.chomp.split("\n").map(&:to_f)
-  return df, ddf
+  return [
+    File.open($df, "r").read.chomp.split("\n").map(&:to_f),
+    File.open($ddf, "r").read.chomp.split("\n").map(&:to_f)
+  ]
 end
 
+# This function performs a newton step in the form:
+#
+# x[+] = x - (hessian(f)(x))^(-1) gradient(f)(x)
 def newton_step(x, df, ddf)
+  dt = 1.0 / (ddf[0] * ddf[3] - ddf[1] * ddf[2])
   return [
-    x[0] - (ddf[0] * df[0] + ddf[1] * df[1]),
-    x[1] - (ddf[2] * df[0] + ddf[3] * df[1])
+    x[0] - dt * (ddf[3] * df[0] - ddf[2] * df[1]),
+    x[1] - dt * (ddf[0] * df[1] - ddf[1] * df[0])
   ]
 end
+
+# Evaluates the norm two of 2D vectors
+def norm(a, b)
+  return ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** (0.5)
+end
+
+if __FILE__ == $0
+  # Maximum iteration
+  ITERMAX = 1000
+  # Tolerance. If ||x[+] - x|| < tol we stop execution
+  # of iterative algorithm.
+  TOL = 1e-5
+
+  # Starting Point
+  x = z = [0.0, 0.0]
+  # Parameter configuration and write configuration
+  write_p([1.0, 100.0])
+
+  # Newton Algorithm
+  i = 0
+  for i in (0..ITERMAX)
+    write(z) # write will evaluate again the function
+    puts(z.map(&:to_s).join(","))
+    df, ddf = read()
+    z = newton_step(x, df, ddf)
+    break if norm(z, x) < TOL
+  end
+
+  # Print on screen results
+  puts "Itermax reached" if i == ITERMAX
+  puts "Solution (after #{i} iterations):"
+  puts "  f(#{z.map(&:to_s).join(", ")}) = #{File.open($f, "r").read.chomp}"
+end # __FILE__ == $0