MainPage.xaml.cs
// http://sydney.edu.au/engineering/it/~comp4302/ann4-6s.pdf
// Fyfe C.: Artificial neural network, pages: 38-40

using System;
using System.Collections.Generic;
using System.Linq;
using System.Net;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Animation;
using System.Windows.Shapes;
using System.Windows.Media.Imaging;     // for WriteableBitmap

namespace TrainedSin
{
    public partial class MainPage : UserControl
    {
        public MainPage()
        {
            InitializeComponent();
        }

        // we train function myFunction
        // code taken from John Bullinaria at http://www.cs.bham.ac.uk/~jxb/NN/nn.html
        private void button1_Click(object sender, RoutedEventArgs e)
        {
            float error, eta = 0.4f, alpha = 0f;
            int p, j, i, k, m, n, np, op, epoch;
            int numPattern, numHidden, numInput, numOutput;
            int[] ranpat;
            float[,] sumH, sumO;
            float[,] input, hidden, output, target;
            float[,] weightIH, weightHO;
            float[,] deltaWeightIH, deltaWeightHO;
            float[] deltaO, deltaH;
            float[] sumTemp;
            Random RandomNumber;
            float[] rinput, rsumH, rhidden, rsumO, routput;
            WriteableBitmap m_RenderTarget;
            float v;

            m_RenderTarget = new WriteableBitmap(400, 200);
            image1.Source = m_RenderTarget;

            numPattern = 81;
            numInput = 1;
            numHidden = Convert.ToInt32(number1.Text);
            numOutput = 1;

            input = new float[numPattern + 1, numInput + 1];
            hidden = new float[numPattern + 1, numHidden + 1];
            output = new float[numPattern + 1, numOutput + 1];
            target = new float[numPattern + 1, numOutput + 1];
            sumH = new float[numPattern + 1, numHidden + 1];
            sumO = new float[numPattern + 1, numOutput + 1];
            weightIH = new float[numInput + 1, numHidden + 1];
            weightHO = new float[numHidden + 1, numOutput + 1];
            deltaWeightIH = new float[numInput + 1, numHidden + 1];
            deltaWeightHO = new float[numHidden + 1, numOutput + 1];
            deltaO = new float[numOutput + 1];
            deltaH = new float[numHidden + 1];
            sumTemp = new float[numHidden + 1];
            ranpat = new int[numPattern + 1];
            RandomNumber = new Random();

            rinput = new float[numInput + 1];
            rsumH = new float[numHidden + 1];
            rhidden = new float[numHidden + 1];
            rsumO = new float[numOutput + 1];
            routput = new float[numOutput + 1];

            v = -8.0f;
            for (i = 1; i <= numPattern; i++)        // patterns
            {
                input[i, 1] = v;
                target[i, 1] = myFunction(input[i, 1]);
                v += 0.2f;
            }

            for (j = 1; j <= numHidden; j++)    /* initialize WeightIH and DeltaWeightIH */
            {
                for (i = 0; i <= numInput; i++)
                {
                    deltaWeightIH[i,j] = 0.0f;
                    weightIH[i, j] = 1.0f * ((float) RandomNumber.NextDouble() - 0.5f);
                }
            }

            for (k = 1; k <= numOutput; k++)    /* initialize WeightHO and DeltaWeightHO */
            {
                for (j = 0; j <= numHidden; j++)
                {
                    deltaWeightHO[j, k] = 0.0f;
                    weightHO[j, k] = 1.0f * ((float)RandomNumber.NextDouble() - 0.5f);
                }
            }

            for (epoch = 1; epoch < 10000; epoch++)
            {
                for (p = 1; p <= numPattern; p++)
                {    /* randomize order of individuals */
                    ranpat[p] = p;
                }
                for (p = 1; p <= numPattern; p++)
                {
                    np = p + (int)(RandomNumber.NextDouble() * (numPattern + 1 - p));
                    op = ranpat[p]; ranpat[p] = ranpat[np]; ranpat[np] = op;
                }

                for (np = 1; np <= numPattern; np++)
                {
                    p = np;
                    for (j = 1; j <= numHidden; j++)    /* compute hidden unit activations */
                    {
                        sumH[p, j] = weightIH[0, j];
                        for (i = 1; i <= numInput; i++)
                        {
                            sumH[p, j] += input[p, i] * weightIH[i, j];
                        }
                        hidden[p, j] = 1.0f / (1.0f + (float)Math.Pow(Math.E, -sumH[p, j]));
                    }
                    for (k = 1; k <= numOutput; k++)    /* compute output unit activations and errors */
                    {
                        sumO[p, k] = weightHO[0, k];
                        for (j = 1; j <= numHidden; j++)
                        {
                            sumO[p, k] += hidden[p, j] * weightHO[j, k];
                        }
                        output[p, k] = 1.0f / (1.0f + (float)Math.Pow(Math.E, -sumO[p, k]));
                        deltaO[k] = (target[p, k] - output[p, k]) * output[p, k] * (1 - output[p, k]);
                    }
                    for (j = 1; j <= numHidden; j++)    /* 'back-propagate' errors to hidden layer */
                    {
                        sumTemp[j] = 0.0f;
                        for (k = 1; k <= numOutput; k++)
                        {
                            sumTemp[j] += weightHO[j, k] * deltaO[k];
                        }
                        deltaH[j] = sumTemp[j] * hidden[p, j] * (1 - hidden[p, j]);
                    }
                    for (j = 1; j <= numHidden; j++)    /* update weights WeightIH */
                    {
                        deltaWeightIH[0, j] = eta * deltaH[j] + alpha * deltaWeightIH[0, j];
                        weightIH[0, j] += deltaWeightIH[0, j];
                        for (i = 1; i <= numInput; i++)
                        {
                            deltaWeightIH[i, j] = eta * input[p, i] * deltaH[j] + alpha * deltaWeightIH[i, j];
                            weightIH[i, j] += deltaWeightIH[i, j];
                        }
                    }
                    for (k = 1; k <= numOutput; k++)    /* update weights WeightHO */
                    {
                        deltaWeightHO[0, k] = eta * deltaO[k] + alpha * deltaWeightHO[0, k];
                        weightHO[0, k] += deltaWeightHO[0, k];
                        for (j = 1; j <= numHidden; j++)
                        {
                            deltaWeightHO[j, k] = eta * hidden[p, j] * deltaO[k] + alpha * deltaWeightHO[j, k];
                            weightHO[j, k] += deltaWeightHO[j, k];
                        }
                    }
                }

                k = 1;
                error = 0.0f;
                for (p = 1; p <= numPattern; p++)
                {
                    error += 0.5f * (target[p, k] - output[p, k]) * (target[p, k] - output[p, k]);
                }

                if (error < Convert.ToDouble(number2.Text)) break;
            }

            for (m = 0; m < 400; m++)       // background of the plot
            {
                for (n = 0; n < 200; n++)
                {
                    m_RenderTarget.Pixels[m * 200 + n] = RGBA((byte)(200), (byte)(200), (byte)(0), 255);
                }
            }

            for (m = 0; m < 400; m++)       // calculating the result
            {
                rinput[1] = 1.0f / 25.0f * m - 8;

                for (j = 1; j <= numHidden; j++)
                {
                    rsumH[j] = weightIH[0, j];
                    for (i = 1; i <= numInput; i++)
                    {
                        rsumH[j] += rinput[i] * weightIH[i, j];
                    }
                    rhidden[j] = 1.0f / (1.0f + (float)Math.Pow(Math.E, -rsumH[j]));
                }
                for (k = 1; k <= numOutput; k++)
                {
                    rsumO[k] = weightHO[0, k];
                    for (j = 1; j <= numHidden; j++)
                    {
                        rsumO[k] += rhidden[j] * weightHO[j, k];
                    }
                    routput[k] = 1.0f / (1.0f + (float)Math.Pow(Math.E, -rsumO[k]));
                }

                n = (int) ((1.0f - routput[1]) * 100.0f);
                if (n < 200)
                {
                    m_RenderTarget.Pixels[n * 400 + m] = RGBA((byte)(200), (byte)(0), (byte)(0), 255);
                }

                n = (int)((1.0f - myFunction(rinput[1])) * 100.0f);
                if (n < 200)
                {
                    m_RenderTarget.Pixels[n * 400 + m] = RGBA((byte)(100), (byte)(100), (byte)(100), 255);
                }

            }

            epochstring.Text = "epochs: " + Convert.ToString(epoch);

        }

        float myFunction(float x)
        {
            return (1 + (float)Math.Sin(x)) / 2;
        }

        static public int RGBA(byte byRed, byte byGreen, byte byBlue, byte alpha)
        {
            return (alpha * 16777216) + (byRed * 65536) + (byGreen * 256) + byBlue;
        }
    }
}