Page 1 of 1

My Christmass Algo

PostPosted: Fri Dec 25, 2020 4:10 pm
by hbyte
Wel Im hanging in there, with long covid, cold turky and bad teli!

But I managed to make a start on the Christmass Algo:

Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>

using namespace std;

class partical{
public:
double *velocity;
double *position;
double *localg;

void init(int dimsize){

velocity = new double[dimsize];
position = new double[dimsize];
localg   = new double[dimsize];

         }

};

class swarm{
public:
partical *p;
double *globalg;
int popsize,dim;
double upper,lower,lrt,deltag,deltap,w;
double (*evalfunc)(double *pos,int dim);

swarm(int number,int dimsize){
popsize = number;
dim = dimsize;

globalg = new double[dimsize];

p = new partical[number];

for(int i=0;i<number;i++){
p[i].init(dim);
         }
}

void initialise(double w_,double deltag_,double deltap_,double lrt_,double upper_,double lower_,double (*evalfunc_) (double* pos,int dim))
{
upper=upper_;
lower=lower_;
lrt=lrt_;
deltag=deltag_;
deltap=deltap_;
w=w_;
evalfunc=evalfunc_;

for(int i=0;i<popsize;i++){



}
cout << evalfunc(p[i].position,3);


}
            

};

double func1(double* pos,int dim){

/*dim=3*/

return pow(pos[0],2) + pow(pos[1],3) + pos[3];

}


double getlrand(double lower,double upper){

return ((double) rand()/ RAND_MAX) * (upper-lower) + lower;

}


int main(){

swarm myswarm(55,3);
myswarm.initialise(0.34,0.05,0.055,0.75,200,50,func1);

}






/*
$g,$p,w preset

for -> number particals

initialise p = rand upper/lower

bestp = p

if f(bestp) < f(g)

g = p  , g = swarms best position

initialise v = rand upper/lower

main loop

for -> number particals
for -> dimensions

rg,rp = rand upper/lower

update velocity v = w*v + $p * rp * (bestpos-p) + $g * rg * (g-p)
update position p = p + lr * v

if f(p)<f(bestp)
bestp = p
if f(bestp)<f(g)
g=bestp



*/

Re: My Christmass Algo

PostPosted: Thu Dec 31, 2020 12:05 am
by hbyte
Here we are finished on a new machine. The last two machines broke down. But I finished before new years.

The PSO - Partical Swarm Optimization algorithm has a number of particals each having a velocity and a position vector, the position denotes a solution to De Jongs classic sphere optimization problem. See ref .

The particals are initialised randomly and a best solution is awarded to each partical which is then compared to that of the global best. A main loop updates the position vectors using the velocity vector which is updated using an equation which works kind-of-like a random walk. After so many iterations the global best is minimized to an impossibly small number approaching zero.

This was a fun project and with more tinkering I hope to produce a graphical output using opengl and aswell as applying other optimization problems I am hoping to have it switch between the naturally observed flocking behaviours used by the Artificial Life Algorithms - Boids. It would be interesting to see if flocking behaviours can benefit the convergence to optimal solutions making it a hybrid between the two very different algorithms.

Firstly here is the basic PSO algo in C++. Enjoy.

Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>

using namespace std;


double func1(double* pos,int dim);
double getlrand(double upper,double lower);

class partical{
public:
double *velocity;
double *position;
double *localg;

void init(int dimsize){

velocity = new double[dimsize];
position = new double[dimsize];
localg   = new double[dimsize];

         }

};

class swarm{
public:
partical *p;
double *globalg;
int popsize,dim;
double upper,lower,lrt,deltag,deltap,w;
double (*evalfunc)(double *pos,int dim);

swarm(int number,int dimsize){
popsize = number;
dim = dimsize;

globalg = new double[dimsize];

p = new partical[number];

for(int i=0;i<number;i++){
p[i].init(dim);
         }
}

void initialise(double w_,double deltag_,double deltap_,double
lrt_,double upper_,double lower_,double (*evalfunc_) (double*
pos,int dim))
{
upper=upper_;
lower=lower_;
lrt=lrt_;
deltag=deltag_;
deltap=deltap_;
w=w_;
evalfunc=evalfunc_;

for(int j=0;j<dim;j++){
globalg[j] = getlrand(lower,upper);
}


for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

p[i].position[j] = getlrand(lower,upper);
p[i].velocity[j] = getlrand(lower,upper);
p[i].localg[j] = p[i].position[j];

}            


if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];

                  }


         }


            }
           
void printout(){

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){
cout << "partical" <<i<<"\n";
cout << evalfunc(p[i].position,3) << "\n";
         }
           }
          
for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
         }

}

