#include"stdio.h"
#include <iostream.h>
#include<math.h>
#define MAXLENGTH 10

//向量定義
typedef struct{             
    int tag;                //行列向量標志。行向量為0，列向量為1。
    int dimension;          //向量的維數
    double elem[MAXLENGTH];    //向量的元素
}vector;

vector vecCreat(int tag, int n){
//建立維數為n的向量
    vector x;
    x.tag=tag;
    x.dimension=n;
    for(int i=0;i<n;++i){
        cout<<"請輸入第"<<i+1<<"分量:";
        cin>>x.elem[i];
    }
    return x;
}

double vecMultiply(vector a, vector b){
//行向量與列向量乘法運算
    if((a.tag!=b.tag)&&(a.dimension==b.dimension)){//相乘條件
        double c=0;
        for(int i=0;i<a.dimension;++i)
            c+=a.elem[i]*b.elem[i];
        return c;
    }
}

vector vecAdd(vector a, vector b){
//向量加法運算
    if((a.tag==b.tag)&&(a.dimension==b.dimension)){//相加條件
        vector c;
        c.dimension=a.dimension;
        c.tag=a.tag;
        for(int i=0;i<c.dimension;++i)
            c.elem[i]=a.elem[i]+b.elem[i];
        return c;
    }
}

vector vecConvert(vector a){
//向量轉置
    if(a.tag==0)a.tag=1;
    if(a.tag==1)a.tag=0;
    return a;
}

double vecMole(vector a){
//求向量模
    double sum=0;
    for(int i=0;i<a.dimension;++i)
        sum+=a.elem[i]*a.elem[i];
    sum=sqrt(sum);    
    return sum;
}

vector numMultiply(double n,vector a){
//數乘向量
    for(int i=0;i<a.dimension;++i)
        a.elem[i]=n*a.elem[i];
    return a;
}

void showPoint(vector x, double f){
//顯示點，解.
    cout<<"--- x=(";
    for(int i=0;i<x.dimension;++i){
        cout<<x.elem[i];
        if(i!=x.dimension-1)cout<<",";
    }
    cout<<") --- f(x)="<<f<<endl;
    cout<<endl;    
}
/////////////////////////////////////////////////////////
////////////// function.h頭文件 /////////////////////////
//////// 函數及其導數的設置均在此文件完成////////////////
/////////////////////////////////////////////////////////
// * 無 約 束 問 題 *                                  //
// 目標函數--在vecFun(vector vx)中修改,x1改成x[1]      //
//              x2改成x[2],依此類推...                 //
/////////////////////////////////////////////////////////


//#include "vector.h"
#define SIZE 10
#define MAX 1e300

double min(double a, double b){
//求兩個數最小
    return a<b?a:b;
}

double max(double a, double b){
//求兩個數最大
    return a>b?a:b;
}

vector vecGrad(double (*pf)(vector x),vector pointx){
//求解向量的梯度
//采用理查德外推算法計算函數pf在pointx點的梯度grad;
    vector grad;
    grad.tag=1; grad.dimension=pointx.dimension;//初始化偏導向量
    vector tempPnt1,tempPnt2;           //臨時向量
    double h=1e-3;
    for(int i=0;i<pointx.dimension;++i){
        tempPnt1=tempPnt2=pointx;
        tempPnt1.elem[i]+=0.5*h;
        tempPnt2.elem[i]-=0.5*h;
        grad.elem[i]=(4*(pf(tempPnt1)-pf(tempPnt2)))/(3*h);
        tempPnt1.elem[i]+=0.5*h;
        tempPnt2.elem[i]-=0.5*h;
        grad.elem[i]-=(pf(tempPnt1)-pf(tempPnt2))/(6*h);
    }
    return grad;
}

double vecFun(vector vx){
//最優目標多元函數
    double x[SIZE];
    for(int i=0;i<vx.dimension;++i)
        x[i+1]=vx.elem[i];
    //----------約束目標函數--------------
    //return x[1]*x[1]+x[2]*x[2];
    //----------無約束正定函數--------------
    return x[1]*x[1]+4*x[2]*x[2]+9*x[3]*x[3]-2*x[1]+18*x[3];//例一
    //----------無約束非二次函數--------------
    //return (1-x[1])*(1-x[1])+(1-x[4])*(1-x[4])+(x[1]*x[1]-x[2])*(x[1]*x[1]-x[2])+(x[2]*x[2]-x[3])*(x[2]*x[2]-x[3])+(x[3]*x[3]-x[4])*(x[3]*x[3]-x[4]);//例二

}

