algorithm - Nelder-Mead optimization in C++ - Code Review Stack Exchange

PHOTO

Mon Jan 30 2023 15:22:55 GMT+0000 (Coordinated Universal Time)

Saved by @yoonchoi #cpp

#include "Eigen/Dense"
#include <vector>
#include <limits>
#include <cstring>
#include <iostream>

using namespace Eigen;

namespace Optimization
{
    template<int nDims>
    class NelderMead {
    public:
        typedef Array<double, 1, nDims> FunctionParameters;
        typedef std::function<double(FunctionParameters)> ErrorFunction;

        NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
                   double minError, double initialEdgeLength,
                   double shrinkCoeff = 1, double contractionCoeff = 0.5,
                   double reflectionCoeff = 1, double expansionCoeff = 1)
                : errorFunction(errorFunction), minError(minError),
                shrinkCoeff(shrinkCoeff), contractionCoeff(contractionCoeff),
                reflectionCoeff(reflectionCoeff), expansionCoeff(expansionCoeff),
                worstValueId(-1), secondWorstValueId(-1), bestValueId(-1)
        {
            this->errors = std::vector(nDims + 1, std::numeric_limits<double>::max());

            const double b = initialEdgeLength / (nDims * SQRT2) * (sqrt(nDims + 1) - 1);
            const double a = initialEdgeLength / SQRT2;

            this->values = initial.replicate(nDims + 1, 1);

            for (int i = 0; i < nDims; i++)
            {
                FunctionParameters simplexRow;
                simplexRow.setConstant(b);
                simplexRow(0, i) = a;
                simplexRow += initial;

                this->values.row(i+1) = simplexRow;
            }
        }

        void optimize()
        {
            for (int i = 0; i < nDims+1; i++)
            {
                this->errors.at(i) = this->errorFunction(this->values.row(i));
            }

            this->invalidateIdsCache();

            while (this->errors.at(this->bestValueId) > this->minError)
            {
                step();
                auto bestError = this->errorFunction(this->best());
                auto worstError = this->errorFunction(this->worst());

                std::cout << "Best error " << std::to_string(bestError) << " with : " << this->best();
                std::cout << " Worst error " << std::to_string(worstError) << " with : " << this->worst();
                std::cout << '\n';
            }
        }

        void step()
        {
            auto meanWithoutWorst = this->getMeanWithoutWorst();

            auto reflectionOfWorst = this->getReflectionOfWorst(meanWithoutWorst);
            auto reflectionError = this->errorFunction(reflectionOfWorst);

            FunctionParameters newValue = reflectionOfWorst;
            double newError = reflectionError;

            bool shrink = false;

            if (reflectionError < this->errors.at(this->bestValueId))
            {
                auto expansionValue = this->expansion(meanWithoutWorst, reflectionOfWorst);
                double expansionError = this->errorFunction(expansionValue);

                if (expansionError < this->errors.at(this->bestValueId))
                {
                    newValue = expansionValue;
                    newError = expansionError;
                }
            }
            else if (reflectionError > this->errors.at(this->worstValueId))
            {
                newValue = this->insideContraction(meanWithoutWorst);
                newError = this->errorFunction(newValue);

                if (newError > this->errors.at(this->worstValueId)) { shrink = true; }
            }
            else if (reflectionError > this->errors.at(this->secondWorstValueId))
            {
                newValue = this->outsideContraction(meanWithoutWorst);
                newError = this->errorFunction(newValue);

                if (newError > reflectionError) { shrink = true; }
            }
            else
            {
                newValue = reflectionOfWorst;
                newError = reflectionError;
            }

            if (shrink)
            {
                this->shrink();
                this->invalidateIdsCache();
                return;
            }

            this->values.row(this->worstValueId) = newValue;
            this->errors.at(this->worstValueId) = newError;
            this->invalidateIdsCache();
        }

        inline FunctionParameters worst() { return this->values.row(this->worstValueId); }
        inline FunctionParameters best() { return this->values.row(this->bestValueId); }

    private:

