/*
The BackPropagation procedure included in the second
volume "Finance et Econometrie" of the book
"Introduction a la programmation sous GAUSS" doesn't
work with the Windows 95/NT version of GAUSS. There
are two problems. First, the procedure uses the medit
procedure of the DOS version. Second, it uses
special DOS characters.

The following BackPropagation procedure takes account of these
problems. If you have a windows 95/NT version of
GAUSS, you may replace the backp.src file with
the following backp.src file.

Th. Roncalli   Thierry.Roncalli@montesquieu.u-bordeaux.fr
28 Nov 1996

---> Backp.src
*/


proc _BackP_Fonction(x);
    local y;
    if _BackP_Options[1] == 3;
        y = tanh(x);
    elseif _BackP_Options[1] == 2;
        y = (1-exp(-x))./(1+exp(-x));
    else;
        _BackP_Options[1] = 1;
        y = 1./(1+exp(-x));
    endif;
    retp(y);
endp;

proc (4) = _BackP_Feedforward(x_n,a,b);
    local u_n,y_n,v_n,z_n;

    u_n = a'(1Ýx_n);
    y_n = _BackP_Fonction(u_n);

    if _BackP_Options[4] == 1;
        v_n = b'(1Ýy_nÝx_n);
    else;
        v_n = b'(1Ýy_n);
    endif;
    z_n = _BackP_Fonction(v_n);

    retp(u_n,y_n,v_n,z_n);
endp;

proc (1) = _BackP_Learning_Laws_1(epsilon,d_m);
    local c_m;
    c_m = -epsilon*d_m;
    retp(c_m);
endp;

proc (1) = _BackP_Learning_Laws_2(epsilon,mu,d_m,c_mm);
    local c_m;
    c_m = mu*c_mm-(1-mu)*epsilon*d_m;
    retp(c_m);
endp;

proc (1) = _BackP_Learning_Laws_3(kappa,phi,theta,d_m,f_mm,e_mm);
    local g_m,e_m,c_m,f_m;
    g_m = (d_m.*f_mm).>0;
    e_m = g_m.*(e_mm+kappa)+(1-g_m).*(phi*e_mm);
    c_m = -e_m.*d_m;
    f_m = theta*f_mm+(1-theta)*d_m;
    retp(c_m);
endp;

proc (1) = _BackP_Learning_Laws_4(kappa,phi,theta,mu,d_m,f_mm,e_mm,c_mm);
    local g_m,e_m,c_m,f_m;
    g_m = (d_m.*f_mm).>0;
    e_m = g_m.*(e_mm+kappa)+(1-g_m).*(phi*e_mm);
    c_m = mu*c_mm-(1-mu)*e_m.*D_m;
    f_m = theta*f_mm+(1-theta)*d_m;
    retp(c_m);
endp;

proc _BackP_Deriv_Fonction(x);
    local y;
    if _BackP_Options[1] == 3;
        y = (2./(exp(x)+exp(-x)))^2;
    elseif _BackP_Options[1] == 2;
        y = 0.5*(1-x).*(1+x);
    else;
        _BackP_Options[1] = 1;
        y = x.*(1-x);
    endif;
    retp(y);
endp;

proc (2) = _BackP_Error_Rules(x_n,a,b,u_n,y_n,v_n,z_n,t_n);
    local Deriv_z_v,Deriv_y_u;
    local Deriv_E_z,p_n,Db_n;
    local Deriv_E_y,q_n,Da_n;

    if _Backp_Options[1] == 3;
        Deriv_z_v = _BackP_Deriv_Fonction(v_n);
        Deriv_y_u = _BackP_Deriv_Fonction(u_n);
    elseif _Backp_Options[1] == 2;
        Deriv_z_v = _BackP_Deriv_Fonction(z_n);
        Deriv_y_u = _BackP_Deriv_Fonction(u_n);
    else;
        Deriv_z_v = _BackP_Deriv_Fonction(z_n);
        Deriv_y_u = _BackP_Deriv_Fonction(y_n);
    endif;

    Deriv_E_z = z_n-t_n;
    p_n = Deriv_E_z.*Deriv_z_v;
    if _BackP_Options[4] == 1;
        Db_n = (1Ýy_nÝx_n)*p_n';
    else;
        Db_n = (1Ýy_n)*p_n';
    endif;

    if _BackP_Options[4] == 1;
        Deriv_E_y = trimr(b,1,rows(x_n))*p_n;
    else;
        Deriv_E_y = trimr(b,1,0)*p_n;
    endif;
    q_n = Deriv_E_y.*Deriv_y_u;
    Da_n = (1Ýx_n)*q_n';

    retp(Da_n,Db_n);
