Tor Olav Kristensen wrote:

> Therefore I'm working on macros to implement the Nelder-
> Mead simplex method to locate the minima of the "error field".

And here they are.

I'll try to explain later (if its necessary).

Tor Olav

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Copyright 2000 by Tor Olav Kristensen
// mailto:tor### [at] hotmailcom
// http://www.crosswinds.net/~tok/tokrays.html
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#version 3.1;

global_settings {
  #max_trace_level 1
  ambient_light 2
}

#include "colors.inc"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// The spline points to put a circle around

#declare SplinePts =
array[6] {
  <1.0, 0, 1.9>,
  <2.0, 0, 1.1>,
  <1.7, 0, 1.5>,
  <2.9, 0, 1.7>,
  <3.1, 0, 3.0>,
  <3.1, 0, 3.4>
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Just to show the spline points

#macro PrintVector(Vector)

  #debug "< "
  #debug str(Vector.x, 0, -1)
  #debug ", "
  #debug str(Vector.y, 0, -1)
  #debug ", "
  #debug str(Vector.z, 0, -1)
  #debug " >"

#end // macro PrintVector

#macro PtPlot(PointArray, SphereRadius)

  #local NrOfPoints = dimension_size(PointArray, 1);

  union {
    #local PtCnt = 0;
    #while(PtCnt < NrOfPoints)
      sphere { PointArray[PtCnt], SphereRadius }
      #local PtCnt = PtCnt + 1;
    #end // while
  }

#end // macro PtPlot


object {
  PtPlot(SplinePts, 0.06)
  pigment { color Cyan*2 }
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Functions used to evalute fitness of a point

#macro SumSqdErrors(Point, PointArray, Radius)

  #local NrOfPoints = dimension_size(PointArray, 1);
  #local ErrorSum = 0;

  #local PtCnt = 0;
  #while(PtCnt < NrOfPoints)
    #local ErrorSum = ErrorSum +
           pow(vlength(Point - PointArray[PtCnt]) - Radius, 2);
    #local PtCnt = PtCnt + 1;
  #end // while

  ErrorSum

#end // macro SumSqdErrors


#macro MeanDist(Point, PointArray)

  #local NrOfPoints = dimension_size(PointArray, 1);
  #local DistSum = 0;

  #local PtCnt = 0;
  #while(PtCnt < NrOfPoints)
    #local ThisDist = vlength(Point - PointArray[PtCnt]);
    #local DistSum = DistSum + ThisDist;
    #local PtCnt = PtCnt + 1;
  #end // while

  (DistSum/NrOfPoints)

#end // macro MeanDist


#macro UserFunction(Point)

  #local BestRadius = MeanDist(Point, SplinePts);
  #local MinimumError = SumSqdErrors(Point, SplinePts, BestRadius);

  MinimumError

#end // macro UserFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Nelder-Mead simplex method macros
// Made to work for several dimensions
// (as many as POV-Ray can handle in standard vectors: 4 or 5?)
// but I have only tested them for 2D points/vectors

#macro RankVertices(SimplexArray,
                    BestIndex, BestValue,
                    NextWorstIndex, NextWorstValue,
                    WorstIndex, WorstValue)

  #local NrOfVertices = dimension_size(SimplexArray, 1);
  #declare BestIndex      = 0;
  #declare WorstIndex     = 0;
  #declare NextWorstIndex = 0;
  #local FValue = UserFunction(SimplexArray[0]);
  #declare BestValue      = FValue;
  #declare WorstValue     = FValue;
  #declare NextWorstValue = FValue;

  #local IndexCnt = 1;
  #while (IndexCnt < NrOfVertices)
    #local FValue = UserFunction(SimplexArray[IndexCnt]);
    #if (FValue <= BestValue)
      #declare BestIndex = IndexCnt;
      #declare BestValue = FValue;
    #else
      #if (FValue >= WorstValue)
        #declare WorstIndex = IndexCnt;
        #declare WorstValue = FValue;
      #else
        #if (FValue >= NextWorstValue)
          #declare NextWorstIndex = IndexCnt;
          #declare NextWorstValue = FValue;
        #end // if
      #end // if
    #end // if
    #local IndexCnt = IndexCnt + 1;
  #end // while

#end // macro RankVertices


#macro FindCentroid(SimplexArray, WorstVertex)

  #local Dimensions = dimension_size(SimplexArray, 1) - 1;
  #local Centroid = -WorstVertex;

  #local IndexCnt = 0;
  #while (IndexCnt <= Dimensions)
    #local Centroid = Centroid + SimplexArray[IndexCnt];
    #local IndexCnt = IndexCnt + 1;
  #end // while

  (Centroid/Dimensions)

#end // macro FindCentroid


#macro ProjVertex(Vertex, Centroid, ScaleCoeff)

  #local vDisplace = Centroid - Vertex;

  (Centroid + ScaleCoeff*vDisplace)

#end // macro ProjVertex


#macro ShrinkSimplex(SimplexArray, Vertex)

  #local NrOfVertices = dimension_size(SimplexArray, 1);

  #local IndexCnt = 0;
  #while (IndexCnt < NrOfVertices)
    #declare SimplexArray[IndexCnt] =
             (SimplexArray[IndexCnt] + Vertex)/2;
    #local IndexCnt = IndexCnt + 1;
  #end // while

#end // macro ShrinkSimplex


#macro MoveSimplex(SimplexArray)

  #local BestIndex      = 0;
  #local BestValue      = 0;
  #local WorstIndex     = 0;
  #local WorstValue     = 0;
  #local NextWorstIndex = 0;
  #local NextWorstValue = 0;

  #local Crefl  =  1.0; // Reflection  coefficient
  #local Cexp   =  2.0; // Expansion   coefficient
  #local Ccontr = -0.5; // Contraction coefficient

  RankVertices(SimplexArray,
               BestIndex, BestValue,
               NextWorstIndex, NextWorstValue,
               WorstIndex, WorstValue)

  #local WorstVertex = SimplexArray[WorstIndex];
  #local Centroid = FindCentroid(SimplexArray, WorstVertex);
  #local NextIndex = WorstIndex;
  #local ReflPoint = ProjVertex(WorstVertex, Centroid, Crefl);
  #local ReflValue = UserFunction(ReflPoint);
  #if (ReflValue < BestValue)
    #if (Debug) #debug "Trying Expansion " #end // if
    #local ExpPoint = ProjVertex(WorstVertex, Centroid, Cexp);
    #local ExpValue = UserFunction(ExpPoint);
    #if (ExpValue < BestValue)
      #declare SimplexArray[NextIndex] = ExpPoint;
      #if (Debug) #debug "-> Expansion" #end // if
    #else
      #declare SimplexArray[NextIndex] = ReflPoint;
      #if (Debug) #debug "-> Reflection" #end // if
    #end
  #else
    #if (ReflValue < NextWorstValue)
      #declare SimplexArray[NextIndex] = ReflPoint;
      #if (Debug) #debug "Reflection" #end // if
    #else
      #if (Debug) #debug "Trying Contraction " #end // if
      #local ContrPoint = ProjVertex(WorstVertex, Centroid, Ccontr);
      #local ContrValue = UserFunction(ContrPoint);
      #if (ContrValue < WorstValue)
        #declare SimplexArray[NextIndex] = ContrPoint;
        #if (Debug) #debug "-> Contraction" #end // if
      #else
        #local BestVertex = SimplexArray[BestIndex];
        ShrinkSimplex(SimplexArray, BestVertex)
        #declare NextIndex = BestIndex;
        #if (Debug) #debug "-> Shrinking" #end // if
      #end // if
    #end // if
  #end // if

  NextIndex

#end // macro MoveSimplex


#macro FindMinima(SimplexArray, Epsilon, MaxIterations, Point, F_Value)

  #local LastFvalue = -1;
  #local Converged = false;

  #local IterCnt = 0;
  #while (IterCnt < MaxIterations & !Converged)
    #if (Debug)
      #debug concat("\nIteration: ", str(IterCnt, 0, -1), " ")
    #end // if
    #local PointIndex = MoveSimplex(SimplexArray);
    #declare Point = SimplexArray[PointIndex];
    #declare F_value = UserFunction(Point);
    #if (Debug)
      #debug "\n                                                        "
    #end // if
    #if (Debug) PrintVector(Point + y*F_value) #end // if
    #local Diff = LastFvalue - F_value;
    #if ((Diff >= 0) & (Diff < Epsilon))
      #if (Debug) #debug "\nConverged!\n" #end // if
      #local Converged = true;
    #end // if
    #local LastFvalue = F_value;
    #local IterCnt = IterCnt + 1;
  #end // while
  #if (Debug) #debug "\n" #end // if

  Converged

#end // macro FindMinima

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Use of Nelder-Mead simplex method macros in 2D

#macro PointsCombine2D(PointArray, Epsilon, MaxIterations)

  #local NrOfPts = dimension_size(PointArray, 1);
  #local BestPoint = PointArray[0];
  #local LeastValue = UserFunction(PointArray[0]);

  #local Dimensions = 2;
  #local SimplexArray = array[Dimensions + 1]
  #local PtCnt1 = 0;
  #while (PtCnt1 < NrOfPts)
    #local SimplexArray[0] = PointArray[PtCnt1];
    #local PtCnt2 = PtCnt1 + 1;
    #while (PtCnt2 < NrOfPts)
    #local SimplexArray[1] = PointArray[PtCnt2];
      #local PtCnt3 = PtCnt2 + 1;
      #while (PtCnt3 < NrOfPts)
        #local SimplexArray[2] = PointArray[PtCnt3];
        #local GoodPoint = BestPoint; // Just to init it
        #local SmallValue = 0; // Just to init it
        #if (FindMinima(SimplexArray, Epsilon, MaxIterations,
                        GoodPoint, SmallValue))
          #local SmallValue = UserFunction(GoodPoint);
          // The above line should not be necessary (bug?)
        #end // if
          #if (Debug)
            #debug "                           "
            #debug "                           "
            PrintVector(GoodPoint)
            #debug "\n"
          #end // if
          #if (SmallValue < LeastValue)
            #local BestPoint = GoodPoint;
            #local LeastValue = SmallValue;
            #if (Debug)
              #debug "> ==== > ==== > ==== > ===="
              #debug "> ==== > ==== > ==== >     "
              #debug "Found a better point!\n"
            #end // if
          #end // if
        #local PtCnt3 = PtCnt3 + 1;
      #end // while
      #local PtCnt2 = PtCnt2 + 1;
    #end // while
    #local PtCnt1 = PtCnt1 + 1;
  #end // while
  #if (Debug) #debug "\n" #end // if

  BestPoint

#end // macro PointsCombine2D

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#declare Debug = false;

#declare PerfectPoint = PointsCombine2D(SplinePts, 0.0001, 50);
sphere { PerfectPoint, 0.2 pigment { color Red*2 } }

camera {
  location PerfectPoint + 4*y
  look_at PerfectPoint
}

background { color Blue/2 }

light_source { <10,  100, 10>, 1 }
light_source { <10, -100, 10>, 1 }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

Post a reply to this message