        void shrink()
        {
            auto bestVertex = this->values.row(this->bestValueId);

            for (int i = 0; i < nDims + 1; i++)
            {
                if (i == this->bestValueId) { continue; }

                this->values.row(i) = bestVertex + this->shrinkCoeff * (this->values.row(i) - bestVertex);
                this->errors.at(i) = this->errorFunction(this->values.row(i));
            }
        }

        inline FunctionParameters expansion(FunctionParameters meanWithoutWorst, FunctionParameters reflection)
        {
            return reflection + this->expansionCoeff * (reflection - meanWithoutWorst);
        }

        inline FunctionParameters insideContraction(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst - this->contractionCoeff * (meanWithoutWorst - this->worst());
        }

        inline FunctionParameters outsideContraction(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst + this->contractionCoeff * (meanWithoutWorst - this->worst());
        }

        FunctionParameters getReflectionOfWorst(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst + this->reflectionCoeff * (meanWithoutWorst - this->worst());
        }

        FunctionParameters getMeanWithoutWorst()
        {
            FunctionParameters mean(0);
            for (int i = 0; i < nDims + 1; i++)
            {
                if (i == this->worstValueId) { continue; }

                mean += this->values.row(i);
            }

            // Not divided by nDims+1 because there's one ignored value.
            mean /= nDims;

            return mean;
        }

        void invalidateIdsCache() {
            double worstError = std::numeric_limits<double>::min();
            int worstId = -1;
            double secondWorstError = std::numeric_limits<double>::max();
            int secondWorstId = -1;
            double bestError = std::numeric_limits<double>::max();
            int bestId = -1;

            for (int i = 0; i < nDims + 1; i++)
            {
                auto error = this->errors.at(i);

                if (error > worstError)
                {
                    secondWorstError = worstError;
                    secondWorstId = worstId;
                    worstError = error;
                    worstId = i;
                }
                else if (error > secondWorstError)
                {
                    secondWorstError = error;
                    secondWorstId = i;
                }

                if (error < bestError)
                {
                    bestError = error;
                    bestId = i;
                }
            }

            // If we deal with a problem in 1D, it won't be set.
            if (secondWorstId == -1)
            {
                secondWorstId = worstId;
            }

            this->bestValueId = bestId;
            this->worstValueId = worstId;
            this->secondWorstValueId = secondWorstId;
        }

        ErrorFunction errorFunction;
        Array<double, nDims + 1, nDims> values;
        std::vector<double> errors;

        int worstValueId;
        int secondWorstValueId;
        int bestValueId;

        double minError;

        double shrinkCoeff;
        double expansionCoeff;
        double contractionCoeff;
        double reflectionCoeff;

        const double SQRT2 = sqrt(2);
    };
}

https://codereview.stackexchange.com/questions/272951/nelder-mead-optimization-in-c

Save snippets that work from anywhere online with our extensions

Available in the Chrome Web Store

Get Firefox Add-on

Get VS Code extension

Comments

More like this

int intvar = char - '0';

What is the meaning of a C++ Wrapper Class? - Stack Overflow

class int_ptr_wrapper
{
public:
    int_ptr_wrapper(int value = 0) :
    mInt(new int(value))
    {}

    // note! needs copy-constructor and copy-assignment operator!

    ~int_ptr_wrapper()
    {
        delete mInt;
    }

private:
    int* mInt;
};

Read file line by line using ifstream in C++ - Stack Overflow

#include <fstream>

std::ifstream file(FILENAME);
if (file.is_open()) {
    std::string line;
    while (std::getline(file, line)) {
        // using printf() in all tests for consistency
        printf("%s", line.c_str());
    }
    file.close();
}

Eigen Tensor C++, multiply each row by its transpose

#include <unsupported/Eigen/CXX11/Tensor>
int main(){
    // Define a tensor
    Eigen::Tensor<double,3> A(55,50,2);
    A.setRandom();

    // Define the pairs of indices to multiply (i.e. index 1 on A twice)
    std::array<Eigen::IndexPair<long>, 1> idx = { Eigen::IndexPair<long> {1, 1} };
    
    // Define a parallel device 
    int num_threads = 4;
    Eigen::ThreadPool       tpl (num_threads);
    Eigen::ThreadPoolDevice dev (&tpl, num_threads);

    // Initialize result buffer B. B will have dimensions [55,2,55,2]
    Eigen::Tensor<double,4> B(55,2,55,2);

    // Contract using the parallel device
    B.device(dev) = A.contract(A, idx);
}