endp;


proc (3) = _BackP_Network(X,T,a,b,NN);
    local n,x_n,t_n,Da,Db,u_n,y_n,v_n,z_n,Da_n,Db_n,r;

    Da = 0; Db = 0;
    r=zeros(rows(T),NN);

    n=1;
    do until n>NN;

         x_n = x[.,n];
         {u_n,y_n,v_n,z_n} = _BackP_Feedforward(x_n,a,b);
         r[.,n] = z_n;

         t_n = t[.,n];
         {Da_n,Db_n} =
         _BackP_Error_Rules(x_n,a,b,u_n,y_n,v_n,z_n,t_n);

         Da = Da+Da_n;
         Db = Db+Db_n;

         n=n+1;
    endo;

    retp(Da,Db,r);
endp;


proc (23) = _BackP_Initialisation(I,J,K);
    local a,b,mat;
    local da_m,da_mm,fa_m,fa_mm,ca_m,ca_mm,ea_m,ea_mm;
    local db_m,db_mm,fb_m,fb_mm,cb_m,cb_mm,eb_m,eb_mm;
    local phi,theta,mu,kappa,epsilon;

    if _BackP_a /= 0;
        a = _Backp_a;
    else;
        a = .2*(rndu(I+1,J)-.5);
    endif;

    if _BackP_b /= 0;
        b = _BackP_b;
    else;
        if _BackP_Options[4] == 1;
            b = ( seqa(0,1,1+J+I) % 2 ).==0;
            b = b.*ones(1+J+I,K)-(1-b).*ones(1+J+I,K);
        else;
            b = ( seqa(0,1,1+J) % 2 ).==0;
            b = b.*ones(1+J,K)-(1-b).*ones(1+J,K);
        endif;
    endif;

    phi=_BackP_Parametres[1];
    theta=_BackP_Parametres[2];
    mu=_BackP_Parametres[3];
    kappa=_BackP_Parametres[4];
    epsilon=_BackP_Parametres[5];

    mat = zeros(1+I,J);
    da_m = mat; da_mm = mat; fa_m = mat; fa_mm = mat; ca_m = mat; ca_mm = mat;
    ea_m = kappa*ones(1+I,J); ea_mm = ea_m;

    if _BackP_Options[4] == 1;
        mat = zeros(1+J+I,K);
        eb_m = kappa*ones(1+J+I,K); eb_mm = eb_m;
    else;
        mat = zeros(1+J,K);
        eb_m = kappa*ones(1+J,K); eb_mm = eb_m;
    endif;
    db_m = mat; db_mm = mat; fb_m = mat; fb_mm = mat; cb_m = mat; cb_mm = mat;

    retp(a,b,
         da_m,da_mm,fa_m,fa_mm,ca_m,ca_mm,ea_m,ea_mm,
         db_m,db_mm,fb_m,fb_mm,cb_m,cb_mm,eb_m,eb_mm,
         phi,theta,mu,kappa,epsilon);
endp;

