#include <stdio.h>
#include <stdlib.h> 
#include <time.h>
#include <math.h>
#include <fcntl.h>
#define NUMPAT 4
#define NUMIN  2
#define NUMHID 9
#define NUMOUT 1
#define randn() ((double)rand()/((double)RAND_MAX+1))

main(a) {
    int i, j, k, p, np, op, ranpat[NUMPAT], epoch;
    int NumPattern = NUMPAT, NumInput = NUMIN, NumHidden = NUMHID, NumOutput = NUMOUT;
    double Input[NUMPAT][NUMIN] = { {0, 0},  {1, 0},  {0, 1},  {1, 1} };
    double Target[NUMPAT][NUMOUT] = { {0},  {1},  {1},  {0} };
    double SumH[NUMPAT][NUMHID], WeightIH[NUMIN][NUMHID], Hidden[NUMPAT][NUMHID];
    double SumO[NUMPAT][NUMOUT], WeightHO[NUMHID][NUMOUT], Output[NUMPAT][NUMOUT];
    double DeltaO[NUMOUT], SumDOW[NUMHID], DeltaH[NUMHID];
    double DeltaWeightIH[NUMIN][NUMHID], DeltaWeightHO[NUMHID][NUMOUT1];
    double Error, eta = 0.5, alpha = 0.9, smallwt = 0.5;
  
    for( j = 0 ; j < NumHidden ; j++ ) {
        for( i = 0 ; i < NumInput ; i++ ) { 
            DeltaWeightIH[i][j] = 0.0 ;
            WeightIH[i][j] = 2.0 * ( randn() - 0.5 ) * smallwt ;
        }
    }
    for( k = 0 ; k < NumOutput ; k ++ ) {
        for( j = 0 ; j < NumHidden ; j++ ) {
            DeltaWeightHO[j][k] = 0.0 ;              
            WeightHO[j][k] = 2.0 * ( randn() - 0.5 ) * smallwt ;
        }
    }
     
    for( epoch = 0 ; epoch < 100000 ; epoch++) {
        for( p = 0 ; p < NumPattern ; p++ ) {
            ranpat[p] = p ;
        }
        for( p = 0 ; p < NumPattern ; p++) {
            np = p + randn() * ( NumPattern - p ) ;
            op = ranpat[p] ; ranpat[p] = ranpat[np] ; ranpat[np] = op ;
        }
        Error = 0.0 ;
        for( np = 0 ; np < NumPattern ; np++ ) {
            p = ranpat[np];
            for( j = 0 ; j < NumHidden ; j++ ) {
                SumH[p][j] = WeightIH[0][j] ;
                for( i = 0 ; i < NumInput ; i++ ) {
                    SumH[p][j] += Input[p][i] * WeightIH[i][j] ;
                }
                Hidden[p][j] = 1.0/(1.0 + exp(-SumH[p][j])) ;
            }
            for( k = 0 ; k < NumOutput ; k++ ) {
                SumO[p][k] = WeightHO[0][k] ;
                for( j = 0 ; j < NumHidden ; j++ ) {
                    SumO[p][k] += Hidden[p][j] * WeightHO[j][k] ;
                }
                Output[p][k] = 1.0/(1.0 + exp(-SumO[p][k])) ;
                Error += 0.5 * (Target[p][k] - Output[p][k]) * (Target[p][k] - Output[p][k]) ;
                DeltaO[k] = (Target[p][k] - Output[p][k]) * Output[p][k] * (1.0 - Output[p][k]) ;
            }
            for( j = 0 ; j < NumHidden ; j++ ) {
                SumDOW[j] = 0.0 ;
                for( k = 0 ; k < NumOutput ; k++ ) {
                    SumDOW[j] += WeightHO[j][k] * DeltaO[k] ;
                }
                DeltaH[j] = SumDOW[j] * Hidden[p][j] * (1.0 - Hidden[p][j]) ;
            }
            for( j = 0 ; j < NumHidden ; j++ ) {
                DeltaWeightIH[0][j] = eta * DeltaH[j] + alpha * DeltaWeightIH[0][j] ;
                WeightIH[0][j] += DeltaWeightIH[0][j] ;
                for( i = 0 ; i < NumInput ; i++ ) { 
                    DeltaWeightIH[i][j] = eta * Input[p][i] * DeltaH[j] + alpha * DeltaWeightIH[i][j];
                    WeightIH[i][j] += DeltaWeightIH[i][j] ;
                }
            }
            for( k = 0 ; k < NumOutput ; k ++ ) {
                DeltaWeightHO[0][k] = eta * DeltaO[k] + alpha * DeltaWeightHO[0][k] ;
                WeightHO[0][k] += DeltaWeightHO[0][k] ;
                for( j = 0 ; j < NumHidden ; j++ ) {
                    DeltaWeightHO[j][k] = eta * Hidden[p][j] * DeltaO[k] + alpha * DeltaWeightHO[j][k] ;
                    WeightHO[j][k] += DeltaWeightHO[j][k] ;
                }
            }
        }
        if( epoch%100 == 0 ) fprintf(stdout, "\nEpoch %-5d :   Error = %f", epoch, Error) ;
        if( Error < 0.0004 ) break ;
    }
    
    fprintf(stdout, "\n\nNETWORK DATA - EPOCH %d\n\nPat\t", epoch) ;
    for( i = 0 ; i < NumInput ; i++ ) {
        fprintf(stdout, "Input%-4d\t", i) ;
    }
    for( k = 0 ; k < NumOutput ; k++ ) {
        fprintf(stdout, "Target%-4d\tOutput%-4d\t", k, k) ;
    }
    for( p = 0 ; p < NumPattern ; p++ ) {        
    fprintf(stdout, "\n%d\t", p) ;
        for( i = 0 ; i < NumInput ; i++ ) {
            fprintf(stdout, "%f\t", Input[p][i]) ;
        }
        for( k = 0 ; k < NumOutput ; k++ ) {
            fprintf(stdout, "%f\t%f\t", Target[p][k], Output[p][k]) ;
        }
    }
    return 1 ;
}