c++ - Passing a shared pointer by reference or by value as parameter to a class - Stack Overflow

class Foo {
public:
    Foo(std::shared_ptr<Boo> boo) 
        : m_Boo(std::move(boo)) {}
private:
    std::shared_ptr<Boo> m_Boo;
};

c++ - using stringstream to print a rounded floating point number - Stack Overflow

#include <sstream>
#include <iostream>
#include <iomanip>

int main()
{
  double d = 12.3456789;
  std::stringstream ss;

  ss << std::fixed << std::setprecision( 4 ) << d;

  std::cout << ss.str() << std::endl;
}

* versus ^ allocation opperators

// pointer to new std::string object -> memory is not garbage-collected
std::string* strPtr = new std::string;

// pointer to System::String object -> memory is garbage-collected
System::String^ manStr = gcnew System::String;

Hashtable^ tempHash = gcnew Hashtable(iterators_);

IDictionaryEnumerator^ enumerator = tempHash->GetEnumerator();

.NET pointer comparison

// here normal pointer
P* ptr = new P; // usual pointer allocated on heap
P& nat = *ptr; // object on heap bind to native object

//.. here CLI managed 
MO^ mngd = gcnew MO; // allocate on CLI heap
MO% rr = *mngd; // object on CLI heap reference to gc-lvalue

Tumpukan / Stack (C++) | Rahmat Subekti

#include <iostream>
#define MAX 10
using namespace std;

//Deklarasi struct tumpukan
struct Stack {
	int top, data[MAX];
}Tumpukan;

void init(){
	Tumpukan.top = -1;
}

bool isEmpty() {
  return Tumpukan.top == -1;
}

bool isFull() {
	return Tumpukan.top == MAX-1;
}
void push() {
   if (isFull()) {
		cout << "\nTumpukan penuh"<<endl;
	}
	else {
    Tumpukan.top++;
    cout << "\nMasukkan data = "; cin >> Tumpukan.data[Tumpukan.top];
    cout << "Data " << Tumpukan.data[Tumpukan.top] << " masuk ke stack"<<endl;
	}
}
void pop() {
	if (isEmpty()) {
		cout << "\nData kosong\n"<<endl;
	}
	else {
    cout << "\nData "<<Tumpukan.data[Tumpukan.top]<<" sudah terambil"<<endl;
    Tumpukan.top--;
	}
}
void printStack() {
	if (isEmpty()) {
		cout << "Tumpukan kosong";
	}
	else {
    cout << "\nTumpukan : ";
		for (int i = Tumpukan.top; i >= 0; i--)
			cout << Tumpukan.data[i] << ((i == 0) ? "" : ",");
	}
}
int main() {
	int pilihan;
	init();
	do {
    printStack();
		cout << "\n1. Input (Push)\n"
        <<"2. Hapus (Pop)\n"
        <<"3. Keluar\n"
        <<"Masukkan Pilihan: ";
		cin >> pilihan;
		switch (pilihan)
		{
		case 1:
			push();
			break;
		case 2:
			pop();
			break;
		default:
      cout << "Pilihan tidak tersedia" << endl;
			break;
		}
	} while (pilihan!=3);
}

#include<iostream>
using namespace std;
#define n 100
//stack class
class stack{
    public:
        int *arr;
        int top;
        //constructor for allocating memory
        stack()
        {
            arr = new int[n];
            top = -1;
        }

        //push function to insert an element at the top of stack
        void push(int data)
        {
            if(this->top==n-1)
            {
                cout << "no space for the new element " << endl;
                return;
            }

            ++top;
            arr[top] = data;
            

            return;
        }
        //Function for removing the element from the top of the stack 
        void pop()
        {
            if(top==-1)
            {
                cout << "Empty Stack " << endl;
                return;
            }
            top--;
        }

        //Function to get the top element value

        int Top()
        {
            if(top==-1)
            {
                cout << "Empty Stack " << endl;
                return -1;
            }
            if(top>n-1)
            {
                cout << "stack overflow " << endl;
                return -1;
            }
            return arr[top];
        }

