text:
/*
 *  maxiMFCC.cpp
 *  mfccs
 *
 *  Created by Chris on 08/03/2011.
 *  Copyright 2011 Goldsmiths Creative Computing. All rights reserved.
 *
 */

#include "maxiMFCC.h"

#ifdef __APPLE_CC__
template <>
void maxiMFCCAnalyser<double>::dct(double *mfccs) {
	vDSP_mmulD(melBands, 1, dctMatrix, 1, mfccs, 1, 1, numCoeffs, numFilters);
	double n = (double) numCoeffs;
	vDSP_vsdivD(mfccs, 1, &n, mfccs, 1, numCoeffs);
}

template <>
void maxiMFCCAnalyser<float>::dct(float *mfccs) {
	vDSP_mmul(melBands, 1, dctMatrix, 1, mfccs, 1, 1, numCoeffs, numFilters);
	float n = (float) numCoeffs;
	vDSP_vsdiv(mfccs, 1, &n, mfccs, 1, numCoeffs);
}
#endif

template <>
void maxiMFCCAnalyser<double>::melFilterAndLogSquare(float* powerSpectrum) {
#ifdef __APPLE_CC__
	//conv to double
	vDSP_vspdp(powerSpectrum, 1, doubleSpec, 1, numBins);
	//	vDSP_mmulD(doubleSpec, 1, melFilters, 1, melBands, 1, 1, 42, 512);
	vDSP_mmulD(doubleSpec, 1, melFilters, 1, melBands, 1, 1, numFilters, numBins);
#endif
	melFilterAndLogSq_Part2(powerSpectrum);
}

template <>
void maxiMFCCAnalyser<float>::melFilterAndLogSquare(float* powerSpectrum) {
#ifdef __APPLE_CC__
	vDSP_mmul(powerSpectrum, 1, melFilters, 1, melBands, 1, 1, numFilters, numBins);
#endif
	melFilterAndLogSq_Part2(powerSpectrum);
}

template <class T>
void maxiMFCCAnalyser<T>::melFilterAndLogSq_Part2(float *powerSpectrum) {
#ifdef __APPLE_CC__
#else
	for (unsigned int filter = 0;filter < numFilters;filter++) {
		melBands[filter] = 0.0;
		for (unsigned int bin=0;bin<numBins;bin++) {
			//			int idx = (numBins * filter) + bin;
			int idx = filter + (bin * numFilters);
			//			melBands[filter] += (melFilters[filter][bin] * powerSpectrum[bin]);
			melBands[filter] += (melFilters[idx] * powerSpectrum[bin]);
		}
	}
#endif
	for(unsigned int filter=0; filter < numFilters; filter++) {
		// log the square
		melBands[filter] = melBands[filter] > 0.000001 ? log(melBands[filter] * melBands[filter]) : 0.0;
	}
}


text:
/*
 *  maxiMFCC.h
 *  mfccs
 *
 *  Created by Chris on 08/03/2011.
 *  Copyright 2011 Goldsmiths Creative Computing. All rights reserved.
 *

 Based on Matthew Yee-King's MFCCMYK java class
 */

#pragma once
#pragma pack(16)

#include "maxiFFT.h"
#include <math.h>
#include <iostream>
#include <cstdlib>
#ifdef __APPLE_CC__
#include <Accelerate/Accelerate.h>
#endif

using namespace std;

// implements this formula:
// mel = 2595 log10(Hz/700 + 1)
inline double hzToMel(double hz){
	return 2595.0 * (log10(hz/700.0 + 1.0));
}

// implements this formula
// Hz = 700 (10^(mel/2595) - 1)
inline double melToHz(double mel){
	return 700.0 * (pow(10, mel/2595.0) - 1.0);
}

template <class T>
class maxiMFCCAnalyser {
public:
	T *melBands;
	maxiMFCCAnalyser():melFilters(NULL),dctMatrix(NULL), melBands(NULL){};
	~maxiMFCCAnalyser() {
		if (melFilters) {
			delete[] melFilters;
			delete[] melBands;
			delete[] dctMatrix;
#ifdef __APPLE_CC__
			delete doubleSpec;
#endif
		}
	}

	void setup(unsigned int numBins, unsigned int numFilters, unsigned int numCoeffs, double minFreq, double maxFreq, unsigned int sampleRate)
	{
		this->numFilters = numFilters;
		this->numCoeffs = numCoeffs;
		this->minFreq = minFreq;
		this->maxFreq = maxFreq;
		this->sampleRate = sampleRate;
		this->numBins = numBins;
		melFilters = NULL;
		melBands = (T*) malloc(sizeof(T) * numFilters);
#ifdef __APPLE_CC__
		doubleSpec = (T*)malloc(sizeof(T) * numBins);
#endif
		//create new matrix
		dctMatrix = (T*)malloc(sizeof(T) * numCoeffs * numFilters);
		calcMelFilterBank(sampleRate, numBins);
		createDCTCoeffs();
	}
	void mfcc(float* powerSpectrum, T *mfccs) {
		melFilterAndLogSquare(powerSpectrum);
		dct(mfccs);
	}

private:
	unsigned int numFilters, numCoeffs;
	double minFreq, maxFreq;
	unsigned int sampleRate;
	T *melFilters;
	unsigned int numBins;
	T *dctMatrix;
#ifdef __APPLE_CC__
	T *doubleSpec;
#endif

#ifdef __APPLE_CC__
	void dct(T *mfccs); //define later
#else
	void dct(T *mfccs) {
		for(int i=0; i < numCoeffs; i++) {
			mfccs[i] = 0.0;
		}
		for(int i=0; i < numCoeffs; i++ ) {
			for(int j=0; j < numFilters; j++) {
				int idx = i + (j * numCoeffs);
				mfccs[i] += (dctMatrix[idx] * melBands[j]);
			}
		}
		for(int i=0; i < numCoeffs; i++) {
			mfccs[i] /= numCoeffs;
		}
	}
#endif

	void melFilterAndLogSquare(float* powerSpectrum);
	void melFilterAndLogSq_Part2(float *powerSpectrum);

	void calcMelFilterBank(double sampleRate, int numBins) {

		double mel, dMel, maxMel, minMel, nyquist, binFreq, start, end, thisF, nextF, prevF;
		int numValidBins;

		// ignore bins over nyquist
		numValidBins = numBins;

		nyquist = sampleRate/2;
		if (maxFreq > nyquist) {
			maxFreq = nyquist;
		}

		maxMel = hzToMel(maxFreq);
		minMel = hzToMel(minFreq);

		dMel = (maxMel - minMel) / (numFilters + 2 - 1);

		T *filtPos = (T*) malloc(sizeof(double) * (numFilters + 2));

		// first generate an array of start and end freqs for each triangle
		mel = minMel;
		for (int i=0;i<numFilters + 2;i++) {
			// start of the triangle
			filtPos[i] = melToHz(mel);
			//		std::cout << "[" << i << "] MFCC: centre is at " <<filtPos[i]<<"hz "<<mel<<" mels" << endl;
			mel += dMel;
		}
		// now generate the coefficients for the mag spectrum
		melFilters = (T*) malloc(sizeof(T) * numFilters * numValidBins);

		for (int filter = 1; filter < numFilters; filter++) {
			for (int bin=0;bin<numValidBins;bin++) {
				// frequency this bin represents
				binFreq = (T) sampleRate / (T) numValidBins * (T) bin;
				thisF = filtPos[filter];
				nextF = filtPos[filter+1];
				prevF = filtPos[filter-1];
				int idx = filter + (bin * numFilters);
				if (binFreq > nextF || binFreq < prevF) {
					// outside this filter
					melFilters[idx] = 0;
					//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
				}
				else {
					T height = 2.0 / (nextF - prevF);

					if (binFreq < thisF) {
						// up
						start = prevF;
						end = thisF;
						melFilters[idx] = (binFreq - start) * (height / (thisF - start));
					}
					else {
						// down
						start = thisF;
						end = nextF;
						melFilters[idx] = height + ((binFreq - thisF) * (-height /(nextF - thisF)));
					}
					//				cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[filter][bin] << endl;
					//cout << "MFCCMYK: filter at " <<thisF << " bin at " <<binFreq <<" coeff " <<melFilters[idx] << endl;
				}
			}
		}
	}
	void createDCTCoeffs() {
		T k = 3.14159265358979323846/numFilters;
		T w1 = 1.0/(sqrt(numFilters));
		T w2 = sqrt(2.0/numFilters);

		//generate dct matrix
		for(int i = 0; i < numCoeffs; i++)
		{
			for(int j = 0; j < numFilters; j++)
			{
				int idx = i + (j * numCoeffs);
				if(i == 0)
					dctMatrix[idx]= w1 * cos(k * (i+1) * (j + 0.5));
				else
					dctMatrix[idx] = w2 * cos(k * (i+1) * (j + 0.5));
			}
		}

	}

};

typedef maxiMFCCAnalyser<double> maxiMFCC;
//typedef maxiMFCCAnalyser<float> maxiFloatMFCC;

repo_name: fedden/RenderMan

mime_type: text/x-c++


text:
/*
 *  maximilian
 *  platform independent synthesis library using portaudio or rtaudio
 *
 *  Created by Mick Grierson on 29/12/2009.
 *  Copyright 2009 Mick Grierson & Strangeloop Limited. All rights reserved.
 *	Thanks to the Goldsmiths Creative Computing Team.
 *	Special thanks to Arturo Castro for the PortAudio implementation.
 *
 *	Permission is hereby granted, free of charge, to any person
 *	obtaining a copy of this software and associated documentation
 *	files (the "Software"), to deal in the Software without
 *	restriction, including without limitation the rights to use,
 *	copy, modify, merge, publish, distribute, sublicense, and/or sell
 *	copies of the Software, and to permit persons to whom the
 *	Software is furnished to do so, subject to the following
 *	conditions:
 *
 *	The above copyright notice and this permission notice shall be
 *	included in all copies or substantial portions of the Software.
 *
 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 *	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 *	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 *	NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 *	HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 *	WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 *	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 *	OTHER DEALINGS IN THE SOFTWARE.
 *
 */
/*

 fft.cpp

 Based on K+R Numerical recipes in C and some other stuff hacked about.

 */

#include "fft.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

int **gFFTBitTable = NULL;
const int MaxFastBits = 16;

int IsPowerOfTwo(int x)
{
	if (x < 2)
		return false;

	if (x & (x - 1))
		return false;

	return true;
}

int NumberOfBitsNeeded(int PowerOfTwo)
{
	int i;

	if (PowerOfTwo < 2) {
		fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
		exit(1);
	}

	for (i = 0;; i++)
		if (PowerOfTwo & (1 << i))
			return i;
}

int ReverseBits(int index, int NumBits)
{
	int i, rev;

	for (i = rev = 0; i < NumBits; i++) {
		rev = (rev << 1) | (index & 1);
		index >>= 1;
	}

	return rev;
}

void InitFFT()
{
    //	gFFTBitTable = new int *[MaxFastBits];
	//use malloc for 16 byte alignment
	gFFTBitTable = (int**) malloc(MaxFastBits * sizeof(int*));

	int len = 2;
	for (int b = 1; b <= MaxFastBits; b++) {

        //		gFFTBitTable[b - 1] = new int[len];
		gFFTBitTable[b - 1] = (int*) malloc(len * sizeof(int));

		for (int i = 0; i < len; i++)
			gFFTBitTable[b - 1][i] = ReverseBits(i, b);

		len <<= 1;
	}
}

inline int FastReverseBits(int i, int NumBits)
{
	if (NumBits <= MaxFastBits)
		return gFFTBitTable[NumBits - 1][i];
	else
		return ReverseBits(i, NumBits);
}

/*
 * Complex Fast Fourier Transform
 */