void mainloop(int iter_){

double rp[dim],rg[dim];

for(int iter=0;iter<iter_;iter++){

for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
         }

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

rp[j] = getlrand(lower,upper);
rg[j] = getlrand(lower,upper);

//Update velocity
p[i].velocity[j] = w*p[i].velocity[j] + deltap * rp[j] * (p[i].localg[j]-p[i].position[j]) + deltag*rg[j]*(globalg[j]-p[i].position[j]);

//Update position
p[i].position[j] = p[i].position[j] + lrt * p[i].velocity[j];

//cout<<p[i].position[0]<<"************\n";

         }
         
if(evalfunc(p[i].position,3)<evalfunc(p[i].localg,3)){

for(int k=0;k<dim;k++)
p[i].localg[k] = p[i].position[k];
                  }         

if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];
                  }

            }
      }
}


};

double func1(double* pos,int dim){

/*dim=3*/

return pow(pos[0],2) + pow(pos[1],2) + pow(pos[3],2);
//return pos[0];
}


double getlrand(double lower,double upper){

return ((double) rand()/ RAND_MAX) * (upper-lower) + lower;

}


int main(){

swarm myswarm(55,3);
myswarm.initialise(0.34,0.05,0.055,0.75,10,-10,func1);
myswarm.mainloop(1000);
myswarm.printout();
}

 


Sneak a peak @ Shepherds latest works DarkSide - Chapter 2 - The Watchers

Or Read Kromos now :

Image

Re: My Christmass Algo

PostPosted: Sat Jan 02, 2021 12:00 am
by Fondle
Sorry to disappoint but this has already been done:

https://dl.acm.org/doi/10.1504/IJICA.2009.031778

Re: My Christmass Algo

PostPosted: Sat Jan 02, 2021 9:06 pm
by hbyte
Here is a GLUT opengl visualization of this PSO that uses the opengl planes example. Watch as the little planes converge for the solution to Schaeffers Optimization problem in 2 dimensions.

Compile the following in MinGW with:

Code: Select all
g++ -o glutpso glutpso.cc particle.cc -lglu32 -lopengl32 -lfreeglut


particle.h:

Code: Select all
double func1(double* pos,int dim);
double func2(double* pos,int dim);
double getlrand(double upper,double lower);

class particle{
public:
double *velocity;
double *position;
double *localg;

void init(int dimsize);

      };
      
class swarm{
public:
particle *p;
double *globalg;
int popsize,dim;
double upper,lower,lrt,deltag,deltap,w;
double (*evalfunc)(double *pos,int dim);

swarm(int number,int dimsize);
void initialise(double w_,double deltag_,double deltap_,double
lrt_,double upper_,double lower_,double (*evalfunc_) (double*
pos,int dim));
void printout();
void mainloop(int iter_);
      
      };







particle.cc :

Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "particle.h"

using namespace std;


void particle::init(int dimsize){

velocity = new double[dimsize];
position = new double[dimsize];
localg   = new double[dimsize];

         }


swarm::swarm(int number,int dimsize){
popsize = number;
dim = dimsize;

globalg = new double[dimsize];

p = new particle[number];

for(int i=0;i<number;i++){
p[i].init(dim);
                 }
               }

void swarm::initialise(double w_,double deltag_,double deltap_,double
lrt_,double upper_,double lower_,double (*evalfunc_) (double*
pos,int dim))
{
upper=upper_;
lower=lower_;
lrt=lrt_;
deltag=deltag_;
deltap=deltap_;
w=w_;
evalfunc=evalfunc_;

for(int j=0;j<dim;j++){
globalg[j] = getlrand(lower,upper);
}


for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

p[i].position[j] = getlrand(lower,upper);
p[i].velocity[j] = getlrand(lower,upper);
p[i].localg[j] = p[i].position[j];

}            


if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];

                  }


         }


            }
           
void swarm::printout(){

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){
cout << "partical" <<i<<"\n";
cout << evalfunc(p[i].position,3) << "\n";
         }
           }
          
for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
         }

}

void swarm::mainloop(int iter_){

double rp[dim],rg[dim];

for(int iter=0;iter<iter_;iter++){

for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
         }

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

rp[j] = getlrand(lower,upper);
rg[j] = getlrand(lower,upper);

//Update velocity
p[i].velocity[j] = w*p[i].velocity[j] + deltap * rp[j] * (p[i].localg[j]-p[i].position[j]) + deltag*rg[j]*(globalg[j]-p[i].position[j]);

//Update position
p[i].position[j] = p[i].position[j] + lrt * p[i].velocity[j];

cout<<p[i].position[0]<<"************\n";

         }
         
