void testApp::FoxLiIntegral()
{
vector src_field;
src_field.resize(m_field.size());
memcpy(&src_field[0], &m_field[0], m_field.size()*sizeof(ofVec4f));
float k = (2*PI/m_wave_length)*1e3;
parallel_for(0,(int)src_field.size(),1,[&](int i){
if (m_mask[i] == 0)
return;
m_field[i].z = m_field[i].w = 0;
for (int si = 0; si < src_field.size(); ++si)
{
if (m_mask[si] == 0)
continue;
float dx = m_field[i].x - src_field[si].x;
float dy = m_field[i].y - src_field[si].y;
float dist = sqrt(dx*dx + dy*dy + m_resonator_length*m_resonator_length);
float kr = k*dist;
float sinkr = sin(kr);
float coskr = cos(kr);
float c = (1 + m_resonator_length/dist)/dist;
m_field[i].z += c*(src_field[si].z*sinkr - src_field[si].w*coskr);
m_field[i].w += c*(src_field[si].z*coskr + src_field[si].w*sinkr);
}
});
NormalizeField();
CheckResidual(m_field, src_field);
}
void testApp::UpdateResult()
{
m_amp.resize(m_field.size());
m_phase.resize(m_field.size());
for (int i = 0; i < m_field.size(); i++)
{
if (m_mask[i] == 0)continue;
ofVec4f& vec = m_field[i];
m_amp[i] = sqrt(vec.z*vec.z + vec.w*vec.w);
m_phase[i] = atan2(vec.w, vec.z)/PI;
}
}
void testApp::DrawDistribution(vector& strength)
{
bool mesh_break = true;
for (int y = 0; y < m_rows - 1; ++y)
{
int y_offset = y*m_cols;
for (int x = 0; x < m_cols; ++x)
{
int index = y_offset + x;
if (m_mask[index] == 0 || m_mask[index + m_cols] == 0)
{
if (!mesh_break)
{
//If there's no point here, the mesh should break.
mesh_break = true;
glEnd();
}
continue;
}
ofVec4f& space_point1 = m_field[index];
ofVec4f& space_point2 = m_field[index + m_cols];
if (mesh_break)
{
//Start connecting points to form mesh.
glBegin(GL_TRIANGLE_STRIP);
mesh_break = false;
}
float amp_1 = strength[index];
float amp_2 = strength[index + m_cols];
//Draw the point and set its normal.
glColor3f(amp_1, 0.2f, 1 - amp_1);
//glNormal3f(dv21.x, dv21.y, dv21.z);
glVertex3f(space_point1.x, amp_1, space_point1.y);
//Draw the point below the prior one to form a triangle.
glColor3f(amp_2, 0.2f, 1 - amp_2);
glVertex3f(space_point2.x, amp_2, space_point2.y);
}
if (!mesh_break)
{
//We break the mesh at the end of the line,.
glEnd();
mesh_break = true;
}
}
}
void testApp::CreateGui()
{
m_gui = new ofxUICanvas(MENU_WIDTH,480);
m_gui->setFont("\GUI\NewMedia Fett.ttf");
m_gui->addLabel("Wave Length (um):", OFX_UI_FONT_SMALL);
m_ui_wavelength = m_gui->addTextInput("WaveLength","0.555",80,26);
m_gui->addLabel("Resonator Length (mm):", OFX_UI_FONT_SMALL);
m_ui_reslength = m_gui->addTextInput("ResLength","1000",80,26);
m_gui->addSpacer();
m_gui->addLabel("Mirror Type:", OFX_UI_FONT_SMALL);
vector<string> vnames; vnames.push_back("One Dimension"); vnames.push_back("Rectangle"); vnames.push_back("Circle");
m_ui_mirrortype = m_gui->addRadio("MirrorType", vnames, OFX_UI_ORIENTATION_VERTICAL);
m_ui_mirrortype->activateToggle("Rectangle");
m_gui->addLabel("Mirror Width/Diameter (mm):", OFX_UI_FONT_SMALL);
m_ui_width = m_gui->addTextInput("Width","1",50,26);
m_gui->addLabel("Mirror Height (mm):", OFX_UI_FONT_SMALL);
m_ui_height = m_gui->addTextInput("Height","1",50,26);
m_gui->addSpacer();
m_gui->addLabel("Iterate Count:", OFX_UI_FONT_SMALL);
m_ui_itcount = m_gui->addTextInput("ItCount","10",80,26);
m_gui->addLabel("Minimum Variance:", OFX_UI_FONT_SMALL);
m_ui_minerror = m_gui->addTextInput("MinError","0.01",80,26);
m_gui->addSpacer();
ofxUILabelButton* btn_start = new ofxUILabelButton("Start",false,80,26);
ofxUILabelButton* btn_stop = new ofxUILabelButton("Stop",false,80,26);
m_gui->addWidgetPosition(btn_start, OFX_UI_WIDGET_POSITION_DOWN, OFX_UI_ALIGN_RIGHT);
m_gui->addWidgetPosition(btn_stop, OFX_UI_WIDGET_POSITION_DOWN, OFX_UI_ALIGN_RIGHT);
m_gui->setColorBack(ofColor(255,100));
m_gui->setWidgetColor(OFX_UI_WIDGET_COLOR_BACK, ofColor(255,100));
ofAddListener(m_gui->newGUIEvent,this,&testApp::guiEvent);
}
class IterationThread : public ofThread
{
protected:
void threadedFunction()
{
int iterate_count = 0;
float residual = FLT_MAX;
while(isThreadRunning() && iterate_count < m_iterate_count && residual > m_min_residual)
{
++iterate_count;
if (Iterate)
residual = Iterate();
if (ProgressCallback)
ProgressCallback(iterate_count, residual);
}
if (FinishedCallback)FinishedCallback();
}
public:
IterationThread():Iterate(NULL),FinishedCallback(NULL),ProgressCallback(NULL){}
void SetIterationFunction(float (*iterate)(), int iterate_count, float min_residual)
{
Iterate = iterate;
m_iterate_count = iterate_count;
m_min_residual = min_residual;
}
void SetFinishCallback(void (*finished)())
{
FinishedCallback = finished;
}
void SetProgressCallback(void (*progress)(int,float))
{
ProgressCallback = progress;
}
private:
int m_iterate_count;
float m_min_residual;
float(*Iterate)();
void(*FinishedCallback)();
void(*ProgressCallback)(int,float);
};
https://zhidao.baidu.com/question/2139883596015832908.html
#include <iostream>
#include <math.h>
#include <time.h>
using namespace std;
#define M 50 //群体数目50
#define N 4 //每个粒子的维数4
//测试类
class TestFunction
{
public:
double resen(double x1,double x2,double x3,double x4)
{
double s=0;
s=100*(x2-x1*x1)*(x2-x1*x1)+(1-x1)*(1-x1)+s;
s=100*(x3-x2*x2)*(x3-x2*x2)+(1-x2)*(1-x2)+s;
s=100*(x4-x3*x3)*(x4-x3*x3)+(1-x3)*(1-x3)+s;
return s;
}
};
class CQPSO
{
private:
double delta;
double (*w)[N];// = new double[50][4]; //总体粒子
double *f;//=new double[M];//适应度值
double *ff;//=new double[M];//相对f的比较值
double (*p)[N];//=new double[M][N];
double *g;//=new double[N];
double *c;//=new double[N];
TestFunction *tf;// = new TestFunction;
double random()
{
double s;
s=(abs(rand())%10000+10000)/10000.0-1.0;
return s;
}
public:
CQPSO(double delta)
{
int i,j;
this->delta=delta;
w=new double[M][N];
f=new double[M];
ff=new double[M];
p=new double[M][N];
g=new double[N];
c=new double[N];
tf=new TestFunction;
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
{
w[i][j]=random();
}
}
}
void CQPSOmethod(int count)
{
int i,j;
if(count==1)
{
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
{
p[i][j]=w[i][j];
}
f[i]=tf->resen(w[i][0],w[i][1],w[i][2],w[i][3]);
}
cqpso_p(f);//得出全局最优
}
if(count>1)
{
cqpso_update( w );
for(i=0;i<M;i++)
{
ff[i]=tf->resen(w[i][0],w[i][1],w[i][2],w[i][3]);
if(ff[i]<f[i])
{
f[i]=ff[i];
for(j=0;j<N;j++) p[i][j]=w[i][j];
}
}
cqpso_p(f);
}
cout<<(tf->resen(g[0],g[1],g[2],g[3]))<<endl;
}
void cqpso_p(double *f)//得到个体最优中最小值——全局最优
{
double temp=f[0];
int i,j;
for(i=1;i<M;i++)
{
if(f[i]<temp)
{
temp=f[i];
}
}
for(i=0;i<M;i++)
{
if(temp==f[i])
{
for(j=0;j<N;j++)
{
g[j]=p[i][j];
}
break;
}
}
}
void cqpso_c(double (*p)[N])
{
int i,j;
for(i=0;i<N;i++) c[i]=0;
for(i=0;i<N;i++)
{
for(j=0;j<M;j++)
{
c[i]=c[i]+p[j][i];
}
}
for(i=0;i<N;i++) c[i]=c[i]/M;
}
void cqpso_update( double (*w)[N] )
{
int i,j;
double *fai=new double[N];
double (*u)[N]=new double[M][N];
double (*pp)[N]=new double[M][N];
for(i=0;i<N;i++)
{
fai[i]=random();
}
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
u[i][j]=random();
}
cqpso_c( p );
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
pp[i][j]=fai[j]*p[i][j]+g[j]*(1-fai[j]);
}
for(i=0;i<M;i++)
{
for(j=0;j<N;j++)
w[i][j]=pp[i][j]+delta*(abs(c[j]-w[i][j]))*log(1/u[i][j]);
}
}
};
int main()
{
int i;
srand((unsigned)time(0));
CQPSO *qo = new CQPSO(0.5);
//qo->w=new double[M][N];
for(i=1;i<100;i++)
qo->CQPSOmethod(i);
}