void FFT(int NumSamples,
         bool InverseTransform,
         float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
{
	int NumBits;                 /* Number of bits needed to store indices */
	int i, j, k, n;
	int BlockSize, BlockEnd;

	double angle_numerator = 2.0 * M_PI;
	float tr, ti;                /* temp real, temp imaginary */

	if (!IsPowerOfTwo(NumSamples)) {
		fprintf(stderr, "%d is not a power of two\n", NumSamples);
		exit(1);
	}

	if (!gFFTBitTable)
		InitFFT();

	if (InverseTransform)
		angle_numerator = -angle_numerator;

	NumBits = NumberOfBitsNeeded(NumSamples);

	/*
	 **   Do simultaneous data copy and bit-reversal ordering into outputs...
	 */

	for (i = 0; i < NumSamples; i++) {
		j = FastReverseBits(i, NumBits);
		RealOut[j] = RealIn[i];
		ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
	}

	/*
	 **   Do the FFT itself...
	 */

	BlockEnd = 1;
	for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {

		double delta_angle = angle_numerator / (double) BlockSize;

		float sm2 = sin(-2 * delta_angle);
		float sm1 = sin(-delta_angle);
		float cm2 = cos(-2 * delta_angle);
		float cm1 = cos(-delta_angle);
		float w = 2 * cm1;
		float ar0, ar1, ar2, ai0, ai1, ai2;

		for (i = 0; i < NumSamples; i += BlockSize) {
			ar2 = cm2;
			ar1 = cm1;

			ai2 = sm2;
			ai1 = sm1;

			for (j = i, n = 0; n < BlockEnd; j++, n++) {
				ar0 = w * ar1 - ar2;
				ar2 = ar1;
				ar1 = ar0;

				ai0 = w * ai1 - ai2;
				ai2 = ai1;
				ai1 = ai0;

				k = j + BlockEnd;
				tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
				ti = ar0 * ImagOut[k] + ai0 * RealOut[k];

				RealOut[k] = RealOut[j] - tr;
				ImagOut[k] = ImagOut[j] - ti;

				RealOut[j] += tr;
				ImagOut[j] += ti;
			}
		}

		BlockEnd = BlockSize;
	}

	/*
	 **   Need to normalize if inverse transform...
	 */

	if (InverseTransform) {
		float denom = (float) NumSamples;

		for (i = 0; i < NumSamples; i++) {
			RealOut[i] /= denom;
			ImagOut[i] /= denom;
		}
	}
}

/*
 * Real Fast Fourier Transform
 *
 * This function was based on the code in Numerical Recipes in C.
 * In Num. Rec., the inner loop is based on a single 1-based array
 * of interleaved real and imaginary numbers.  Because we have two
 * separate zero-based arrays, our indices are quite different.
 * Here is the correspondence between Num. Rec. indices and our indices:
 *
 * i1  <->  real[i]
 * i2  <->  imag[i]
 * i3  <->  real[n/2-i]
 * i4  <->  imag[n/2-i]
 */

void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
{
	int Half = NumSamples / 2;
	int i;

	float theta = M_PI / Half;

	float *tmpReal = (float*) malloc(Half * sizeof(float));
	float *tmpImag = (float*) malloc(Half * sizeof(float));

	for (i = 0; i < Half; i++) {
		tmpReal[i] = RealIn[2 * i];
		tmpImag[i] = RealIn[2 * i + 1];
	}

	FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);

	float wtemp = float (sin(0.5 * theta));

	float wpr = -2.0 * wtemp * wtemp;
	float wpi = float (sin(theta));
	float wr = 1.0 + wpr;
	float wi = wpi;

	int i3;

	float h1r, h1i, h2r, h2i;

	for (i = 1; i < Half / 2; i++) {

		i3 = Half - i;

		h1r = 0.5 * (RealOut[i] + RealOut[i3]);
		h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
		h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
		h2i = -0.5 * (RealOut[i] - RealOut[i3]);

		RealOut[i] = h1r + wr * h2r - wi * h2i;
		ImagOut[i] = h1i + wr * h2i + wi * h2r;
		RealOut[i3] = h1r - wr * h2r + wi * h2i;
		ImagOut[i3] = -h1i + wr * h2i + wi * h2r;

		wr = (wtemp = wr) * wpr - wi * wpi + wr;
		wi = wi * wpr + wtemp * wpi + wi;
	}

	RealOut[0] = (h1r = RealOut[0]) + ImagOut[0];
	ImagOut[0] = h1r - ImagOut[0];

	free(tmpReal);
	free(tmpImag);
}

/*
 * PowerSpectrum
 *
 * This function computes the same as RealFFT, above, but
 * adds the squares of the real and imaginary part of each
 * coefficient, extracting the power and throwing away the
 * phase.
 *
 * For speed, it does not call RealFFT, but duplicates some
 * of its code.
 */

void PowerSpectrum(int NumSamples, float *In, float *Out)
{
	int Half = NumSamples / 2;
	int i;

	float theta = M_PI / Half;

	float *tmpReal = new float[Half];
	float *tmpImag = new float[Half];
	float *RealOut = new float[Half];
	float *ImagOut = new float[Half];

	for (i = 0; i < Half; i++) {
		tmpReal[i] = In[2 * i];
		tmpImag[i] = In[2 * i + 1];
	}

	FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);

	float wtemp = float (sin(0.5 * theta));

	float wpr = -2.0 * wtemp * wtemp;
	float wpi = float (sin(theta));
	float wr = 1.0 + wpr;
	float wi = wpi;

	int i3;

	float h1r, h1i, h2r, h2i, rt, it;
	//float total=0;

	for (i = 1; i < Half / 2; i++) {

		i3 = Half - i;

		h1r = 0.5 * (RealOut[i] + RealOut[i3]);
		h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
		h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
		h2i = -0.5 * (RealOut[i] - RealOut[i3]);

		rt = h1r + wr * h2r - wi * h2i; //printf("Realout%i = %f",i,rt);total+=fabs(rt);
		it = h1i + wr * h2i + wi * h2r; // printf("  Imageout%i = %f\n",i,it);

		Out[i] = rt * rt + it * it;

		rt = h1r - wr * h2r + wi * h2i;
		it = -h1i + wr * h2i + wi * h2r;

		Out[i3] = rt * rt + it * it;

		wr = (wtemp = wr) * wpr - wi * wpi + wr;
		wi = wi * wpr + wtemp * wpi + wi;
	}
	//printf("total = %f\n",total);
	rt = (h1r = RealOut[0]) + ImagOut[0];
	it = h1r - ImagOut[0];
	Out[0] = rt * rt + it * it;

	rt = RealOut[Half / 2];
	it = ImagOut[Half / 2];
	Out[Half / 2] = rt * rt + it * it;

	delete[]tmpReal;
	delete[]tmpImag;
	delete[]RealOut;
	delete[]ImagOut;
}

void WindowFunc(int whichFunction, int NumSamples, float *in)
{
	int i;

	if (whichFunction == 1) {
		// Bartlett (triangular) window
		for (i = 0; i < NumSamples / 2; i++) {
			in[i] *= (i / (float) (NumSamples / 2));
			in[i + (NumSamples / 2)] *=
			(1.0 - (i / (float) (NumSamples / 2)));
		}
	}

	if (whichFunction == 2) {
		// Hamming
		for (i = 0; i < NumSamples; i++)
			in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
	}

	if (whichFunction == 3) {
		// Hanning
		for (i = 0; i < NumSamples; i++)
			in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
	}
}

void fft::genWindow(int whichFunction, int NumSamples, float *window)
{
	int i;

	if (whichFunction == 1) {
		// Bartlett (triangular) window
		for (i = 0; i < NumSamples / 2; i++) {
			window[i] = (i / (float) (NumSamples / 2));
			window[i + (NumSamples / 2)] =
			(1.0 - (i / (float) (NumSamples / 2)));
		}
	}

	if (whichFunction == 2) {
		// Hamming
		for (i = 0; i < NumSamples; i++)
			window[i] = 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
	}

	if (whichFunction == 3) {
		// Hanning
		for (i = 0; i < NumSamples; i++)
			window[i] = 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
	}
}

/* constructor */
fft::fft(int fftSize) {
	n = fftSize;
	half = fftSize / 2;
	//use malloc for 16 byte alignment
	in_real = (float *) malloc(n * sizeof(float));
	in_img = (float *) malloc(n * sizeof(float));
	out_real = (float *) malloc(n * sizeof(float));
	out_img = (float *) malloc(n * sizeof(float));

#ifdef __APPLE_CC__
	log2n = log2(n);
	A.realp = (float *) malloc(half * sizeof(float));
    A.imagp = (float *) malloc(half * sizeof(float));
	setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);
    if (setupReal == NULL) {
        printf("\nFFT_Setup failed to allocate enough memory  for"
               "the real FFT.\n");
    }
	polar = (float *) malloc(n * sizeof(float));
#endif
}

/* destructor */
fft::~fft() {
	delete[] in_real, out_real, in_img, out_img;
#ifdef __APPLE_CC__
	vDSP_destroy_fftsetup(setupReal);
    delete[] A.realp;
    delete[] A.imagp;
	delete[] polar;
#endif

}

/* Calculate the power spectrum */
void fft::powerSpectrum(int start, float *data, float *window, float *magnitude,float *phase) {
	int i;

	//windowing
	for (i = 0; i < n; i++) {
		in_real[i] = data[start + i] * window[i];
	}

	RealFFT(n, in_real, out_real, out_img);

	for (i = 0; i < half; i++) {
		/* compute power */
		float power = out_real[i]*out_real[i] + out_img[i]*out_img[i];
		/* compute magnitude and phase */
		magnitude[i] = sqrt(power);
		phase[i] = atan2(out_img[i],out_real[i]);

        //		if (magnitude[i] < 0.000001){ // less than 0.1 nV
        //			magnitude[i] = 0; // out of range
        //		} else {
        //			magnitude[i] = 20.0*log10(magnitude[i] + 1);  // get to to db scale
        //		}
	}

}

void fft::convToDB(float *in, float *out) {
	for (int i = 0; i < half; i++) {
		if (in[i] < 0.000001){ // less than 0.1 nV
			out[i] = 0; // out of range
		} else {
			out[i] = 20.0*log10(in[i] + 1);  // get to to db scale
		}
	}
}

#ifdef __APPLE_CC__

/* Calculate the power spectrum */
void fft::powerSpectrum_vdsp(int start, float *data, float *window, float *magnitude,float *phase) {

    uint32_t        i;

	//multiply by window
	vDSP_vmul(data, 1, window, 1, in_real, 1, n);

	//convert to split complex format - evens and odds
	vDSP_ctoz((COMPLEX *) in_real, 2, &A, 1, half);

	//calc fft
	vDSP_fft_zrip(setupReal, &A, 1, log2n, FFT_FORWARD);

	//scale by 2 (see vDSP docs)
	static float scale=0.5 ;
	vDSP_vsmul(A.realp, 1, &scale, A.realp, 1, half);
	vDSP_vsmul(A.imagp, 1, &scale, A.imagp, 1, half);

	//back to split complex format
	vDSP_ztoc(&A, 1, (COMPLEX*) out_real, 2, half);

	//convert to polar
	vDSP_polar(out_real, 2, polar, 2, half);

	for (i = 0; i < half; i++) {
		magnitude[i]=polar[2*i];
		phase[i]=polar[2*i + 1];
	}

}

void fft::convToDB_vdsp(float *in, float *out) {
	float ref = 1.0;
	vDSP_vdbcon(in, 1, &ref, out, 1, half, 1);
	//get rid of any -infs
	float vmin=0.0;
	float vmax=9999999.0;
	vDSP_vclip(out, 1, &vmin, &vmax, out, 1, half);
}

#endif

void fft::inversePowerSpectrum(int start, float *finalOut, float *window, float *magnitude,float *phase) {
	int i;

	/* get real and imag part */
	for (i = 0; i < half; i++) {
        //		float mag = pow(10.0, magnitude[i] / 20.0) - 1.0;
        //		in_real[i] = mag *cos(phase[i]);
        //		in_img[i]  = mag *sin(phase[i]);
		in_real[i] = magnitude[i] *cos(phase[i]);
		in_img[i]  = magnitude[i] *sin(phase[i]);
	}

	/* zero negative frequencies */
	memset(in_real+half, half, 0.0);
	memset(in_img+half, half, 0.0);

	FFT(n, 1, in_real, in_img, out_real, out_img); // second parameter indicates inverse transform

	for (i = 0; i < n; i++) {
		finalOut[start + i] += out_real[i] * window[i]
		;
	}

}

#ifdef __APPLE_CC__
void fft::inversePowerSpectrum_vdsp(int start, float *finalOut, float *window, float *magnitude,float *phase) {
	uint32_t i;

	for (i = 0; i < half; i++) {
        //		polar[2*i] = pow(10.0, magnitude[i] / 20.0) - 1.0;
		polar[2*i] = magnitude[i];
		polar[2*i + 1] = phase[i];
	}

	vDSP_rect(polar, 2, in_real, 2, half);

	vDSP_ctoz((COMPLEX*) in_real, 2, &A, 1, half);
	vDSP_fft_zrip(setupReal, &A, 1, log2n, FFT_INVERSE);
	vDSP_ztoc(&A, 1, (COMPLEX*) out_real, 2, half);

    float scale = 1./n;
	vDSP_vsmul(out_real, 1, &scale, out_real, 1, n);

	//multiply by window
	vDSP_vmul(out_real, 1, window, 1, finalOut, 1, n);

}

