IntPoly.cpp File Reference

#include <stdio.h>
#include "Poly.h"
#include "Vector.h"
#include "tools.h"
#include "KeepBests.h"
#include "IntPoly.h"

Include dependency graph for IntPoly.cpp:

Go to the source code of this file.

Functions
Vector	LAGMAXModified (Vector G, Matrix H, double rho, double &VMAX)
int	findK (double *ValuesF, int n, ObjectiveFunction *of, Vector *points)
int	calculateNParallelJob (int n, double *vf, Vector *cp, ObjectiveFunction *of, int *notsuccess)

Function Documentation

calculateNParallelJob()

int calculateNParallelJob	(	int	n,
		double *	vf,
		Vector *	cp,
		ObjectiveFunction *	of,
		int *	notsuccess
	)

Definition at line 43 of file parallel.cpp.

{
    int i,r,nsuccess=0;
    for (i=0; i<n; i++)
    {
        r=0;
        vf[i]=of->eval(cp[i],&r);
        notsuccess[i]=r;
        if (!r) nsuccess++;
    }
    return nsuccess;
    
}

findK()

int findK	(	double *	ValuesF,
		int	n,
		ObjectiveFunction *	of,
		Vector *	points
	)

Definition at line 214 of file IntPoly.cpp.

{
    // if (of->isConstrained) return 0;
    // find index k of the best value of the function 
    double minimumValueF=INF;
    int i,k=-1; 
    for (i=0; i<n; i++)
        if ((ValuesF[i]<minimumValueF)
             &&((of==NULL)||(of->isFeasible(points[i])))) 
        { k=i; minimumValueF=ValuesF[i]; }
    if (k==-1) k=0;
    return k;
}

LAGMAXModified()

Vector LAGMAXModified	(	Vector	G,
		Matrix	H,
		double	rho,
		double &	VMAX
	)

Definition at line 40 of file MSSolver.cpp.

