My Christmass Algo

Put your Vanilla code here

My Christmass Algo

Postby hbyte » Fri Dec 25, 2020 4:10 pm

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



*/
hbyte
Site Admin
 
Posts: 80
Joined: Thu Aug 13, 2020 6:11 pm

Re: My Christmass Algo

Postby hbyte » Thu Dec 31, 2020 12:05 am

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
hbyte
Site Admin
 
Posts: 80
Joined: Thu Aug 13, 2020 6:11 pm

Re: My Christmass Algo

Postby Fondle » Sat Jan 02, 2021 12:00 am

Sorry to disappoint but this has already been done:

https://dl.acm.org/doi/10.1504/IJICA.2009.031778
Fondle
 
Posts: 5
Joined: Fri Aug 14, 2020 10:33 am

Re: My Christmass Algo

Postby hbyte » Sat Jan 02, 2021 9:06 pm

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.
hbyte
Site Admin
 
Posts: 80
Joined: Thu Aug 13, 2020 6:11 pm


Return to Python and ML

Who is online

Users browsing this forum: No registered users and 2 guests

cron