#endif


text:
/*
 *  maximilian
 *  platform independent synthesis library using portaudio or rtaudio
 *
 *  Created by Mick Grierson on 29/12/2009.
 *  Copyright 2009 Mick Grierson & Strangeloop Limited. All rights reserved.
 *	Thanks to the Goldsmiths Creative Computing Team.
 *	Special thanks to Arturo Castro for the PortAudio implementation.
 *
 *	Permission is hereby granted, free of charge, to any person
 *	obtaining a copy of this software and associated documentation
 *	files (the "Software"), to deal in the Software without
 *	restriction, including without limitation the rights to use,
 *	copy, modify, merge, publish, distribute, sublicense, and/or sell
 *	copies of the Software, and to permit persons to whom the
 *	Software is furnished to do so, subject to the following
 *	conditions:
 *
 *	The above copyright notice and this permission notice shall be
 *	included in all copies or substantial portions of the Software.
 *
 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 *	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 *	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 *	NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 *	HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 *	WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 *	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 *	OTHER DEALINGS IN THE SOFTWARE.
 *
 */

#include "maxiFFT.h"

#if 0 //#ifdef __APPLE_CC__
	#include "maximilian.h"
#else
	#include "../maximilian.h"
#endif

#include <iostream>
#include "math.h"

using namespace std;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//F F T
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void maxiFFT::setup(int _fftSize, int _windowSize, int _hopSize) {
	_fft = new fft(_fftSize);
	fftSize = _fftSize;
	windowSize = _windowSize;
	bins = fftSize / 2;
	hopSize = _hopSize;
	buffer = (float *) malloc(fftSize * sizeof(float));
	magnitudes = (float *) malloc(bins * sizeof(float));
	magnitudesDB = (float *) malloc(bins * sizeof(float));
	phases = (float *) malloc(bins * sizeof(float));
	avgPower = new float;
	memset(buffer, 0, fftSize * sizeof(float));
	memset(magnitudes, 0, bins * sizeof(float));
	memset(magnitudesDB, 0, bins * sizeof(float));
	memset(phases, 0, bins * sizeof(float));
	*avgPower = 0;
	pos =windowSize - hopSize;
	newFFT = 0;
	window = (float*) malloc(fftSize * sizeof(float));
	memset(window, 0, fftSize * sizeof(float));
	fft::genWindow(3, windowSize, window);
}

bool maxiFFT::process(float value) {
	//add value to buffer at current pos
	buffer[pos++] = value;
	//if buffer full, run fft
	newFFT = pos == windowSize;
	if (newFFT) {
#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
		_fft->powerSpectrum_vdsp(0, buffer, window, magnitudes, phases);
#else
		_fft->powerSpectrum(0, buffer, window, magnitudes, phases);
#endif
		//shift buffer back by one hop size
		memcpy(buffer, buffer + hopSize, (windowSize - hopSize) * sizeof(float));
		//set the new data to zero (is this necessary?, it will get overwritten)
		memset(buffer + windowSize - hopSize, 0, hopSize * sizeof(float));
		//reset pos to start of hop
		pos= windowSize - hopSize;
	}
	return newFFT;
}

float* maxiFFT::magsToDB() {
#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
	_fft->convToDB_vdsp(magnitudes, magnitudesDB);
#else
	_fft->convToDB(magnitudes, magnitudesDB);
#endif
	return magnitudesDB;
}

float maxiFFT::spectralFlatness() {
	float geometricMean=0, arithmaticMean=0;
	for(int i=0; i < bins; i++) {
		if (magnitudes[i] != 0)
			geometricMean += logf(magnitudes[i]);
		arithmaticMean += magnitudes[i];
	}
	geometricMean = expf(geometricMean / (float)bins);
	arithmaticMean /= (float)bins;
	return arithmaticMean !=0 ?  geometricMean / arithmaticMean : 0;
}

float maxiFFT::spectralCentroid() {
	float x=0, y=0;
	for(int i=0; i < bins; i++) {
		x += fabs(magnitudes[i]) * i;
		y += fabs(magnitudes[i]);
	}
	return y != 0 ? x / y * ((float) maxiSettings::sampleRate / fftSize) : 0;
}