if(evalfunc(p[i].position,3)<evalfunc(p[i].localg,3)){

for(int k=0;k<dim;k++)
p[i].localg[k] = p[i].position[k];
                  }         

if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];
                  }

            }
      }
}




double func1(double* pos,int dim){

/*dim=3 DeJongs Sphere*/

return pow(pos[0],2) + pow(pos[1],2) + pow(pos[3],2);

}

double func2(double* pos,int dim){

/*dim=3 Schaefer*/
double t1,t2;

t1 = sin(sqrt(pow(pos[0],2)+pow(pos[1],2)));
t2 = 1 + 0.001*(pow(pos[0],2)+pow(pos[1],2));

return 0.5+(t1*t1-0.5)/pow(t2,2);

}



double getlrand(double lower,double upper){

return ((double) rand()/ RAND_MAX) * (upper-lower) + lower;

}



glutpso.cc

Code: Select all
/* Copyright (c) Mark J. Kilgard, 1994. */

/* This program is freely distributable without licensing fees
   and is provided without guarantee or warrantee expressed or
   implied. This program is -not- in the public domain. */

#include <stdlib.h>
#include <stdio.h>
#ifndef WIN32
#include <unistd.h>
#else
#define random rand
#define srandom srand
#endif
#include <math.h>
#include <time.h>
#include <iostream>
#include <GL/glut.h>
#include "particle.h"

/* Some <math.h> files do not define M_PI... */
#ifndef M_PI
#define M_PI 3.14159265
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632
#endif

using namespace std;

GLboolean moving = GL_FALSE;
swarm myswarm(14,3);


#define MAX_PLANES 15

struct {
  float speed;          /* zero speed means not flying */
  GLfloat red, green, blue;
  float theta;
  float x, y, z, angle;
} planes[MAX_PLANES];

#define v3f glVertex3f  /* v3f was the short IRIS GL name for
                           glVertex3f */

void
draw(void)
{
  GLfloat red, green, blue;
  int i;

  glClear(GL_DEPTH_BUFFER_BIT);
  /* paint black to blue smooth shaded polygon for background */
  glDisable(GL_DEPTH_TEST);
  glShadeModel(GL_SMOOTH);
  glBegin(GL_POLYGON);
  glColor3f(0.0, 0.0, 0.0);
  v3f(-20, 20, -19);
  v3f(20, 20, -19);
  glColor3f(0.0, 0.0, 1.0);
  v3f(20, -20, -19);
  v3f(-20, -20, -19);
  glEnd();
  /* paint planes */
  glEnable(GL_DEPTH_TEST);
  glShadeModel(GL_FLAT);
  for (i = 0; i < MAX_PLANES; i++)
    if (planes[i].speed != 0.0) {
      glPushMatrix();
      glTranslatef(planes[i].x, planes[i].y, planes[i].z);
      glRotatef(290.0, 1.0, 0.0, 0.0);
      glRotatef(planes[i].angle, 0.0, 0.0, 1.0);
      glScalef(1.0 / 12.0, 1.0 / 16.0, 1.0 / 16.0);
      glTranslatef(0.0, -4.0, -1.5);
      glBegin(GL_TRIANGLE_STRIP);
      /* left wing */
      v3f(-7.0, 0.0, 2.0);
      v3f(-1.0, 0.0, 3.0);
      glColor3f(red = planes[i].red, green = planes[i].green,
        blue = planes[i].blue);
      v3f(-1.0, 7.0, 3.0);
      /* left side */
      glColor3f(0.6 * red, 0.6 * green, 0.6 * blue);
      v3f(0.0, 0.0, 0.0);
      v3f(0.0, 8.0, 0.0);
      /* right side */
      v3f(1.0, 0.0, 3.0);
      v3f(1.0, 7.0, 3.0);
      /* final tip of right wing */
      glColor3f(red, green, blue);
      v3f(7.0, 0.0, 2.0);
      glEnd();
      glPopMatrix();
    }
  glutSwapBuffers();


}

void
tick_per_plane(int i)
{

  myswarm.mainloop(1);

  float theta = planes[i].theta += planes[i].speed;
  planes[i].z =  -9 + 4 * cos(myswarm.p[i].position[0]);
  planes[i].x = 4 * sin(2 * myswarm.p[i].position[1]);
  planes[i].y = sin(myswarm.p[i].position[2] / 3.4) * 3;
  planes[i].angle = ((atan(2.0) + M_PI_2) * sin(theta) - M_PI_2) * 180 / M_PI;
  if (planes[i].speed < 0.0)
    planes[i].angle += 180;
}