proc (1) = _BackP_Erreur(X,T,a,b,I,J,K);
    local message;

    if cols(X) /=  cols(T);
        ERRORLOG "Les dimensions des matrices X et T ne sont pas conformes";
        goto Erreur;
    endif;

    if ismiss(X) or ismiss(T);
        ERRORLOG "Les matrices X et T contiennent des donnees manquantes";
        goto Erreur;
    endif;

    if (rows(a)/=(1+I)) or (cols(a)/=J);
        ERRORLOG "Mauvaise dimension de la matrices A";
        message = "A doit etre de dimension "
                   $+ ftos(1+I,"%lf",3,0) $+ ftos(J," * %lf.",3,0);
        ERRORLOG message;
        goto Erreur;
    endif;

    if _BackP_Options[4] == 1;
        if (rows(b)/=1+J+I) or (cols(b)/=K);
            ERRORLOG "Mauvaise dimension de la matrices B";
            message = "B doit etre de dimension "$+ftos(1+J+I,"%lf",3,0)
                     $+ ftos(K," * %lf.",3,0);
            ERRORLOG message;
            goto Erreur;
        endif;
    else;
        if (rows(b)/=1+J) or (cols(b)/=K);
            ERRORLOG "Mauvaise dimension de la matrice B";
            message = "B doit etre de dimension "$+ftos(1+J,"%lf",3,0)
                       $+ftos(K," * %lf.",3,0);
            ERRORLOG message;
            goto Erreur;
         endif;
    endif;

    retp(0);

    Erreur:
        retp(1);
endp;

proc (3) = BackPropagation(X,T,J);
    local I,K,NN,r,m;
    local u_n,y_n,v_n,z_n,n,x_n;
    local a,b,da_m,da_mm,fa_m,fa_mm,ca_m,ca_mm,ea_m,ea_mm;
    local
db_m,db_mm,fb_m,fb_mm,cb_m,cb_mm,eb_m,eb_mm,erreur,phi,theta,mu,kappa,epsilon;
    local old_cursor,varnames,fh,test,choix,st;
    local E,epoque,EQM,LS_epsilon,LS_E,LS_a,LS_b,LS_r,l;

    if (J/=trunc(J)) or (J<0);
        ERRORLOG "Mauvaise initialisation du nombre de noeuds caches";
        retp(error(0),error(0),error(0));
    endif;

    /* Initialisation */
    I=rows(X); K=rows(T); NN=cols(X); r=zeros(K,NN);
    EQM = {};

    {a,b,da_m,da_mm,fa_m,fa_mm,ca_m,ca_mm,ea_m,ea_mm,
     db_m,db_mm,fb_m,fb_mm,cb_m,cb_mm,eb_m,eb_mm,phi,theta,mu,kappa,epsilon}
        = _BackP_Initialisation(I,J,K);

    old_cursor=csrtype(0);

    /* Gestion des erreurs */

    erreur = _BackP_Erreur(X,T,a,b,I,J,K);
    if erreur == 1;
        retp(error(0),error(0),error(0));
    endif;

    /* Creation du fichier de sauvegarde */

    if type(_Backp_Sauvegarde) == 13;
        varnames = { R } ;
        create fh = ^_BackP_sauvegarde with varnames,(NN*K),8;
    endif;

    if _print==1;
        gosub PRESENTATION;
    else;
        cls;
        _print=0;
    endif;


    m=1;
    do until m>_BackP_Options[3];

        /* Reseau */

        {da_m,db_m,r} = _BackP_Network(X,T,a,b,NN);

        /* Erreur Quadratique Moyenne */

        E = 0.5*sumc(sumc((R-T)^2))/(K*NN);
        EQM = EQMÝ(m~E);

        if type(_BackP_Sauvegarde) == 13;
          gosub SAUVEGARDE;
        endif;

        /* Modification des poids */

        if _BackP_Options[2] == 5;

            LS_E = zeros(rows(_BackP_Line_Search),1);

            l=1;
            do until l>rows(_BackP_Line_Search);

                LS_epsilon = _BackP_Line_Search[l];
                ca_m = _BackP_Learning_Laws_1(LS_epsilon,da_m);
                cb_m = _BackP_Learning_Laws_1(LS_epsilon,db_m);

                LS_a = a+ca_m;
                LS_b = b+cb_m;

                Ls_r=zeros(K,NN);

                n=1;
                do until n>NN;
                    x_n = x[.,n];
                    {u_n,y_n,v_n,z_n} = _BackP_Feedforward(x_n,a,b);
                    Ls_r[.,n] = z_n;
                    n=n+1;
                endo;

                LS_E[l] = 0.5*sumc(sumc((LS_r-t)^2))/(K*N);
                if _print == 1;
                    locate 22,60; print ftos(LS_epsilon," %lf ---",3,2);
                endif;

                l=l+1;

            endo;

            Ls_epsilon=_BackP_Line_Search[minindc(Ls_E)];
            ca_m = _BackP_Learning_Laws_1(LS_epsilon,da_m);
            cb_m = _BackP_Learning_Laws_1(LS_epsilon,db_m);

            if _print == 1;
                locate 22,70; print ftos(LS_epsilon," %lf",3,2);
             endif;

        elseif _BackP_Options[2] == 4;

            ca_m = _BackP_Learning_Laws_1(epsilon,da_m);
            cb_m = _BackP_Learning_Laws_1(epsilon,db_m);

        elseif _BackP_Options[2] == 3;

            ca_m = _BackP_Learning_Laws_2(epsilon,mu,da_m,ca_mm);
            cb_m = _BackP_Learning_Laws_2(epsilon,mu,db_m,cb_mm);

        elseif _BackP_Options[2] == 2;

            ca_m = _BackP_Learning_Laws_3(kappa,phi,theta,da_m,fa_mm,ea_mm);
            cb_m = _BackP_Learning_Laws_3(kappa,phi,theta,db_m,fb_mm,eb_mm);

        else;
            _BackP_Options[2] = 1;

            ca_m =