maxiFFT::~maxiFFT() {
	delete _fft;
	if (buffer)
		delete[] buffer,magnitudes,phases,window, avgPower, magnitudesDB;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//I N V E R S E  F F T
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void maxiIFFT::setup(int _fftSize, int _windowSize, int _hopSize) {
	_fft = new fft(_fftSize);
	fftSize = _fftSize;
	windowSize = _windowSize;
	bins = fftSize / 2;
	hopSize = _hopSize;
	buffer = (float *) malloc(fftSize * sizeof(float));
	ifftOut = (float *) malloc(fftSize * sizeof(float));
	memset(buffer, 0, windowSize * sizeof(float));
	pos =0;
	window = (float*) malloc(fftSize * sizeof(float));
	memset(window, 0, fftSize * sizeof(float));
	fft::genWindow(3, windowSize, window);
}

float maxiIFFT::process(float *magnitudes, float *phases) {
	if (0==pos) {
		//do ifft
		memset(ifftOut, 0, fftSize * sizeof(float));
#if defined(__APPLE_CC__) && !defined(_NO_VDSP)
		_fft->inversePowerSpectrum_vdsp(0, ifftOut, window, magnitudes, phases);
#else
		_fft->inversePowerSpectrum(0, ifftOut, window, magnitudes, phases);
#endif
		//add to output
		//shift back by one hop
		memcpy(buffer, buffer+hopSize, (fftSize - hopSize) * sizeof(float));
		//clear the end chunk
		memset(buffer + (fftSize - hopSize), 0, hopSize * sizeof(float));
		//merge new output
		for(int i=0; i < fftSize; i++) {
			buffer[i] += ifftOut[i];
		}
	}

	nextValue = buffer[pos];
	//limit the values, this alg seems to spike occasionally (and break the audio drivers)
	if (nextValue > 0.99999f) nextValue = 0.99999f;
	if (nextValue < -0.99999f) nextValue = -0.99999f;
	if (hopSize == ++pos ) {
		pos=0;
	}

	return nextValue;
}

maxiIFFT::~maxiIFFT() {
	delete _fft;
	delete[] ifftOut, buffer, window;
};

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//O C T A V E  A N A L Y S E R
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void maxiFFTOctaveAnalyzer::setup(float samplingRate, int nBandsInTheFFT, int nAveragesPerOctave){

    samplingRate = samplingRate;
    nSpectrum = nBandsInTheFFT;
    spectrumFrequencySpan = (samplingRate / 2.0f) / (float)(nSpectrum);
    nAverages = nBandsInTheFFT;
    // fe:  2f for octave bands, sqrt(2) for half-octave bands, cuberoot(2) for third-octave bands, etc
    if (nAveragesPerOctave==0) // um, wtf?
		nAveragesPerOctave = 1;
//    nAveragesPerOctave = nAveragesPerOctave;
    averageFrequencyIncrement = pow(2.0f, 1.0f/(float)(nAveragesPerOctave));
    // this isn't currently configurable (used once here then no effect), but here's some reasoning:
    // 43 is a good value if you want to approximate "computer" octaves: 44100/2/2/2/2/2/2/2/2/2/2
    // 55 is a good value if you'd rather approximate A-440 octaves: 440/2/2/2
    // 65 is a good value if you'd rather approximate "middle-C" octaves:  ~262/2/2
    // you could easily double it if you felt the lowest band was just rumble noise (as it probably is)
    // but don't go much smaller unless you have a huge fft window size (see below for more why)
    // keep in mind, if you change it, that the number of actual bands may change +/-1, and
    // for some values, the last averaging band may not be very useful (may extend above nyquist)
    firstOctaveFrequency = 55.0f;
    // for each spectrum[] bin, calculate the mapping into the appropriate average[] bin.
    // this gives us roughly log-sized averaging bins, subject to how "fine" the spectrum bins are.
    // with more spectrum bins, you can better map into the averaging bins (especially at low
    // frequencies) or use more averaging bins per octave.  with an fft window size of 2048,
    // sampling rate of 44100, and first octave around 55, that's about enough to do half-octave
    // analysis.  if you don't have enough spectrum bins to map adequately into averaging bins
    // at the requested number per octave then you'll end up with "empty" averaging bins, where
    // there is no spectrum available to map into it.  (so... if you have "nonreactive" averages,
    // either increase fft buffer size, or decrease number of averages per octave, etc)
    spe2avg = new int[nSpectrum];
    int avgidx = 0;
    float averageFreq = firstOctaveFrequency; // the "top" of the first averaging bin
    // we're looking for the "top" of the first spectrum bin, and i'm just sort of
    // guessing that this is where it is (or possibly spectrumFrequencySpan/2?)
    // ... either way it's probably close enough for these purposes
    float spectrumFreq = spectrumFrequencySpan;
    for (int speidx=0; speidx < nSpectrum; speidx++) {
		while (spectrumFreq > averageFreq) {
			avgidx++;
			averageFreq *= averageFrequencyIncrement;
		}
		spe2avg[speidx] = avgidx;
		spectrumFreq += spectrumFrequencySpan;
    }
    nAverages = avgidx;
    averages = new float[nAverages];
    peaks = new float[nAverages];
    peakHoldTimes = new int[nAverages];
    peakHoldTime = 0; // arbitrary
    peakDecayRate = 0.9f; // arbitrary
    linearEQIntercept = 1.0f; // unity -- no eq by default
    linearEQSlope = 0.0f; // unity -- no eq by default
}

void maxiFFTOctaveAnalyzer::calculate(float * fftData){

	int last_avgidx = 0; // tracks when we've crossed into a new averaging bin, so store current average
    float sum = 0.0f; // running total of spectrum data
    int count = 0; // count of spectrums accumulated (for averaging)
    for (int speidx=0; speidx < nSpectrum; speidx++) {
		count++;
		sum += fftData[speidx] * (linearEQIntercept + (float)(speidx) * linearEQSlope);
		int avgidx = spe2avg[speidx];
		if (avgidx != last_avgidx) {

			for (int j = last_avgidx; j < avgidx; j++){
				averages[j] = sum / (float)(count);
			}
			count = 0;
			sum = 0.0f;
		}
		last_avgidx = avgidx;
    }
    // the last average was probably not calculated...
    if ((count > 0) && (last_avgidx < nAverages)){
		averages[last_avgidx] = sum / (float)(count);
	}

    // update the peaks separately
    for (int i=0; i < nAverages; i++) {
		if (averages[i] >= peaks[i]) {
			// save new peak level, also reset the hold timer
			peaks[i] = averages[i];
			peakHoldTimes[i] = peakHoldTime;
		} else {
			// current average does not exceed peak, so hold or decay the peak
			if (peakHoldTimes[i] > 0) {
				peakHoldTimes[i]--;
			} else {
				peaks[i] *= peakDecayRate;
			}
		}
    }
}


text:
/*
 *  maximilian
 *  platform independent synthesis library using portaudio or rtaudio
 *
 *  Created by Mick Grierson on 29/12/2009.
 *  Copyright 2009 Mick Grierson & Strangeloop Limited. All rights reserved.
 *	Thanks to the Goldsmiths Creative Computing Team.
 *	Special thanks to Arturo Castro for the PortAudio implementation.
 *
 *	Permission is hereby granted, free of charge, to any person
 *	obtaining a copy of this software and associated documentation
 *	files (the "Software"), to deal in the Software without
 *	restriction, including without limitation the rights to use,
 *	copy, modify, merge, publish, distribute, sublicense, and/or sell
 *	copies of the Software, and to permit persons to whom the
 *	Software is furnished to do so, subject to the following
 *	conditions:
 *
 *	The above copyright notice and this permission notice shall be
 *	included in all copies or substantial portions of the Software.
 *
 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 *	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 *	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 *	NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 *	HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 *	WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 *	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 *	OTHER DEALINGS IN THE SOFTWARE.
 *
 */

#ifndef _MAXI_FFT
#define _MAXI_FFT

//#define _NO_VDSP  //set this if you don't want to use apple's vDSP fft functions

#include "fft.h"
#include "stddef.h"

class maxiFFT {

public:
	maxiFFT(){
		_fft = NULL;
		buffer = magnitudes = phases  = window = avgPower = NULL;
	};
	~maxiFFT();
	void setup(int fftSize, int windowSize, int hopSize);
	bool process(float value);
	float* magsToDB();
	float *magnitudes, *phases, *magnitudesDB;
	float *avgPower;
	int windowSize;
	int hopSize;
	int bins;

	//features
	float spectralFlatness();
	float spectralCentroid();

private:
	float *buffer, *window;
	int pos;
	float nextValue;
	int fftSize;
	fft *_fft;
	bool newFFT;

};

class maxiIFFT {

public:
	maxiIFFT(){
		_fft=0;
	};
	~maxiIFFT();
	void setup(int fftSize, int windowSize, int hopSize);
	float process(float *magnitudes, float *phases);

private:
	float *ifftOut, *buffer, *window;
	int windowSize;
	int bins;
	int hopSize;
	int pos;
	float nextValue;
	int fftSize;
	fft *_fft;
};

class maxiFFTOctaveAnalyzer {
    /*based on code by David Bollinger, http://www.davebollinger.com/
     */
public:

	float samplingRate; // sampling rate in Hz (needed to calculate frequency spans)
	int nSpectrum; // number of spectrum bins in the fft
	int nAverages; // number of averaging bins here
	int nAveragesPerOctave; // number of averages per octave as requested by user
	float spectrumFrequencySpan; // the "width" of an fft spectrum bin in Hz
	float firstOctaveFrequency; // the "top" of the first averaging bin here in Hz
	float averageFrequencyIncrement; // the root-of-two multiplier between averaging bin frequencies
	float * averages; // the actual averages
	float * peaks; // peaks of the averages, aka "maxAverages" in other implementations
	int * peakHoldTimes; // how long to hold THIS peak meter?  decay if == 0
	int peakHoldTime; // how long do we hold peaks? (in fft frames)
	float peakDecayRate; // how quickly the peaks decay:  0f=instantly .. 1f=not at all
	int * spe2avg; // the mapping between spectrum[] indices and averages[] indices
	// the fft's log equalizer() is no longer of any use (it would be nonsense to log scale
	// the spectrum values into log-sized average bins) so here's a quick-and-dirty linear
	// equalizer instead:
	float linearEQSlope; // the rate of linear eq
	float linearEQIntercept; // the base linear scaling used at the first averaging bin
	// the formula is:  spectrum[i] * (linearEQIntercept + i * linearEQSlope)
	// so.. note that clever use of it can also provide a "gain" control of sorts
	// (fe: set intercept to 2f and slope to 0f to double gain)

	void setup(float samplingRate, int nBandsInTheFFT, int nAveragesPerOctave);

	void calculate(float * fftData);

};

#endif


text:
/*
 *  maximilian
 *  platform independent synthesis library using portaudio or rtaudio
 *
 *  Created by Mick Grierson on 29/12/2009.
 *  Copyright 2009 Mick Grierson & Strangeloop Limited. All rights reserved.
 *	Thanks to the Goldsmiths Creative Computing Team.
 *	Special thanks to Arturo Castro for the PortAudio implementation.
 *
 *	Permission is hereby granted, free of charge, to any person
 *	obtaining a copy of this software and associated documentation
 *	files (the "Software"), to deal in the Software without
 *	restriction, including without limitation the rights to use,
 *	copy, modify, merge, publish, distribute, sublicense, and/or sell
 *	copies of the Software, and to permit persons to whom the
 *	Software is furnished to do so, subject to the following
 *	conditions:
 *
 *	The above copyright notice and this permission notice shall be
 *	included in all copies or substantial portions of the Software.
 *
 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 *	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 *	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 *	NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 *	HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 *	WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 *	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 *	OTHER DEALINGS IN THE SOFTWARE.
 *
 */

#ifndef _FFT
#define _FFT

#ifndef M_PI
#define	M_PI		3.14159265358979323846  /* pi */
#endif

#ifdef __APPLE_CC__
#include <Accelerate/Accelerate.h>
#endif

class fft {

public:

	fft(int fftSize);
	~fft();

	int n; //fftSize
	int half; //halfFFTSize

	float			*in_real, *out_real, *in_img, *out_img;

#ifdef __APPLE_CC__
	int log2n; //log2(n);
    FFTSetup        setupReal;
    COMPLEX_SPLIT   A;
	float *polar;
	void powerSpectrum_vdsp(int start, float *data, float *window, float *magnitude,float *phase);
	void inversePowerSpectrum_vdsp(int start, float *finalOut, float *window, float *magnitude,float *phase);
	void convToDB_vdsp(float *in, float *out);
#endif

	/* Calculate the power spectrum */
	void powerSpectrum(int start, float *data, float *window, float *magnitude, float *phase);
	/* ... the inverse */
	void inversePowerSpectrum(int start, float *finalOut, float *window, float *magnitude,float *phase);
	void convToDB(float *in, float *out);

	static void genWindow(int whichFunction, int NumSamples, float *window);

};

#endif


text:
[yapf]
based_on_style = pep8
column_limit = 69
coalesce_brackets = false
dedent_closing_brackets = false
each_dict_entry_on_separate_line = false
indent_dictionary_value = true
space_between_ending_comma_and_closing_bracket = false
split_all_comma_separated_values = false
split_all_top_level_comma_separated_values = false
split_before_first_argument = false
split_before_closing_bracket = false
split_before_expression_after_opening_paren = false
split_before_named_assigns = false
split_complex_comprehension = false
split_penalty_after_opening_bracket = 3000
split_penalty_for_added_line_split = 300

repo_name: lesteve/scikit-learn-tutorial

repo_language: Jupyter Notebook


text:
name: scikit-learn-tutorial

dependencies:
  - python
  - scikit-learn
  - pandas
  - conda-forge::pandas-profiling
  - matplotlib
  - seaborn
  - jupyter


text:
PYTHON_SCRIPTS_DIR = python_scripts
NOTEBOOKS_DIR = notebooks
RENDERED_NOTEBOOKS_DIR = rendered_notebooks
JUPYTER_KERNEL := python3
MINIMAL_NOTEBOOK_FILES = $(shell ls $(PYTHON_SCRIPTS_DIR)/*.py | perl -pe "s@$(PYTHON_SCRIPTS_DIR)@$(NOTEBOOKS_DIR)@" | perl -pe "s@\.py@.ipynb@")
MINIMAL_RENDERED_NOTEBOOK_FILES = $(shell ls $(PYTHON_SCRIPTS_DIR)/*.py | perl -pe "s@$(PYTHON_SCRIPTS_DIR)@$(RENDERED_NOTEBOOKS_DIR)@" | perl -pe "s@\.py@.ipynb@")

all: $(RENDERED_NOTEBOOKS_DIR)

.PHONY: $(NOTEBOOKS_DIR) $(RENDERED_NOTEBOOKS_DIR) sanity_check_$(PYTHON_SCRIPTS_DIR) sanity_check_$(NOTEBOOKS_DIR) sanity_check_$(RENDERED_NOTEBOOKS_DIR) all

$(NOTEBOOKS_DIR): sanity_check_$(PYTHON_SCRIPTS_DIR) $(MINIMAL_NOTEBOOK_FILES) sanity_check_$(NOTEBOOKS_DIR)

$(RENDERED_NOTEBOOKS_DIR): $(MINIMAL_RENDERED_NOTEBOOK_FILES) sanity_check_$(RENDERED_NOTEBOOKS_DIR)

$(NOTEBOOKS_DIR)/%.ipynb: $(PYTHON_SCRIPTS_DIR)/%.py
	jupytext --set-formats $(PYTHON_SCRIPTS_DIR)//py:percent,$(NOTEBOOKS_DIR)//ipynb $<
	jupytext --sync $<

$(RENDERED_NOTEBOOKS_DIR)/%.ipynb: $(NOTEBOOKS_DIR)/%.ipynb
	cp $< $@
	jupyter nbconvert --to notebook --execute --ExecutePreprocessor.timeout=None \
                    --ExecutePreprocessor.kernel_name=$(JUPYTER_KERNEL) --inplace $@

sanity_check_$(PYTHON_SCRIPTS_DIR):
	python build_scripts/check-python-scripts.py $(PYTHON_SCRIPTS_DIR)
	yapf --recursive --in-place --parallel $(PYTHON_SCRIPTS_DIR)

sanity_check_$(NOTEBOOKS_DIR):
	python build_scripts/sanity-check.py $(PYTHON_SCRIPTS_DIR) $(NOTEBOOKS_DIR)

sanity_check_$(RENDERED_NOTEBOOKS_DIR):
	python build_scripts/sanity-check.py $(NOTEBOOKS_DIR) $(RENDERED_NOTEBOOKS_DIR)


text:
scikit-learn
pandas
pandas-profiling
matplotlib
seaborn
plotly
jupyter


text:
Note: to edit on HackMD https://hackmd.io/6fFy3_y6SOWonOYeGaNiHw?edit

# Intro module: What is Machine Learning
 (only slides / videos i.e. no code)
## Learning objectives:

- familiarity with scope of most common ML applications in business or science
- understand difference between memorization and generalization
- familiarity with main machine learning concepts and vocabulary

## Content
Why and when? Example applications.
- Mention the iris example, pitch it as historical, but also as a
  botanical and agriculture problem.
  The benefit of this example is that it forces to think about
  measurement, but also because it has one class that is easy to separate
  only one feature
    - Aurelie: real irises for the video?
- The "adult" dataset
    - Maybe looking at it with excel, to be in an environment familiar to
      people
    - Mention the importance of data visualization: intuitions about the
      data can be very helpful

Learning from data vs expertly engineered decision rules
- One the iris example, show that cutting on one specific feature
  separates well one class
- How do we automate this? How do we achieve this on more complex data
  such as the census dataset?

Descriptive vs predictive analysis
- Generalization (Out of sample properties)
- An example of where it makes a difference: if the data has redundant
  variables, such as expressing the education level as the name of the
  degree or the corresponding number of years of education

Generalization vs memorization: the need for a train / test split
- The nearest neighbors example to illustrate this

Features and samples
- The data matrix
  - Build the data matrix of Iris

Supervised vs Unsupervised
- Formalize supervised learning (define "X" and "y")
- Introduce unsupervised learning, for instance dimensionality reduction
  (and go back to the example of redundant variables: if we have many of
  these, we should be able to reduce the problem without even looking at
  y

Regression vs Classification
- In the adult data: it would make more sense to do a continuous
  prediction
- In the iris example, it is naturally a classification problem

A few words about the style and scope of this course: it is centered
around code, though we strive to keep it simple

## Quizz:
Given a case study (e.g. pricing apartments based on a real estate website database) and sample toy dataset: say whether it’s an application of supervised vs unsupervised, classification vs regression, what are the features, what is the target variable, what is a record.

Propose a hand engineer decision rule that can be used as a baseline

Propose a quantitative evaluation of the success of this decision rule.

# The Predictive Modeling Pipeline

## Notebook module #1: exploratory analysis
### Learning objectives:

- load tabular data with pandas
- visualize marginal distribution with histograms
- visualize pairwise interactions with scatter plots
- identify outlier and dynamic range of each column

### Content
Defining a predictive task that relates to the business or scientific case

Pandas read_csv

Simple exploratory data analysis with pandas and matplotlib

## Notebook module #2: basic preprocessing for minimal model fit

### Learning objectives:
- Know the difference between a numerical and a categorical variable
- use a scaler
- convert category labels to dummy variables
- combine feature preprocessing and model with pipeline
- evaluate generalization of model with cross-validation

### Content

Prepare a train  / test split

Basic model on numerical features only

Basic processing: missing values and scaling

Use a pipeline to evaluate model with cross-validation with and without scaling

Handling categorical variables with one-hot encoding

Use the column transformer to build pipeline with heterogeneous dtype

Model fitting and performance evaluation with cross-validation

- Gael thinks that we could use a video here for cross-validation
  (in particular, the "plot_cv_indices" in the notebook gets a bit in the
  way of being accessible and didactic

## Notebook module #3: basic parameter tuning and final test score evaluation

### Learning objectives:
- Learn to no trust blindly the default parameters of scikit-learn estimators

### Content
Parameter tuning with Grid and Random hyperparameter search
Nested cross-validation

Confirmation of performance with final test set

# Supervised learning
## Learning objectives:
Understand decision rules for a few important algorithms
Know how to diagnose model generalization errors (overfitting especially)
How to use variable selection and generalization to fight overfitting
Feature engineering to limit underfitting

## Olivier: Overfitting/Underfitting validation curves, learning curves, regularisation with linear models

- Video about overfitting?

## Loïc: Trees in depth + ensembles

## Guillaume: Evaluation of supervised learning models:

Confusion matrix for classifiers / precision / recall / ROC AUC curve (Mention imbalanced classes)
Predict vs True plot for regressors

## Olivier: Linear models in depth
Logistic Regression, linear regression, classification vs regression, multi-class, linear separability. Pros and cons
        L1 and L2 penalty for linear models
        Learning curves and validation curves (video: how to read curves)

## Baselines: majority class classifier (already in second module) and k-nearest neighbors

## Feature engineering to augment the expressivity of linear models:

Binning / Polynomial feature extraction / Nystroem method

Feature selection to combat overfitting and speed-up models

## Univariate feature selection
Show catastrophic example where feature selection is done on the whole dataset rather than only on train

## Evaluating the feature importance with permutations
Failure Mode : cardinality bias of overfitting random forest feature importances

## Looking at the decision function with partial dependence plots

Gael thinks that explaining the difference between conditional and
marginal interpretation is important.

Stability of hyperparameter during cross-validation


text: # This repo has moved to: https://github.com/INRIA/scikit-learn-mooc/


text:
name: CI

on:
  push:
    branches:
    - master

jobs:
  build:

    runs-on: ubuntu-latest
    steps:
    - uses: actions/checkout@v1
    - name: Setup Node
      uses: actions/setup-node@v1
      with:
        node-version: '12.x'
    - uses: actions/setup-python@v1
      with:
        python-version: '3.x'
        architecture: 'x64'

    - name: Install python dependencies
      working-directory: ./slides
      run: pip install -r requirements.txt

    - name: Install npm dependencies
      run: npm install -g bellbind/remarkjs-pdf

    - name: Export html
      working-directory: ./slides
      run: make && mkdir ../dist && mv *.html ../dist/
    - name: Export pdf
      working-directory: ./dist
      run: for f in *.html ; do remarkjs-pdf "$f"; done

    - name: Upload distributed files
      uses: actions/upload-artifact@v1
      with:
        name: dist
        path: dist/


text:
"""
Simple styling used for matplotlib figures
"""

from matplotlib import pyplot as plt

# Configuration settings to help visibility on small screen / prints
plt.rcParams['xtick.labelsize'] = 20
plt.rcParams['ytick.labelsize'] = 20
plt.rcParams['figure.titlesize'] = 15
plt.rcParams['font.size'] = 20
plt.rcParams['axes.labelsize'] = 20
plt.rcParams['axes.facecolor'] = 'none'
plt.rcParams['legend.fontsize'] = 18
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['figure.figsize'] = [.8 * 6.4, .8 * 4.8]
plt.rcParams['legend.frameon'] = False
plt.rcParams['legend.columnspacing'] = 1.8
plt.rcParams['legend.handlelength'] = 1.5
plt.rcParams['legend.handletextpad'] = 0.5

# Utility functions
def light_axis():
    "Hide the top and right spines"
    ax = plt.gca()
    for s in ('top', 'right'):
        ax.spines[s].set_visible(False)
    plt.xticks(())
    plt.yticks(())
    plt.subplots_adjust(left=.01, bottom=.01, top=.99, right=.99)

def no_axis():
    plt.axis('off')
    plt.subplots_adjust(left=.0, bottom=.0, top=1, right=1)


text:
"""
Simple example of overfit with splines
"""
import numpy as np
from matplotlib import pyplot as plt
import style_figs

from sklearn import datasets, linear_model

# Load the diabetes dataset
diabetes = datasets.load_diabetes()

# Use only one feature
diabetes_X = diabetes.data[:, np.newaxis]
diabetes_X_temp = diabetes_X[:, :, 2]

# Split the data into training/testing sets
diabetes_X_train = diabetes_X_temp[:-200:3]
diabetes_X_test = diabetes_X_temp[-200:].T

# Split the targets into training/testing sets
diabetes_y_train = diabetes.target[:-200:3]
diabetes_y_test = diabetes.target[-200:]

# Sort the data and remove duplicates (for interpolation)
order = np.argsort(diabetes_X_train.ravel())
X_train = diabetes_X_train.ravel()[order]
y_train = diabetes_y_train[order]
# Avoid duplicates
y_train_ = list()
for this_x in np.unique(X_train):
    y_train_.append(np.mean(y_train[X_train == this_x]))
X_train = np.unique(X_train)

y_train = np.array(y_train_)

# Create linear regression object
regr = linear_model.LinearRegression()

# Train the model using the training sets
regr.fit(X_train.reshape((-1, 1)), y_train)

plt.figure(1, figsize=(.8*4, .8*3), facecolor='none')
# Plot with test data
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

plt.scatter(X_train, y_train,  color='k', s=9)

plt.plot([-.08, .12], regr.predict([[-.08, ], [.12, ]]),
        linewidth=3)

plt.axis('tight')
ymin, ymax = plt.ylim()
style_figs.light_axis()
plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('ols_simple.svg', facecolor='none', edgecolor='none')

plt.scatter(diabetes_X_test, diabetes_y_test,  color='C1', s=9)
plt.ylim(ymin, ymax)
plt.xlim(-.08, .12)

plt.savefig('ols_test.svg', facecolor='none', edgecolor='none')

# Plot cubic splines
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

from scipy import interpolate
f = interpolate.interp1d(X_train, y_train,
                         kind="quadratic",
                         bounds_error=False, fill_value="extrapolate")
plt.scatter(X_train, y_train,  color='k', s=9, zorder=20)
x_spline = np.linspace(-.08, .12, 600)
y_spline = f(x_spline)
plt.plot(x_spline, y_spline, linewidth=3)

plt.axis('tight')
plt.xlim(-.08, .12)
plt.ylim(ymin, ymax)

style_figs.light_axis()

plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('splines_cubic.svg', facecolor='none', edgecolor='none')

plt.scatter(diabetes_X_test, diabetes_y_test,  color='C1', s=9)
plt.savefig('splines_test.svg', facecolor='none', edgecolor='none')

plt.show()


text:
"""
Some simple visualizations on the iris data.
"""

import numpy as np
from sklearn import datasets
from matplotlib import pyplot as plt
import style_figs

iris = datasets.load_iris()

# Plot the histograms of each class for each feature

X = iris.data
y = iris.target
for x, feature_name in zip(X.T, iris.feature_names):
    plt.figure(figsize=(2.5, 2))
    patches = list()
    for this_y, target_name in enumerate(iris.target_names):
        patch = plt.hist(x[y == this_y],
                         bins=np.linspace(x.min(), x.max(), 20),
                         label=target_name)
        patches.append(patch[-1][0])
    style_figs.light_axis()
    feature_name = feature_name.replace(' ', '_')
    feature_name = feature_name.replace('(', '')
    feature_name = feature_name.replace(')', '')
    plt.savefig('iris_{}_hist.svg'.format(feature_name))

plt.figure(figsize=(6, .25))
plt.legend(patches, iris.target_names, ncol=3, loc=(0, -.37),
           borderaxespad=0)
style_figs.no_axis()
plt.savefig('legend_irises.svg')


text:
import numpy as np
from matplotlib import pyplot as plt

# Set up figures look and feel
import style_figs

# %%
# Our data-generating process

def f(t):
    return 1.2 * t ** 2 + .1 * t ** 3 - .4 * t ** 5 - .5 * t ** 9

N_SAMPLES = 50

rng = np.random.RandomState(0)
x = 2 * rng.rand(N_SAMPLES) - 1

y = f(x) + .4 * rng.normal(size=N_SAMPLES)

x_test = 2 * rng.rand(N_SAMPLES) - 1

y_test = f(x_test) + .4 * rng.normal(size=N_SAMPLES)

plt.figure()
plt.scatter(x, y, s=20, color='k')

style_figs.no_axis()
plt.subplots_adjust(top=.96)
plt.xlim(-1.1, 1.1)
plt.ylim(-.74, 2.1)
plt.savefig('polynomial_overfit_0.svg', facecolor='none', edgecolor='none')

# %%
# Our model (polynomial regression)

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# %%
# Fit model with various complexity in the polynomial degree

plt.figure()
plt.scatter(x, y, s=20, color='k')

t = np.linspace(-1, 1, 100)

for d in (1, 2, 5, 9):
    model = make_pipeline(PolynomialFeatures(degree=d), LinearRegression())
    model.fit(x.reshape(-1, 1), y)
    plt.plot(t, model.predict(t.reshape(-1, 1)), label='Degree %d' % d,
             linewidth=4)

    style_figs.no_axis()
    plt.legend(loc='upper center', borderaxespad=0, borderpad=0)
    plt.subplots_adjust(top=.96)
    plt.ylim(-.74, 2.1)

    plt.savefig('polynomial_overfit_%d.svg' % d, facecolor='none',
                edgecolor='none')

plt.plot(t, f(t), 'k--', label='Truth')

style_figs.no_axis()
plt.legend(loc='upper center', borderaxespad=0, borderpad=0)
plt.ylim(-.74, 2.1)
plt.subplots_adjust(top=.96)

plt.savefig('polynomial_overfit.svg', facecolor='none', edgecolor='none')

# %%
# A figure with the true model and the estimated one

plt.figure(figsize=[.5 * 6.4, .5 * 4.9])
plt.scatter(x, y, s=20, color='k')
plt.plot(t, model.predict(t.reshape(-1, 1)), color='C3',
         label='$\hat{f}$')

plt.plot(t, f(t), 'k--', label='$f^{\star}$')
style_figs.no_axis()
plt.ylim(-1.25, 2.5)
plt.legend(loc='upper center', borderaxespad=0, borderpad=0,
           labelspacing=.4, fontsize=26)
plt.subplots_adjust(top=1)

plt.savefig('polynomial_overfit_simple.svg', facecolor='none', edgecolor='none')

# %%
# A figure with the true model and the estimated one

plt.figure(figsize=[.5 * 6.4, .5 * 4.9])
plt.scatter(x, y, s=20, color='k')
plt.plot(t, model.predict(t.reshape(-1, 1)), color='C3',
         label='Fitted model')

plt.plot(t, f(t), 'k--', label='Best possible fit\n$\\approx$generative process')
style_figs.no_axis()
plt.ylim(-1.25, 2.5)
plt.legend(loc='upper right', borderaxespad=0, borderpad=0,
           labelspacing=.6, fontsize=16)
plt.subplots_adjust(top=1)

# More detailed legend
plt.savefig('polynomial_overfit_simple_legend.svg', facecolor='none', edgecolor='none')

# %%
# Underfit settings

model = make_pipeline(PolynomialFeatures(degree=1), LinearRegression())
model.fit(x.reshape(-1, 1), y)

plt.figure(figsize=[.5 * 6.4, .5 * 4.9])
plt.scatter(x, y, s=20, color='k')
plt.plot(t, model.predict(t.reshape(-1, 1)), color='C0',
         label='Fitted model\n$\\approx$best possible fit')

plt.plot(t, f(t), 'k--', label='Generative process')
style_figs.no_axis()
plt.ylim(-1.25, 2.5)
plt.legend(loc='upper right', borderaxespad=0, borderpad=0,
           labelspacing=.4, fontsize=16)
plt.subplots_adjust(top=1)

plt.savefig('polynomial_overfit_assymptotic.svg', facecolor='none', edgecolor='none')

# %%
# Train and test set with various complexity in the polynomial degree

plt.figure()
plt.scatter(x, y, s=20, color='k')
plt.scatter(x_test, y_test, s=20, color='C1')

t = np.linspace(-1, 1, 100)

for d in (1, 2, 5, 9):
    model = make_pipeline(PolynomialFeatures(degree=d), LinearRegression())
    model.fit(x.reshape(-1, 1), y)
    plt.plot(t, model.predict(t.reshape(-1, 1)), label='Degree %d' % d,
             linewidth=4)

    style_figs.no_axis()
    plt.legend(loc='upper center', borderaxespad=0, borderpad=0)
    plt.subplots_adjust(top=.96)
    plt.ylim(-.74, 2.1)

    plt.savefig('polynomial_overfit_test_%d.svg' % d, facecolor='none',
                edgecolor='none')

# %%
# Validation curves
from sklearn import model_selection

N_SAMPLES = 150

rng = np.random.RandomState(0)
x = 2 * rng.rand(N_SAMPLES) - 1

y = f(x) + .4 * rng.normal(size=N_SAMPLES)

param_range = np.arange(1, 15)

train_scores, test_scores = model_selection.validation_curve(
    model, x[::2].reshape((-1, 1)), y[::2],
    param_name='polynomialfeatures__degree',
    param_range=param_range,
    cv=model_selection.ShuffleSplit(n_splits=20, test_size=.5,
                                    random_state=1),
    scoring='r2')

plotted_degrees = [1, 2, 5, 9, 15]
for i, degree in enumerate(plotted_degrees):
    plt.figure(figsize=(4.5, 3))
    if degree > 1:
        symbol_train = '--'
        symbol_test = ''
    else:
        symbol_train = 'o'
        symbol_test = 'o'
    plt.plot(param_range[:degree],
             -np.mean(test_scores, axis=1)[:degree],
             symbol_test, color='C1',
             label='Test\nerror')
    plt.plot(param_range[:degree],
             -np.mean(train_scores, axis=1)[:degree],
             symbol_train, color='k',
             label='Train\nerror')

    ax = plt.gca()
    for s in ('top', 'right'):
        ax.spines[s].set_visible(False)

    plt.ylim(ymin=-.8, ymax=.5)
    plt.xlim(0.5, 15)
    plt.legend(loc='upper right', labelspacing=1.5)

    plt.yticks(())
    plt.xticks(plotted_degrees[:i+1])
    plt.ylabel('Error')
    plt.xlabel('Polynomial degree             ')
    plt.subplots_adjust(left=.08, bottom=.23, top=.99, right=1.04)

    plt.savefig('polynomial_validation_curve_%i.svg' % degree,
                facecolor='none', edgecolor='none')

# %%
# Learning curves
rng = np.random.RandomState(0)
x = 2 * rng.rand(100 * N_SAMPLES) - 1

y = f(x) + .4 * rng.normal(size=100 * N_SAMPLES)

X = x.reshape((-1, 1))

np.random.seed(42)

def savefig(name):
    " Ajust layout, and then save"
    ax = plt.gca()
    for s in ('top', 'right'):
        ax.spines[s].set_visible(False)
    plt.ylim(-.65, .0)
    plt.xlim(.5*train_sizes.min(), train_sizes.max())
    plt.xticks((100, 1000), ('100', '1000'), size=13)
    plt.yticks(())

    plt.ylabel('Error')
    plt.xlabel('Number of samples      ')
    plt.subplots_adjust(left=.07, bottom=.22, top=.99, right=.99)
    plt.savefig(name, edgecolor='none', facecolor='none')

# Degree 9
model = make_pipeline(PolynomialFeatures(degree=9), LinearRegression())
train_sizes, train_scores, test_scores = model_selection.learning_curve(
    model, X, y, cv=model_selection.ShuffleSplit(n_splits=20),
    train_sizes=np.logspace(-2.5, -.3, 30))

idx_to_plot = [0, 7, 19, 29]

for i in idx_to_plot:
    n_train = train_sizes[i]
    if i > 0:
        symbol_train = '--'
        symbol_test = ''
    else:
        symbol_train = 'o'
        symbol_test = 'o'
    plt.figure(figsize=(4.5, 3))
    test_plot = plt.semilogx(train_sizes[:i+1],
                             -np.mean(test_scores, axis=1)[:i+1],
                             symbol_test,
                             color='C1')
    train_plot = plt.semilogx(train_sizes[:i+1],
                              -np.mean(train_scores, axis=1)[:i+1],
                              symbol_train, color='k')
    leg1 = plt.legend(['Test error', 'Train error'],
                      loc='upper right', borderaxespad=-.2)
    savefig('polynomial_learning_curve_%i.svg' % n_train)

# %%
# Training with varying sample size
t = np.linspace(-1, 1, 100)

d = 9
model = make_pipeline(PolynomialFeatures(degree=d), LinearRegression())
for i in idx_to_plot:
    plt.figure()
    n_train = train_sizes[i]
    plt.scatter(x[::2], y[::2], marker='.', s=20, color='C1', alpha=.1)
    plt.scatter(x[:min(n_train, 3000)], y[:min(n_train, 3000)], s=20, color='k')

    model.fit(x[:n_train].reshape(-1, 1), y[:n_train])
    plt.plot(t, model.predict(t.reshape(-1, 1)), label='Degree %d' % d,
             linewidth=4, color='C3')

    style_figs.no_axis()
    plt.legend(loc='upper center', borderaxespad=0, borderpad=0)
    plt.subplots_adjust(top=.96)
    plt.ylim(-.74, 2.1)

    plt.savefig('polynomial_overfit_ntrain_%i.svg' % n_train,
                facecolor='none',
                edgecolor='none')

# %%
# Various notions of model complexity
N_SAMPLES = 50

rng = np.random.RandomState(0)
x = 2 * rng.rand(N_SAMPLES) - 1

y = f(x) + .4 * rng.normal(size=N_SAMPLES)

from sklearn import tree

for degree in (4, 16):
    plt.figure(figsize=(.8*4, .8*3), facecolor='none')
    plt.clf()
    ax = plt.axes([.1, .1, .9, .9])

    poly = make_pipeline(PolynomialFeatures(degree=degree),
                                            LinearRegression())
    decision_tree = tree.DecisionTreeRegressor(
        max_depth=int(np.log2(degree)))

    poly.fit(x.reshape((-1, 1)), y)
    decision_tree.fit(x.reshape((-1, 1)), y)

    plt.scatter(x, y,  color='k', s=9, alpha=.8)

    plt.plot(t, poly.predict(t.reshape((-1, 1))),
            color='C3', linewidth=3, label='Polynome')
    plt.plot(t, decision_tree.predict(t.reshape((-1, 1))),
            color='C0', linewidth=3, label='Decision Tree')

    plt.axis('tight')
    plt.ylim(-.74, 2.1)
    style_figs.no_axis()
    plt.legend(loc='upper center', borderaxespad=0, borderpad=0,
               handlelength=1)

    plt.savefig('different_models_complex_%i.svg' % degree,
                facecolor='none', edgecolor='none')

# %%
# Simple figure to demo overfit: with our polynomial data

N_SAMPLES = 50

rng = np.random.RandomState(0)
x = 2 * rng.rand(N_SAMPLES) - 1

y = f(x) + .4 * rng.normal(size=N_SAMPLES)

x_test = 2 * rng.rand(10*N_SAMPLES) - 1

y_test = f(x_test) + .4 * rng.normal(size=10*N_SAMPLES)

plt.figure(figsize=(.8*4, .8*3), facecolor='none')
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

from sklearn import linear_model
# Create linear regression object
regr = linear_model.LinearRegression()
regr.fit(x.reshape((-1, 1)), y)

plt.scatter(x, y,  color='k', s=9)

plt.plot([-1, 1], regr.predict([[-1, ], [1, ]]),
        linewidth=3)

plt.axis('tight')
plt.xlim(-1, 1)
ymin, ymax = plt.ylim()
style_figs.light_axis()
plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('ols_simple.svg', facecolor='none', edgecolor='none')

plt.scatter(x_test, y_test,  color='C1', s=9, zorder=20, alpha=.4)
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

plt.savefig('ols_simple_test.svg', facecolor='none', edgecolor='none')

# %%
# Plot cubic splines
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

from scipy import interpolate
f = interpolate.interp1d(x, y,
                         kind="quadratic",
                         bounds_error=False, fill_value="extrapolate")
plt.scatter(x, y,  color='k', s=9, zorder=20)
x_spline = np.linspace(-1, 1, 600)
y_spline = f(x_spline)
plt.plot(x_spline, y_spline, linewidth=3)

plt.axis('tight')
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

style_figs.light_axis()

plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('splines_cubic.svg', facecolor='none', edgecolor='none')

plt.scatter(x_test, y_test,  color='C1', s=9, zorder=20, alpha=.4)
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

plt.savefig('splines_cubic_test.svg', facecolor='none', edgecolor='none')

# %%
# Simple figure to demo overfit: with linearly-generated data

N_SAMPLES = 50

rng = np.random.RandomState(0)
x = 2 * rng.rand(N_SAMPLES) - 1

y = regr.coef_ * x + regr.intercept_ + .4 * rng.normal(size=N_SAMPLES)

x_test = 2 * rng.rand(10*N_SAMPLES) - 1

y_test = regr.coef_ * x_test + regr.intercept_ + .4 * rng.normal(size=10*N_SAMPLES)

plt.figure(figsize=(.8*4, .8*3), facecolor='none')
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

from sklearn import linear_model
# Create linear regression object
regr = linear_model.LinearRegression()
regr.fit(x.reshape((-1, 1)), y)

plt.scatter(x, y,  color='k', s=9)

plt.plot([-1, 1], regr.predict([[-1, ], [1, ]]),
        linewidth=3)

plt.axis('tight')
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)
style_figs.light_axis()
plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('linear_ols.svg', facecolor='none', edgecolor='none')

plt.scatter(x_test, y_test,  color='C1', s=9, zorder=20, alpha=.4)
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

plt.savefig('linear_ols_test.svg', facecolor='none', edgecolor='none')

# %%
# Plot cubic splines
plt.clf()
ax = plt.axes([.1, .1, .9, .9])

from scipy import interpolate
f = interpolate.interp1d(x, y,
                         kind="quadratic",
                         bounds_error=False, fill_value="extrapolate")
plt.scatter(x, y,  color='k', s=9, zorder=20)
x_spline = np.linspace(-1, 1, 600)
y_spline = f(x_spline)
plt.plot(x_spline, y_spline, linewidth=3)

plt.axis('tight')
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

style_figs.light_axis()

plt.ylabel('y', size=16, weight=600)
plt.xlabel('x', size=16, weight=600)

plt.savefig('linear_splines.svg', facecolor='none', edgecolor='none')

plt.scatter(x_test, y_test,  color='C1', s=9, zorder=20, alpha=.4)
plt.xlim(-1, 1)
plt.ylim(ymin, ymax)

plt.savefig('linear_splines_test.svg', facecolor='none', edgecolor='none')


text: This directory contains didactic figures and scripts that generate them.


text:
class: center, middle

# MOOC Machine learning with scikit-learn

???

Hi, welcome to the MOOC on machine learning with scikit-learn

---

class: titlepage

.header[MOOC Machine learning with scikit-learn]

# What is machine learning?

A first, introductory lesson, focusing on general concepts rather than
coding or maths.

???

In this lesson, we will discover what is machine learning and its
general underlying concepts. This lesson is an introduction and focuses
on general concepts, rather than coding or maths.

---
# What is machine learning?

Machine learning deals with building predictive models

???

So, what is machine learning?

In a nutshell, machine learning is about building predictive models. I
will explain what we mean by predictive models later on.

---
class: center, middle

# Why and when?

Some examples of machine learning

???

First of all, I would like to discuss : Why and when would we use machine learning?

Let me start with a couple examples.

---

# Which iris is that?

???

*Pick up a iris*,

Consider flowers, for instance. I am not very good at recognizing them.
This is an iris. But what type?

--

.tight[
* Setosa
* Versicolor
* Virginica
]

???

There are three classic types of iris,

> Setosa - Versicolor - Virginica

Maybe I can use machine learning to build mathematical rules to set them
apart?

For this, I will describe the irises with numbers. I can measure the
dimensions of their sepal and petal. *Pick up a ruler, and pretend to
measure*

--

| Sepal length | Sepal width | Petal length | Petal width | Iris type  |
| ------------ | ----------- | ------------ | ----------- | ---------- |
| 6cm          | 3.4cm       | 4.5cm        | 1.6cm       | versicolor |
| 5.7cm        | 3.8cm       | 1.7cm        | 0.3cm       | setosa     |
| 6.5cm        | 3.2cm       | 5.1cm        | 2cm         | virginica  |
| 5cm          | 3.cm        | 1.6cm        | 0.2cm       | setosa     |

???

Now I have a set of numbers describing irises. Can I build decision rules
from these numbers? Maybe. See, for instance, it seems that setosa irises have a
very small petal.

---

# Is this person rich?

???

Let us consider a problem closer to a business case; that of estimating
the income of someone.

--

US Census data:

.very-small[

| Age | Workclass | Education    | Marital-status     | Occupation         | Relationship | Race  | Sex  | Capital-gain | Hours-per-week | Native-country | Class |
| --- | --------- | ------------ | ------------------ | ------------------ | ------------ | ----- | ---- | ------------ | -------------- | -------------- | ----- |
| 25  | Private   | 11th         | Never-married      | Machine-op-inspct  | Own-child    | Black | Male | 0            | 40             | United-States  | <=50K |
| 38  | Private   | HS-grad      | Married-civ-spouse | Farming-fishing    | Husband     | White  | Male | 0            | 50             | United-States   | <=50K |
| 28  | Local-gov | Assoc-acdm   | Married-civ-spouse | Protective-serv    | Husband      | White | Male | 0            | 40             | United-States   | >50K  |
| 44  | Private   | Some-college | Married-civ-spouse | Machine-op-inspct  | Husband      | Black | Male | 7688         | 40             | United-States   | >50K  |

]

???

For this, we can use data from the US census.

As you see, it gives is fairly diverse information about individuals : *age, workclass, education, marital-statut, occupation, relationship, and more information*, as well as whether they earn more or less than 50 thousand dollars a year.
We can use machine learning to build rules that will predict this income
information from the rest of the demographic information.

You can also see that, even in a simple example like this, having
intuitions about the information available on many individuals can be a
bit challenging. Data visualization will come in handy.

---
# Engineering rules: data versus experts

Expert knowledge: setosa irises have small petals

???

So, machine learning is about engineering decision rules from the data.
But experts can also engineer decision rules from their knowledge of the
problem. For instance, flower experts know that setosa irises have small petals.

--

	Sepal length
	Sepal width
	Petal length
	Petal width

This rule can be inferred from the data

???

This rule can also be inferred by looking at the data. If we look at the
distribution of the sepal and petal measurements for the different iris
species, we can see that petal length and width single out well setosa.

*But how can we deal with more complex, more numerous data? *

**The benefit of machine-learning is that it automates the creation of rules from the data, including their details, such here as where exactly
to set the threshold on the petal length.**

*In concrete terms, machine learning builds the prediction rules from the data.*

---
class: center, middle

# Predictive analysis

Beyond classic statistical tools

???
We will now get to the heart of the matter by focusing on predictive models.

To build predictive models we are going to rely on statistical
properties of the data. However, the questions and the tools that we
will use differ a bit from those traditionally used in statistics.

---
# Generalizing

Concluding on new instances

???

In machine learning, we want to conclude on new instances. In the
example on the census, I want to be able to predict the income of new
individuals, with a combination of jobs and demographics that I have
never seen.

--

Many sources of variability:

.tight.columns[
* age
* workclass
* education
* marital-status
* occupation
* relationship
* race
* sex
* capital-gain
* hours-per-week
* native-country
]

???

The challenge is that there are many ways an individual can
vary, even within the limited description given by our data.

--

*+* Noise: unexplainable variance

???

An additional difficulty is the noise in the data. By noise, we mean the
aspects that cannot by explained solely from the data. For instance,
an individual's income may have been influenced by the mood of his manager
during his annual review, which is not in our database.

---
# Memorizing

* store all known individual .small[(the census)]
* given a new individual, predict the income of its closest match in our database

???

One possibility for prediction is to memorize the whole available data.

Given the problem of predicting income, we can store all known
individuals, the census. Then, given a new individual, we predict the
income of its closest match in our database.

This strategy is known is machine learning as a "nearest neighbor"
predictor.

--

Trying out this strategy on the data we have, the census, *what error
rate do we expect?*

???

If we try this strategy on the data we have, the census, what error rate
do we expect?

--

> **0 errors**

???

Every individual for which we are asking a prediction is in our database.
Thus, its closest match will be itself, and as a consequence, we will
have zero prediction error.

--

.red[Yet, we will make errors on **new** data]

???

However, if we try our strategy on unseen data, it will not be able to
find exact matches. Hence, it is likely to make some errors.

---

.center.large.reversed[**Generalizing  ≠  Memorizing**  ]

???

What we are seeing here is that generalizing is very different from
memorizing. This is a fundamental challenge of machine learning.

--

.center["test" data  ≠  "train" data]

.pull-left[Data on which the predictive model is applied]

.pull-right[Data used to learn the predictive model]

* Different sampling of noise
* Unobserved combination of features

???

The reason is that the data on which we apply the predictive model, known
as the test data, is different from the data used to build the predictive
model, known as the train data.

They are different because there might be different noise, but also
because there might be individuals with new configuration of features
that we have not observed: different combinations of occupation, age, or
marital status.

---

# The machine-learning workflow

???

The typical workflow in machine-learning is to use a given dataset to
learn a predictive model, for instance predicting the type of iris, and
then to apply it to new data, or what we will call a "test" set, to put
the model in production, or to check its validity.

---

class: center, middle

# Some vocabulary

???

To go further, it is useful to define a bit the concepts of machine
learning

---

# The data matrix

We deal with a table of data:

* Rows are different observations, or **samples**
* Columns are different descriptors, or **features**

.small[
| Sepal length | Sepal width | Petal length | Petal width | Iris type  |
| ------------ | ----------- | ------------ | ----------- | ---------- |
| 6cm          | 3.4cm       | 4.5cm        | 1.6cm       | versicolor |
| 5.7cm        | 3.8cm       | 1.7cm        | 0.3cm       | setosa     |
| 6.5cm        | 3.2cm       | 5.1cm        | 2cm         | virginica  |
| 5cm          | 3.cm        | 1.6cm        | 0.2cm       | setosa     |
]

???

All the data that we will consider will be in what we call a "data
matrix". It can be seen as describing our problem in a 2D table:

* the different rows of this table are different observations, for
  instance different irises. We call these "samples"

* the columns of this table give the different measures or descriptors
  that we have for these samples. We call them "features".

---

# Supervised machine learning

* A data matrix *X* with *n* observations

* A **target** *y*: a property of each observation

The goal is to **predict** *y*

???

In supervised machine learning,
the data that we have are annotated. In other words, they are associated with a label or a target class.
For Irises, each data point is associated with a type of Iris: our target class.
The purpose of supervised learning is to predict this target (here the type of iris) on new data without annotation: entering only petal and sepal dimensions.

In mathematical terms, for supervised machine learning, we are given

* a data matrix, that we shall denote *X*, with *n* observations

* and a **target** *y*, which gives a characteristic for each observation

The goal of supervised learning is to predict y from *X*.

---

# Unsupervised machine learning

* A data matrix *X* with *n* observations

The goal is to extract from *X* structure that generalizes.

Very wide variety of different problems.

???

In unsupervised learning, we are given the data matrix *X*, but we have
no available target.

The goal is then to extract some form of structure from *X* that
generalize to new data.

If we take the example of Irises, in unsupervised learning, the input data would not include the type of Iris because the data is not annotated.
The goal could therefore be to find similarities and structures within the data, or to group together observations that share common characteristics.

Unsupervised learning covers a very wide variety of different problems. We
will not cover them for now.

---

# Regression and classification

Supervised learning: predicting a target *y*

* **Classification**: *y* is discrete, made of different classes

  *eg* irises: Setosa, Versicolor, Virginica

* **Regression**: *y* is continuous, a numerical quantity

  *eg* wage prediction

???

Going back to supervised learning, we have to predict a target *y*. This
target is a property of our data.

It may be discrete, describing different classes of the data. For
instance with irises, we are trying to predict the type of iris. In such
a situation, we say that it is a classification task.

It may be continuous, describing a numerical property of our observation.
For instance, when trying to predict from the census who is rich, it
would be interesting to predict the income, in dollars. In such a
sitation, we say that it is a regression task.

---

.center[
# Take home messages
]

* Machine learning is about extracting from *data* rules that *generalize*
  to new observations

* We work with:

  - a data matrix *X* *n_samples* x *n_features*
  - a target *y* of length *n_samples*, made of numbers for *regression*
    and of discrete classes for *classification*

???

To summarize, machine learning is about extracting from data rules that
generalize to new observations.

In practice, we will work with a data matrix, that we will call *X*, with
*n_samples* rows times *n_features* columns. For supervised learning,
we have a target vector *y*, of length *n_samples*, which is made of
numbers characterizing each observation for regression problems and of
discrete classes for classification problems.


text:
class: titlepage

.header[MOOC Machine learning with scikit-learn]

# Overfit and underfit

This lesson covers overfit and underfit. These are concepts important to
understand why a model generalizes well or not to new data

???

This lesson covers overfit and underfit, important concepts to understand
why a model generalizes well or not to new data

---
# Which data fit do you prefer?

.shift-left.pull-left[]

.pull-right[]

???

Let me start with a simple question: given the following data, which of
the two models do you prefer? I'll give you a moment to think about it.

Most people reply that they prefer the one with the straight line.
However, the one with the wiggly line fits the data perfectly, doesn't it?
So why might we prefer the straight line?

---
# Which data fit do you prefer?

.shift-left.pull-left[]
.pull-right[]
.centered.reversed[**On new data**]

???

Answering this question might be hard. However, in the context of machine
learning we aim for models that generalize. Hence, the good way to frame
the question is: how will the model perform on new data?

---
# Which data fit do you prefer?

.shift-left.pull-left[]
.pull-right[]
.centered[A harder example]

???

How about a slightly harder example? Which one should we choose?
This is a difficult question.

In this lesson, we will study useful concepts to understand these
tradeoffs.

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

???

In the latest example, we have generated the data so that y is a
9th-degree polynomial function of X, with some additional noise.

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

* Fit polynomials of various degrees

???

What we will now do is fit to this data polymonials of various degrees.
We'll start with a polynomial of degree 1: a simple linear regression of
y on X. Clearly, this model does not explain well the data.

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

* Fit polynomials of various degrees

???

If we fit a polynomial of degree 2, the fit is better.

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

* Fit polynomials of various degrees

???

Polynomial of degree 5: it's hard to tell whether it explains the data
better or not

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

* Fit polynomials of various degrees

???

And if we go all the way up to degree 9, the fit actually looks really
bad

---
# Varying model complexity

.polynomial[]

* Data generated with 9th-degree polynomial + noise

* Fit polynomials of various degrees

???

The actual function that was used to generate the data looks like this,
though we added observational noise.

---
# Overfit: model too complex

.pull-left.shift-left[]

.pull-right.width50.shift-left[.shift-left[Model too complex for the data:]

* Its best possible fit would approximate well the generative process

* But its flexibility captures noise
]

???

In the case of the polynomial of degree 9, the problem that we face is
that the model that we use is too complex for the data at hand. This
problem is know as overfit in machine learning. With such a rich model,
its best possible fit would approximate well the data-generating process.
Indeed, here we are fitting a polynomial of degree 9 on data generated
with a polynomial of degree 9. However, due to limited data,
the model fit captures noise because it is too flexible.

--

.reversed[**Not enough data**     **Too much noise**]

???

This problem is typically encountered when there is not enough data, or
too much noise.

---
# Overfit: variance

.pull-left.shift-left[]

.pull-right.width50.shift-left[]

???

The problem of overfit is one of variance: on average, the predictions
are not really off, but each tends to fall far from the target. This can
be seen by their large spread around the best possible prediction. A useful
mental picture is that of the spread of arrows on a target.

---
# Underfit: model too simple

.pull-left.shift-left[]

.pull-right.width50.shift-left[Model too simple for the data:

* Its best fit does not approximate well the generative process

* Yet it captures little noise
]

???

At the opposite end of the spectrum, when we are fitting a polynomial of
degree 1, the model is too simple for the data at hand. We say that it
underfits. Its best possible fit cannot approximate well the
data-generating process. On the positive side, it captures little noise,
As a consequence even with limited data, the empirical fit is close to
the best possible fit on an infinite amount of data.

--

.reversed[**Plenty of data**     **Low noise**]

???

Underfit is more common when there is plenty of data compared to the
complexity of the model, or in low-noise situations.

---
# Underfit: bias

.pull-left.shift-left[]

.pull-right.width50.shift-left[]

???

Unlike overfit, underfit leads to systematic biases:
the predictions cannot be on target on average, because the model that we
use to predict is systematically off the data-generating process.

---
# Underfit versus overfit

.shift-left.pull-left[]

.pull-right[]

.shift-up.pull-left.shift-left[.centered.reversed[Bias]]
.shift-up.pull-right.width50[.centered.reversed[Variance]]

???

So we have these two opposit behaviors:
* underfit, with systematic bias
* and overfit, with large variance

The challenge is to find the right tradeoff between the two.

---
# Underfit versus overfit

.pull-left.width50[]
.pull-right.width50.shift-left[]

.shift-up.pull-left.shift-left[.centered.reversed[Bias]]
.shift-up.pull-right.width50[.centered.reversed[Variance]]

???

This bias-variance tradeoff is classic in statistics. Often, adding a
little bit of bias helps reducing the variance. For instance, as with
throwing darts at a target, where throwing the darts less strong might
lead to being below the target on average, but with less scatter.

---
class: center, middle

# How to choose this tradeoff

Between underfit and overfit

???

So, now that we understand that there is tradeoff between underfit and
overfit, the question is: from a practical standpoint, how do we choose
it?

---

# Train vs test error

.shift-left.pull-left[]
.width50.pull-right[
Measure:

* errors on test data (generalization)

* errors on the train data
]

???

To probe the tradeoff between underfit and overfit, our central
tool will be to measure both the generalization error, on unseen
test data, and the error on the data used to train the model

---

# Train vs test error: increasing complexity

.shift-left.pull-left[]
.width50.pull-right[]

???

We can look at these errors while varying model complexity.

If we start with a very simple model, the training error is similar to
the testing error: the model does not have enough capacity to capture
noise.

---

# Train vs test error: increasing complexity

.shift-left.pull-left[]
.width50.pull-right[]

???

As we increase model complexity, both the train and the test errors go
down: the model explains well the data, but does not capture noise yet.

---

# Train vs test error: increasing complexity

.shift-left.pull-left[]
.width50.pull-right[]

???

Even more complex models fit the training data better, but they start
capturing noise in their model fit. As a result, their error on the test
data is larger.

---

# Train vs test error: increasing complexity

.shift-left.pull-left[]
.width50.pull-right[]

???

As we keep increasing model complexity, the train error keeps going down,
but the test increases sharply. The model is overfitting.

---

# Train vs test error: increasing complexity

.centered[.shift-left[]]

.shift-up.centered[
.reversed[Underfit]
.reversed[Sweet spot]
.reversed[Overfit]
]

???

The big picture is that models that are too simple have similar train and
test error, while models that are too complex have a small train error
but a very large test error. There is a sweet spot in the middle, and
this is where good machine-learning models lie.

---
class: center, middle

# With varying sample size

Learning curves

???

Another useful way to look at the tradeoff between underfit and overfit
is with varying sample size. Such an analysis is often known as a
learning curve.

---

# Train vs test error: varying sample size

.shift-left.pull-left[]
.width50.pull-right[]

.reversed[Overfit]

???

If we fit a polynomial of degree 9 on a small dataset, we will not have
enough data and we will be in an overfit situation. As a result the train
error will be low, but the test error will be high.

---

# Train vs test error: varying sample size

.shift-left.pull-left[]
.width50.pull-right[]

???

As we increase the sample size, the test error decreases while
the train error increases: the models are less overfitting.

---

# Train vs test error: varying sample size

.shift-left.pull-left[]
.width50.pull-right[]

.reversed[Sweet spot?]

???

With enough data, the train and test errors converge: the model is no
longer overfit.

---

# Train vs test error: varying sample size

.shift-left.pull-left[]
.width50.pull-right[]

.reversed[Diminishing returns]

???

After a while, we reach diminishing returns: the test and the train error
are no longer changing: more data does not bring any benefits.

--

Try more complex models?

???

Given that more data are not improving model prediction, it may be useful
to try more complex models.

---

# Bayes rate: no improvement possible

.shift-left.pull-left[]
.width50.pull-right[
Here, the data-generating process is a degree-9 polynomial

We cannot do better

Predictions limited by noise

]

???

In general, when adding more data does not provide improvements,
it is useful to try more complex models.

However, in the specific case of our example, the data-generating process
in a degree-9 polynomial. As a result, more complex models will not
improve the prediction: the present model has already captured all the
non-random link between X and y.

The corresponding prediction is imperfect. Yet, we cannot do better: the
prediction is limited by the intrinsic randomness of the link between X
and y.

---
class: middle

.center[
# Model families
]

.centered[
Crucial to match

.tight[
* statistical model

* data-generating process
]
]

.centered[So far: polynomes for both]

???

The excellent performance that we just observed for degree-9 polynomes
when there is plenty of data comes from the perfect match between the
statistical model used to analyse the data and the data-generating process. So far, we have used
polynomes for both of them. But in practice, given data, we seldom know
a simple form of model in which the data is drawn.

For this reason, the choice of family of model is crucial.

---

# Different model families

.shift-left.pull-left[]
.pull-right.width[
* Different "inductive bias"

* Different notion of "complexity"
]

???

Different model families come with different form of bias, which we call
inductive bias, and of complexity.

For instance, polynomial regression tend to be smooth, and their
complexity is controlled by the degree of the polynome used to fit the
data. On the other hand, decision tree are locally constant. There
complexity is controlled by the number of constant regions that they fit.

---

# Different model families

.shift-left.pull-left[]
.pull-right[]

.shift-up.pull-left[.reversed[Simple variant]]
.shift-up.pull-right[.reversed[Complex variant]]

.centered[⇦ regularization ⇦]

???

Varying model complexity for these different families of model explores
different underfit and overfit behavior.

In general forcing a model towards a simpler fit is called
"regularization". In scikit-learn, and in machine learning in general,
almost every model comes with a specific regularization strategy. Part of
machine-learning practice is about understanding the inductive biases and
the regularization strategy to choose the right model. We will cover this
later, as we explore various models.

---

.center[
# Take home messages
]

Models too complex for the data **overfit**:

.tight[
- they explain too well the data that they have seen
- they do not generalize
]

Models too simple for the data **underfit**:

.tight[
- they capture no noise
- but are limited by their expressivity
]

Comparing train and test error reveals this tradeoff

Different models = different complexity and bias

???

Let me conclude this lesson by recapitulating the take home messages.

When the models are too complex for the data at hand, they overfit. This
means that they explain too well the data that they have seen as they
capture noise, and thus do not generalize to new data.

On the opposite, when models are too simple for the data at hand, they
underfit. This means that they capture no noise, but their ability to
generalize is then limited by their expressivity.

Comparing train and test error, for instance varying model complexity, or
amount of data, can reveal this tradeoff.

The notion of complexity, or the bias imposed by a model, are both
specific to a model family. We will later study different model families.


text:
<!DOCTYPE html>
<html>
<head>
  <title>Presentation</title>
  <meta charset="utf-8">
  <link rel="stylesheet" type="text/css" href="custom.css">
</head>
<body>
<textarea id="source"></textarea>
<script src="https://remarkjs.com/downloads/remark-latest.min.js">
</script>
<script>
  var slideshow = remark.create({
    ratio: '16:9',
    slideNumberFormat: '(%current%/%total%)',
    countIncrementalSlides: false,
    highlightLines: true,
    sourceUrl: getUrlParameter('file') ? getUrlParameter('file') : 'slides.md'
  });

  function getUrlParameter(name) {
    name = name.replace(/[\[]/, '\\[').replace(/[\]]/, '\\]');
    var regex = new RegExp('[\\?&]' + name + '=([^&#]*)');
    var results = regex.exec(location.search);
    return results === null ? '' : decodeURIComponent(results[1].replace(/\+/g, ' '));
  };
</script>
</body>
</html>


text:
.remark-slide-number {
    display: none
}

div.remark-slide {
    font-family: "Lucida Grande", "Trebuchet MS", Verdana, sans-serif;
}

div.remark-slide-content.middle {
    padding-top: 0em;
    padding-bottom 3em;
}

div.remark-slide-content.center {
    padding-left: 4em;
    padding-right 4em;
}

div.remark-slide p, div.remark-slide li, .remark-slide-content {
    font-size: 35px;
}

div.remark-slide-content {
    padding-left: 2.5em;
}

div.remark-slide-content:not(.center) h1 {
    margin-right: -2em
}

div.remark-slide-content ul {
    margin-top: .5ex;
}

div.remark-slide-content p {
    margin-bottom: .5ex;
}

div.remark-slide-content li p {
    margin-top: .5ex;
}

div.remark-slide-content li {
    margin-top: 1ex;
}

div.remark-slide-content li p:first-child {
    margin-top: 1.75ex;
}

div.remark-slide-content li li {
    margin-top: .5ex;
}

thead {
    background: #daf2ff;
}

tbody tr:nth-child(even){
  background-color: #e9e9e9;
}

body.remark-container {
    background-color: #eee;
    /* scikit-learn ish gradient */
    background: linear-gradient(160deg, rgb(150, 215, 249) 0%, rgb(255, 243, 211) 100%);
}

/* Enlarge the presenter's note display */

.remark-notes-area {
    left: 50%;
}

.remark-container.remark-presenter-mode .remark-slides-area,
.remark-container.remark-presenter-mode .remark-preview-area {
    left: 1%;
}

.remark-notes-area .remark-bottom-area .remark-notes-current-area {
    height: 80%;
}

.remark-notes-area .remark-bottom-area {
  top: 50px;
}

/* Custom css for layout in slides */
div.remark-slide-content.titlepage h1 {
    font-size: 65px;
}

div.remark-slide-content.titlepage img {
    bottom: 4%;
    left: 4%;
    position: absolute;
    width: 66%;
    opacity: 20%;
}

.header {
    margin-top: 0px;
    font-weight: bold;
    font-size: 30px;
    color: #00517b;
}

div.remark-slide .very-small {
    font-size: 20px;
}

div.remark-slide .small {
    font-size: 28px;
}

div.remark-slide .large {
    font-size: 45px;
}

.larger {
    transform: scale(1.2);
    transform-origin: left;
}

.pull-left {
    float: left;
    width: 45%;
}

.pull-right {
    float: right;
    width: 45%;
}

.pull-right p:first-child {
    margin-top: 0px;
}

.pull-right + p, .pull-right + span {
    clear: both;
}

.clear, p + span.clear {
    clear: both;
}

.width30 {
    width: 30%;
}

.width40 {
    width: 40%;
}

.width50 {
    width: 50%;
}

.width60 {
    width: 60%;
}

.shift-left {
    margin-left: -1ex;
    margin-right: -1ex;
    display: block;
}

.shift-up {
    margin-top: -3ex;
}

.shift-up-more {
    margin-top: -5ex;
}

img.shift-up {
    margin-top: -5ex;
}

.red {
    color: #911;
}

span.reversed {
    color: white;
    background-color: black;
    padding: .5ex 1ex;
    border-radius: .4ex;
    display: inline-block;
}

.centered, span.centered {
    margin-left: auto;
    margin-right: auto;
    display: table;
}

/*.columms ul li {
    display: inline-block;
    width: 45%;
}*/

.columns ul li {
    float: left;
    min-width: 10ex;
    width: 28%;
    margin-right: 1.25em;
    vertical-align: text-top;
}

div.tight li {
    margin-top: .2ex;
}

.columns ul:after {
     visibility: hidden;
     display: block;
     font-size: 0;
     content: " ";
     clear: both;
     height: 0;
}

/* Adhoc classes */
.polynomial {
    float: left;
    width: 65%;
    margin-top: -1ex;
    margin-left: -3ex;
}


text:
# Compilation is done via remarker and htmlark, both pip instalable

all: ml_concepts.html overfit.html

%.html: %.md custom.css
	# HTMLArk is to embed images and css
	remarker $< -c custom.css | htmlark --ignore-js > $@


text:
Hi,

Welcome to the Machine-learning with scikit-learn course. The goal of
this course is to teach you practical aspects of machine learning. It
focuses on tabular data, given that such data is often encountered in the
industry. The course is light on maths, and focuses on practical aspects,
not only about pure machine learning, but also about the basics of data
preparation and visualization for machine learning. Most of the content is
centered on executable Python code that teaches how to analyse the data,
with tools such as scikit-learn.

Our goal is to be didactic. If you know Python programming and basic
numerics, you should be able to follow along. We hope that this course
will help introduce more people to machine learning.


text:
remarker
htmlark[parsers]


text:
# Preview

To preview slides `https://lesteve.github.io/scikit-learn-tutorial/slides/?file=[FILENAME].md`

Example : https://lesteve.github.io/scikit-learn-tutorial/slides/?file=ml_concepts.md

# Export

To install packages needed to generate the slides:

```
pip install -r requirements.txt
```

Note: for some reason if you `pip install` `htmlark` and not `html[parsers]`, you'll get a blank HTML page.

Then use `make` to export html files.


text:
`cps_85_wages.csv` is available at https://www.openml.org/d/534
`adult-census.csv` is available at https://www.openml.org/d/15950


text:
import sys
import os

# TODO: we could get the list from .gitignore
IGNORE_LIST = [
    '.ipynb_checkpoints',
]

folder = sys.argv[1]
non_python_files = [
    fn
    for fn in os.listdir(folder)
    if fn not in IGNORE_LIST and not fn.endswith('.py')
]

if non_python_files:
    raise RuntimeError(f'Looks like you have non python files in {folder}:\n'
                       f'{non_python_files}')


text:
import os
import sys
import difflib

# TODO: we could get the list from .gitignore
IGNORE_LIST = [
    '.ipynb_checkpoints',
]

folder1, folder2 = sys.argv[1:3]

def get_basename(folder):
    contents = []
    for fn in os.listdir(folder):
        content = os.path.splitext(os.path.basename(fn))[0]
        if content not in IGNORE_LIST:
            contents.append(content)
    return contents

basenames1 = sorted(get_basename(folder1))
basenames2 = sorted(get_basename(folder2))

if basenames1 != basenames2:
    only_in_folder1 = set(basenames1) - set(basenames2)
    only_in_folder2 = set(basenames2) - set(basenames1)

    raise RuntimeError(
        f'Inconsistency between folder {folder1} and {folder2}\n'
        f'Only in folder {folder1}: {only_in_folder1}\n'
        f'Only in folder {folder2}: {only_in_folder2}')