        //Function to check is stack is empty or not
        bool empty()
        {
            return top == -1;
        }
//Function for printing the stack
void print()
{
    if(top==-1)
    {
        cout << "Empty Stack " << endl;
        return;
    }
    for (int i = 0; i <= top;i++)
    {
        cout << arr[i] << "  ";
    }
    return;
}

~stack()
{
    delete arr;
}

};


//main function
int main()
{
    //making object of stack class
    stack object;
    object.push(1);
    object.push(2);
    object.push(3);
    object.push(4);
    cout << "1st print " << endl;
    object.print();
    object.pop();
    cout << endl;
    cout << "2nd print" << endl;
    object.print();
    cout << endl;
   cout<< object.Top();
   cout << endl;
   //checking if the stack is Empty
   cout<<object.empty();
   cout << endl;

   return 0;
}

c++ - What is definition of reference type? - Stack Overflow

{ // Block scope
     Foo fooVal = makeFoo(); // Say makeFoo() returns a (temporary, unnamed) Foo
     // Here the temporary Foo is dead (fooVal is a copy).

     // Foo &fooRef = makeFoo(); // Error, reference is non-const
     Foo const &fooCRef = makeFoo(); // All good

     // ...

     // The second temporary is still alive
     fooCRef.doSomethingFunny(); // Works like a charm !

} // The second temporary dies with fooRef

//------------------------------------------

Foo *fooPtr = new Foo; // Here is a Foo
Foo &fooRef = *fooPtr; // Here is an alias for that Foo

delete fooPtr; // Here is the end of that Foo's life

fooRef.doSomethingFunny(); // Here comes trouble...

How do I convert a double into a string in C++? - Stack Overflow

std::ostringstream strs;
strs << dbl;
std::string str = strs.str();

Gauss Jordan Method Using C++ with Output

/*
Program: Gauss Jordan Method
All array indexes are assumed to start from 1
*/

#include<iostream>
#include<iomanip>
#include<math.h>
#include<stdlib.h>

#define   SIZE   10

using namespace std;

int main()
{
	 float a[SIZE][SIZE], x[SIZE], ratio;
	 int i,j,k,n;

     /* Setting precision and writing floating point values in fixed-point notation. */
     cout<< setprecision(3)<< fixed;

	 /* Inputs */
	 /* 1. Reading number of unknowns */
	 cout<<"Enter number of unknowns: ";
	 cin>>n;

	 /* 2. Reading Augmented Matrix */
	 cout<<"Enter Coefficients of Augmented Matrix: "<< endl;
	 for(i=1;i<=n;i++)
	 {
		  for(j=1;j<=n+1;j++)
		  {
			   cout<<"a["<< i<<"]"<< j<<"]= ";
			   cin>>a[i][j];
		  }
	 }
    /* Applying Gauss Jordan Elimination */
     for(i=1;i<=n;i++)
     {
          if(a[i][i] == 0.0)
          {
               cout<<"Mathematical Error!";
               exit(0);
          }
          for(j=1;j<=n;j++)
          {
               if(i!=j)
               {
                    ratio = a[j][i]/a[i][i];
                    for(k=1;k<=n+1;k++)
                    {
                        a[j][k] = a[j][k] - ratio*a[i][k];
                    }
               }
          }
     }
     /* Obtaining Solution */
     for(i=1;i<=n;i++)
     {
        x[i] = a[i][n+1]/a[i][i];
     }

	 /* Displaying Solution */
	 cout<< endl<<"Solution: "<< endl;
	 for(i=1;i<=n;i++)
	 {
	  	cout<<"x["<< i<<"] = "<< x[i]<< endl;
	 }

	 return(0);
}

parsing - Parse (split) a string in C++ using string delimiter (standard C++) - Stack Overflow

std::string s = "scott>=tiger>=mushroom";
std::string delimiter = ">=";

size_t pos = 0;
std::string token;
while ((pos = s.find(delimiter)) != std::string::npos) {
    token = s.substr(0, pos);
    std::cout << token << std::endl;
    s.erase(0, pos + delimiter.length());
}
std::cout << s << std::endl;

Browse more snippets >>