Phương pháp đơn hình 2 pha Code C++ (Bài toán quy hoạch tuyến tính)

Solve:

Để biết thêm về phương pháp đơn hình: Link 

Tài liệu của thầy Nguyễn Văn Hiệu: Link 

Code C++11:

#include<bits/stdc++.h>

using namespace std;

vector<int> v;

int soAnThiet;

int pha= 0;

 

vector<int> HeSo(vector<int> c, vector< vector<double>> a){

    vector<int> cb;

    set<pair<int, int>> sp;             //dung de tim he so theo thu tu hang 1 -> hang m

   

   

    for(int j=0; j<a[0].size(); j++){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii>=0) sp.insert({ii, j});      //neu la 1 cot trong ma tran don vi

    }

   

    for(auto p : sp)

        cb.push_back(c[p.second]);

    return cb;

}

vector<double> TinhDelta(vector<int> cb, vector<int> c, vector< vector<double>> a){

    //deltai = sum(aij * cbi) - cj

    vector<double> v;

   

    for(int j=0; j<a[0].size(); j++){

        double s= 0;

        for(int i=0; i<a.size(); i++)

            s+= a[i][j] * cb[i];

        v.push_back(s-c[j]);

    }

    return v;

}

double TinhFx(vector<int> cb, vector<double> b){

    double s= 0;

    int i= 0;

    for(auto val_b: b)

        s += val_b * cb[i++];

    return s;

}

bool VoNghiem(vector<double>delta, vector< vector<double>> a){

    //neu ton tai deltaj > 0 ma aij <= 0 (voi moi i=1:m) thi vn

   

    for(int j=0; j<delta.size(); j++)

    if (delta[j] > 0){

       

        bool check= true;

        for(int i=0; i<a.size(); i++)

            if (a[i][j] > 0) check= false;

       

        if (check) return true;

    }

    return false;

}

bool PATU(vector<double>delta){

    for(auto val:delta)

        if (val > 0) return false;

    return true;

}

vector<double> Nghiem(vector< vector<double>> a, vector<double> b){

    vector<double> x;

   

    for(int j=0; j<a[0].size(); j++){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii>=0) x.push_back(b[ii]);  //neu la 1 cot trong ma tran don vi

            else x.push_back(0);        

    }

   

    return x;

}

int TimCotVao(vector<double>delta){

    double valMax= delta[0];

    int index= 0;

   

    for(int i=1; i<delta.size(); i++){

        if (valMax < delta[i]){

            valMax= delta[i];

            index= i;

        }

    }

    return index;

}

int TimCotRa(int jVao, vector<vector<double>>a, vector<double> b, int &iVao){

    double valMin= (a[0][jVao]<=0 ? 100111000 : b[0]*1.0 / a[0][jVao]);

   

    for(int i=1; i<a.size(); i++){

       

        if (a[i][jVao] <= 0) continue;

        double tiSoKoAm= b[i]*1.0 / a[i][jVao];  

       

       

        if (valMin > tiSoKoAm){

           

            valMin= tiSoKoAm;

            iVao= i;

        }

    }

   

    for(int j=0; j<a[0].size(); j++)

    if (j != jVao){

               

        int ii= -1;

        for(int i=0; i<a.size(); i++){

            if (a[i][j] == 1 && ii == -1){

                ii= i;  

            }

            else

                if (a[i][j]) ii= -2;    //co the them break de giam dpt

        }

        if (ii == iVao) return j;  

    }

    return -1;      //error

}

void output(vector<int> c, vector<int> cb, vector<vector<double>> a, vector<double> delta, double fx, vector<double> b){

    cout<<"\n\n-----------------------BANG DON HINH--------------------------------\n\n";

    cout<<setw(10)<<""; for(auto val: c) cout<<setw(10)<<val; cout<<setw(10)<<" <-- [c]";

    cout<<"\n";

    int i= 0;

    for(auto v : a){

        cout<<"\n";

       

        cout<<setw(10)<<cb[i];

       

        for(auto val : v)

            cout<<setw(10)<<fixed<<setprecision(3)<<val;

       

        cout<<setw(10)<<b[i++];

       

    }

    cout<<"\n";

    cout<<setw(10)<<"[cb]";

    cout<<"\n";

    cout <<setw(10)<<"[delta"; for(auto val:delta) cout<<setw(10)<<fixed<<setprecision(3)<<val; cout<<"]"<<setw(10)<<fixed<<setprecision(3)<<fx<<"=fx\n";

    cout<<"\n\n-------------------------------------------------------------------------\n\n";

}