{
    //not tested in depth but running already quite good

//    SUBROUTINE LAGMAX (N,G,H,RHO,D,V,VMAX)
//    IMPLICIT REAL*8 (A-H,O-Z)
//    DIMENSION G(*),H(N,*),D(*),V(*)
//
//    N is the number of variables of a quadratic objective function, Q say.
//    G is the gradient of Q at the origin.
//    H is the symmetric Hessian matrix of Q. Only the upper triangular and
//      diagonal parts need be set.
//    RHO is the trust region radius, and has to be positive.
//    D will be set to the calculated vector of variables.
//    The array V will be used for working space.
//    VMAX will be set to |Q(0)-Q(D)|.
//
//     Calculating the D that maximizes |Q(0)-Q(D)| subject to ||D|| .EQ. RHO
//    requires of order N**3 operations, but sometimes it is adequate if
//    |Q(0)-Q(D)| is within about 0.9 of its greatest possible value. This
//    subroutine provides such a solution in only of order N**2 operations,
//    where the claim of accuracy has been tested by numerical experiments. 
/*
    int n=G.sz();
    Vector D(n), V(n);
    lagmax_(&n, (double *)G, *((double**)H), &rho, 
                    (double*)D, (double*)V, &VMAX);
    return D;
*/
    int i,n=G.sz();
    Vector D(n);

    Vector V=H.getMaxColumn();
    D=H.multiply(V);
    double vv=V.square(),
           dd=D.square(),
           vd=V.scalarProduct(D), 
           dhd=D.scalarProduct(H.multiply(D)),
           *d=D, *v=V, *g=G;

//
//     Set D to a vector in the subspace spanned by V and HV that maximizes
//     |(D,HD)|/(D,D), except that we set D=HV if V and HV are nearly parallel.
//
    if (sqr(vd)<0.9999*vv*dd)
    {
        double a=dhd*vd-dd*dd,
               b=.5*(dhd*vv-dd*vd),
               c=dd*vv-vd*vd,
               tmp1=-b/a;
        if (b*b>a*c)
        {
            double tmp2=sqrt(b*b-a*c)/a, dd1, dd2, dhd1, dhd2;
            Vector D1=D.clone();
            D1.multiply(tmp1+tmp2);
            D1+=V;
            dd1=D1.square();
            dhd1=D1.scalarProduct(H.multiply(D1));

            Vector D2=D.clone();
            D2.multiply(tmp1-tmp2);
            D2+=V;
            dd2=D2.square();
            dhd2=D2.scalarProduct(H.multiply(D2));

            if (condorAbs(dhd1/dd1)>condorAbs(dhd2/dd2)) { D=D1; dd=dd1; dhd=dhd1; } 
            else { D=D2; dd=dd2; dhd=dhd2; }
            d=(double*)D;
        }
    };


//
//     We now turn our attention to the subspace spanned by G and D. A multiple
//     of the current D is returned if that choice seems to be adequate.
//
    double gg=G.square(),
           normG=sqrt(gg),
           gd=G.scalarProduct(D), 
           temp=gd/gg,
           scale=sign(rho/sqrt(dd), gd*dhd);

    i=n; while (i--) v[i]=d[i]-temp*g[i];
    vv=V.square();

    if ((normG*dd)<(0.5-2*rho*condorAbs(dhd))||(vv/dd<1e-4))
    {
        D.multiply(scale);
        VMAX=condorAbs(scale*(gd+0.5*scale*dhd));
        return D;
    }

//
//     G and V are now orthogonal in the subspace spanned by G and D. Hence
//     we generate an orthonormal basis of this subspace such that (D,HV) is
//     negligible or zero, where D and V will be the basis vectors.
//
    H.multiply(D,G); //  D=HG;
    double ghg=G.scalarProduct(D),
           vhg=V.scalarProduct(D),
           vhv=V.scalarProduct(H.multiply(V));
    double theta, cosTheta, sinTheta;

    if (condorAbs(vhg)<0.01*mmax(condorAbs(vhv),condorAbs(ghg)))
    {
        cosTheta=1.0;
        sinTheta=0.0;
    } else
    {
        theta=0.5*atan(0.5*vhg/(vhv-ghg));
        cosTheta=cos(theta);
        sinTheta=sin(theta);
    }
    i=n;
    while(i--)
    {
        d[i]= cosTheta*g[i]+ sinTheta*v[i];
        v[i]=-sinTheta*g[i]+ cosTheta*v[i];
    };

//
//     The final D is a multiple of the current D, V, D+V or D-V. We make the
//     choice from these possibilities that is optimal.
//

    double norm=rho/D.euclidianNorm();
    D.multiply(norm);
    dhd=(ghg*sqr(cosTheta)+vhv*sqr(sinTheta))*sqr(norm);

    norm=rho/V.euclidianNorm();
    V.multiply(norm);
    vhv=(ghg*sqr(sinTheta)+vhv*sqr(cosTheta)*sqr(norm));

    double halfRootTwo=sqrt(0.5),   // =sqrt(2)/2=cos(PI/4)
           t1=normG*cosTheta*rho,   // t1=condorAbs(D.scalarProduct(G));
           t2=normG*sinTheta*rho,   // t2=condorAbs(V.scalarProduct(G));
           at1=condorAbs(t1),
           at2=condorAbs(t2),
           t3=0.25*(dhd+vhv),
           q1=condorAbs(at1+0.5*dhd),
           q2=condorAbs(at2+0.5*vhv),
           q3=condorAbs(halfRootTwo*(at1+at2)+t3),
           q4=condorAbs(halfRootTwo*(at1-at2)+t3);
    if ((q4>q3)&&(q4>q2)&&(q4>q1))
    {
        double st1=sign(t1*t3), st2=sign(t2*t3);
        i=n; while (i--) d[i]=halfRootTwo*(st1*d[i]-st2*v[i]);
        VMAX=q4;
        return D;
    }
    if ((q3>q2)&&(q3>q1))
    {
        double st1=sign(t1*t3), st2=sign(t2*t3);
        i=n; while (i--) d[i]=halfRootTwo*(st1*d[i]+st2*v[i]);
        VMAX=q3;
        return D;
    }
    if (q2>q1)
    {
        if (t2*vhv<0) V.multiply(-1);
        VMAX=q2;
        return V;
    }
    if (t1*dhd<0) D.multiply(-1);
    VMAX=q1;
    return D;
}