1. Lev Panov
  2. mo

Commits

Lev Panov  committed 5c07568

2d_min update

  • Participants
  • Parent commits 6bfce02
  • Branches master

Comments (0)

Files changed (3)

File labs/2d_min/2d_min.pro

View file
+SOURCES += \
+    min_2d.cpp
+
+OTHER_FILES += \
+    min_2d.c

File labs/2d_min/min_2d.c

View file
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include <stdlib.h>
-
-double Function( double x1, double x2 )
-{
-  double x1_2 = x1 * x1, x2_2 = x2 * x2;
-
-  return x1_2 * x1_2 + 2 * x2_2 * x2_2 + x1_2 * x2_2 + 2 * x1 + x2;
-}
-
-void Gradient( double x[], double y[] )
-{
-  y[0] = 4 * x[0] * x[0] * x[0] + 2 * x[0] * x[1] * x[1] + 2;
-  y[1] = 8 * x[1] * x[1] * x[1] + 2 * x[1] * x[0] * x[0] + 1;
-}
-
-void Jacobian( double x[], double J[2][2] )
-{
-  J[0][0] = 12 * x[0] * x[0] + 2 * x[1] * x[1];
-  J[0][1] = J[1][0] = 4 * x[0] * x[1];
-  J[1][1] = 24 * x[1] * x[1] + 2 * x[0] * x[0];
-}
-
-void Gradient2( double x[], double y[], int recalculate )
-{
-  static double H[2][2], invH[2][2], grad[2], det;
-
-  if (recalculate)
-  {
-    Jacobian(x, H);
-    det = H[0][0] * H[1][1] - H[0][1] * H[1][0];
-    if (det < 1e-5)
-    {
-      y[0] = y[1] = 0;
-      return;
-    }
-    invH[0][0] =  H[1][1] / det;
-    invH[0][1] = -H[0][1] / det;
-    invH[1][0] = -H[1][0] / det;
-    invH[1][1] =  H[0][0] / det;
-  }
-  Gradient(x, grad);
-  y[0] = (invH[0][0] * grad[0] + invH[0][1] * grad[1]);
-  y[1] = (invH[1][0] * grad[0] + invH[1][1] * grad[1]);
-}
-  
-double eps = 1e-6;
-
-double Step_subdivide( double x[], double y[], double *f, double norm )
-{
-  double step = 1, fnew, epsilon = eps;
-
-  if (norm < eps)
-    return 0;
-
-  do
-  {
-    step *= 0.5;
-    fnew = Function(x[0] - step * y[0], x[1] - step * y[1]);
-  } while (fnew - *f > -epsilon * step * norm); 
-
-  x[0] -= step * y[0];
-  x[1] -= step * y[1];
-  *f = fnew;
-  return step;
-}
-
-double Step_Newton( double x[], double y[], double *f, double norm )
-{
-  x[0] -= y[0];
-  x[1] -= y[1];
-  *f = Function(x[0], x[1]);
-  return 1;
-}
-
-double lipsh = 9.2;
-
-double Step_fixed( double x[], double y[], double *f, double norm )
-{
-  double step = (1 - eps) / lipsh;
-
-  x[0] -= step * y[0];
-  x[1] -= step * y[1];
-  *f = Function(x[0], x[1]);
-
-  return step;
-}
-
-double Step_optimal( double x[], double y[], double *f, double norm )
-{
-  double alpha = (sqrt(5) - 1) / 2;
-  double a = 0, b = 1;
-  double l, m, fl, fm;
-
-  l = a + (1 - alpha) * (b - a);
-  m = a + b - l;
-  fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
-  fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
- 
-  while (fabs(m - l) > .001)
-  {
-    if (fl < fm)
-    {
-      b = m;
-      m = l;
-      fm = fl;
-      l = a + (1 - alpha) * (b - a);
-      fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
-    }
-    else 
-    {
-      a = l;
-      l = m;
-      fl = fm;
-      m = a + alpha * (b - a);
-      fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
-    }
-  }
-
-  x[0] -= l * y[0];
-  x[1] -= l * y[1];
-  *f = Function(x[0], x[1]);
-  return l;
-}
-
-int period;
-
-int Test( double (*Step)(double[], double[], double *, double), int i )
-{
-  int    step = 0;
-  double x[2] = {-1, -1}, x_prev[2], delta[2], old_delta[2], y[2],
-         f = Function(x[0], x[1]), norm, mod, prev_mod;
-  static double solution[2];
-
-  do
-  {
-    Gradient(x, y);
-    if (i == 3)
-      Gradient2(x, y, 1);
-    if (i == 4)
-      Gradient2(x, y, (step % period) == 0);
-    norm = sqrt(y[0] * y[0] + y[1] * y[1]);
-    printf("x = (% 7.4lf, % 7.4lf) y = (% 7.4lf, % 7.4lf) norm = % 7.4lf,"
-           " f=%7.4lf ", x[0], x[1], y[0], y[1], norm, f);
-    memcpy(x_prev, x, 2 * sizeof(double));
-    printf("step = %6.4lf", Step(x, y, &f, norm));
-
-    if (i == 1)
-    {
-      delta[0] = x[0] - solution[0];
-      delta[1] = x[1] - solution[1];
-      mod = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
-      if (step != 0)
-        printf(", |x - x*| / |x_prev - x*| = %lf", mod / prev_mod);
-    }
-    if (i == 2)
-    {
-      delta[0] = x[0] - x_prev[0];
-      delta[1] = x[1] - x_prev[1];
-      printf(", delta = (%8.5lf, %8.5lf)", delta[0], delta[1]);
-      if (step != 0)
-        printf(", delta * old_delta = %.5lf", delta[0] * old_delta[0] +
-	                                      delta[1] * old_delta[1]);
-    }
-    printf("\n");
-      
-    step++;
-    old_delta[0] = delta[0];
-    old_delta[1] = delta[1];
-    prev_mod = mod;
-  } while (norm > eps);
-
-  if (i == 0)
-     memcpy(solution, x, 2 * sizeof(double));
-  printf("Done in %d steps\n", step);
-  return step;
-}
-
-double FindLipshitzConstant( double left, double right, double bottom, double top )
-{
-  double R, R_max = 0;
-  double t, x1, x2, y1, y2, dx, dy, df, denom;
-
-  while (1)
-  {
-    t = (double)(rand() % 100) / 100.0;
-    x1 = left + (1 - t) + right * t;
-    t = (double)(rand() % 100) / 100.0;
-    x2 = left + (1 - t) + right * t;
-    t = (double)(rand() % 100) / 100.0;
-    y1 = bottom + (1 - t) + top * t;
-    t = (double)(rand() % 100) / 100.0;
-    y2 = bottom + (1 - t) + top * t;
-
-    //printf("%lf, %lf, %lf, %lf\n", x1, x2, y1, y2);
-    dx = x2 - x1;
-    dy = y2 - y1;
-    df = fabs(Function(x2, y2) - Function(x1, y1));
-    
-    denom = sqrt(dx * dx + dy * dy);
-    if (denom < 1e-5)
-      continue;
-    R = df / denom;
-    if (R > R_max)
-      R_max = R;
-    printf("%lf\n", R_max);
-  }
-}
-
-int main( void) 
-{
-  int i;
-  //FindLipshitzConstant(-1, 0, -1, 0);
-  printf("\n2-dimensional minimization example\n");
-  printf("\nFixed step:\n");
-  Test(Step_fixed, 0);
-  printf("\nStep subdivision:\n");
-  Test(Step_subdivide, 1);
-  printf("\nOptimal step:\n");
-  Test(Step_optimal, 2);
-
-  printf("\nSecond order gradient:\n");
-  Test(Step_Newton, 3);
-  period = 2;
-  for (i = 0; i < 4; i++)
-  {
-    printf("\nSecond order gradient, with calculating invH once at %d steps:\n",
-	   period);
-    Test(Step_Newton, 4);
-    period *= 2;
-  }
-  //scanf("%*c");
-  return 0;
-}
+/*************************************************************
+ * FILENAME   : min_2d.c
+ * AUTHOR     : Zhidkov Evgueny, 3057/2
+ * PURPOSE    : 2D gradient methods
+ * LAST UPDATE: 14.03.2004
+ * NOTE       : none
+ *************************************************************/
+
+//////////////////////////////////////////////////////////////
+// LOCAL INCLUDES
+//////////////////////////////////////////////////////////////
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include <conio.h>
+
+//////////////////////////////////////////////////////////////
+// LOCAL DEFINES
+//////////////////////////////////////////////////////////////
+
+#define EPSILON 1.e-6
+
+//////////////////////////////////////////////////////////////
+// MATH FUNCTION DECLARATION
+//////////////////////////////////////////////////////////////
+
+double Function( double x1, double x2 )
+{
+  return exp(x1*x1 + 2 *x2*x2) + sqrt(1.0 + 2*x1 * x1 + x2 * x2);
+}
+
+void Gradient( double x[], double y[] )
+{
+  double sqrt_1 = sqrt(1 + x[0] * x[0] + x[1] * x[1]);
+
+  y[0] = 2 * x[0] * exp(x[0]*x[0] + 2 *x[1]*x[1])  + 2*x[0] / sqrt(1.0 + 2*x[0]*x[0] + x[1]*x[1]);
+  y[1] = 4 * x[1] * exp(x[0]*x[0] + 2 *x[1]*x[1])  + x[1] / sqrt(1.0 + 2*x[0]*x[0] + x[1]*x[1]);
+}
+
+void Hessian( double x[], double H[2][2] )
+{
+  double sqr_1 = 1 + x[0] * x[0] + x[1] * x[1];
+  double sqrt_1 = sqrt(sqr_1);
+
+  H[0][0] = 4 * (sqrt_1 - x[0] * x[0] / sqrt_1) / sqr_1;
+  H[0][1] = H[1][0] = -4 * x[0] * x[1] / (sqr_1 * sqrt_1);
+  H[1][1] = 4 * (sqrt_1 - x[1] * x[1] / sqrt_1) / sqr_1;
+}
+
+void Gradient2( double x[], double y[], int recalculate )
+{
+  static double H[2][2], invH[2][2], grad[2], det;
+
+  if (recalculate)
+  {
+    Hessian(x, H);
+    det = H[0][0] * H[1][1] - H[0][1] * H[1][0];
+    if (det < 1e-12)
+    {
+      y[0] = y[1] = 0;
+      return;
+    }
+    invH[0][0] =  H[1][1] / det;
+    invH[0][1] = -H[0][1] / det;
+    invH[1][0] = -H[1][0] / det;
+    invH[1][1] =  H[0][0] / det;
+  }
+  Gradient(x, grad);
+  y[0] = (invH[0][0] * grad[0] + invH[0][1] * grad[1]);
+  y[1] = (invH[1][0] * grad[0] + invH[1][1] * grad[1]);
+}
+  
+//////////////////////////////////////////////////////////////
+// OPTIMIZATIONS METHODS
+//////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////
+// STEP SUBDIVISION
+//////////////////////////////////////////////////////////////
+double Step_subdivide( double x[], double y[], double *f, double norm )
+{
+  double step = 2, fnew;
+
+  if (norm < EPSILON)
+    return 0;
+
+  do
+  {
+    step *= 0.5;
+    fnew = Function(x[0] - step * y[0], x[1] - step * y[1]);
+  } while (fnew - *f > -EPSILON * step * norm * norm); 
+
+  x[0] -= step * y[0];
+  x[1] -= step * y[1];
+  *f = fnew;
+  return step;
+}
+
+//////////////////////////////////////////////////////////////
+// STEP NEWTON
+//////////////////////////////////////////////////////////////
+double Step_Newton( double x[], double y[], double *f, double norm )
+{
+  x[0] -= y[0];
+  x[1] -= y[1];
+  *f = Function(x[0], x[1]);
+  return 1;
+}
+
+//////////////////////////////////////////////////////////////
+// STEP FIXED LENGTH
+//////////////////////////////////////////////////////////////
+double lipsh = 4.0;
+double m = 1e12, M = -1e12;
+double Step_fixed( double x[], double y[], double *f, double norm )
+{
+  //double step = (1 - EPSILON) / lipsh;
+  double step = 2.0 / (m + M);
+
+  x[0] -= step * y[0];
+  x[1] -= step * y[1];
+  *f = Function(x[0], x[1]);
+
+  return step;
+}
+
+//////////////////////////////////////////////////////////////
+// STEP OPTIMAL
+//////////////////////////////////////////////////////////////
+double Step_optimal( double x[], double y[], double *f, double norm )
+{
+  double alpha = (sqrt(5) - 1) / 2;
+  double a = 0, b = 1;
+  double l, m, fl, fm;
+
+  l = a + (1 - alpha) * (b - a);
+  m = a + b - l;
+  fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
+  fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
+ 
+  while (fabs(m - l) > .000001)
+  {
+    if (fl < fm)
+    {
+      b = m;
+      m = l;
+      fm = fl;
+      l = a + (1 - alpha) * (b - a);
+      fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
+    }
+    else 
+    {
+      a = l;
+      l = m;
+      fl = fm;
+      m = a + alpha * (b - a);
+      fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
+    }
+  }
+
+  x[0] -= l * y[0];
+  x[1] -= l * y[1];
+  *f = Function(x[0], x[1]);
+  return l;
+}
+
+//////////////////////////////////////////////////////////////
+// STEP FLETCHER-REEVES
+//////////////////////////////////////////////////////////////
+int    fp_loc_iter;
+double Step_FletcherReeves( double x[], double y[], double *f, double norm )
+{
+  double p_cur[2];
+  double a = 0, step = .5;
+
+  static double grad_prev[2] = {0.0, 0.0}, p_prev[2] = {0.0, 0.0};
+
+  if (fp_loc_iter % 2 == 0)
+  {
+    p_cur[0] = -y[0];
+    p_cur[1] = -y[1];
+  }
+  else
+  {
+    double betta;
+    betta = (y[0] * y[0] + y[1] * y[1]) / (grad_prev[0] * grad_prev[0] + grad_prev[1] * grad_prev[1]);
+    p_cur[0] = -y[0] + betta * p_prev[0];
+    p_cur[1] = -y[1] + betta * p_prev[1];
+  }
+  grad_prev[0] = y[0];
+  grad_prev[1] = y[1];
+  p_prev[0] = p_cur[0];
+  p_prev[1] = p_cur[1];
+  fp_loc_iter++;
+
+  while (step > 0.000001)
+  {
+    double f_prev, f_cur;
+    int    i = 2;
+
+    f_prev = Function(x[0] + (a + step) * p_cur[0], x[1] + (a + step) * p_cur[1]);
+    while (i < 100)
+    {
+      f_cur = Function(x[0] + (a + step * i) * p_cur[0], x[1] + (a + step * i) * p_cur[1]);
+      if (f_cur > f_prev) {
+        break;
+      }
+      f_prev = f_cur;
+      i++;
+    }
+    a += step * (i - 2);
+    step *= 0.1;
+  }
+  a += step;
+
+  x[0] += a * p_cur[0];
+  x[1] += a * p_cur[1];
+  *f = Function(x[0], x[1]);
+
+  return a;
+}
+
+//////////////////////////////////////////////////////////////
+// TEST METHODS ROUTINE
+//////////////////////////////////////////////////////////////
+int period;
+int Test( double (*Step)(double[], double[], double *, double), int i )
+{
+  int    step = 0;
+  double x[2] = {-1, -1}, x_prev[2], delta[2], old_delta[2], y[2],
+         f = Function(x[0], x[1]), norm, mod, prev_mod;
+
+  static double solution[2];
+
+  fp_loc_iter = 0;
+  do
+  {
+    Gradient(x, y);
+    norm = sqrt(y[0] * y[0] + y[1] * y[1]);
+    printf("x=(%06.4lf,%06.4lf) y=(%06.4lf, %06.4lf) norm=%08.6lf,"
+           " f=%07.5lf", x[0], x[1], y[0], y[1], norm, f);
+    if (norm <= EPSILON)
+      break;
+    if (i == 3)
+      Gradient2(x, y, 1);
+    if (i == 4)
+      Gradient2(x, y, (step % period) == 0);
+    memcpy(x_prev, x, 2 * sizeof(double));
+    printf(" step=%6.5lf\n", Step(x, y, &f, norm));
+
+    if (i == 1)
+    {
+      delta[0] = x[0] - solution[0];
+      delta[1] = x[1] - solution[1];
+      mod = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
+      if (step != 0)
+        printf("|x-x*|/|x_prev-x*| = %8.7lf\n", mod / prev_mod);
+    }
+    if (i == 2)
+    {
+      delta[0] = x[0] - x_prev[0];
+      delta[1] = x[1] - x_prev[1];
+      printf("delta=(%08.6lf, %08.6lf)", delta[0], delta[1]);
+      if (step != 0)
+        printf(", <delta,old_delta> = %.7lf", delta[0] * old_delta[0] +
+	                                            delta[1] * old_delta[1]);
+      printf("\n");
+    }
+      
+    step++;
+    old_delta[0] = delta[0];
+    old_delta[1] = delta[1];
+    prev_mod = mod;
+  } while (norm > EPSILON);
+
+  if (i == 0)
+     memcpy(solution, x, 2 * sizeof(double));
+  printf("\nDone in %d steps\n", step);
+  getch();
+  return step;
+}
+
+//////////////////////////////////////////////////////////////
+// FIND LIPSHITZ CONSTANT
+//////////////////////////////////////////////////////////////
+double FindLipshitzConstant( double left, double right, double bottom, double top )
+{
+  double R, R_max = 0, x1[2], x2[2], y1[2], y2[2];
+  double t, dx, dy, df, denom;
+
+  while (1)
+  {
+    t = (double)(rand() % 100) / 100.0;
+    x1[0] = left + (1 - t) + right * t;
+    t = (double)(rand() % 100) / 100.0;
+    x1[1] = left + (1 - t) + right * t;
+    t = (double)(rand() % 100) / 100.0;
+    x2[0] = bottom + (1 - t) + top * t;
+    t = (double)(rand() % 100) / 100.0;
+    x2[1] = bottom + (1 - t) + top * t;
+
+    dx = x2[0] - x1[0];
+    dy = x2[1] - x1[1];
+    Gradient(x1, y1);
+    Gradient(x2, y2);
+    df = sqrt((y2[0] - y1[0]) * (y2[0] - y1[0]) + (y2[1] - y1[1]) * (y2[1] - y1[1]));
+    
+    denom = sqrt(dx * dx + dy * dy);
+    if (denom < EPSILON)
+      continue;
+    R = df / denom;
+    if (R > R_max)
+    {
+      R_max = R;
+      printf("%lf\n", R_max);
+    }
+  }
+}
+
+//////////////////////////////////////////////////////////////
+// FIND EIGEN VALUES (m and M)
+//////////////////////////////////////////////////////////////
+void FindEigenVals( double left, double right, double bottom, double top )
+{
+  int    i, j, num = 100;
+  double x[2], H[2][2];
+  double p, q, descr, x1, x2;
+
+  for (i = 0; i <= num; i++)
+  {
+    x[0] = left + ((right - left) * i) / num;
+    for (j = 0; j <= num; j++)
+    {
+      x[1] = bottom + ((top - bottom) * j) / num;
+      Hessian(x, H);
+      p = -H[0][0] - H[1][1];
+      q = H[0][0] * H[1][1] - H[0][1] * H[1][0];
+      descr = p * p - 4 * q;
+      if (descr < 0)
+      {
+        printf("SHIT\n");
+        continue;
+      }
+      descr = sqrt(descr);
+      x1 = -0.5 * p - 0.5 * descr;
+      x2 = x1 + descr;
+      if (x1 < m)
+      {
+        m = x1;
+      }
+      if (x2 > M)
+      {
+        M = x2;
+      }
+    }
+  }
+  printf("m=%08.6lf; M=%08.6lf\n", m, M);
+  getch();
+}
+
+//////////////////////////////////////////////////////////////
+// MAIN FUNCTION
+//////////////////////////////////////////////////////////////
+int main( void) 
+{
+  int i;
+//  FindLipshitzConstant(-1, 0, -1, 0);
+  FindEigenVals(-1, 0, -1, 0);
+  printf("\n2-dimensional minimization example\n");
+  printf("\nFixed step:\n");
+  getch();
+  Test(Step_fixed, 0);
+//  printf("\nStep subdivision:\n");
+//  getch();
+//  Test(Step_subdivide, 1);
+  printf("\nOptimal step:\n");
+  getch();
+  Test(Step_optimal, 2);
+  printf("\nFletcher-Reeves step:\n");
+  getch();
+  Test(Step_FletcherReeves, 5);
+
+  printf("\nSecond order gradient:\n");
+  getch();
+  Test(Step_Newton, 3);
+  period = 4;
+  for (i = 0; i < 2; i++)
+  {
+    printf("\nSecond order gradient, with calculating invH once at %d steps:\n", period);
+     getch();
+    Test(Step_Newton, 4);
+    period *= 2;
+  }
+  return 0;
+}