vector< vector<double> > XoaCotGia(vector< vector<double> > a){

    vector< vector<double> > res;

   

    for(auto v:a){

        vector<double> vd;

        for(int i=0; i<soAnThiet; i++)

            vd.push_back(v[i]);

        res.push_back(vd);

    }

    return res;

}

vector<double> process(int n, vector<int> c, int m, vector< vector<double>> a, vector<double> b){

    cout<<"\n\n\n\n@@@@@@@@@@@@@@@ PHA "<<++pha<<" @@@@@@@@@@@@@@@\n\n\n\n";

       

    vector<int> cb= HeSo(c, a);

    vector<double> delta = TinhDelta(cb, c, a);

    double fx = TinhFx(cb, b);

   

    int cnt= 0;

    while(true){            output(c, cb, a, delta, fx, b);

        if (VoNghiem(delta, a)){

            cout<<"Bai toan vo nghiem!";

            exit(0);

        }

        if (PATU(delta)){

            if  (pha == 1){

                vector<double> nghiem = Nghiem(a, b);

               

                double s= 0;

                for(int i=soAnThiet; i<nghiem.size(); i++)

                    s += nghiem[i];

                if (s > 0.00000001){

                    cout<<"Vo nghiem do g(t*) > 0";

                    exit(0);

                }

            }

            return (pha == 1 ? process(soAnThiet, v, m, XoaCotGia(a), b) : Nghiem(a, b));

        }

        int jv= TimCotVao(delta);

        int iv= 0;

        int jr= TimCotRa(jv, a, b, iv);

       

        //bien doi ma tran a va vector b

        double soNhan= 1 / a[iv][jv];

        for(int j=0; j<a[0].size(); j++) a[iv][j] *= soNhan;

        b[iv] *= soNhan;

       

        for(int i=0; i<a.size(); i++)

        if (i != iv){

           

            soNhan = -a[i][jv] / a[iv][jv];    

            for(int j=0; j<a[i].size(); j++)

                a[i][j] += a[iv][j]*soNhan;

           

            b[i] += b[iv] * soNhan;

        }

       

        //bien doi vector delta

        soNhan= -delta[jv] / a[iv][jv];

        for(int j= 0; j<a[0].size(); j++)

            delta[j] += a[iv][j] * soNhan;

        fx += b[iv] * soNhan;

       

       

    }

    return vector<double>();                                                                                              

}  

       

   

void input(){

    freopen("donhinh2pha.txt","r",stdin);

    //cout<<"nhap n (so an x (co so + ko co so)) : ";

    int n;  cin >> n;

   

    vector<int> c(n);  

    for(int i=0; i<n; i++){

        //cout<<"nhap c"<<i<<" : ";

        cin >> c[i];    

    }

   

    //cout<<"nhap m (so an x (co so)): ";

    int m;  cin >> m;

   

    vector< vector<double> > a(m, vector<double>(n));

    vector<double> b(m);

   

    for(int i=0; i<a.size(); i++){

        for(int j=0; j<a[i].size(); j++){

       

            //cout<<"nhap a"<<i<<j<<" : ";

            cin >> a[i][j];

        }

        //cout<<" sum(Aij * xj) = bi, nhap b"<<i<<" : ";

        cin >> b[i];

    }

   

    //cout<<"Nhap so an thiet:";

    cin >> soAnThiet;

    for(int i=0; i<soAnThiet; i++){

        int x;

        cin >> x;

        v.push_back(x);

    }

    vector<double> x= process(n, c, m, a, b);

    cout<<"VAY NGHIEM TOI UU LAN LUOT LA: \n";

    for(int i=0; i<x.size(); i++)

        cout<<"x"<<i<<" = "<<x[i]<<"\n";

}

int main(){

   

    input();

   

}

Test: (bỏ freopen nếu nhập từ console)

5

0 0 0 1 1

2

4 3 4 1 0 4

4 1 6 0 1 5

3

-4 3 1

(don hinh 2 pha vi du cua thay)

5

0 0 0 1 1

2

4 3 4 1 0 4

4 1 -3 0 1 5

3

-4 3 1 

(vi du 3 bai toan m - vo nghiem)

5

0 0 0 1 1

2

4 3 4 1 0 4

4 1 -3 0 1 4

3

-4 3 1

(vi du 2 bai toan m)

7

