Commit 13506dea authored by Vivien Kraus's avatar Vivien Kraus
Browse files

Merge branch '1-implement-document' into 'master'

Resolve "Implement & document"

Closes #1

See merge request !1
parents 430f7f2b d3936ce8
Pipeline #187749 passed with stages
in 9 minutes and 30 seconds
# julien -- train a multi-layer perceptron
# Copyright (C) 2019 Vivien Kraus
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
image: debian
stages:
- publish
pages:
stage: publish
script: "./publish"
artifacts:
paths:
- public
expire_in: 1 week
[submodule ".gnulib"]
path = .gnulib
url = git://git.sv.gnu.org/gnulib.git
.gnulib @ a4740cc6
Subproject commit a4740cc6fa277fb41aab06563657458e03ab0346
# julien -- train a multi-layer perceptron
# Copyright (C) 2019 Vivien Kraus
# Licensed under the Apache License, Version 2.0 (the "License"); you
# may not use this file except in compliance with the License. You
# may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
# implied. See the License for the specific language governing
# permissions and limitations under the License.
AM_CPPFLAGS += -I $(srcdir)/%reldir% -I %reldir%
noinst_LTLIBRARIES += %reldir%/libjulien.la
%canon_reldir%_libjulien_la_SOURCES = \
%reldir%/julien.h \
%reldir%/julien.c
%canon_reldir%_libjulien_la_LIBADD = \
lib/libgnu.la
/*
* julien -- train a multi-layer perceptron
* Copyright (C) 2019 Vivien Kraus
* Licensed under the Apache License, Version 2.0 (the "License"); you
* may not use this file except in compliance with the License. You
* may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
* implied. See the License for the specific language governing
* permissions and limitations under the License.
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /* HAVE_CONFIG_H */
#include <stdlib.h>
#include "julien.h"
#include <assert.h>
#include <stdio.h>
static inline void *
xmalloc (size_t sz)
{
void *ret = malloc (sz);
if (ret == NULL)
{
abort ();
}
return ret;
}
struct julien
{
size_t dim_input;
size_t dim_output;
size_t n_hidden;
size_t *restrict dim_hidden;
double *restrict parameters;
double *restrict forward;
double *restrict backward;
};
static size_t
count_parameters (const struct julien *julien)
{
size_t i;
size_t ret = 0;
if (julien->n_hidden == 0)
{
return julien->dim_input * julien->dim_output;
}
ret += julien->dim_input * julien->dim_hidden[0];
ret += julien->dim_hidden[julien->n_hidden - 1] * julien->dim_output;
for (i = 1; i < julien->n_hidden; i++)
{
ret += julien->dim_hidden[i - 1] * julien->dim_hidden[i];
}
return ret;
}
static size_t
count_hidden_neurons (const struct julien *julien)
{
size_t ret = 0;
size_t i_layer;
for (i_layer = 0; i_layer < julien->n_hidden; i_layer++)
{
ret += julien->dim_hidden[i_layer];
}
return ret;
}
struct julien *
julien_alloc (size_t dim_input, size_t dim_output, size_t n_hidden,
const size_t *dim_hidden)
{
struct julien *ret = xmalloc (sizeof (struct julien));
size_t i, n_parameters;
double initial_value;
size_t n_activations = dim_input + dim_output;
ret->dim_input = dim_input;
ret->dim_output = dim_output;
ret->n_hidden = n_hidden;
ret->dim_hidden = xmalloc (n_hidden * sizeof (size_t));
for (i = 0; i < n_hidden; i++)
{
ret->dim_hidden[i] = dim_hidden[i];
}
n_parameters = count_parameters (ret);
ret->parameters = xmalloc (n_parameters * sizeof (double));
initial_value = 1.0 / n_parameters;
for (i = 0; i < n_parameters; i++)
{
double ampl = rand () / (RAND_MAX + 1.0);
ret->parameters[i] = (ampl * 2 - 1) * initial_value;
}
n_activations += count_hidden_neurons (ret);
ret->forward = xmalloc (n_activations * sizeof (double));
ret->backward = xmalloc (n_activations * sizeof (double));
return ret;
}
void
julien_free (struct julien *julien)
{
free (julien->parameters);
free (julien->dim_hidden);
free (julien->forward);
free (julien->backward);
free (julien);
}
const double *
julien_data (const struct julien *julien, size_t *n_parameters)
{
if (n_parameters)
{
*n_parameters = count_parameters (julien);
}
return julien->parameters;
}
void
julien_set_data (struct julien *julien, size_t start, size_t n,
const double *restrict parameters)
{
size_t n_total = count_parameters (julien);
size_t i;
for (i = 0; i < n && i + start < n_total; i++)
{
julien->parameters[i + start] = parameters[i];
}
}
static void
forward_layer (size_t n_inputs, const double *restrict inputs,
size_t n_outputs, double *restrict outputs,
const double *restrict weights, int first)
{
size_t i_output;
size_t i_weight = 0;
for (i_output = 0; i_output < n_outputs; i_output++)
{
size_t i_link;
outputs[i_output] = 0;
for (i_link = 0; i_link < n_inputs; i_link++)
{
if (first || inputs[i_link] >= 0)
{
outputs[i_output] += weights[i_weight] * inputs[i_link];
}
i_weight++;
}
}
assert (i_weight == n_outputs * n_inputs);
}
static void
backward_layer (size_t n_inputs, double *restrict inputs,
size_t n_outputs, const double *restrict outputs,
const double *restrict weights,
const double *restrict activations, int relu)
{
size_t i_input, i_output;
size_t i_weight = 0;
for (i_input = 0; i_input < n_inputs; i_input++)
{
inputs[i_input] = 0;
}
for (i_output = 0; i_output < n_outputs; i_output++)
{
for (i_input = 0; i_input < n_inputs; i_input++)
{
const double weight = weights[i_weight++];
if (!relu || activations[i_input] > 0)
{
inputs[i_input] += weight * outputs[i_output];
}
}
}
}
static void
gradient_layer (size_t n_inputs, const double *restrict inputs,
size_t n_outputs, const double *restrict outputs,
double *restrict weights, double learning_rate,
double *restrict sum_updates, double *restrict n_updates,
int relu_in)
{
size_t i_input, i_output;
size_t i_weight = 0;
for (i_output = 0; i_output < n_outputs; i_output++)
{
for (i_input = 0; i_input < n_inputs; i_input++)
{
double in_value = inputs[i_input];
double error_value = outputs[i_output];
double g, abs_g;
if (relu_in && in_value < 0)
{
in_value = 0;
}
g = in_value * error_value;
abs_g = g;
if (abs_g < 0)
{
abs_g = -g;
}
weights[i_weight++] -= learning_rate * g;
*sum_updates += abs_g;
*n_updates += 1.0;
}
}
}
static inline void
forward_net_aux (struct julien *julien, size_t *restrict i_neuron,
size_t *restrict i_parameter, size_t layer_dim_input,
size_t layer_dim_output, double *restrict activations,
int *restrict first)
{
forward_layer (layer_dim_input, &(activations[*i_neuron]), layer_dim_output,
&(activations[*i_neuron + layer_dim_input]),
&(julien->parameters[*i_parameter]), *first);
*first = 0;
*i_neuron += layer_dim_input;
*i_parameter += (layer_dim_input) * (layer_dim_output);
}
static inline size_t
forward_net (struct julien *julien, size_t dim_input,
const double *restrict input)
{
size_t i_neuron, i_parameter, i_layer;
size_t layer_dim_input = julien->dim_input;
size_t layer_dim_output;
int first = 1;
for (i_neuron = 0; i_neuron < dim_input && i_neuron < julien->dim_input;
i_neuron++)
{
julien->forward[i_neuron] = input[i_neuron];
}
for (i_neuron = dim_input; i_neuron < julien->dim_input; i_neuron++)
{
julien->forward[i_neuron] = 0;
}
i_neuron = 0;
i_parameter = 0;
for (i_layer = 0; i_layer < julien->n_hidden; i_layer++)
{
layer_dim_output = julien->dim_hidden[i_layer];
forward_net_aux (julien, &i_neuron, &i_parameter, layer_dim_input,
layer_dim_output, julien->forward, &first);
layer_dim_input = layer_dim_output;
}
layer_dim_output = julien->dim_output;
forward_net_aux (julien, &i_neuron, &i_parameter, layer_dim_input,
layer_dim_output, julien->forward, &first);
assert (i_parameter == count_parameters (julien));
return i_neuron + layer_dim_output;
}
static inline void
backward_net_aux (struct julien *julien, size_t *restrict i_neuron,
size_t *restrict i_parameter, size_t layer_dim_input,
size_t layer_dim_output, double *restrict activations,
int relu)
{
*i_neuron -= layer_dim_input;
*i_parameter -= layer_dim_input * layer_dim_output;
backward_layer (layer_dim_input, &(activations[*i_neuron]),
layer_dim_output,
&(activations[*i_neuron + layer_dim_input]),
&(julien->parameters[*i_parameter]),
&(julien->forward[*i_neuron]), relu);
}
static inline size_t
backward_net (struct julien *julien, size_t dim_output,
const double *restrict output)
{
size_t n_activations =
julien->dim_input + julien->dim_output + count_hidden_neurons (julien);
size_t i_neuron = n_activations, i_parameter =
count_parameters (julien), j, i_layer;
size_t layer_dim_input;
size_t layer_dim_output = julien->dim_output;
/* The errors for the last layer is easy, since we are doing least
* squares regression: we simply subtract the target */
i_neuron -= layer_dim_output;
for (j = 0; j < julien->dim_output; j++)
{
julien->backward[n_activations - j - 1] =
julien->forward[n_activations - j - 1];
}
for (j = 0; j < dim_output && j < julien->dim_output; j++)
{
julien->backward[i_neuron + j] -= output[j];
}
for (i_layer = julien->n_hidden; i_layer != 0; i_layer--)
{
layer_dim_input = julien->dim_hidden[i_layer - 1];
backward_net_aux (julien, &i_neuron, &i_parameter, layer_dim_input,
layer_dim_output, julien->backward, 1);
layer_dim_output = layer_dim_input;
}
layer_dim_input = julien->dim_input;
backward_net_aux (julien, &i_neuron, &i_parameter, layer_dim_input,
layer_dim_output, julien->backward, 0);
assert (i_parameter == 0);
return i_neuron;
}
static inline double
gradient_net (struct julien *julien, double learning_rate)
{
size_t layer_dim_input = julien->dim_input;
size_t layer_dim_output;
double n_updates = 0;
double sum_updates = 0;
size_t i_layer;
size_t i_neuron = 0, i_parameter = 0;
int relu_in = 0;
for (i_layer = 0; i_layer < julien->n_hidden; i_layer++)
{
layer_dim_output = julien->dim_hidden[i_layer];
gradient_layer (layer_dim_input, &(julien->forward[i_neuron]),
layer_dim_output,
&(julien->backward[i_neuron + layer_dim_input]),
&(julien->parameters[i_parameter]), learning_rate,
&sum_updates, &n_updates, relu_in);
relu_in = 1;
i_neuron += layer_dim_input;
i_parameter += layer_dim_input * layer_dim_output;
layer_dim_input = layer_dim_output;
}
layer_dim_output = julien->dim_output;
gradient_layer (layer_dim_input, &(julien->forward[i_neuron]),
layer_dim_output,
&(julien->backward[i_neuron + layer_dim_input]),
&(julien->parameters[i_parameter]), learning_rate,
&sum_updates, &n_updates, relu_in);
i_neuron += layer_dim_input;
i_parameter += layer_dim_input * layer_dim_output;
assert (i_neuron == count_hidden_neurons (julien) + julien->dim_input);
assert (i_parameter == count_parameters (julien));
return sum_updates / n_updates;
}
double
julien_learn (struct julien *julien, size_t dim_input,
const double *restrict input, size_t dim_output,
const double *restrict output, double learning_rate)
{
double average_update = 0;
size_t n_activations =
julien->dim_input + julien->dim_output + count_hidden_neurons (julien);
if (forward_net (julien, dim_input, input) != n_activations)
{
assert (0);
}
/* activations are progressively replaced with the errors */
if (backward_net (julien, dim_output, output) != 0)
{
assert (0);
}
/* fprintf (stderr, "%s:%d: parameters:", __FILE__, __LINE__); */
/* for (size_t i = 0; i < count_parameters (julien); i++) */
/* { */
/* fprintf (stderr, " %f", julien->parameters[i]); */
/* } */
/* fprintf (stderr, "\n%s:%d: forward:", __FILE__, __LINE__); */
/* for (size_t i = 0; i < n_activations; i++) */
/* { */
/* fprintf (stderr, " %f", julien->forward[i]); */
/* } */
/* fprintf (stderr, "\n%s:%d: backward:", __FILE__, __LINE__); */
/* for (size_t i = 0; i < n_activations; i++) */
/* { */
/* fprintf (stderr, " %f", julien->backward[i]); */
/* } */
/* fprintf (stderr, "\n"); */
return gradient_net (julien, learning_rate);
}
size_t
julien_predict (struct julien *julien, size_t dim_input, const double *input,
size_t max_output, size_t start_output,
double *restrict output)
{
size_t n_activations =
julien->dim_input + julien->dim_output + count_hidden_neurons (julien);
size_t j;
if (forward_net (julien, dim_input, input) != n_activations)
{
assert (0);
}
for (j = 0; j < max_output && j + start_output < julien->dim_output; j++)
{
output[j] =
julien->forward[n_activations - julien->dim_output + start_output +
j];
}
return julien->dim_output;
}
/*
* julien -- train a multi-layer perceptron
* Copyright (C) 2019 Vivien Kraus
* Licensed under the Apache License, Version 2.0 (the "License"); you
* may not use this file except in compliance with the License. You
* may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
* implied. See the License for the specific language governing
* permissions and limitations under the License.
*/
#ifndef H_JULIEN_INCLUDED
#define H_JULIEN_INCLUDED
#include <stddef.h>
#ifdef __cplusplus
extern "C"
{
#endif /* __cplusplus */
struct julien;
/* Create a new perceptron with the given input and output
* dimensions, and with hidden layers having different dimensions.
* Cannot return NULL. */
struct julien *julien_alloc (size_t dim_input,
size_t dim_output,
size_t n_hidden, const size_t *dim_hidden);
/* Reciprocal to julien_alloc. */
void julien_free (struct julien *julien);
/* Get the parameters. */
const double *julien_data (const struct julien *julien,
size_t *n_parameters);
/* Load the parameters. */
void julien_set_data (struct julien *julien,
size_t start, size_t n, const double *parameters);
/* Learn from (input, output) at a given learning rate, and return
* the mean abs update. dim_input and dim_output are the sizes of
* the input and output arrays. They should match the input and
* output dimensions, but if they are less than that, then the
* inputs and outputs are padded with 0s. */
double julien_learn (struct julien *julien,
size_t dim_input,
const double *input,
size_t dim_output,
const double *output, double learning_rate);
/* Fill the output array, discarding the first start_output neurons,
* considering up to max_output neurons, and return the total number
* of output neurons. */
size_t julien_predict (struct julien *julien,
size_t dim_input,
const double *input,
size_t max_output,
size_t start_output, double *output);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* not H_JULIEN_INCLUDED */
#+title: Learn a multi-layer perceptron with Julien
#+author: Vivien Kraus
#+email: vivien@planete-kraus.eu
#+language: en
# Copyright (C) 2019 Vivien Kraus
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="style.css" />
#+OPTIONS: tex:dvipng
* Introduction
A multi-layer perceptron is a machine learning algorithm that is able
to learn some useful functions. In this package, the outputs are
real-valued.
* First example: the SUM function!
Let us take an example! Suppose that we want to learn the /sum/
function. The function has two values, that can be either 0 or 1, and
returns the sum of the two variables. So the dataset is in
table [[tbl-sum-data]].
#+name: tbl-sum-data
#+caption: The "SUM" function
| Input x | Input y | Output |
|---------+---------+--------|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 2 |
So, how do we learn this? Let us write our first julien-enabled
program! Write code [[optimize-sum-code]] to a file, [[file:main_sum.c]].
#+name: optimize-sum-code
#+caption: Learn to sum with julien!
#+BEGIN_SRC C -i -n -r :eval no :tangle main_sum.c
/* Ensure that julien.h is in the current directory */
#include "julien.h"
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
/* Print the model internals as a CSV row */
static void dump_model (const struct julien *jl);
int
main ()
{
/* We have two inputs, and one output. We do not use the other
,* arguments for now. */
struct julien *jl = julien_alloc (2, 1, 0, NULL); // (ref:sum-call-alloc)
/* Keep track of how much we updated our model */
double update;
/* Keep learning! */
do
{
double input_data[2];
update = 0;
/* Iterate over the dataset */
for (input_data[0] = 0; input_data[0] <= 1; input_data[0]++)
{
for (input_data[1] = 0; input_data[1] <= 1; input_data[1]++)
{
/* We want to learn to produce the sum! */
double expected_output = input_data[0] + input_data[1];
/* Call the learning function (see the commentary) */
update += julien_learn (jl, 2, input_data, 1, &expected_output, 1); // (ref:sum-call-learn)
}
}
/* Visualize our progress! */
fprintf (stderr, "%s:%d: mean update is %g\n",
__FILE__, __LINE__, update / 4);
}
while (update > 1e-8);
dump_model (jl);