regression.c

本文在VC++环境下用演化算法实现了非线性参数回归

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <ctype.h>
#include <values.h>
#include <string.h>

/*  Definition of control parameters */
#define POPSIZE 50      /*population size*/
#define NVARS 3         /*number of parameters*/
#define PXOVER 0.9      /*probability of crossover*/
#define PMUTATION 0.15  /*probability of mutation*/

typedef struct   /*genotype (GT), an individual of the population*/
{
    float chrom[NVARS];   /*a string of variables*/
    float lsquare;        /*GT's fitness*/
    float r;              /*GT's relative coefficient*/
    int fitness;          /*relative fitness*/
}individual;

static individual newpop[POPSIZE+1];  /*new population*/
static individual oldpop[POPSIZE+1];  /*old population*/
static int location[POPSIZE+1];       /*location of individuals after sorting*/

static int gen;      /*current generation no.*/
static int maxgen;   /*maximum generation*/
static int change=0; /*no. of generations that fitness has no progress*/

static int ndata;    /*pairs of experimental data*/
static float *datax; /*sampling points*/
static float *datay; /*obtained values*/
static float sumdatay,sumdatax,sum2datax,ymax;
                     /*values for calculating relative coefficient*/
static float a[NVARS],b[NVARS];
                     /*lower bound and upper bound of variables*/

/* Macro definitions of random number generators */
/*#define rdint(i)   (rand%(int)(i))*/
/*#define rdft()     ((float)rdint(16384)/(16383.0))*/
/*#define rnd(a,b)   (rdint((int)(b)-(int)(a)+1)+(int)(a))*/
/*#define flip(prob) (rdft()<=prob?1:0)*/

int rdint(int i)
{
    int tmp=0;

    tmp=rand()%(int)(i);

    return tmp;
}

float rdft(void)
{
    float tmp=0.0;

    tmp=(float)rdint(16384)/(16383.0);

    return tmp;
}

int flip(float prob)
{
    int tmp;

    tmp=rdft()<=prob?1:0;

    return tmp;
}

float rnd(float a,float b)
{
    float tmp;

    tmp=rdint((int)(b)-(int)(a)+1)+(int)(a);

    return tmp;
}

/* Evaluation function */
float evaluate(float *prms)
{
    int i;
    float y,serr=0;

    for(i=0;i<ndata;i++)
    {
        if((prms[2]<=0)||(datax[i]*prms[2]>500))
            return 1.0e10;  /*dealing with overflow*/
        else
            y=prms[0]-prms[1]*pow(prms[2],datax[i]);

        serr+=(y-datay[i])*(y-datay[i]);
    }

    return serr;
}

/* Relative coefficient calculation function */
float cal_coeff(float *prms)
{
    int i;
    float coeff,sigma1,sigma2;
    float fy,sumfy=0,sum2fy=0,sumxfy=0;

    for(i=0;i<ndata;i++)
    {
        if((prms[2]<=0)||(datax[i]*prms[2]>500))
            return 1.0e10;
        else
    	{
            fy=prms[0]-prms[1]*pow(prms[2],datax[i]);
            sumfy+=fy;
            sum2fy+=fy*fy;
            sumxfy+=datax[i]*fy;
    	}
    }

    sigma1=fabs((ndata*sum2datax)-(sumdatax*sumdatax));
    sigma2=fabs((ndata*sum2fy)-(sumfy*sumfy));
    coeff=(ndata*sumxfy-sumdatax*sumfy)/sqrt(sigma1*sigma2);

    return coeff;
}

/* Function to determine the sorting location of an array by quick sort */
void quick_sort(float *item,int left,int right, int *locat)
{
    register int i,j,k;
    float fa,fb;

    i=left;j=right;
    fa=item[(left+right)/2];

    do
    {
        while((item[i]>fa)&&(i<right))i++;
        while((fa>item[j])&&(j>left))j--;

        if(i<=j)
    	{
            fb=item[i];
            item[i]=item[j];
            item[j]=fb;
            k=locat[i];
            locat[i]=locat[j];
            locat[j]=k;
        	i++;
        	j--;
    	}
    }while(i<=j);

    if(left<j)
        quick_sort(item,left,right,locat);
    if(i<right)
        quick_sort(item,left,right,locat);
}

/* Function to determine the sorting location of individual in a population
   and normalize the fitness of individuals for linear ranking selection
   and a minimization problem */
void normalfitness(individual *pop)
{
    float fit[POPSIZE];
    int i;

    for(i=0;i<POPSIZE;i++)
    {
        fit[i]=pop[i].lsquare;
        location[i]=i;
    }

    quick_sort(fit,0,POPSIZE-1,location);

    for(i=0;i<POPSIZE;i++)
        pop[location[i]].fitness=2*(i+1);
}

/* Function to initialize the population */
void initpop(void)
{
    int i,j;

    for(i=0;i<POPSIZE;i++)
    {
        for(j=0;j<NVARS;j++)
            oldpop[i].chrom[j]=a[j]+rdft()*(b[j]-a[j]);

        oldpop[i].lsquare=evaluate(oldpop[i].chrom);
        oldpop[i].fitness=0;
    }

    normalfitness(oldpop);
    oldpop[POPSIZE]=oldpop[location[POPSIZE-1]];
}

/* Selection function: Standard roulette wheel selection */
int select(individual *pop)
{
    float rand1;
    int partsum=0;
    int i=0;

    rand1=rdft()*POPSIZE*(POPSIZE+1);

    for(i=0;(rand1>=partsum);i++)
        partsum+=pop[i].fitness;

    return i-1;
}

