a_fn.cxx

Go to the documentation of this file.

/*
         Copyright 2011-2019 Pierre SMARS (smars@yuntech.edu.tw)
         This program is free software: you can redistribute it and/or modify
         it under the terms of the GNU General Public License as published by
         the Free Software Foundation, either version 2 of the License, or
         (at your option) any later version.
 
         This program is distributed in the hope that it will be useful,
         but WITHOUT ANY WARRANTY; without even the implied warranty of
         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
         GNU General Public License for more details.
 
         You should have received a copy of the GNU General Public License
         along with this program.  If not, see <http://www.gnu.org/licenses/>.
         */
#include "a_fn.h"
#include "a_mat_sq.h"
#include "a_mat_c.h"
#include "a_mat_r.h"
 
// .PORTABILITY : ansi C++
 
//---------------------------------------------------------------------------
a_fn::a_fn(int nv, double dmin, double gmin) : nv_(nv), dmin_(dmin), gmin_(gmin)
{
        x_ = new a_mat_c(nv_);
        x0_ = new a_mat_c(nv_);
}
//---------------------------------------------------------------------------
a_fn::~a_fn()
{
        delete x_;
}
//---------------------------------------------------------------------------
void a_fn::showX(std::ostream& o)
{
        if (nv_==1)
                o << (*x_)(0) << std::endl;
        else
        {
                for (int i = 0; i < nv_-1; i++)
                        o << (*x_)(i) << " ";
                o << (*x_)(nv_-1) << std::endl;
        }
}
//---------------------------------------------------------------------------
void a_fn::setX(const a_mat_c& x)
{
        for (int i = 0; i < nv_; i++)
                (*x_)(i) = x(i);
}
//---------------------------------------------------------------------------
double a_fn::getdF(int ai)
{
        double p0 = (*x_)(ai);
        (*x_)(ai) = 0.;
        (*x_)(ai) = p0-dmin_/2.;
        double Fmin = this->getF();
        (*x_)(ai) = p0+dmin_/2.;
        double Fmax = this->getF();
        (*x_)(ai) = p0;
        return (Fmax-Fmin)/dmin_;
}
//---------------------------------------------------------------------------
void a_fn::getGrad(a_mat_c& grad)
{
        for (int i = 0; i < nv_; i++)
                grad(i) = this->getdF(i);
}
//---------------------------------------------------------------------------
void a_fn::dx(int ai)
{
        double ax = this->getX(ai);
        double dfp = this->getdF(ai);
        double fp = this->getF();
        (*x_)(ai) = ax-fp/dfp;
}
//---------------------------------------------------------------------------
void a_fn::newton(int ai)
{
        double f0 = this->getF();
        double df = this->getdF(ai);
        if (df != 0)
        {
                do
                {
                        this->dx(ai);
                        double f = this->getF();
                        df = fabs(f0-f);
                        f0 = f;
                }
                while (df>dmin_);
        }
}
//---------------------------------------------------------------------------
void a_fn::incr(double alpha, const a_mat_c& dir)
{
        for (int i = 0; i < nv_; i++)
                (*x_)(i) += dir(i)*alpha;
}
//---------------------------------------------------------------------------
void a_fn::go(double alpha, const a_mat_c& dir)
{
        for (int i = 0; i < nv_; i++)
                (*x_)(i) = (*x0_)(i) + dir(i)*alpha;
        /*      if (alpha != alpha0_)
                        {
                        this->incr(alpha-alpha0_,dir);
                        alpha0_ = alpha;
                        }*/
}
//---------------------------------------------------------------------------
double a_fn::gF(double alpha, const a_mat_c& dir)
{
        this->go(alpha,dir);
        return this->getF();
}
//---------------------------------------------------------------------------
void a_fn::minGold0(const a_mat_c& dir, double A1, double B1)
{
        double f0 = this->getF();
        const double p = 0.381966011249;
        double crit0 = (B1-A1)/1.e10;
        bool t1 = true;
        bool t2 = true;
        this->setx0();
        double x1,x2;
        do
        {
                std::cout << ".";
                double f1,f2;
                double L1 = B1-A1;
                double L2 = p*L1;
                if (t1)
                {
                        x1 = A1+L2;
                        f1 = this->gF(x1,dir);
                }
                if (t2)
                {
                        x2 = B1-L2;
                        f2 = this->gF(x2,dir);
                }
                if (f1<f2)
                {
                        B1 = x2;
                        x2 = x1;
                        f2 = f1;
                        t1 = true;      t2 = false;
                }
                else if (f1==f2)
                {
                        A1= x1;
                        B1 = x2;
                        t1 = true;      t2 = true;
                }
                else
                {
                        A1 = x1;
                        x1 = x2;
                        f1 = f2;
                        t1 = false;     t2 = true;
                }
        }
        while (B1-A1>crit0);
        //      this->go(A1,dir);
        double f = this->getF();
        if (f>f0)
        {
                std::cout << "f>f0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
                this->go(0,dir);
        }
}
//---------------------------------------------------------------------------
void a_fn::minGold(const a_mat_c& dir)
{
        double f0 = this->getF();
        this->setx0();
        int k = -1;
        int iter;
        double di0 = 0.1;
        double i,di;
        double prov = 0;
        do
        {
                k++;
                //              this->go(0,dir);
                iter = 0;
                double f0;
                double f = this->gF(0,dir);
                i = 0;
                di = di0*exp(k*log(0.1));
                //              std::cout << "di: " << di << std::endl;
                do
                {
                        //                      std::cout << "ok" << std::endl;
                        iter++;
                        di *= 2.;
                        i += di;
                        f0 = f;
                        f = this->gF(i,dir);
                }
                while (f<=f0);
                if (iter ==1 ) prov++;
        }
        while (iter<4);
        double A1 = i-di-di/2.;
        double B1 = i;
        std::cout << "A1: " << A1 << "    B1: " << B1 << std::endl;
        this->go(A1,dir);
        double f = this->getF();
        if (f>f0)
        {
                std::cout << "minGold----------------------------" << std::endl;
                //              std::cout << "f-f0: " << f-f0 << std::endl;
                //              std::cout << "A1: " << A1 << "    B1: " << B1 << std::endl;
        }
        f0 = f;
        this->minGold0(dir,0,B1-A1);
        f = this->getF();
        if (f>f0)
                std::cout << "minGold0++++++++++++++++++++++++++++" << std::endl;
}
//---------------------------------------------------------------------------
void a_fn::mingrad()
{
        auto dir = new a_mat_c(nv_);
        do
        {
                this->getGrad(*dir);
                dir->normalise();
                *dir *= -1;
                this->minGold(*dir);
                std::cout << "f: " << this->getF() << std::endl;
        }
        while (!this->minver());
        delete dir;
}
//---------------------------------------------------------------------------
void a_fn::minsyst()
{
        auto dir = new a_mat_c(nv_);
        do
        {
                for (int j = 0; j < nv_; j++)
                {
                        if (this->getdF(j)>0)
                                (*dir)(j) = -1;
                        else
                                (*dir)(j) = 1;
                        this->minGold(*dir);
                        (*dir)(j) = 0;
                }
        }
        while (!this->minver());
        delete dir;
}
//---------------------------------------------------------------------------
bool a_fn::minver()
{
        auto dir = new a_mat_c(nv_);
        this->getGrad(*dir);
        double test = sqrt(dir->sumsq()/nv_);
        delete dir;
        return (test<gmin_);
}
//***************************************************************************
//---------------------------------------------------------------------------
class a_mat_sq_DVP : public a_mat_sq
{
        public:
                a_mat_sq_DVP(UI  nv);
};
//---------------------------------------------------------------------------
a_mat_sq_DVP::a_mat_sq_DVP(UI  nv) : a_mat_sq(nv)
{
        for (int i = 0; i < nv; i++) 
                (*this)(i,i) = 1;
}
//***************************************************************************
//---------------------------------------------------------------------------
void a_fn::minDFP()
{
        auto H = new a_mat_sq_DVP(nv_);
        a_mat_c S(nv_);
        a_mat_c g1(nv_);
        a_mat_c g2(nv_);
        a_mat_sq Mi(nv_);
        a_mat_sq Ni(nv_);
        a_mat_r mr(nv_);
        a_mat_sq m0(nv_);
        a_mat_c x0(nv_);
        a_mat_c x1(nv_);
 
        bool stop = false;
        int iter = 0;
        double f = this->getF();
        do
        {
                double f0 = f;
                iter++;
                //              if (iter==500)
                //                      this->showX(std::cout);
 
                this->getGrad(g1);
                S = -(*H)*g1;
                x0 = *x_;
                this->minGold(S);
                f = this->getF();
                //              if (f>f0)
                //              {
                //                      std::cout << "minGold---------------" << std::endl;
                //                      this->showX(std::cout);
                //              }
                std::cout << "i" << iter << "  f: " << f << std::endl;
                stop = this->minver();
                if (!stop)
                {
                        x1 = (*x_)-x0;
                        double d1 = sqrt(x1.sumsq()/S.sumsq());
                        this->getGrad(g2);
                        g1 = g2-g1;
                        mr = S.transpose();
                        Mi = S*mr;
                        Mi *= d1;
                        m0 = g1*mr;
                        d1 = m0(1,1);
                        if ((d1 != 0)&&(iter % 10))
                        {
                                Mi *= 1/d1;
                                //calcul de Ni
                                mr = g1.transpose();
                                g2 = (*H)*g1;
                                m0 = g2*mr;
                                d1 = m0(1,1);
                                mr = g2.transpose();
                                Ni = g2*mr;
                                if (d1 != 0)
                                        Ni *= -1/d1;
                                else
                                        stop = true;
                                //nouveau H
                                (*H) += Mi+Ni;
                        }
                        else
                        {
                                delete H;
                                H = new a_mat_sq_DVP(nv_);
                        }
                }
        }
        while (!stop);
        delete H;
}