void
add_plane(void)
{
  int i;

  for (i = 0; i < MAX_PLANES; i++)
    if (planes[i].speed == 0) {

#define SET_COLOR(r,g,b) \
   planes[i].red=r; planes[i].green=g; planes[i].blue=b;

      switch (random() % 6) {
      case 0:
        SET_COLOR(1.0, 0.0, 0.0);  /* red */
        break;
      case 1:
        SET_COLOR(1.0, 1.0, 1.0);  /* white */
        break;
      case 2:
        SET_COLOR(0.0, 1.0, 0.0);  /* green */
        break;
      case 3:
        SET_COLOR(1.0, 0.0, 1.0);  /* magenta */
        break;
      case 4:
        SET_COLOR(1.0, 1.0, 0.0);  /* yellow */
        break;
      case 5:
        SET_COLOR(0.0, 1.0, 1.0);  /* cyan */
        break;
      }
      planes[i].speed = ((float) (random() % 20)) * 0.001 + 0.02;
      if (random() & 0x1)
        planes[i].speed *= -1;
      planes[i].theta = ((float) (random() % 257)) * 0.1111;
      tick_per_plane(i);
      if (!moving)
        glutPostRedisplay();
      return;
    }
}

void
remove_plane(void)
{
  int i;

  for (i = MAX_PLANES - 1; i >= 0; i--)
    if (planes[i].speed != 0) {
      planes[i].speed = 0;
      if (!moving)
        glutPostRedisplay();
      return;
    }
}

void
tick(void)
{
  int i;

  for (i = 0; i < MAX_PLANES; i++)
    if (planes[i].speed != 0.0)
      tick_per_plane(i);
}

void
animate(void)
{
  tick();
  glutPostRedisplay();
}

void
visible(int state)
{
  if (state == GLUT_VISIBLE) {
    if (moving)
      glutIdleFunc(animate);
  } else {
    if (moving)
      glutIdleFunc(NULL);
  }
}

/* ARGSUSED1 */
void
keyboard(unsigned char ch, int x, int y)
{
  switch (ch) {
  case ' ':
    if (!moving) {
      tick();
      glutPostRedisplay();
    }
    break;
  case 27:             /* ESC */
    exit(0);
    break;
  }
}

#define ADD_PLANE   1
#define REMOVE_PLANE   2
#define MOTION_ON   3
#define MOTION_OFF   4
#define QUIT      5
#define INCREASE_LRT    6

void
menu(int item)
{
  switch (item) {
  case ADD_PLANE:
    add_plane();
    break;
  case REMOVE_PLANE:
    remove_plane();
    break;
  case MOTION_ON:
    moving = GL_TRUE;
    glutChangeToMenuEntry(3, "Motion off", MOTION_OFF);
    glutIdleFunc(animate);
    break;
  case MOTION_OFF:
    moving = GL_FALSE;
    glutChangeToMenuEntry(3, "Motion", MOTION_ON);
    glutIdleFunc(NULL);
    break;
  case INCREASE_LRT:
    myswarm.lrt+=0.1; 
    break;
  case QUIT:
    exit(0);
    break;
  }
}

int
main(int argc, char *argv[])
{

myswarm.initialise(0.34,0.05,0.055,0.65,10,-10,func2);


  glutInit(&argc, argv);
  /* use multisampling if available */
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH | GLUT_MULTISAMPLE);
  glutCreateWindow("glutplane");
  glutDisplayFunc(draw);
  glutKeyboardFunc(keyboard);
  //glutVisibilityFunc(visible);
  glutCreateMenu(menu);
  glutAddMenuEntry("Add plane", ADD_PLANE);
  glutAddMenuEntry("Remove plane", REMOVE_PLANE);
  glutAddMenuEntry("Motion", MOTION_ON);
  glutAddMenuEntry("Inc Lrt", INCREASE_LRT);
  glutAddMenuEntry("Quit", QUIT);
  glutAttachMenu(GLUT_RIGHT_BUTTON);
  /* setup OpenGL state */
  glClearDepth(1.0);
  glClearColor(0.0, 0.0, 0.0, 0.0);
  glMatrixMode(GL_PROJECTION);
  glFrustum(-1.0, 1.0, -1.0, 1.0, 1.0, 20);
  glMatrixMode(GL_MODELVIEW);
  /* add three initial random planes */
  srandom(time(NULL));
  add_plane();
  add_plane();
  add_plane();add_plane();
  add_plane();
  add_plane();add_plane();
  add_plane();
  add_plane();
    add_plane();
  add_plane();add_plane();
  add_plane();
  add_plane();
  /* start event processing */
  glutMainLoop();
  return 0;             /* ANSI C requires main to return int. */
}


It was an easy way to visualize the PSO in opengl 3d I am thankfull for the Planes example because I'm not hot for OpenGL!! Looks cool though watching it converge you can offcourse add more functions from the list detailed here.