/* Mutation: Random unitorm mutation */
float mutate(int i)
{
    return a[i]+rdft()*(b[i]-a[i]);
}

/* Crossover: performs atithmetical crossover of the two
   selected parents and products two offsprings */
void xover(float *parent1,float *parent2,float *child1,float *child2)
{
    float alpha;
    int i;

    for(i=0;i<NVARS;i++)
    {
        alpha=rdft();
        child1[i]=alpha*parent1[i]+(1-alpha)*parent2[i];
        child2[i]=alpha*parent2[i]+(1-alpha)*parent1[i];

        if(rdft()<PMUTATION)
            child1[i]=mutate(i);
        if(rdft()<PMUTATION)
            child2[i]=mutate(i);
    }
}

/* Initialize function: Reports parameters of the algorithm */
void initialize(int argc,char *argv[])
{
    int i;
    FILE *datafp;  /*pointer of experimental data file*/
    srand((unsigned)time(NULL));  /*initialize random number generator*/

    /*printf("")*/

    switch(argc) 
    {
    case 1:
        printf("\n\07Usage:%s datafile\n\n",argv[0]);
        exit(-1);
        break;
    case 2:
        if((datafp==fopen(argv[1],"rt"))==NULL)
    	{
            printf("cannot open data file %s!\n",argv[1]);
            exit(1);
    	}
        break;
    default:
        printf("invalid parameters!\n");
        exit(1);
    }
    fscanf(datafp,"%d",&ndata);

    if((datax==(float *)malloc((ndata+1)*sizeof(float)))==NULL)
    {
        printf("Out of memory!\n");
        exit(1);
    }
    if((datay==(float *)malloc((ndata+1)*sizeof(float)))==NULL)
    {
        printf("Out of memory!\n");
        exit(1);
    }

    for(i=0;i<ndata;i++)
        fscanf(datafp,"%f %f",&datax[i],&datay[i]);

    fclose(datafp);

    for(i=0,sumdatay=0;i<ndata;i++)
        sumdatay+=datay[i];
    for(i=0,sumdatax=0;i<ndata;i++)
        sumdatax+=datax[i];
    for(i=0,sum2datax=0;i<ndata;i++)
        sum2datax+=(datax[i])*(datax[i]);
    for(i=0,ymax=-1.0e10;i<ndata;i++)
    {
        if(ymax<datay[i])
            ymax=datay[i];
    }

    printf("\nEnter max generation:");
    scanf("%d",&maxgen);
    printf("G-Generation u-unchanged generations");
    printf("err-least square error r-relative coefficient.\n\n");

    printf("Enter the lower and upper bound of the parameters:\n");
    for(i=0;i<NVARS;i++)
    {
        printf("\rlower bound of parameter %d=",i);
        scanf("%f",&a[i]);
        printf("\rupper bound of parameter %d=",i);
        scanf("%f",&b[i]);
    }

    location[POPSIZE]=POPSIZE;
    initpop();
}

/* Recombine function: Creats a new population from an old one
   through selection, crossover, mutation. */
void recombine(void)
{
    int i;
    int mate1,mate2;

    while(i<POPSIZE)
    {
        /*select two parents from the old population*/
        mate1=select(oldpop);

        do 
    	{
            mate2=select(oldpop);
        }while(mate1==mate2);

        if(rdft()<PXOVER)
    	{
            xover(oldpop[mate1].chrom,oldpop[mate2].chrom,
                newpop[i].chrom,newpop[i+1].chrom);

            newpop[i].lsquare=evaluate(newpop[i].chrom);
            newpop[i+1].lsquare=evaluate(newpop[i+1].chrom);
    	}
        else
    	{
            newpop[i]=oldpop[mate1];
            newpop[i+1]=oldpop[mate2];
    	}
        i=i+2;
    }
    /* select an individual from the new population randomly,
     * if it is better than the best individual from the old
     * population, then replace it; otherwise leave it unchanged.
     */
    i=rnd(0,POPSIZE-1);
    if(newpop[i].lsquare>oldpop[POPSIZE].lsquare)
        newpop[i]=oldpop[POPSIZE];

    normalfitness(newpop);

    /* if best individual from the new population is better
     * than the best individual from the previous popultaion,
     * copy the best from the new population; else replace the 
     * best individual from the current population with the best
     * one from the previous generation.
     */
    
    if(newpop[location[POPSIZE-1]].lsquare<oldpop[POPSIZE].lsquare)
    {
        newpop[POPSIZE]=newpop[location[POPSIZE-1]];
        change=0;
    }
    else
    {
        newpop[POPSIZE]=oldpop[POPSIZE];
        change++;
    }
    for(i=0;i<POPSIZE;i++)
        newpop[i].r=cal_coeff(newpop[i].chrom);
}

/* Report function */
void report(int gen)
{
    int i;
    printf("\rG=%3d u=%2d",gen,change);

    for(i=0;i<NVARS;i++)
        printf("par%d=%f",i,newpop[POPSIZE].chrom[i]);

    printf("err=%f ",newpop[POPSIZE].lsquare);

    getch();
}

/* Main function */
void main(int argc, char *argv[])
{
    int i;
    long key;

    gen=0;

    initialize(argc,argv);

    for(gen=1;gen<=maxgen;i++)
    {
        recombine();
        report(gen);

        /* copy the new population to the old one */

        for(i=0;i<=POPSIZE;i++)
            oldpop[i]=newpop[i];
    }
}