0 0 0 0 0 1 1

3

-1 4 -2 1 0 0 0 6

1 1 2 0 -1 1 0 6  

2 -1 2 0 0 0 1 4 

5

-1 -2 1 0 0

(vi du 1 bai toan m)

7

0 0 0 0 0 1 1

3

1 -4 2 -5 9 0 0 3

0 1 -3 4 -5 1 0 6

0 1 -1 1 -1 0 1 1

5

2 6 -5 1 4

(don hinh 2 pha tren https://www.youtube.com/watch?v=S0nh0OZos6c)

Gradient descent cho hàm nhiều biến - Code Python

 Problem:

Bài toán thực tế:

Một căn nhà rộng $x_1{\text{m}}^2$, có $x_2$ phòng ngủ và cách trung tâm thành phố $x_3\text{km}$ có giá là bao nhiêu. Giả sử chúng ta đã có số liệu thống kê từ 1000 căn nhà trong thành phố đó, liệu rằng khi có một căn nhà mới với các thông số về diện tích, số phòng ngủ và khoảng cách tới trung tâm, chúng ta có thể dự đoán được giá của căn nhà đó không?

Diễn dãi theo code:

h = so hang

c = so cot

Input:

    h, c

    x[1][1] ... x[h][c]

y[1] ... y[h]

    x[h+1][1] ... x[h+1][c]

Ouput:

    y[h+1]

Solve:

Ta co

θ[1] * x[1][1] + ...  + θ[c] * x[1][c] ~ y[1]

...

...

θ[1] * x[h][1] + ...  + θ[c] * x[h][c] ~ y[h]

=> θ[1]* x[h+1][1] + ...  + θ[c]*x[h+1][c] ~ y[h+1]

=> Can tim cac θ de giai quyet bai toan

Ví dụ đơn giản:

Input: Cho tập điểm (x, y)

Ouput: Tìm 1 đường thằng gần tất cả các điểm nhất (tìm a, b của đường thằng y= ax + b)

Solve:

B1: Dự đoán θ0, θ1 (~ a, b)

B2: Cập nhập θ cho đến khi đạt điều kiện dừng (thường  abs(dHMM) < epsilon, loop = ...)

              θ_new = θ_old  - alpha * dHMM

  Cách tìm θi của bài toán trên 

Code Python:

import numpy as np

import matplotlib.pyplot as plt

np.random.seed(2)

def grad(w):

    N = Xbar.shape[0]

    return 1/N * Xbar.T.dot(Xbar.dot(w) - y)

def l(w):

    N = Xbar.shape[0]

    return .5/N*np.linalg.norm(Xbar.dot(w)-y, 2)**2

def myGradientDescent(w_init, grad, alpha, loop = 1000, esilon = 1e-4):

    w = [w_init]

    for i in range(loop):

        w_new = w[-1] - alpha*grad(w[-1])

        if np.linalg.norm(grad(w_new))/len(w_new) < esilon:

            break

        w.append(w_new)

    return (w, i)

if __name__ == '__main__':

    # dataset

    X = np.random.rand(1000, 1)

    y = 4 + 3 * X + .2 * np.random.randn(1000, 1)  # noise added

    # Building Xbar

    one = np.ones((X.shape[0], 1))

    Xbar = np.concatenate((one, X), axis=1)

    # Gradient descent

    w_init = np.array([[2], [1]])

    (w1, it1) = myGradientDescent(w_init, grad, 1)

    print(' Phương pháp GradientDescent: w = ', w1[-1].T, ',\n after %d iterations.' %(it1+1),',\n l = %f ' % l(w1[-1]))

    # Linear regression

    A = np.dot(Xbar.T, Xbar)

    b = np.dot(Xbar.T, y)

    w_lr = np.dot(np.linalg.pinv(A), b)

    print('Phuong phap nghich dao: w = ', w_lr.T)

    # Display result

    w = w_lr

    w_0 = w[0][0]

    w_1 = w[1][0]

    x0 = np.linspace(0, 1, 2, endpoint=True)

    y0 = w_0 + w_1 * x0

   

    # Draw the fitting line

    plt.plot(X.T, y.T, 'b.')  # data

    plt.plot(x0, y0, 'y', linewidth=2)  # the fitting line

    plt.axis([0, 1, 0, 10])

    plt.show()

Ouput:

Các nguồn tham khảo: Link1, Link2

Xem thêm bài viết cùng chủ đề: Gradient Descent cho hàm 1 biến