File labs/2d_min/min_2d.cpp

View file
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+double Function( double x1, double x2 )
+{
+  double x1_2 = x1 * x1, x2_2 = x2 * x2;
+
+  return x1_2 * x1_2 + 2 * x2_2 * x2_2 + x1_2 * x2_2 + 2 * x1 + x2;
+}
+
+void Gradient( double x[], double y[] )
+{
+  y[0] = 4 * x[0] * x[0] * x[0] + 2 * x[0] * x[1] * x[1] + 2;
+  y[1] = 8 * x[1] * x[1] * x[1] + 2 * x[1] * x[0] * x[0] + 1;
+}
+
+void Jacobian( double x[], double J[2][2] )
+{
+  J[0][0] = 12 * x[0] * x[0] + 2 * x[1] * x[1];
+  J[0][1] = J[1][0] = 4 * x[0] * x[1];
+  J[1][1] = 24 * x[1] * x[1] + 2 * x[0] * x[0];
+}
+
+void Gradient2( double x[], double y[], int recalculate )
+{
+  static double H[2][2], invH[2][2], grad[2], det;
+
+  if (recalculate)
+  {
+    Jacobian(x, H);
+    det = H[0][0] * H[1][1] - H[0][1] * H[1][0];
+    if (det < 1e-5)
+    {
+      y[0] = y[1] = 0;
+      return;
+    }
+    invH[0][0] =  H[1][1] / det;
+    invH[0][1] = -H[0][1] / det;
+    invH[1][0] = -H[1][0] / det;
+    invH[1][1] =  H[0][0] / det;
+  }
+  Gradient(x, grad);
+  y[0] = (invH[0][0] * grad[0] + invH[0][1] * grad[1]);
+  y[1] = (invH[1][0] * grad[0] + invH[1][1] * grad[1]);
+}
+  
+double eps = 1e-6;
+
+double Step_subdivide( double x[], double y[], double *f, double norm )
+{
+  double step = 1, fnew, epsilon = eps;
+
+  if (norm < eps)
+    return 0;
+
+  do
+  {
+    step *= 0.5;
+    fnew = Function(x[0] - step * y[0], x[1] - step * y[1]);
+  } while (fnew - *f > -epsilon * step * norm); 
+
+  x[0] -= step * y[0];
+  x[1] -= step * y[1];
+  *f = fnew;
+  return step;
+}
+
+double Step_Newton( double x[], double y[], double *f, double norm )
+{
+  x[0] -= y[0];
+  x[1] -= y[1];
+  *f = Function(x[0], x[1]);
+  return 1;
+}
+
+double lipsh = 9.2;
+
+double Step_fixed( double x[], double y[], double *f, double norm )
+{
+  double step = (1 - eps) / lipsh;
+
+  x[0] -= step * y[0];
+  x[1] -= step * y[1];
+  *f = Function(x[0], x[1]);
+
+  return step;
+}
+
+double Step_optimal( double x[], double y[], double *f, double norm )
+{
+  double alpha = (sqrt(5) - 1) / 2;
+  double a = 0, b = 1;
+  double l, m, fl, fm;
+
+  l = a + (1 - alpha) * (b - a);
+  m = a + b - l;
+  fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
+  fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
+ 
+  while (fabs(m - l) > .001)
+  {
+    if (fl < fm)
+    {
+      b = m;
+      m = l;
+      fm = fl;
+      l = a + (1 - alpha) * (b - a);
+      fl = Function(x[0] - y[0] * l, x[1] - y[1] * l);
+    }
+    else 
+    {
+      a = l;
+      l = m;
+      fl = fm;
+      m = a + alpha * (b - a);
+      fm = Function(x[0] - y[0] * m, x[1] - y[1] * m);
+    }
+  }
+
+  x[0] -= l * y[0];
+  x[1] -= l * y[1];
+  *f = Function(x[0], x[1]);
+  return l;
+}
+
+int period;
+
+int Test( double (*Step)(double[], double[], double *, double), int i )
+{
+  int    step = 0;
+  double x[2] = {-1, -1}, x_prev[2], delta[2], old_delta[2], y[2],
+         f = Function(x[0], x[1]), norm, mod, prev_mod;
+  static double solution[2];
+
+  do
+  {
+    Gradient(x, y);
+    if (i == 3)
+      Gradient2(x, y, 1);
+    if (i == 4)
+      Gradient2(x, y, (step % period) == 0);
+    norm = sqrt(y[0] * y[0] + y[1] * y[1]);
+    printf("x = (% 7.4lf, % 7.4lf) y = (% 7.4lf, % 7.4lf) norm = % 7.4lf,"
+           " f=%7.4lf ", x[0], x[1], y[0], y[1], norm, f);
+    memcpy(x_prev, x, 2 * sizeof(double));
+    printf("step = %6.4lf", Step(x, y, &f, norm));
+
+    if (i == 1)
+    {
+      delta[0] = x[0] - solution[0];
+      delta[1] = x[1] - solution[1];
+      mod = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
+      if (step != 0)
+        printf(", |x - x*| / |x_prev - x*| = %lf", mod / prev_mod);
+    }
+    if (i == 2)
+    {
+      delta[0] = x[0] - x_prev[0];
+      delta[1] = x[1] - x_prev[1];
+      printf(", delta = (%8.5lf, %8.5lf)", delta[0], delta[1]);
+      if (step != 0)
+        printf(", delta * old_delta = %.5lf", delta[0] * old_delta[0] +
+	                                      delta[1] * old_delta[1]);
+    }
+    printf("\n");
+      
+    step++;
+    old_delta[0] = delta[0];
+    old_delta[1] = delta[1];
+    prev_mod = mod;
+  } while (norm > eps);
+
+  if (i == 0)
+     memcpy(solution, x, 2 * sizeof(double));
+  printf("Done in %d steps\n", step);
+  return step;
+}
+
+double FindLipshitzConstant( double left, double right, double bottom, double top )
+{
+  double R, R_max = 0;
+  double t, x1, x2, y1, y2, dx, dy, df, denom;
+
+  while (1)
+  {
+    t = (double)(rand() % 100) / 100.0;
+    x1 = left + (1 - t) + right * t;
+    t = (double)(rand() % 100) / 100.0;
+    x2 = left + (1 - t) + right * t;
+    t = (double)(rand() % 100) / 100.0;
+    y1 = bottom + (1 - t) + top * t;
+    t = (double)(rand() % 100) / 100.0;
+    y2 = bottom + (1 - t) + top * t;
+
+    //printf("%lf, %lf, %lf, %lf\n", x1, x2, y1, y2);
+    dx = x2 - x1;
+    dy = y2 - y1;
+    df = fabs(Function(x2, y2) - Function(x1, y1));
+    
+    denom = sqrt(dx * dx + dy * dy);
+    if (denom < 1e-5)
+      continue;
+    R = df / denom;
+    if (R > R_max)
+      R_max = R;
+    printf("%lf\n", R_max);
+  }
+}
+
+int main( void) 
+{
+  int i;
+  //FindLipshitzConstant(-1, 0, -1, 0);
+  printf("\n2-dimensional minimization example\n");
+  printf("\nFixed step:\n");
+  Test(Step_fixed, 0);
+  printf("\nStep subdivision:\n");
+  Test(Step_subdivide, 1);
+  printf("\nOptimal step:\n");
+  Test(Step_optimal, 2);
+
+  printf("\nSecond order gradient:\n");
+  Test(Step_Newton, 3);
+  period = 2;
+  for (i = 0; i < 4; i++)
+  {
+    printf("\nSecond order gradient, with calculating invH once at %d steps:\n",
+	   period);
+    Test(Step_Newton, 4);
+    period *= 2;
+  }
+  //scanf("%*c");
+  return 0;
+}