_BackP_Learning_Laws_4(kappa,phi,theta,mu,da_m,fa_mm,ea_mm,ca_mm);
            cb_m =
_BackP_Learning_Laws_4(kappa,phi,theta,mu,db_m,fb_mm,eb_mm,cb_mm);

        endif;

        a = a+ca_m;
        b = b+cb_m;

        da_mm = da_m; db_mm=db_m; ca_mm = ca_m; cb_m=cb_mm; ea_mm = ea_m;
eb_mm = eb_m;

        if _print == 1;
            locate 17,22; print "EPOQUE :";
            locate 17,50; print ftos(m,"%lf",10,0);
            locate 18,22; print "EQM :";
            locate 18,49; print ftos(E,"%le",10,5);
            locate 22,16;

            if _BackP_Options[2] == 5;
                locate 22,8;
                print "Learning Laws : Steepest Descent + Line Search ---";
            elseif _BackP_Options[2] == 4;
                print "        Learning Laws : Steepest Descent          ";
            elseif _BackP_Options[2] == 3;
                print "            Learning Laws : Momentum              ";
            elseif _BackP_Options[2] == 2;
                print "     Learning Laws : Adaptive Learning Rates      ";
            else;
                print "Learning Laws : Adaptive Learning Rates + Momentum";
            endif;

        endif;

        test=key;
        if test == 1025;
            _print=1-_print;
            test=0;
            if _print == 1;
                gosub PRESENTATION;
            else;
                cls;
            endif;
        endif;

        if test/=0;
            gosub MENU;
            call csrtype(0);
        endif;

        m=m+1;
    endo;

    cls;
    call csrtype(old_cursor);
    closeall fh;
    retp(EQM,a,b);


    PRESENTATION:
        cls;
        locate 1,6;  print chrs(45*ones(1,74));
        l=1;
        do until l>3;
            locate 1+l,6; print chrs(124);
            locate 1+l,79; print chrs(124);
            l=l+1;
        endo;
        locate 5,6; print chrs(45*ones(1,74));
        locate 3,35; print "BACKPROPAGATION";
        locate 6,20;  print chrs(45*ones(1,42));
        l=1;
        do until l>13;
            locate 6+l,20; print chrs(124);
            locate 6+l,61; print chrs(124);
            l=l+1;
        endo;
        locate 19,20;  print chrs(45*ones(1,42));
        locate 7,22; print "Phi:";
        locate 7,50; print ftos(phi,"%lf",10,5);
        locate 8,22; print "Kappa :";
        locate 8,50; print ftos(kappa,"%lf",10,5);
        locate 9,22; print "Mu :  ";
        locate 9,50; print ftos(mu,"%lf",10,5);
        locate 10,22; print "Theta :";
        locate 10,50; print ftos(theta,"%lf",10,5);
        locate 11,22; print "Epsilon :";
        locate 11,50; print ftos(epsilon,"%lf",10,5);
        locate 13,22; print "ENTREES :";
        locate 13,50; print ftos(I,"%lf",10,0);
        locate 14,22; print "CACHES :";
        locate 14,50; print ftos(J,"%lf",10,0);
        locate 15,22; print "SORTIES :";
        locate 15,50; print ftos(K,"%lf",10,0);
        locate 20,6;  print chrs(45*ones(1,74));
        locate 21,6; print chrs(124);
        locate 21,79; print chrs(124);
        locate 22,6; print chrs(124);
        locate 22,79; print chrs(124);
        locate 23,6; print chrs(45*ones(1,74));
        if _BackP_Options[4] == 1;
            locate 21,15; print "Liens directs         --- ";
        else;
            locate 21,12; print "Pas de liens directs  --- ";
        endif;
        if _BackP_Options[1] == 3;
            locate 21,43; print "Fonction : Tangente Hyperbolique";
        elseif _BackP_Options[1] == 2;
            locate 21,43; print "Fonction : Logistique alternative";
        else;
            locate 21,43; print "Fonction : Logistique";
        endif;
    return;

    SAUVEGARDE:
        call writer(fh,vecr(r)');
    return;

    MENU:
        let string st[7,1] = "Changement des parametres" "Graphe (EQM)"
        "Graphe (Cible)" "Affichage des matrices" "Learning Laws"
        "Reprendre l'algorithme" "Quitter";
        choix = _BackP_Affichage(st);

        if choix==1;
            {phi,theta,mu,kappa,epsilon} =
_BackP_A1(phi,theta,mu,kappa,epsilon);
        elseif choix==2;
            call _BackP_A2(EQM);
        elseif choix==3;
            call _BackP_A3(r,T);
        elseif choix==4;
            {a,b} = _BackP_A4(a,b,r,T);
        elseif choix==5;
            {_BackP_Options,_BackP_Line_Search} =
                _BackP_A5(_BackP_Options,_BackP_Line_Search);
            _BackP_a=a; _BackP_b=b;
            {a,b,da_m,da_mm,fa_m,fa_mm,ca_m,ca_mm,ea_m,ea_mm,

db_m,db_mm,fb_m,fb_mm,cb_m,cb_mm,eb_m,eb_mm,phi,theta,mu,kappa,epsilon}
                = _BackP_Initialisation(I,J,K);
        elseif choix==6;
            if _print==1;
                gosub presentation;
            else;
                cls;
            endif;
            return;
        elseif choix==7;
            cls;
            call csrtype(old_cursor);
            closeall fh;

            retp(EQM,a,b);
        endif;

        if choix /= 6;
            goto MENU;
        endif;
    return;

endp;

proc _BackP_Affichage(st);
    local i,choix;

    cls;
    locate 2,15; print chrs(45*ones(1,44));
    i=3;
    do until i>17;
      locate i,15; print chrs(124);
      locate i,58; print chrs(124);
      i=i+1;
    endo;
    locate 18,15; print chrs(45*ones(1,44));
    i=1;
    do until i>rows(st);
        locate (1+2*i),18; print ftos(i,"%lf.",1,0)$+st[i];
        i=i+1;
    endo;
    locate 17,24; print "Votre choix :"; locate 17,38; choix=con(1,1);

    retp(choix);
endp;

proc (5) = _BackP_A1(phi,theta,mu,kappa,epsilon);
    local choix,st,c;
    choix =0;
    let string st[6,1] = "Phi" "Kappa" "Mu" "Theta" "Epsilon"
                         "Retour au menu principal";
    do until choix == 6;
        choix = _BackP_Affichage(st);
        if sumc(choix.==seqa(1,1,5)) == 1;
            locate 20,2; print "Entrez la nouvelle valeur :";
            locate 20,40; c=con(1,1);
        endif;
        if choix==1;
            phi=c;
        elseif choix==2;
            kappa=c;
        elseif choix==3;
            mu=c;
        elseif choix==4;
            theta=c;
        elseif choix==5;
            epsilon=c;
        endif;
      endo;

      retp(phi,theta,mu,kappa,epsilon);
endp;

proc (0) = _BackP_A2(EQM);
    local debut,fin;
    cls;
    locate 3,2; print ftos(rows(EQM),"EPOQUES : 1 - %lf",2,0);
    locate 5,2; print "Entrez la valeur inferieure de l'intervalle :";
    locate 5,60; debut=con(1,1);
    locate 7,2; print "Entrez la valeur superieure de l'intervalle :";
    locate 7,60; fin=con(1,1);
    if debut>fin or debut<1 or fin>rows(EQM);
        retp;
    endif;
    graphset;
    _pdate=""; _pltype=6;
    title("Erreur Quadratique Moyenne");
    xlabel("EPOQUE");
    ylabel("EQM");
    EQM = EQM[debut:fin,.];
    xy(EQM[.,1],EQM[.,2]);
    retp;
endp;

proc (0) = _BackP_A3(r,T);
    local K,s;
    cls;
    K = rows(T);
    s=seqa(1,1,cols(T));

    graphset;
    _pdate="";
    _paxes={0,1};
    xlabel("observation");
    ylabel("");

    _pltype=(6*ones(K,1))Ý(1*ones(K,1));
    title("CIBLE - PREVISION");
    xy(s,(t')~(r'));

    _pltype=(6*ones(K,1));
    title("ERREUR");
    xy(s,(t')-(r'));

    retp;
endp;

proc (2) = _BackP_A4(a,b,r,T);
    local choix,st,m1,m2,m3;
    choix=0;
    let string st[5,1] = "Matrice A" "Matrice B" "Previsions" "Erreurs"
                         "Retour au menu principal";

    do until choix == 5;
        choix = _BackP_Affichage(st);
/*
        if choix==1;
            {m1,m2,m3}=medit(a,1,1);
            if ismiss(m1)==0;
                a=m1;
            endif;
        elseif choix==2;
            {m1,m2,m3}=medit(b,1,1);
            if ismiss(m1)==0;
                b=m1;
            endif;
        elseif choix==3;
            call medit(r,1,1);
        elseif choix==4;
            call medit(r-t,1,1);
        endif;
*/
    endo;

    retp(a,b);
endp;

proc (2) =  _BackP_A5(_BackP_Options,_BackP_Line_Search);
    local choix,st;
    choix=0;
    let string st[6,1] = "Adaptative Learning Rates + Momentum"
    "Adaptative Learning Rates" "Momentum" "Steepest Descent"
    "Steepest Descent + Line Search" "Retour au menu principal";

    do until choix == 6;
        choix = _BackP_Affichage(st);
        if choix == 1;
            _BackP_Options[2] = 1;
        elseif choix == 2;
            _BackP_Options[2] = 2;
        elseif choix == 3;
            _BackP_Options[2] = 3;
        elseif choix == 4;
            _BackP_Options[2] = 4;
        elseif choix == 5;
            _BackP_Options[2] = 5;
            _BackP_Line_Search = _BackP_A5a(_BackP_Line_Search);
        endif;
    endo;

    retp(_BackP_Options,_BackP_Line_Search);
endp;

proc (1) = _BackP_A5a(_BackP_Line_Search);
    local choix,st,c;
    choix =0;
    st = "Visualiser le vecteur" $Ý "Modifier le vecteur";
    st = st $Ý "Retour au menu principal";
    do until choix == 3;
        choix = _BackP_Affichage(st);
        if choix == 1;
/*
            call medit(_BackP_Line_Search,1,1);
*/
        elseif choix == 2;
            locate 20,2; print "Entrez le nombre d'elements :";
            locate 20,40; c=con(1,1);
            cls;
            _BackP_Line_Search = con(c,1);
        endif;
    endo;
    retp(_BackP_Line_Search);
endp;

proc (0) = BackPset;
    _BackP_a = 0;
    _BackP_b = 0;
    _BackP_Sauvegarde = 0;
    _BackP_Options = {1,4,10000,0};
    _BackP_Parametres = {0.5,0.7,0.9,0.1,0.5};
    _BackP_Line_Search = {0.2,0.4,0.6,0.8,1};
  retp;
endp;