double vecFun_Si(vector vx){
//不等式約束函數
    double x[SIZE];
    for(int i=0;i<vx.dimension;++i)
        x[i+1]=vx.elem[i];
    return x[1]+1;//不等式約束函數
}

double vecFun_Hi(vector vx){
//等式約束函數
    double x[SIZE];
    for(int i=0;i<vx.dimension;++i)
        x[i+1]=vx.elem[i];
    return x[1]+x[2]-2;//等式約束函數
}

double vecFun_S(vector x,double v,double l,double u){
//約束問題轉化為無約束問題的增廣目標函數F(x,v,l,u)
    double sum1=0,sum2=0,sum3=0;//分別定義三項的和
    //sum1
    double temp=max(0,v-2*u*vecFun_Si(x));
    sum1=(temp*temp-v*v)/(4*u);
    //sum2
    sum2=l*vecFun_Hi(x);
    //sum3
    sum3=u*vecFun_Hi(x)*vecFun_Hi(x);    
    //F(x,v,l,u)=f(x)+sum1-sum2+sum3
    return vecFun(x)+sum1-sum2+sum3;
}

vector vecFunD_S(vector x,double v,double l,double u){//利用重載函數實現目標函數F(x,v,l,u)的導數
//約束問題轉化為無約束問題的增廣目標函數F(x,v,l,u)的導函數
    vector sum1,sum2,sum3;//分別定義三項導數的和
    //sum1
    sum1.dimension=x.dimension;    sum1.tag=1;
    sum2.dimension=x.dimension;    sum2.tag=1;
    sum3.dimension=x.dimension;    sum3.tag=1;
    double temp=max(0,v-2*u*vecFun_Si(x));
    if(temp==0){
        for(int i=0;i<sum1.dimension;++i)
            sum1.elem[i]=0;
    }
    else{
        sum1=numMultiply(-(v-2*u*vecFun_Si(x)),vecGrad(vecFun_Si,x));
    }
    //-sum2
    sum2=numMultiply(-l,vecGrad(vecFun_Hi,x));
    //sum3
    sum3=numMultiply(2*u*vecFun_Hi(x),vecGrad(vecFun_Hi,x));    
    //F=f(x)+sum1-sum2+sum3
    return vecAdd(vecAdd(vecGrad(vecFun,x),sum1),vecAdd(sum2,sum3));
}
///////////////////////////////////////////////////
/////            nonrestrict.h頭文件          /////
/////  包含無約束問題的求解函數: 直線搜索     /////
/////         共軛梯度法,H終止準則            /////
///////////////////////////////////////////////////
//include "restrict.h"

vector lineSearch(vector x, vector p ,double t){
//從點x沿直線p方向對目標函數f(x)作直線搜索得到極小點x2,t為初始步長。
    vector x2;
    double l=0.1,n=0.4;         //條件1和2的參數  
    double a=0,b=MAX;           //區間 
    double f1=0,f2=0,g1=0,g2=0;
    int i=0;                    //迭代次數
    do{
        if(f2-f1>l*t*g1){b=t;t=(a+b)/2;}         //改變b,t循環
        do{
            if(g2<n*g1){a=t;t=min((a+b)/2,2*a);} //改變a,t循環
            f1=vecFun(x);                                   //f1=f(x)
            g1=vecMultiply(vecConvert(vecGrad(vecFun,x)),p);//g1=∨f(x)*p
            x2=vecAdd(numMultiply(t,p),x);                  //x2=x+t*p
            if(i++ && vecFun(x2)==f2){return x2;}           //【直線搜索進入無限跌代，則此次跌代結束。返回當前最優點】
            f2=vecFun(x2);                                  //f2=f(x2)
            g2=vecMultiply(vecConvert(vecGrad(vecFun,x2)),p);      //g2=∨f(x2)*p    
            //cout<<"("<<x2.elem[0]<<","<<x2.elem[1]<<","<<x2.elem[2]<<","<<x2.elem[3]<<");"; //輸出本次結果
            //cout<<"t="<<t<<";";
            //cout<<"f(x)="<<f2<<endl;
        }
        while(g2<n*g1);      //不滿足條件ii,則改變a,t循環
    }
    while(f2-f1>l*t*g1);     //不滿足條件i,則改變b,t循環
    return x2;
}


