Speech/music discrimination (not used for anything yet)
authorJean-Marc Valin <jmvalin@jmvalin.ca>
Thu, 17 Nov 2011 11:21:07 +0000 (19:21 +0800)
committerJean-Marc Valin <jmvalin@jmvalin.ca>
Fri, 13 Jul 2012 18:50:34 +0000 (14:50 -0400)
Also, reducing the VBR rate on panned mono

celt/celt.c
celt/celt.h
src/analysis.c
src/mlp.c [new file with mode: 0644]
src/mlp.h [new file with mode: 0644]
src/mlp_data.c [new file with mode: 0644]
src/mlp_train.c [new file with mode: 0644]
src/mlp_train.h [new file with mode: 0644]
src/tansig_table.h [new file with mode: 0644]

index 6d50a27..8c5e56e 100644 (file)
@@ -791,7 +791,8 @@ static void init_caps(const CELTMode *m,int *cap,int LM,int C)
 }
 
 static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
-      const opus_val16 *bandLogE, int end, int LM, int C, int N0, AnalysisInfo *analysis)
+      const opus_val16 *bandLogE, int end, int LM, int C, int N0,
+      AnalysisInfo *analysis, opus_val16 *stereo_saving)
 {
    int i;
    opus_val32 diff=0;
@@ -819,6 +820,12 @@ static int alloc_trim_analysis(const CELTMode *m, const celt_norm *X,
          trim_index-=2;
       else if (sum > QCONST16(.8f,10))
          trim_index-=1;
+#ifndef FIXED_POINT
+      *stereo_saving = -.5*log2(1.01-sum*sum);
+      /*printf("%f\n", *stereo_saving);*/
+#else
+      *stereo_saving = 0;
+#endif
    }
 
    /* Estimate spectral tilt */
@@ -944,6 +951,7 @@ int celt_encode_with_ec(CELTEncoder * restrict st, const opus_val16 * pcm, int f
    int anti_collapse_on=0;
    int silence=0;
    opus_val16 tf_estimate=0;
+   opus_val16 stereo_saving = 0;
    ALLOC_STACK;
 
    if (nbCompressedBytes<2 || pcm==NULL)
@@ -1409,7 +1417,7 @@ int celt_encode_with_ec(CELTEncoder * restrict st, const opus_val16 * pcm, int f
    if (tell+(6<<BITRES) <= total_bits - total_boost)
    {
       alloc_trim = alloc_trim_analysis(st->mode, X, bandLogE,
-            st->end, LM, C, N, &st->analysis);
+            st->end, LM, C, N, &st->analysis, &stereo_saving);
       ec_enc_icdf(enc, alloc_trim, trim_icdf, 7);
       tell = ec_tell_frac(enc);
    }
@@ -1470,6 +1478,8 @@ int celt_encode_with_ec(CELTEncoder * restrict st, const opus_val16 * pcm, int f
 #ifndef FIXED_POINT
      if (st->analysis.valid && st->analysis.activity<.4)
         target -= (coded_bins<<BITRES)*2*(.4-st->analysis.activity);
+
+     target -= MIN32(target/3, stereo_saving*(st->mode->eBands[intensity]<<LM<<BITRES));
 #endif
 
 #ifdef FIXED_POINT
@@ -1548,6 +1558,7 @@ int celt_encode_with_ec(CELTEncoder * restrict st, const opus_val16 * pcm, int f
         /*printf ("+%d\n", adjust);*/
      }
      nbCompressedBytes = IMIN(nbCompressedBytes,nbAvailableBytes+nbFilledBytes);
+     /*printf("%d\n", nbCompressedBytes*50*8);*/
      /* This moves the raw bits to take into account the new compressed size */
      ec_enc_shrink(enc, nbCompressedBytes);
    }
index 54bca44..8948401 100644 (file)
@@ -57,6 +57,7 @@ typedef struct {
    opus_val16 activity;
    int boost_band[2];
    opus_val16 boost_amount[2];
+   opus_val16 music_prob;
 }AnalysisInfo;
 
 #define __celt_check_mode_ptr_ptr(ptr) ((ptr) + ((ptr) - (const CELTMode**)(ptr)))
index 6c1db2f..b39defe 100644 (file)
 #include "arch.h"
 #include "quant_bands.h"
 #include <stdio.h>
+#ifndef FIXED_POINT
+#include "mlp.c"
+#include "mlp_data.c"
+#endif
 
 #ifndef M_PI
 #define M_PI 3.141592653
@@ -103,6 +107,7 @@ void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEnc
     float slope=0;
     float frame_stationarity;
     float relativeE;
+    float frame_prob;
     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
 
     kfft = mode->mdct.kfft[0];
@@ -294,6 +299,26 @@ void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEnc
     features[25] = info->activity;
     features[26] = frame_stationarity;
 
+#ifndef FIXED_POINT
+    mlp_process(&net, features, &frame_prob);
+    frame_prob = .5*(frame_prob+1);
+    /*frame_prob = .45*frame_prob + .55*frame_prob*frame_prob*frame_prob;*/
+    /*printf("%f\n", frame_prob);*/
+    {
+       float alpha, beta;
+       float p0, p1;
+       alpha = .01;
+       beta = .2;
+       p0 = (1-info->music_prob)*(1-alpha) +    info->music_prob *alpha;
+       p1 =    info->music_prob *(1-alpha) + (1-info->music_prob)*alpha;
+       p0 *= pow(1-frame_prob, beta);
+       p1 *= pow(frame_prob, beta);
+       info->music_prob = p1/(p0+p1);
+       /*printf("%f\n", info->music_prob);*/
+    }
+#else
+    info->music_prob = 0;
+#endif
     /*for (i=0;i<27;i++)
        printf("%f ", features[i]);
     printf("\n");*/
diff --git a/src/mlp.c b/src/mlp.c
new file mode 100644 (file)
index 0000000..200f1cb
--- /dev/null
+++ b/src/mlp.c
@@ -0,0 +1,103 @@
+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#include <math.h>
+#include "mlp.h"
+#include "arch.h"
+#include "tansig_table.h"
+#define MAX_NEURONS 100
+
+#ifdef FIXED_POINT
+extern const opus_val16 tansig_table[501];
+static inline opus_val16 tansig_approx(opus_val32 _x) /* Q19 */
+{
+       int i;
+       opus_val16 xx; /* Q11 */
+       /*double x, y;*/
+       opus_val16 dy, yy; /* Q14 */
+       /*x = 1.9073e-06*_x;*/
+       if (_x>=QCONST32(10,19))
+               return QCONST32(1.,14);
+       if (_x<=-QCONST32(10,19))
+               return -QCONST32(1.,14);
+       xx = EXTRACT16(SHR32(_x, 8));
+       /*i = lrint(25*x);*/
+       i = SHR32(ADD32(1024,MULT16_16(25, xx)),11);
+       /*x -= .04*i;*/
+       xx -= EXTRACT16(SHR32(MULT16_16(20972,i),8));
+       /*x = xx*(1./2048);*/
+       /*y = tansig_table[250+i];*/
+       yy = tansig_table[250+i];
+       /*y = yy*(1./16384);*/
+       dy = 16384-MULT16_16_Q14(yy,yy);
+       yy = yy + MULT16_16_Q14(MULT16_16_Q11(xx,dy),(16384 - MULT16_16_Q11(yy,xx)));
+       return yy;
+}
+#else
+/*extern const float tansig_table[501];*/
+static inline double tansig_approx(double x)
+{
+       int i;
+       double y, dy;
+       if (x>=10)
+               return 1;
+       if (x<=-10)
+               return -1;
+       i = lrint(25*x);
+       x -= .04*i;
+       y = tansig_table[250+i];
+       dy = 1-y*y;
+       y = y + x*dy*(1 - y*x);
+       return y;
+}
+#endif
+
+void mlp_process(const MLP *m, const opus_val16 *in, opus_val16 *out)
+{
+       int j;
+       opus_val16 hidden[MAX_NEURONS];
+       const opus_val16 *W = m->weights;
+       /* Copy to tmp_in */
+       for (j=0;j<m->topo[1];j++)
+       {
+               int k;
+               opus_val32 sum = SHL32(EXTEND32(*W++),8);
+               for (k=0;k<m->topo[0];k++)
+                       sum = MAC16_16(sum, in[k],*W++);
+               hidden[j] = tansig_approx(sum);
+       }
+       for (j=0;j<m->topo[2];j++)
+       {
+               int k;
+               opus_val32 sum = SHL32(EXTEND32(*W++),14);
+               for (k=0;k<m->topo[1];k++)
+                       sum = MAC16_16(sum, hidden[k], *W++);
+               out[j] = tansig_approx(EXTRACT16(PSHR32(sum,17)));
+       }
+}
+
diff --git a/src/mlp.h b/src/mlp.h
new file mode 100644 (file)
index 0000000..68ff68d
--- /dev/null
+++ b/src/mlp.h
@@ -0,0 +1,41 @@
+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef _MLP_H_
+#define _MLP_H_
+
+#include "arch.h"
+
+typedef struct {
+       int layers;
+       const int *topo;
+       const opus_val16 *weights;
+} MLP;
+
+void mlp_process(const MLP *m, const opus_val16 *in, opus_val16 *out);
+
+#endif /* _MLP_H_ */
diff --git a/src/mlp_data.c b/src/mlp_data.c
new file mode 100644 (file)
index 0000000..fdd32db
--- /dev/null
@@ -0,0 +1,77 @@
+#include "mlp.h"
+
+/* RMS error was 0.196871, seed was 1321340808 */
+
+static const float weights[291] = {
+
+/* hidden layer */
+1.93994, 0.00636575, -0.0838112, 0.188811, -0.157845, 
+-0.122662, 0.296779, -0.066386, -0.0764464, -0.00372055, 
+-0.0397377, -0.000976218, -0.03931, 0.0111525, -0.0377797, 
+-0.003592, 0.00213057, 0.115952, -0.0864595, 0.170621, 
+-0.139312, -0.125683, 0.226746, -0.148058, -3.11536, 
+-5.7119, -0.325896, -4.37802, 2.1242, -0.119952, 
+0.0232531, -0.0998321, -0.0909719, -0.164338, 0.0370311, 
+0.0196689, 0.0495503, -0.267277, -0.15925, -0.129835, 
+-0.171845, -0.0672326, -0.0319364, -0.0960325, 0.132835, 
+0.0978292, -0.0204049, -0.128357, -0.0582566, -0.21682, 
+0.00496659, -0.0224912, -2.3249, 2.10627, -5.06275, 
+0.300689, -1.05938, 0.111387, 0.100606, 0.122446, 
+-0.0175274, 0.0107236, -0.030947, -0.0712338, -0.0456196, 
+0.0158188, 0.0139863, 0.0122389, 0.0426144, 0.00963211, 
+0.00741379, 0.014572, -0.0365356, 0.0780221, 0.0835844, 
+0.101463, -0.0194471, 0.016752, -0.0360326, -0.0671933, 
+5.35889, -6.06707, 1.35677, -1.90924, 0.0347801, 
+0.0122876, 0.00258179, -0.0217294, 0.0827611, 0.0859281, 
+-0.00417207, -0.109872, -0.238913, -0.288535, -0.0319008, 
+0.156671, -0.00911369, -0.0351284, 0.0355504, 0.101236, 
+-0.140194, -0.128439, 0.0275677, -0.0507381, 0.106048, 
+0.0672367, 0.00438842, -0.0925318, 5.68238, -3.47798, 
+0.246634, 0.0970976, -1.33011, 0.0498353, 0.179046, 
+0.0162675, -0.102764, -0.227255, 0.234701, -0.00777973, 
+0.0767733, -0.00420136, 0.0344874, -0.0332389, 0.062122, 
+-0.0360523, 0.0461029, 0.0861842, 0.0136479, 0.0133092, 
+0.165541, -0.0573712, -0.0694408, -0.196571, 0.222621, 
+0.0197353, 3.42359, 5.23165, -1.10221, 3.66079, 
+0.40144, -0.493484, 0.217106, -0.0143906, 0.295599, 
+-0.614104, 0.596788, 0.956514, 0.107316, -0.172138, 
+-0.111201, 0.0162694, -0.136564, 0.0567972, -0.107051, 
+-0.0578785, 0.0597572, -0.592051, 0.11802, -0.0846178, 
+0.144399, -0.386859, 0.429763, 0.763419, 8.40166, 
+4.25269, -3.25962, 2.04492, -1.54948, 0.0286627, 
+0.0855541, -0.128902, 0.0428149, 0.147296, -0.178688, 
+0.582621, -0.0423034, -0.168806, -0.0930681, -0.0505222, 
+-0.059881, 0.0344017, -0.0538223, -0.0095173, 0.044275, 
+0.178126, 0.0321441, -0.192936, -0.0359919, 0.0449504, 
+-0.255187, 0.330503, 14.3362, -12.7585, 2.10511, 
+1.00446, -1.5146, 0.00315578, -0.0189675, 0.0506854, 
+-0.0306224, -0.0343434, -0.0222091, 0.00040356, -0.179946, 
+-0.213007, -0.046152, 0.0122855, 0.0335543, -0.0172102, 
+0.0236597, 0.088535, -0.0980871, -0.129909, -0.019153, 
+0.0544563, -0.0272701, -0.00304803, -0.00145721, 0.0190295, 
+-6.75401, 2.83619, 2.38708, -0.904901, 0.670252, 
+-0.0809205, -0.077534, -0.0347895, -0.0143415, -0.00527138, 
+0.0400907, 0.041551, 0.00823289, -0.00772847, -0.0172196, 
+-0.0125943, -0.0285652, -0.00141913, -0.010938, -0.0154068, 
+0.0149916, -0.0577316, -0.0750255, -0.019028, -0.0175507, 
+-0.00248046, 0.0350994, 0.0396102, -0.334886, 1.32123, 
+-0.363775, 0.0925417, 7.5025, 0.76236, 0.489961, 
+0.514362, 0.350457, 0.321636, -0.000131804, 0.0942301, 
+0.506788, -0.325235, 0.162356, -0.147705, 0.155451, 
+-0.111074, 0.120173, 0.0586432, 0.407685, 0.374031, 
+0.510908, 0.25445, 0.285288, 0.184939, 0.0386202, 
+0.089713, 12.6662, 2.54239, -14.9728, 7.46559, 
+
+/* output layer */
+-4.63633, -1.25936, -1.33365, 4.91614, 1.1609, 
+1.30642, -0.780207, 1.09432, -1.46686, 8.41454, 
+1.55149, };
+
+static const int topo[3] = {27, 10, 1};
+
+const MLP net = {
+       3,
+       topo,
+       weights
+};
+
diff --git a/src/mlp_train.c b/src/mlp_train.c
new file mode 100644 (file)
index 0000000..a9d548b
--- /dev/null
@@ -0,0 +1,495 @@
+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#include "mlp_train.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <semaphore.h>
+#include <pthread.h>
+#include <time.h>
+#include <signal.h>
+
+int stopped = 0;
+
+void handler(int sig)
+{
+       stopped = 1;
+       signal(sig, handler);
+}
+
+MLPTrain * mlp_init(int *topo, int nbLayers, float *inputs, float *outputs, int nbSamples)
+{
+       int i, j, k;
+       MLPTrain *net;
+       int inDim, outDim;
+       net = malloc(sizeof(*net));
+       net->topo = malloc(nbLayers*sizeof(net->topo[0]));
+       for (i=0;i<nbLayers;i++)
+               net->topo[i] = topo[i];
+       inDim = topo[0];
+       outDim = topo[nbLayers-1];
+       net->in_rate = malloc((inDim+1)*sizeof(net->in_rate[0]));
+       net->weights = malloc((nbLayers-1)*sizeof(net->weights));
+       net->best_weights = malloc((nbLayers-1)*sizeof(net->weights));
+       for (i=0;i<nbLayers-1;i++)
+       {
+               net->weights[i] = malloc((topo[i]+1)*topo[i+1]*sizeof(net->weights[0][0]));
+               net->best_weights[i] = malloc((topo[i]+1)*topo[i+1]*sizeof(net->weights[0][0]));
+       }
+       double inMean[inDim];
+       for (j=0;j<inDim;j++)
+       {
+               double std=0;
+               inMean[j] = 0;
+               for (i=0;i<nbSamples;i++)
+               {
+                       inMean[j] += inputs[i*inDim+j];
+                       std += inputs[i*inDim+j]*inputs[i*inDim+j];
+               }
+               inMean[j] /= nbSamples;
+               std /= nbSamples;
+               net->in_rate[1+j] = .5/(.0001+std);
+               std = std-inMean[j]*inMean[j];
+               if (std<.001)
+                       std = .001;
+               std = 1/sqrt(inDim*std);
+               for (k=0;k<topo[1];k++)
+                       net->weights[0][k*(topo[0]+1)+j+1] = randn(std);
+       }
+       net->in_rate[0] = 1;
+       for (j=0;j<topo[1];j++)
+       {
+               double sum = 0;
+               for (k=0;k<inDim;k++)
+                       sum += inMean[k]*net->weights[0][j*(topo[0]+1)+k+1];
+               net->weights[0][j*(topo[0]+1)] = -sum;
+       }
+       for (j=0;j<outDim;j++)
+       {
+               double mean = 0;
+               double std;
+               for (i=0;i<nbSamples;i++)
+                       mean += outputs[i*outDim+j];
+               mean /= nbSamples;
+               std = 1/sqrt(topo[nbLayers-2]);
+               net->weights[nbLayers-2][j*(topo[nbLayers-2]+1)] = mean;
+               for (k=0;k<topo[nbLayers-2];k++)
+                       net->weights[nbLayers-2][j*(topo[nbLayers-2]+1)+k+1] = randn(std);
+       }
+       return net;
+}
+
+#define MAX_NEURONS 100
+
+double compute_gradient(MLPTrain *net, float *inputs, float *outputs, int nbSamples, double *W0_grad, double *W1_grad, double *error_rate)
+{
+       int i,j;
+       int s;
+       int inDim, outDim, hiddenDim;
+       int *topo;
+       double *W0, *W1;
+       double rms=0;
+       int W0_size, W1_size;
+       double hidden[MAX_NEURONS];
+       double netOut[MAX_NEURONS];
+       double error[MAX_NEURONS];
+
+        *error_rate = 0;
+       topo = net->topo;
+       inDim = net->topo[0];
+       hiddenDim = net->topo[1];
+       outDim = net->topo[2];
+       W0_size = (topo[0]+1)*topo[1];
+       W1_size = (topo[1]+1)*topo[2];
+       W0 = net->weights[0];
+       W1 = net->weights[1];
+       memset(W0_grad, 0, W0_size*sizeof(double));
+       memset(W1_grad, 0, W1_size*sizeof(double));
+       for (i=0;i<outDim;i++)
+               netOut[i] = outputs[i];
+       for (s=0;s<nbSamples;s++)
+       {
+               float *in, *out;
+               in = inputs+s*inDim;
+               out = outputs + s*outDim;
+               for (i=0;i<hiddenDim;i++)
+               {
+                       double sum = W0[i*(inDim+1)];
+                       for (j=0;j<inDim;j++)
+                               sum += W0[i*(inDim+1)+j+1]*in[j];
+                       hidden[i] = tansig_approx(sum);
+               }
+               for (i=0;i<outDim;i++)
+               {
+                       double sum = W1[i*(hiddenDim+1)];
+                       for (j=0;j<hiddenDim;j++)
+                               sum += W1[i*(hiddenDim+1)+j+1]*hidden[j];
+                       netOut[i] = tansig_approx(sum);
+                       error[i] = out[i] - netOut[i];
+                       rms += error[i]*error[i];
+                       *error_rate += fabs(error[i])>1;
+               }
+               /* Back-propagate error */
+               for (i=0;i<outDim;i++)
+               {
+                        float grad = 1-netOut[i]*netOut[i];
+                       W1_grad[i*(hiddenDim+1)] += error[i]*grad;
+                       for (j=0;j<hiddenDim;j++)
+                               W1_grad[i*(hiddenDim+1)+j+1] += grad*error[i]*hidden[j];
+               }
+               for (i=0;i<hiddenDim;i++)
+               {
+                       double grad;
+                       grad = 0;
+                       for (j=0;j<outDim;j++)
+                               grad += error[j]*W1[j*(hiddenDim+1)+i+1];
+                       grad *= 1-hidden[i]*hidden[i];
+                       W0_grad[i*(inDim+1)] += grad;
+                       for (j=0;j<inDim;j++)
+                               W0_grad[i*(inDim+1)+j+1] += grad*in[j];
+               }
+       }
+       return rms;
+}
+
+#define NB_THREADS 8
+
+sem_t sem_begin[NB_THREADS];
+sem_t sem_end[NB_THREADS];
+
+struct GradientArg {
+       int id;
+       int done;
+       MLPTrain *net;
+       float *inputs;
+       float *outputs;
+       int nbSamples;
+       double *W0_grad;
+       double *W1_grad;
+       double rms;
+       double error_rate;
+};
+
+void *gradient_thread_process(void *_arg)
+{
+       int W0_size, W1_size;
+       struct GradientArg *arg = _arg;
+       int *topo = arg->net->topo;
+       W0_size = (topo[0]+1)*topo[1];
+       W1_size = (topo[1]+1)*topo[2];
+       double W0_grad[W0_size];
+       double W1_grad[W1_size];
+       arg->W0_grad = W0_grad;
+       arg->W1_grad = W1_grad;
+       while (1)
+       {
+               sem_wait(&sem_begin[arg->id]);
+               if (arg->done)
+                       break;
+               arg->rms = compute_gradient(arg->net, arg->inputs, arg->outputs, arg->nbSamples, arg->W0_grad, arg->W1_grad, &arg->error_rate);
+               sem_post(&sem_end[arg->id]);
+       }
+       fprintf(stderr, "done\n");
+       return NULL;
+}
+
+float mlp_train_backprop(MLPTrain *net, float *inputs, float *outputs, int nbSamples, int nbEpoch, float rate)
+{
+       int i, j;
+       int e;
+       float best_rms = 1e10;
+       int inDim, outDim, hiddenDim;
+       int *topo;
+       double *W0, *W1, *best_W0, *best_W1;
+       double *W0_old, *W1_old;
+       double *W0_old2, *W1_old2;
+       double *W0_grad, *W1_grad;
+       double *W0_oldgrad, *W1_oldgrad;
+       double *W0_rate, *W1_rate;
+       double *best_W0_rate, *best_W1_rate;
+       int W0_size, W1_size;
+       topo = net->topo;
+       W0_size = (topo[0]+1)*topo[1];
+       W1_size = (topo[1]+1)*topo[2];
+       struct GradientArg args[NB_THREADS];
+       pthread_t thread[NB_THREADS];
+       int samplePerPart = nbSamples/NB_THREADS;
+       int count_worse=0;
+       int count_retries=0;
+
+       topo = net->topo;
+       inDim = net->topo[0];
+       hiddenDim = net->topo[1];
+       outDim = net->topo[2];
+       W0 = net->weights[0];
+       W1 = net->weights[1];
+       best_W0 = net->best_weights[0];
+       best_W1 = net->best_weights[1];
+       W0_old = malloc(W0_size*sizeof(double));
+       W1_old = malloc(W1_size*sizeof(double));
+       W0_old2 = malloc(W0_size*sizeof(double));
+       W1_old2 = malloc(W1_size*sizeof(double));
+       W0_grad = malloc(W0_size*sizeof(double));
+       W1_grad = malloc(W1_size*sizeof(double));
+       W0_oldgrad = malloc(W0_size*sizeof(double));
+       W1_oldgrad = malloc(W1_size*sizeof(double));
+       W0_rate = malloc(W0_size*sizeof(double));
+       W1_rate = malloc(W1_size*sizeof(double));
+       best_W0_rate = malloc(W0_size*sizeof(double));
+       best_W1_rate = malloc(W1_size*sizeof(double));
+       memcpy(W0_old, W0, W0_size*sizeof(double));
+       memcpy(W0_old2, W0, W0_size*sizeof(double));
+       memset(W0_grad, 0, W0_size*sizeof(double));
+       memset(W0_oldgrad, 0, W0_size*sizeof(double));
+       memcpy(W1_old, W1, W1_size*sizeof(double));
+       memcpy(W1_old2, W1, W1_size*sizeof(double));
+       memset(W1_grad, 0, W1_size*sizeof(double));
+       memset(W1_oldgrad, 0, W1_size*sizeof(double));
+       
+       rate /= nbSamples;
+       for (i=0;i<hiddenDim;i++)
+               for (j=0;j<inDim+1;j++)
+                       W0_rate[i*(inDim+1)+j] = rate*net->in_rate[j];
+       for (i=0;i<W1_size;i++)
+               W1_rate[i] = rate;
+       
+       for (i=0;i<NB_THREADS;i++)
+       {
+               args[i].net = net;
+               args[i].inputs = inputs+i*samplePerPart*inDim;
+               args[i].outputs = outputs+i*samplePerPart*outDim;
+               args[i].nbSamples = samplePerPart;
+               args[i].id = i;
+               args[i].done = 0;
+               sem_init(&sem_begin[i], 0, 0);
+               sem_init(&sem_end[i], 0, 0);
+               pthread_create(&thread[i], NULL, gradient_thread_process, &args[i]);
+       }
+       for (e=0;e<nbEpoch;e++)
+       {
+               double rms=0;
+                double error_rate = 0;
+               for (i=0;i<NB_THREADS;i++)
+               {
+                       sem_post(&sem_begin[i]);
+               }
+               memset(W0_grad, 0, W0_size*sizeof(double));
+               memset(W1_grad, 0, W1_size*sizeof(double));
+               for (i=0;i<NB_THREADS;i++)
+               {
+                       sem_wait(&sem_end[i]);
+                       rms += args[i].rms;
+                       error_rate += args[i].error_rate;
+                       for (j=0;j<W0_size;j++)
+                               W0_grad[j] += args[i].W0_grad[j];
+                       for (j=0;j<W1_size;j++)
+                               W1_grad[j] += args[i].W1_grad[j];
+               }
+
+               float mean_rate = 0, min_rate = 1e10;
+               rms = (rms/(outDim*nbSamples));
+               error_rate = (error_rate/(outDim*nbSamples));
+               fprintf (stderr, "%f (%f %f) ", error_rate, rms, best_rms);
+               if (rms < best_rms)
+               {
+                       best_rms = rms;
+                       for (i=0;i<W0_size;i++)
+                       {
+                               best_W0[i] = W0[i];
+                               best_W0_rate[i] = W0_rate[i];
+                       }
+                       for (i=0;i<W1_size;i++)
+                       {
+                               best_W1[i] = W1[i];
+                               best_W1_rate[i] = W1_rate[i];
+                       }
+                       count_worse=0;
+                       count_retries=0;
+               } else {
+                       count_worse++;
+                       if (count_worse>30)
+                       {
+                           count_retries++;
+                               count_worse=0;
+                               for (i=0;i<W0_size;i++)
+                               {
+                                       W0[i] = best_W0[i];
+                                       best_W0_rate[i] *= .7;
+                                       if (best_W0_rate[i]<1e-15) best_W0_rate[i]=1e-15;
+                                       W0_rate[i] = best_W0_rate[i];
+                                       W0_grad[i] = 0;
+                               }
+                               for (i=0;i<W1_size;i++)
+                               {
+                                       W1[i] = best_W1[i];
+                                       best_W1_rate[i] *= .8;
+                                       if (best_W1_rate[i]<1e-15) best_W1_rate[i]=1e-15;
+                                       W1_rate[i] = best_W1_rate[i];
+                                       W1_grad[i] = 0;
+                               }
+                       }
+               }
+               if (count_retries>10)
+                   break;
+               for (i=0;i<W0_size;i++)
+               {
+                       if (W0_oldgrad[i]*W0_grad[i] > 0)
+                               W0_rate[i] *= 1.01;
+                       else if (W0_oldgrad[i]*W0_grad[i] < 0)
+                               W0_rate[i] *= .9;
+                       mean_rate += W0_rate[i];
+                       if (W0_rate[i] < min_rate)
+                               min_rate = W0_rate[i];
+                       if (W0_rate[i] < 1e-15)
+                               W0_rate[i] = 1e-15;
+                       /*if (W0_rate[i] > .01)
+                               W0_rate[i] = .01;*/
+                       W0_oldgrad[i] = W0_grad[i];
+                       W0_old2[i] = W0_old[i];
+                       W0_old[i] = W0[i];
+                       W0[i] += W0_grad[i]*W0_rate[i];
+               }
+               for (i=0;i<W1_size;i++)
+               {
+                       if (W1_oldgrad[i]*W1_grad[i] > 0)
+                               W1_rate[i] *= 1.01;
+                       else if (W1_oldgrad[i]*W1_grad[i] < 0)
+                               W1_rate[i] *= .9;
+                       mean_rate += W1_rate[i];
+                       if (W1_rate[i] < min_rate)
+                               min_rate = W1_rate[i];
+                       if (W1_rate[i] < 1e-15)
+                               W1_rate[i] = 1e-15;
+                       W1_oldgrad[i] = W1_grad[i];
+                       W1_old2[i] = W1_old[i];
+                       W1_old[i] = W1[i];
+                       W1[i] += W1_grad[i]*W1_rate[i];
+               }
+               mean_rate /= (topo[0]+1)*topo[1] + (topo[1]+1)*topo[2];
+               fprintf (stderr, "%g %d", mean_rate, e);
+               if (count_retries)
+                   fprintf(stderr, " %d", count_retries);
+               fprintf(stderr, "\n");
+               if (stopped)
+                       break;
+       }
+       for (i=0;i<NB_THREADS;i++)
+       {
+               args[i].done = 1;
+               sem_post(&sem_begin[i]);
+               pthread_join(thread[i], NULL);
+               fprintf (stderr, "joined %d\n", i);
+       }
+       free(W0_old);
+       free(W1_old);
+       free(W0_grad);
+       free(W1_grad);
+       free(W0_rate);
+       free(W1_rate);
+       return best_rms;
+}
+
+int main(int argc, char **argv)
+{
+       int i, j;
+       int nbInputs;
+       int nbOutputs;
+       int nbHidden;
+       int nbSamples;
+       int nbEpoch;
+       int nbRealInputs;
+       unsigned int seed;
+       int ret;
+       float rms;
+       float *inputs;
+       float *outputs;
+       if (argc!=6)
+       {
+               fprintf (stderr, "usage: mlp_train <inputs> <hidden> <outputs> <nb samples> <nb epoch>\n");
+               return 1;
+       }
+       nbInputs = atoi(argv[1]);
+       nbHidden = atoi(argv[2]);
+       nbOutputs = atoi(argv[3]);
+       nbSamples = atoi(argv[4]);
+       nbEpoch = atoi(argv[5]);
+       nbRealInputs = nbInputs;
+       inputs = malloc(nbInputs*nbSamples*sizeof(*inputs));
+       outputs = malloc(nbOutputs*nbSamples*sizeof(*outputs));
+       
+       seed = time(NULL);
+       fprintf (stderr, "Seed is %u\n", seed);
+       srand(seed);
+       build_tansig_table();
+       signal(SIGTERM, handler);
+       signal(SIGINT, handler);
+       signal(SIGHUP, handler);
+       for (i=0;i<nbSamples;i++)
+       {
+               for (j=0;j<nbRealInputs;j++)
+                       ret = scanf(" %f", &inputs[i*nbInputs+j]);
+               for (j=0;j<nbOutputs;j++)
+                       ret = scanf(" %f", &outputs[i*nbOutputs+j]);
+               if (feof(stdin))
+               {
+                       nbSamples = i;
+                       break;
+               }
+       }
+       int topo[3] = {nbInputs, nbHidden, nbOutputs};
+       MLPTrain *net;
+
+       fprintf (stderr, "Got %d samples\n", nbSamples);
+       net = mlp_init(topo, 3, inputs, outputs, nbSamples);
+       rms = mlp_train_backprop(net, inputs, outputs, nbSamples, nbEpoch, 1);
+       printf ("#include \"mlp.h\"\n\n");
+       printf ("/* RMS error was %f, seed was %u */\n\n", rms, seed);
+       printf ("static const float weights[%d] = {\n", (topo[0]+1)*topo[1] + (topo[1]+1)*topo[2]);
+       printf ("\n/* hidden layer */\n");
+       for (i=0;i<(topo[0]+1)*topo[1];i++)
+       {
+               printf ("%g, ", net->weights[0][i]);
+               if (i%5==4)
+                       printf("\n");
+       }
+       printf ("\n/* output layer */\n");
+       for (i=0;i<(topo[1]+1)*topo[2];i++)
+       {
+               printf ("%g, ", net->weights[1][i]);
+               if (i%5==4)
+                       printf("\n");
+       }
+       printf ("};\n\n");
+       printf ("static const int topo[3] = {%d, %d, %d};\n\n", topo[0], topo[1], topo[2]);
+       printf ("const MLP net = {\n");
+       printf ("\t3,\n");
+       printf ("\ttopo,\n");
+       printf ("\tweights\n};\n");
+       return 0;
+}
diff --git a/src/mlp_train.h b/src/mlp_train.h
new file mode 100644 (file)
index 0000000..1857f64
--- /dev/null
@@ -0,0 +1,86 @@
+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef _MLP_TRAIN_H_
+#define _MLP_TRAIN_H_
+
+#include <math.h>
+#include <stdlib.h>
+
+double tansig_table[501];
+static inline double tansig_double(double x) 
+{
+       return 2./(1.+exp(-2.*x)) - 1.;
+}
+static inline void build_tansig_table()
+{
+       int i;
+       for (i=0;i<501;i++)
+               tansig_table[i] = tansig_double(.04*(i-250));
+}
+
+static inline double tansig_approx(double x)
+{
+       int i;
+       double y, dy;
+       if (x>=10)
+               return 1;
+       if (x<=-10)
+               return -1;
+       i = lrint(25*x);
+       x -= .04*i;
+       y = tansig_table[250+i];
+       dy = 1-y*y;
+       y = y + x*dy*(1 - y*x);
+       return y;
+}
+
+inline float randn(float sd)
+{
+   float U1, U2, S, x;
+   do {
+      U1 = ((float)rand())/RAND_MAX;
+      U2 = ((float)rand())/RAND_MAX;
+      U1 = 2*U1-1;
+      U2 = 2*U2-1;
+      S = U1*U1 + U2*U2;
+   } while (S >= 1 || S == 0.0f);
+   x = sd*sqrt(-2 * log(S) / S) * U1;
+   return x;
+}
+
+
+typedef struct {
+       int layers;
+       int *topo;
+       double **weights;
+       double **best_weights;
+       double *in_rate;
+} MLPTrain;
+
+
+#endif /* _MLP_TRAIN_H_ */
diff --git a/src/tansig_table.h b/src/tansig_table.h
new file mode 100644 (file)
index 0000000..ccf43da
--- /dev/null
@@ -0,0 +1,105 @@
+/* This file is auto-generated by gen_tables */
+
+static const opus_val16 tansig_table[501] = {
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -0.999999, -0.999999, -0.999999, -0.999999, 
+-0.999999, -0.999999, -0.999999, -0.999999, -0.999999, 
+-0.999999, -0.999999, -0.999999, -0.999999, -0.999998, 
+-0.999998, -0.999998, -0.999998, -0.999998, -0.999998, 
+-0.999997, -0.999997, -0.999997, -0.999997, -0.999997, 
+-0.999996, -0.999996, -0.999996, -0.999995, -0.999995, 
+-0.999994, -0.999994, -0.999994, -0.999993, -0.999992, 
+-0.999992, -0.999991, -0.999990, -0.999990, -0.999989, 
+-0.999988, -0.999987, -0.999986, -0.999984, -0.999983, 
+-0.999982, -0.999980, -0.999978, -0.999977, -0.999975, 
+-0.999973, -0.999970, -0.999968, -0.999965, -0.999962, 
+-0.999959, -0.999956, -0.999952, -0.999948, -0.999944, 
+-0.999939, -0.999934, -0.999929, -0.999923, -0.999916, 
+-0.999909, -0.999902, -0.999893, -0.999885, -0.999875, 
+-0.999865, -0.999853, -0.999841, -0.999828, -0.999813, 
+-0.999798, -0.999781, -0.999763, -0.999743, -0.999722, 
+-0.999699, -0.999673, -0.999646, -0.999617, -0.999585, 
+-0.999550, -0.999513, -0.999472, -0.999428, -0.999381, 
+-0.999329, -0.999273, -0.999213, -0.999147, -0.999076, 
+-0.999000, -0.998916, -0.998826, -0.998728, -0.998623, 
+-0.998508, -0.998384, -0.998249, -0.998104, -0.997946, 
+-0.997775, -0.997590, -0.997389, -0.997172, -0.996937, 
+-0.996682, -0.996407, -0.996108, -0.995784, -0.995434, 
+-0.995055, -0.994644, -0.994199, -0.993718, -0.993196, 
+-0.992631, -0.992020, -0.991359, -0.990642, -0.989867, 
+-0.989027, -0.988119, -0.987136, -0.986072, -0.984921, 
+-0.983675, -0.982327, -0.980869, -0.979293, -0.977587, 
+-0.975743, -0.973749, -0.971594, -0.969265, -0.966747, 
+-0.964028, -0.961090, -0.957917, -0.954492, -0.950795, 
+-0.946806, -0.942503, -0.937863, -0.932862, -0.927473, 
+-0.921669, -0.915420, -0.908698, -0.901468, -0.893698, 
+-0.885352, -0.876393, -0.866784, -0.856485, -0.845456, 
+-0.833655, -0.821040, -0.807569, -0.793199, -0.777888, 
+-0.761594, -0.744277, -0.725897, -0.706419, -0.685809, 
+-0.664037, -0.641077, -0.616909, -0.591519, -0.564900, 
+-0.537050, -0.507977, -0.477700, -0.446244, -0.413644, 
+-0.379949, -0.345214, -0.309507, -0.272905, -0.235496, 
+-0.197375, -0.158649, -0.119427, -0.079830, -0.039979, 
+0.000000, 0.039979, 0.079830, 0.119427, 0.158649, 
+0.197375, 0.235496, 0.272905, 0.309507, 0.345214, 
+0.379949, 0.413644, 0.446244, 0.477700, 0.507977, 
+0.537050, 0.564900, 0.591519, 0.616909, 0.641077, 
+0.664037, 0.685809, 0.706419, 0.725897, 0.744277, 
+0.761594, 0.777888, 0.793199, 0.807569, 0.821040, 
+0.833655, 0.845456, 0.856485, 0.866784, 0.876393, 
+0.885352, 0.893698, 0.901468, 0.908698, 0.915420, 
+0.921669, 0.927473, 0.932862, 0.937863, 0.942503, 
+0.946806, 0.950795, 0.954492, 0.957917, 0.961090, 
+0.964028, 0.966747, 0.969265, 0.971594, 0.973749, 
+0.975743, 0.977587, 0.979293, 0.980869, 0.982327, 
+0.983675, 0.984921, 0.986072, 0.987136, 0.988119, 
+0.989027, 0.989867, 0.990642, 0.991359, 0.992020, 
+0.992631, 0.993196, 0.993718, 0.994199, 0.994644, 
+0.995055, 0.995434, 0.995784, 0.996108, 0.996407, 
+0.996682, 0.996937, 0.997172, 0.997389, 0.997590, 
+0.997775, 0.997946, 0.998104, 0.998249, 0.998384, 
+0.998508, 0.998623, 0.998728, 0.998826, 0.998916, 
+0.999000, 0.999076, 0.999147, 0.999213, 0.999273, 
+0.999329, 0.999381, 0.999428, 0.999472, 0.999513, 
+0.999550, 0.999585, 0.999617, 0.999646, 0.999673, 
+0.999699, 0.999722, 0.999743, 0.999763, 0.999781, 
+0.999798, 0.999813, 0.999828, 0.999841, 0.999853, 
+0.999865, 0.999875, 0.999885, 0.999893, 0.999902, 
+0.999909, 0.999916, 0.999923, 0.999929, 0.999934, 
+0.999939, 0.999944, 0.999948, 0.999952, 0.999956, 
+0.999959, 0.999962, 0.999965, 0.999968, 0.999970, 
+0.999973, 0.999975, 0.999977, 0.999978, 0.999980, 
+0.999982, 0.999983, 0.999984, 0.999986, 0.999987, 
+0.999988, 0.999989, 0.999990, 0.999990, 0.999991, 
+0.999992, 0.999992, 0.999993, 0.999994, 0.999994, 
+0.999994, 0.999995, 0.999995, 0.999996, 0.999996, 
+0.999996, 0.999997, 0.999997, 0.999997, 0.999997, 
+0.999997, 0.999998, 0.999998, 0.999998, 0.999998, 
+0.999998, 0.999998, 0.999999, 0.999999, 0.999999, 
+0.999999, 0.999999, 0.999999, 0.999999, 0.999999, 
+0.999999, 0.999999, 0.999999, 0.999999, 0.999999, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 
+};