void Fletcher_Reeves(vector xx,vector &minx, double &minf){
//Fletcher-Reeves共軛梯度法.
//給定初始點xx.對vecFun函數求得最優點x和最優解f,分別用minx和minf返回
    int k=0,j=0;       //迭代次數,j為總迭代次數
    double c=10e-1;    //終止限c
    int n=xx.dimension;//問題的維數
    double f0,f,a;     //函數值f(x),a為使p方向向量共軛的系數
    vector g0,g;       //梯度g0,g
    vector p0,p;       //搜索方向P0,p
    vector x,x0;       //最優點和初始點
    double t=1;        //直線搜索的初始步長=1
    x=xx;                         //x0
    f=vecFun(x);                  //f0=f(x0)
    g=vecGrad(vecFun,x);          //g0
    p0=numMultiply(-1,g);         //p0=-g0,初始搜索方向為負梯度方向
    do{
        x0=x;f0=f;g0=g;              
        x=lineSearch(x0,p0,t);       //從點x出發,沿p方向作直線搜索
        f=vecFun(x);                 //f=f(x)
        g=vecGrad(vecFun,x);                //g=g(x)
        if(Himmulblau(x0,x,f0,f)==1){//滿足H終止準則,返回最優點及最優解。
            cout<<"* 總迭代次數："<<j<<" ----"<<endl;
            minx=x; minf=f;
            break;
        }
        else{                        //不滿足H準則
            cout<<"* 第"<<j<<"次迭代";
            showPoint(x,f);          //顯示中間結果x,f    

            if(k==n){                //若進行了n+1次迭代
                k=0;                 //重置k=0,p0=-g
                p0=numMultiply(-1,g);
                t=vecMole(vecAdd(x,numMultiply(-1,x0)))/vecMole(p0);//下次直線搜索步長
                continue;            //進入下一次循環,重新直線搜索
            }                        
            else{                     //否則,計算共軛向量p
                a=vecMole(g)*vecMole(g)/(vecMole(g0)*vecMole(g0));
                p=vecAdd(numMultiply(-1,g),numMultiply(a,p0));
            }
        }
        if(fabs(vecMultiply(vecConvert(p),g))<=c){//共軛梯度方向下降或上升量很小
            p0=numMultiply(-1,g);//p0=-g,k=0
            k=0;
        }
        else if(vecMultiply(vecConvert(p),g)<-c){//共軛梯度方向下降并且下降量保證
            p0=p;                //p0=p,k=k+1
            ++k;
        }
        else{                                    //共軛梯度方向上升并且下降量保證
            p0=numMultiply(-1,p);//p0=-p,k=k+1
            ++k;
        }
        t=vecMole(vecAdd(x,numMultiply(-1,x0)))/vecMole(p0);//下次直線搜索步長

    }
    while(++j);
}

///////////////////////////////////////////////////
/////             main.h頭文件                /////
/////  主程序文件，控制整個程序的流程         /////
///////////////////////////////////////////////////
//#include "nonrestrict.h"

void printSelect(){
    cout<<"************** 最優化方法 ***************"<<endl;
    cout<<"*****************************************"<<endl;
    cout<<"***  請選擇：                         ***"<<endl;
    cout<<"***  1. 無約束最小化問題(共軛梯度法)  ***"<<endl;
    cout<<"***  3. 退出(exit)                    ***"<<endl;
    cout<<"*****************************************"<<endl;
}

void sub1(){
    char key;
    cout<<"--------------------共軛梯度法求無約束最小化問題----------------"<<endl;
    cout<<"步驟:"<<endl;
    cout<<"◆ 1. 輸入f(x)及其導數.(參考function.h文件說明)"<<endl;
    cout<<"-----確認是否已輸入<目標函數>及其<導數>? (Y/N)";
    cin>>key;
    if(key=='N' || key=='n')return; 
    else if(key=='Y' || key=='y')
	{
        vector x0; //起始點
        int m;
        cout<<"◆ 2. 起始點X0初始化"<<endl<<"-----請輸入X0的維數:";
        cin>>m;
        x0=vecCreat(1,m);
        vector minx;//求得的極小點
        double minf;//求得的極小值
        Fletcher_Reeves(x0,minx,minf);
        cout<<"◆ 最優解及最優值:";
        showPoint(minx,minf);
    }
}



void main(){
    int mark=1;
    while(mark){
        printSelect();
        int sel;
        cin>>sel;
        switch(sel){
            case 1:
                sub1();
                break;
        
            case 3:
                mark=0;
                break;
        };
    }
}