From news bursts to story graphs, with local LLMs

This post traces how tweets or news stories spread through the media and whether we can visualise the spread of news or items visually. Typically, items start with small bursts and then become full-blown stories as more people connect and discuss it online.

The plan was simple to state:

posts → people, places, organisations and claims → a graph that changes over time → “what’s happening?”

Originally, I wanted to use X f(ormerly Twitter) and listen to tweets live, but the cost of the X stream is prohibitive, so I switched to Bluesky whose whole public stream is free to read.

The the idea was to identify news/event bursts and then see how the event or news item actually becomes a full-blown story. The interesting part is what happens over the next day. So the project was created in two halves 1) finding bursts, and 2) following each burst into a story.

The setup

Data: Bluesky’s public stream through its free Jetstream feed, English posts only. This was done for four days from 18 to 22 October 2026, and then frozen so every number below comes from the same recording: about 8 million original posts and 7 million reposts.
Everything runs offline on the recording. Recording happens once; every experiment replays the saved files. That made it possible to change a threshold and compare results on identical data.

Flowchart illustrating data capture and replay processes for Bluesky Jetstream, detailing steps for polling news feeds, saving JSONL files, freezing data collection, and replaying files for analysis.

Job	Model (all local)	Why this one
Find names in every post	GLiNER, a small entity model	fast enough for every post, about 140 posts a second
Judge whether posts say the same thing	MiniLM sentence embeddings	about 7,000 posts a second
Extract facts and claims for a story	GLM-4.7-Flash, via Ollama	the fastest model
Tamil news, merging names, captions	Gemma 4 31B, via Ollama	best quality , and it reads Tamil
Store the graph	Neo4j Community edition	free and local

How it works

Record the Bluesky stream to hourly files.
Find the names in each post: people, organisations, places, facilities, products and events.
Detect bursts: a set of names that suddenly many more people are writing about than usual.
Follow the story: collect its posts from 3 hours before the burst to 24 hours after.
Extract facts and claims at 12 points in time with a local model, each tied to the posts that support it.
Merge duplicate names, store everything in Neo4j, and draw the story as an interactive page that plays the 12 steps in order.

Flowchart illustrating data sources, capture and replay methods, and entity processing in a system. It includes information about Bluesky Jetstream, news RSS feeds, data recording, and entity cache operations.

Flowchart illustrating a process for tracking and analyzing story bursts, including steps for identifying co-mention graphs, tracking event writers and reposts, managing facts and claims, and generating knowledge graphs with outputs.

Part 1: Bursts, or finding the moment a story breaks

1. Counting hashtags finds bots

The obvious first step is to count hashtags and words and flag whatever grows fastest. In the first morning of data, the top trends were:

Trend	Posts	Authors
roadworks (a UK traffic bot)	564	1
#kickban	197	1
#greatwhitemike	158	1
nebraska	43	32

Three of the top four were single accounts. Everything from here on counts distinct people, not posts, and nearly identical texts count once, because spam networks post the same template from dozens of accounts (one sent the same advert from 82).

2. Names, not hashtags

Only 11% of Bluesky posts carry a hashtag, names were needed inside the text. Gliner and Spacy entity extractors were compared on 200 posts, using Gemma 4 31B as the reference. Gemma is too slow for every post, at 0.3 posts a second, but it is accurate.

Extractor	Precision	Recall	Types right	Speed
spaCy (small)	0.77	0.67	0.59	about 460 posts a second
GLiNER	0.48	0.81	0.72	about 160 posts a second
GLiNER with a junk filter	0.73	0.72	0.75	about 160 posts a second

Unfiltered, GLiNER labelled “I’m 49”, “Dinner” and “He” as entities. A small filter (no pronouns, no dates, no common phrases in lower case) fixed most of it. Over a full hour, 53% of posts named nothing, 24% named one thing, and only 23% named two or more. Short posts mostly mention one part of a story.

3. One story, seven events

The initial detector linked names that people mention together. On the morning of 18 October, Ukrainian drones struck targets in St Petersburg: an oil terminal, a VK data centre, an Ozon warehouse and a bridge. The detector reported seven separate events. One post said “Ozon”, another “VK data center”, a third “oil terminal”, and almost nobody named two of them in one post.

4. Need for embeddings

Names alone could not tell which posts belonged to the same story, because most posts name only one part of it. I needed a way to compare what posts say, not just which names they share.

An embedding turns a post into a list of numbers, so that posts with similar meaning get similar numbers even when they share no words. “Drones hit an oil terminal in St Petersburg” and “Ozon warehouse on fire after the attack” come out close together; a post about a football score does not. A small local model (MiniLM) does this at about 7,000 posts a second, so it can be applied freely. Reading every post with a language model would understand meaning better, but Gemma manages 0.3 posts a second while 16 to 44 arrive every second.

I ended up using embeddings in four places:

Where	What went wrong with names alone	What embeddings did
Grouping a burst into one event	The strike came out as 7 events, because “Ozon”, “VK data center” and “oil terminal” were rarely named together	Names whose posts read alike were linked, and the strike became one event with 8 names
Following an event from one window to the next	When later posts used new names (“Baltic oil terminal”, “Saint Petersburg”), the strike restarted as a new event	The new group was matched to the existing event by meaning, and the strike stayed one event for 105 minutes
Choosing the posts to read for an event	A post that only said “Nintendo” or “Korea” was pulled into unrelated stories: a 22 year old game, a Korean sports team	A post sharing only one name with the event must also read like the event
Following a story for a day	Every passing mention of a name would have joined the story	Posts naming one of the story’s names must also be close in meaning to the story’s own posts

Two limits showed up. The threshold matters: at a similarity of 0.45, unrelated political stories merged, because political posts all read alike to a small model; at 0.65 the strike split again; 0.55 worked. And embeddings cannot rescue a story whose names are too general, which comes up in Part 2.

Each change was tested on the same three hours of data:

Run	Change	Events found	What happened
1	names mentioned together	173	the strike came out as 7 events, a college football game as 7
2	linking posts that read alike	35	the strike disappeared altogether (see the next section)
4	writers counted per event	60	the strike became one event with 8 names
7	thresholds tuned	61	the strike lasted 105 minutes as one event, and unrelated politics stopped merging

5. The burst formula

Putting the pieces together, this is how a burst is scored. Every 5 minutes, for every name, the detector looks at the last 15 minutes:

authors: the distinct people who wrote or reposted a post naming it
total: the distinct people active at all in those 15 minutes
share: authors divided by total

Each name has its own baseline: its usual share, and how much that share normally varies. Both are running averages that remember roughly the last 6 hours (a half life of 72 steps of 5 minutes):

Each part has a job:

Share, not count. Activity on Bluesky roughly triples between the US night and morning. Counts would make every name look like a burst each morning; shares only move when the conversation does.
The expected term allows for ordinary chance: a name expected 100 times can easily turn up 110 times.
The (total × √var)² term asks more of names that are always busy and always vary, such as Trump.
The + 1 stops a name never seen before from scoring infinity with two people, while a new name with 20 people still scores about 20.
The baseline is updated after scoring, so a burst is never compared with itself.

A name is bursting when at least 8 people named it in the window. Bursting names are grouped into events (by co-mention and by embeddings), and an event must then pass two more tests from the previous section: at least 3 people wrote original posts about it, and for at least one of its names the writers alone score at least 2 on a second baseline that counts writers only. The first 2 hours of any run only build the baselines.

Two examples. At 06:45 IST on 18 October, Ozon was named by 21 people in 15 minutes where about 0.8 would be normal, giving z = 7.6, and VK data center by 11 people where almost none was expected, giving z = 9.1: both bursts. For a constantly busy name, take an illustration: a usual share of 3% of 15,000 active people gives an expected 450; if its share normally varies by 0.4 points, the spread term adds 60, so the denominator is about 64. Then 520 people give z = 1.1, which is ordinary, and it would take about 770 to reach z = 5.

7. The burst is short

Every tech event the detector found on Monday 19 October lasted between 10 and 105 minutes. That is when people suddenly pile in. But the stories themselves ran for a day:

Bar chart showing social media post activity over time for three stories: 'St Petersburg drone strike', 'Prop 40 and Sergey Brin', and 'Nvidia and warehoused GPUs', with indicated burst detection periods.

The shaded band is all the detector saw. For the strike, the biggest wave came six hours later, in the Indian evening, as casualty figures and claims arrived. Prop 40 had a second wave the next evening. A graph built from the burst alone had 2 to 5 checkpoints: a snapshot, not a story.

Part 2: From a burst to a story

1. Following a story for a day

So the burst became a seed. Each detected event gives its main names and the average embedding of its posts. One pass over the recording, from 3 hours before the burst to 24 hours after, collects every post that names two of the story’s names, or names one and reads like the story, along with its repost count. All stories are followed in the same pass:

Story	Posts followed	Reposts of them
Steam Next Fest	2,312	7,763
Nvidia and warehoused GPUs	794	3,398
St Petersburg drone strike	414	3,784
Prop 40 and Sergey Brin	342	4,044
Amazon GPUs	286	2,965

2. Twelve checkpoints, facts and claims

The followed posts are cut into 12 slices with similar numbers of posts, so busy periods get finer steps. For each slice, GLM-4.7-Flash reads the most reposted posts (and the first one) and returns:

Facts as subject → verb → object (“Ukrainian drones → attacked → Ozon warehouse”), each citing the posts that state it. A fact citing no real post is dropped.
Claims, with who makes them: “135 drones were shot down” (the regional authorities), “our services are working normally” (VK), “the warehouse was used by the Russian defence ministry” (Ukraine; Ozon denies it).

Separating facts from attributed claims turned out to be one of the most useful parts. A story is mostly claims: casualty figures from three different sources, accusations, denials.

Infographic showing the impact of a drone attack in St. Petersburg, highlighting targeted infrastructure such as oil terminals, data centers, and warehouses.

3. Many names for one thing

Extraction runs slice by slice, and the same thing comes back under different names: “Vijay”, “CM Vijay”, “C Joseph Vijay” and “TN CM”; “IMD”, “India Meteorological Department” and “weather department”. Three layers fixed most of it:

Prevention: each slice is given the names already used for this story, with an instruction to reuse them, plus a short list of preferred verbs (“attacked”, “located in”, “deployed to” and so on).
Clean up in code: “X says Y” becomes a claim by X, not a fact; a fact linking a thing to itself is dropped; “VK data centre was hit by a drone” is turned round to “drone hit VK data centre”.
One Gemma call per story groups the names that mean the same thing and picks one, while being told that a district is not its city (North Chennai is not Chennai).

On the Chennai rain story (below), this cut 127 different names to 70, and 42 different verbs to 11. Gemma also knew that “Kovai” is Coimbatore and that “the Tamil Nadu Weatherman” is the blogger Pradeep John. It made mistakes too: it merged “weather centre”, which in Tamil reports usually means Chennai’s regional centre, into the national IMD.

4. Captions that say only what’s new

My first pages showed each slice’s summary, and every summary began by restating the whole story. So Gemma now writes one fixed sentence for the story, and one sentence per step with only what is new. For Prop 40, California’s proposed tax on billionaires:

Story: Users debate California’s proposed billionaire wealth tax (Prop 40) and criticise tech billionaires’ ties to Donald Trump and Mar-a-Lago.

Users suggest renaming Mar-a-Lago and discuss Props 41 and 42.

Sergey Brin reportedly spent over $100 million against Prop 40 and dined with Trump at Mar-a-Lago.

Reports that Mark Zuckerberg moved to Florida and Gavin Newsom vacationed with Brin.

Musk, Thiel and Zuckerberg reportedly oppose Prop 40; Bernie Sanders supports it.

A poll suggests Prop 40 faces a “spectacular collapse”.

Brin reportedly moved to Nevada and bought a Miami home; the BBC reportedly demanded Mar-a-Lago’s membership list.

Economist Gabriel Zucman explains Prop 40 as a wealth tax; Mark Ruffalo warns of health insurance losses.

And for the strike, the account firms up over the day: targets first, then the radar site and bridge, then Ozon’s denial, then casualty figures from three different sources (Zelensky, the Russian defence ministry and the city’s mayor, each with different numbers), then “over 1,000 drones, the largest attack on the city”. These are what posts claimed, not checked facts, which is why every claim keeps its source.

A visual map illustrating connections related to California's proposed billionaire wealth tax, featuring Sergey Brin's move to Nevada and a Miami purchase, alongside claims related to Mar-a-Lago and political figures like Donald Trump and Peter Thiel.

The pages play the 12 steps in order, with a speed control (2 to 10 seconds per step) and each step’s caption over the graph. They work on a phone.

5. When following goes wrong

The Nvidia story shows the weakness. Its seed names were “Nvidia”, “GPUs”, “RAM” and “Jensen Huang”, all of which people discuss constantly. So the day’s story became everything about Nvidia: warehoused GPUs, “channel stuffing” accusations, a German AI cloud deal, a new CPU. Every Nvidia post read like the story, so embeddings could not hold it together either. You can see it in the figure above: a flat, steady line rather than a wave. A story can only be followed if its names are specific to it. Steam Next Fest, a festival of hundreds of game demos, had the opposite problem: a few prolific developers filled the evidence with their own game.

A detour to Chennai: a story told in two languages

Heavy northeast monsoon rain hit Chennai on 17 and 18 October while I was recording. Bluesky had almost nothing about it (128 posts in 21 hours, mostly automatic copies of Tamil TV headlines), so I followed 1,621 news articles instead, in English and Tamil, alongside Open-Meteo’s hourly rainfall.

GLM cannot read Tamil. It returned nothing for seven Tamil headlines. Gemma 4 31B read them and wrote every name in English (சென்னை → Chennai, தாம்பரம் → Tambaram), so Tamil and English articles about the same place joined up in the graph.
News arrived the morning after. The heaviest rain fell around 00:30 IST; articles peaked six to nine hours later, and Tamil outlets were out first, from 06:30, with English coverage peaking around midday.

Bar graph showing rainfall (in mm/h) in central Chennai over time with a peak of 16 mm/h. Below, line graphs show the number of Tamil news articles, English news articles, and Bluesky posts about Chennai weather over the same period.

The figures firmed up through the day: 175 mm at Parrys, then “30 cm in Chennai”, then 36.34 cm overnight. The official and critical accounts sat side by side: the Corporation’s “no fatalities” and “drained by dawn, barring Perambur and Velachery”, a minister’s “Chennai cannot withstand 24 hours of continuous rain”, and a resident’s “36 cm fell in 3 hours and the weather centre gave no warning”.

One caveat: the rainfall model under reads badly. It peaked at 16 mm an hour, while the reported totals were ten times higher. It gets the timing right, not the amount.

Play the stories

Each story is an interactive page: press Play to watch it build up step by step, with a caption saying what changed, and use the slider to jump to any step. They work on a phone too. All eight are listed at https://tvganesh.github.io/bluesky-stories/.

St Petersburg drone strike: targets, denials and casualty figures from three sources, over a day.
Prop 40 and Sergey Brin: California’s billionaire tax, Brin’s opposition and his reported move.
Chennai monsoon rain: the story from English and Tamil news, alert by alert.
Nvidia and warehoused GPUs: the story that drifted into everything about Nvidia.

The pages show what the models extracted, not the posts themselves: account names, links to posts and quoted posts have been removed, and posts their authors deleted after I recorded them were dropped.

This was an interesting exercise to see how the stories develop. What would be more interesting is so what cause tweets or stories to go viral. That is definitely food for thought.

Also see

To see all posts click Index of posts

Identifying cricketing shots using AI

Image classification using Deep Learning has been around for almost a decade. In fact, this field with the use of Convolutional Neural Networks (CNN) is quite mature and the algorithms work very well in image classification, object detection, facial recognition and self-driving cars. In this post, I use AI image classification to identify cricketing shots. While the problem falls in a well known domain, the application of image classification in identifying cricketing shots is probably new. I have selected three cricketing shots, namely, the front drive, sweep shot, and the hook shot for this purpose. My purpose was to build a proof-of-concept and not a perfect product. I have kept the dataset deliberately small (for obvious reasons) of just about 14 samples for each cricketing shot, and for a total of about 41 total samples for both training and test data. Anyway, I get a reasonable performance from the AI model.

Included below are some examples of the data set

This post is based on this or on Image classification from Hugging face. Interestingly, this, the model used here is based on Vision Transformers (ViT from Google Brain) and not on Convolutional Neural Networks as is usually done.

The steps are to fine-tune ViT Transformer with the ‘strokes’ dataset are

Install the necessary libraries

! pip install transformers[torch] datasets evaluate accelerate -U
! pip install -U accelerate
! pip install -U transformers

b) Login to Hugging Face account

 from huggingface_hub import notebook_login
notebook_login()

Login successful

c) Load the batting strokes dataset with 41 images

from datasets import load_dataset
df1 = load_dataset("tvganesh/strokes",split='train')
type(df1)
len(df1)

41

df1
Dataset({
    features: ['image', 'label'],
    num_rows: 41
})

d) Create a dictionary that maps the label name to an integer and vice versa. Display the labels

labels = df1.features["label"].names
label2id, id2label = dict(), dict()
for i, label in enumerate(labels):
label2id[label] = str(i)
id2label[str(i)] = label

labels

['front drive', 'hook shot', 'sweep shot']

e) Load ViT image processor. To apply the correct transformations, ImageProcessor is initialised with a configuration that was saved along with the pretrained model

from transformers import AutoImageProcessor

checkpoint = "google/vit-base-patch16-224-in21k"
image_processor = AutoImageProcessor.from_pretrained(checkpoint)

f) Apply image transformations to the images to make the model more robust against overfitting

from torchvision.transforms import RandomResizedCrop, Compose, Normalize, ToTensor

normalize = Normalize(mean=image_processor.image_mean, std=image_processor.image_std)
size = (
    image_processor.size["shortest_edge"]
    if "shortest_edge" in image_processor.size
    else (image_processor.size["height"], image_processor.size["width"])
)
_transforms = Compose([RandomResizedCrop(size), ToTensor(), normalize])

g) Create a preprocessing function to apply the transforms and return pixel_values of the image as the inputs to the model – :

def transforms(examples):
    examples["pixel_values"] = [_transforms(img.convert("RGB")) for img in examples["image"]]
    del examples["image"]
    return examples

h) Apply the preprocessing function over the entire dataset, using Hugging Face Dataset’s ‘with_transform’ method

df1 = df1.with_transform(transforms)

from transformers import DefaultDataCollator
data_collator = DefaultDataCollator()

i) Evaluate model’s performance with evaluate

import evaluate
accuracy = evaluate.load("accuracy")

j) Calculate accuracy by passing in predictions and labels

import numpy as np
def compute_metrics(eval_pred):
    predictions, labels = eval_pred
    predictions = np.argmax(predictions, axis=1)
    return accuracy.compute(predictions=predictions, references=labels)

k) Load ViT by specifying the number of labels along with the number of expected labels, and the label mapping

 from transformers import AutoModelForImageClassification, TrainingArguments, Trainer

model = AutoModelForImageClassification.from_pretrained(
    checkpoint,
    num_labels=len(labels),
    id2label=id2label,
    label2id=label2id,
)

Pass the training arguments to Trainer along with the model, dataset, tokenizer, data collator, and compute_metrics function.
Call train() to finetune your model.

training_args = TrainingArguments(
    output_dir="data_classify",
    remove_unused_columns=False,
    evaluation_strategy="epoch",
    save_strategy="epoch",
    learning_rate=5e-5,
    per_device_train_batch_size=8,
    #gradient_accumulation_steps=4,
    per_device_eval_batch_size=6,
    num_train_epochs=20,
    warmup_ratio=0.1,
    logging_steps=10,
    load_best_model_at_end=True,
    metric_for_best_model="accuracy",
    push_to_hub=True,
)

trainer = Trainer(
    model=model,
    args=training_args,
    data_collator=data_collator,
    train_dataset=train_dataset,
    eval_dataset=test_dataset,
    tokenizer=image_processor,
    compute_metrics=compute_metrics,
)

trainer.train()

Epoch	Training Loss	Validation Loss	Accuracy
1	No log	0.434451	1.000000
2	No log	0.388312	1.000000
3	0.361200	0.409932	0.888889
4	0.361200	0.245226	1.000000
5	0.293400	0.196930	1.000000
6	0.293400	0.167858	1.000000
7	0.293400	0.140349	1.000000
8	0.203000	0.153016	1.000000
9	0.203000	0.116115	1.000000
10	0.150500	0.129171	1.000000
11	0.150500	0.103121	1.000000
12	0.150500	0.108433	1.000000
13	0.138800	0.107799	1.000000
14	0.138800	0.093700	1.000000
15	0.107600	0.100769	1.000000
16	0.107600	0.113148	1.000000
17	0.107600	0.100740	1.000000
18	0.104700	0.177483	0.888889
19	0.104700	0.084438	1.000000
20	0.090200	0.112654	1.000000
TrainOutput(global_step=80, training_loss=0.18118578270077706, metrics={'train_runtime': 176.3834, 'train_samples_per_second': 3.628, 'train_steps_per_second': 0.454, 'total_flos': 4.959531785650176e+16, 'train_loss': 0.18118578270077706, 'epoch': 20.0})

m) Push to Hub

trainer.push_to_hub()

You can try out my fine-tuned model at identify_stroke̱

Here are a couple of trials

As I mentioned before, the model should be reasonably accurate but not perfect, since my training dataset is extremely small. This is just a prototype to show that shot identification in cricket with AI is in the realm of the possible.

References

Do take a look at

To see all posts click Index of posts

GenerativeAI:Using T5 Transformer model to summarise Indian Philosophy

Ever since I started to use ChatGPT, I have been fascinated by its capabilities. To a large extent, the abilities of Large Language Models (LLMs) is quite magical – the way it answers questions, the way it summarises passages, the way it creates poems et cetera. All the LLMs need is a large corpus of data from the internet, articles, wikis, blogs, and so on.

On delving a little deeper into Generative AI, LLMs I learnt that, this is based on the principle of being able to predict the most probable word in a given sequence. It made me wonder whether the world of ideas, language and communication are actually governed by probabilities. Does what we communicate fall within the purview of statistics?

As an aside, just by extending further if we visualise a world in which every human action to a situation is assigned an embedding vector, and if we feed the responses of all humans over time in different situations, to the equivalent of a Transformer of a Large Human Reaction Model (LHRM) ;-), we can envisage the model being capable of predicting the response of human in a given situation. In my opinion, the machine would be fairly right most of the occasions as it could select the most probable choice of action, much like ‘The Machine’ in Person of Interest. However, this does not mean that the machine (AI) is actually more intelligent than humans. All it means is that the choice of humans responses are a part of a finite subset possibilities and The Machine (AI) can compute the possibilities and associated probabilities much quicker than humans. Does it mean that the world is deterministic? Possibly.

In this post, I use the T5 transformer to summarise Indian philosophy. For this task, I have fine-tuned the T5 model with a curated dataset taken from random passages on Hindu philosophy available on the internet. For each passage, I had to and hand-create the corresponding summary. This was a fairly tedious and demanding task but an enlightening one. It was interesting to understand how our ancestors, the Rishis, understood reality, the physical world, senses, the mind, the intellect, consciousness (Atman) and universal consciousness (Brahman). (Incidentally I was only able to curate only about 130 rows of philosophical snippets and manually create the corresponding summaries. Probably this is a very small dataset for fine-tuning but I just wanted to see the performance of the T5 model in a new domain.)

In this post the T5 model is fine-tuned with the curated dataset and the rouge1 and rouge2 scores are used to evaluate the model’s performance.

I have used the Hugging Face Hub for the transformer model, corresponding LLM functions and management of the dataset etc. The Hugging Face ecosystem is simply wow!!

Summarisation with T5-small model

a) Install the necessary libraries

! pip install transformers[torch] datasets evaluate rouge_score accelerate -U
! pip install -U accelerate
! pip install -U transformers

b) Login to Hugging Face account


from huggingface_hub import notebook_login
notebook_login()

Login successful

c) Load the curated dataset on Hindu philosophy

from datasets import load_dataset
df1 = load_dataset("tvganesh/philosophy",split='train')

d) Load a T5 tokenizer to process text and summary

Prefix the input with a prompt so T5 knows this is a summarization task.
Use the keyword text_target argument when tokenizing labels.
Truncate sequences to be no longer than the maximum length set by the max_length parameter. The max_length of the text kept at 220 words and the max_length of the summary is kept at 50 words.
The ‘map’ function of the Huggingface dataset can be used to apply the pre_process function across the entire data.

from transformers import AutoTokenizer

checkpoint = "t5-small"
#checkpoint = "facebook/bart-large-cnn"
tokenizer = AutoTokenizer.from_pretrained(checkpoint)

prefix = "summarize: "

def preprocess_function(passages):
    inputs = [prefix + doc for doc in passages["text"]]
    model_inputs = tokenizer(inputs, max_length=220, truncation=True)

    labels = tokenizer(text_target=passages["summary"], max_length=50, truncation=True)

    model_inputs["labels"] = labels["input_ids"]
    return model_inputs

tokenized_df1 = df1.map(preprocess_function, batched=True)

DataCollatorForSeq2Seq can be used to dynamically pad the sentences to the longest length in a batch during collation, instead of padding the whole dataset to the maximum length.

from transformers import DataCollatorForSeq2Seq
data_collator = DataCollatorForSeq2Seq(tokenizer=tokenizer, model=checkpoint)

e) Evaluate performance of Model

The rouge1,rouge2 metric can be used to evaluate the performance of the model

import evaluate
rouge = evaluate.load("rouge")

f)Create a function compute_metrics that passes your predictions and labels to ‘compute’ to calculate the ROUGE metric:

import numpy as np

def compute_metrics(eval_pred):
    # evaluate predictions and labels
    predictions, labels = eval_pred
    decoded_preds = tokenizer.batch_decode(predictions, skip_special_tokens=True)
    labels = np.where(labels != -100, labels, tokenizer.pad_token_id)
    decoded_labels = tokenizer.batch_decode(labels, skip_special_tokens=True)
     # compute rouge score between the labels and predictions
    result = rouge.compute(predictions=decoded_preds, references=decoded_labels, use_stemmer=True)

    prediction_lens = [np.count_nonzero(pred != tokenizer.pad_token_id) for pred in predictions]
    result["gen_len"] = np.mean(prediction_lens)

    return {k: round(v, 4) for k, v in result.items()}

g) Split the data into training(80%) and test(20%) data set

train_dataset = tokenized_df1.shuffle(seed=42).select(range(100))
test_dataset = tokenized_df1.shuffle(seed=42).select(range(30))

len(train_dataset)

h) Train the model with AutoModelForSeq2SeqLM

from transformers import AutoModelForSeq2SeqLM, Seq2SeqTrainingArguments, Seq2SeqTrainer
model = AutoModelForSeq2SeqLM.from_pretrained(checkpoint)

Set training hyperparameters in Seq2SeqTrainingArguments. The Adam optimization, with learning rate, beta1 & beta2 are used
Pass the training arguments to Seq2SeqTrainer along with the model, dataset, tokenizer, data collator, and compute_metrics function.
Call train() to finetune your model.

training_args = Seq2SeqTrainingArguments(
    output_dir="philosophy_model",
    evaluation_strategy="epoch",
    learning_rate= 5.6e-03,
    adam_beta1=0.9,
    adam_beta2=0.99,
    adam_epsilon=1e-06,
    per_device_train_batch_size=8,
    per_device_eval_batch_size=8,
    weight_decay=0.01,
    save_total_limit=3,
    num_train_epochs=20,
    predict_with_generate=True,
    fp16=True,
    push_to_hub=True,
)

trainer = Seq2SeqTrainer(
    model=model,
    args=training_args,
    train_dataset=train_dataset,
    eval_dataset=test_dataset,
    tokenizer=tokenizer,
    data_collator=data_collator,
    compute_metrics=compute_metrics,
)

trainer.train()

Epoch	Training Loss	Validation Loss	Rouge1	Rouge2	Rougel	Rougelsum	Gen Len
1	No log	2.246223	0.363200	0.146200	0.311400	0.312600	18.333300
2	No log	1.461140	0.459000	0.303900	0.417800	0.417800	18.566700
3	No log	0.832312	0.546500	0.425900	0.524700	0.520800	17.133300
4	No log	0.472341	0.616100	0.517600	0.601000	0.600400	18.366700
5	No log	0.312106	0.681200	0.607800	0.674700	0.671400	18.233300
6	No log	0.154585	0.741800	0.702300	0.733800	0.731300	18.066700
7	No log	0.112100	0.783200	0.763000	0.780200	0.778900	18.500000
8	No log	0.069882	0.801400	0.788200	0.802700	0.800900	18.533300
9	No log	0.045941	0.795800	0.780500	0.794600	0.791700	18.500000
10	No log	0.051655	0.809100	0.795800	0.810500	0.809000	18.466700
11	No log	0.035792	0.799400	0.785200	0.797300	0.794600	18.500000
12	No log	0.041766	0.779900	0.754800	0.774700	0.773200	18.266700
13	No log	0.010703	0.810000	0.800400	0.810700	0.809000	18.500000
14	No log	0.006519	0.807700	0.797100	0.809400	0.807500	18.500000
15	No log	0.017779	0.808000	0.796000	0.809400	0.807500	18.366700
16	No log	0.001681	0.810000	0.800400	0.810700	0.809000	18.500000
17	No log	0.005469	0.810000	0.800400	0.810700	0.809000	18.500000
18	No log	0.002003	0.810000	0.800400	0.810700	0.809000	18.500000
19	No log	0.000638	0.810000	0.800400	0.810700	0.809000	18.500000
20	No log	0.000498	0.810000	0.800400	0.810700	0.809000	18.500000
TrainOutput(global_step=260, training_loss=0.6491916949932391, metrics={'train_runtime': 57.99, 'train_samples_per_second': 34.489, 'train_steps_per_second': 4.484, 'total_flos': 101132046434304.0, 'train_loss': 0.6491916949932391, 'epoch': 20.0})

As we can see the rouge1 to rouge2 scores are fairly good. Anything above 0.5 is considered good. Maybe this is because the T5 model has already been pre-trained on a fairly large philosophical dataset

j) Push to hub

trainer.push_to_hub()

k) Summarise using pipeline

text = "summarize: A seeker who has the necessary qualifications, in order that he may be redeemed from his inner weaknesses, attachments, animalisms and false values is advised to serve with devotion a Teacher who is well- established in the experience of the Self."

from transformers import pipeline

summarizer = pipeline("summarization", model="tvganesh/philosophy_model")
summarizer(text)

[{'summary_text': 'A seeker who has the necessary qualifications will be able to free oneself of sense objects, and one cannot expect this to happen without any mental tossing'}]

l) Summarise using model generate

from transformers import AutoTokenizer

tokenizer = AutoTokenizer.from_pretrained("tvganesh/philosophy_model")
inputs = tokenizer(text, return_tensors="pt").input_ids

from transformers import AutoModelForSeq2SeqLM

model = AutoModelForSeq2SeqLM.from_pretrained("tvganesh/philosophy_model")
outputs = model.generate(inputs, max_new_tokens=70, do_sample=False)

tokenizer.decode(outputs[0], skip_special_tokens=True)

'A seeker who has the necessary qualifications will help in his journey to redeem himself'

l) Number of beams

summary_ids = model.generate(inputs,
                                    num_beams=10,
                                    no_repeat_ngram_size=3,
                                    min_length=20,
                                    max_length=70,
                                    early_stopping=True)
output = tokenizer.decode(summary_ids[0], skip_special_tokens=True)
output

'A seeker who has the necessary qualifications will be able to free himself of sense objects and false values'

I also tried Facebook’s BART Large model but the performance was not good at all.

You can try out the model at the following link philosophy_model

Anyway this was a good learning experience.

References

Also see

To see all posts click Index of posts

T20 Win Probability using CTGANs, synthetic data

This should be my last post on computing T20 Win Probability. In this post I compute Win Probability using Augmented Data with the help of Conditional Tabular Generative Adversarial Networks (CTGANs).

A.Introduction

I started the computation of T20 match Win Probability in my earlier post

a) ‘Computing Win-Probability of T20 matches‘ where I used

vanilla Logistic Regression to get an accuracy of 0.67,
Random Forest with Tidy models gave me an accuracy of 0.737
Deep Learning with Keras also with 0.73.

This was done without player embeddings

b) Next I used player embeddings for batsmen and bowlers in my post Boosting Win Probability accuracy with player embeddings , and my accuracies improved significantly

glmnet : accuracy – 0.728 and roc_auc – 0.81
random forest : accuracy – 0.927 and roc_auc – 0.98
mlp-dnn :accuracy – 0.762 and roc_auc – 0.854

c) Third I tried using Deep Learning with Keras using player embeddings

DL network gave an accuracy of 0.8639

This was lightweight and could be easily deployed in my Shiny GooglyPlusPlus app as opposed to the Tidymodel’s Random Forest, which was bulky and slow.

d) Finally I decided to try and improve the accuracy of my Deep Learning Model using Synthetic data. Towards this end, my explorations led me to Conditional Tabular Generative Adversarial Networks (CTGANs). CTGAN are GAN networks that can be used with Tabular data as GAN models are not useful with tabular data. However, the best performance I got for

DL Keras Model + Synthetic data : accuracy =0.77

The poorer accuracy was because CTGAN requires enormous computing power (GPUs) and RAM. The free version of Colab, Kaggle kept crashing when I tried with even 0.1 % of my 1.2 million dataset size. Finally, I tried with just 0.05% and was able to generate synthetic data. Most likely, it is the small sample size and the smaller number of epochs could be the reason for the poor result. In any case, it was worth trying and this approach would possibly work with sufficient computing resources.

B.Generative Adversarial Networks (GANs)

Generative Adversarial Networks (GANs) was the brain child of Ian Goodfellow who demonstrated it in 2014. GANs are capable of generating synthetic text, tables, images, videos using available data. In Adversarial nets framework, the generative model is pitted against an adversary: a
discriminative model that learns to determine whether a sample is from the model distribution or the
data distribution.

GANs have 2 Deep Neural Networks , the Generator and Discriminator which compete against other

The Generator (Counterfeiter) takes random noise as input and generates fake images, tables, text. The generator learns to generate plausible data. The generated instances become negative training examples for the discriminator.
The Discriminator (Police) which tries to distinguish between the real and fake images, text. The discriminator learns to distinguish the generator’s fake data from real data. The discriminator penalises the generator for producing implausible results.

A pictorial representation of the GAN model can be shown below

Theoretically best performance of GANs are supposed to happen when the network reaches the ‘Nash equilibrium‘, i.e. when the Generator produces near fake images and the Discriminator’s loss is f ~0.5 i.e. the discriminator is unable to distinguish between real and fake images.

Note: Though I have mentioned T20 data in the above GAN model, the T20 tabular data is actually used in CTGAN which is slightly different from the above. See Reference 2) below.

C. Conditional Tabular Generative Adversial Networks (CTGANs)

“Modeling the probability distribution of rows in tabular data and generating realistic synthetic data is a non-trivial task. Tabular data usually contains a mix of discrete and continuous columns. Continuous columns may have multiple modes whereas discrete columns are sometimes imbalanced making the modeling difficult.” CTGANs handle these challenges.

I came upon CTGAN after spending some time exploring GANs via blogs, videos etc. For building the model I use real T20 match data. However, CTGAN requires immense raw computing power and a lot of RAM. My initial attempts on Colab, my Mac (12 core, 32GB RAM), took forever before eventually crashing, I switched to Kaggle and used GPUs. Still I was only able to use only a miniscule part of my T20 dataset. My match data has 1.2 million rows, hoanything > 0.05% resulted in Kaggle crashing. Since I was able to use only a fraction, I executed the CTGAN model over several iterations, each iteration with a random 0.05% sample of the dataset. At the end of each iterations I also generate synthetic dataset. Over 12 iterations, I generate close 360K of ‘synthetic‘ T20 match data.

I then augment the 1.2 million rows of ‘real‘ T20 match data with the generated ‘synthetic T20 match data to run my Deep Learning model

D. Executing the CTGAN model

a. Read the real T20 match data

!pip install ctgan
import pandas as pd
import ctgan
from ctgan import CTGAN
from numpy.random import seed

# Read the T20 match data
df = pd.read_csv('/kaggle/input/cricket1/t20.csv')

# Randomly sample 0.05% of the dataset. Note larger datasets cause the algo to crash
train_dataset = df.sample(frac=0.05)

# Print the real T20 match data
print(train_dataset.head(10))
print(train_dataset.shape)

             batsmanIdx  bowlerIdx  ballNum  ballsRemaining  runs   runRate  \
363695         3333        432      134             119   153  1.285714   
1082839        3881       1180      218              30    93  3.100000   
595799         2366        683      187              65   120  1.846154   
737614         4490       1381      148              87   144  1.655172   
410202          934       1003       19             106    35  1.842105   
525627          921       1711      251               1     8  8.000000   
657669         4718       1602      130             115   145  1.260870   
666461         4309       1989       44              87    38  0.863636   
651229         3336        754       30              92    36  1.200000   
709892         3048        421       97              28   119  1.226804   

            numWickets  runsMomentum  perfIndex  isWinner  
363695            0      0.092437  18.333333         1  
1082839           5      0.200000   4.736842         0  
595799            4      0.107692   9.566667         0  
737614            1      0.114943   9.130435         1  
410202            0      0.103774  20.263158         0  
525627            8      3.000000   3.837209         0  
657669            0      0.095652  19.555556         0  
666461            0      0.126437   9.500000         0  
651229            0      0.119565  13.200000         0  
709892            3      0.285714   9.814433         1  
(59956, 10)

b. Run CTGAN model on the real T20 data

import pandas as pd
import ctgan
from ctgan import CTGAN
from numpy.random import seed
from pickle import TRUE

df = pd.read_csv('/kaggle/input/cricket1/t20.csv')

#Specify the categorical features. batsmanIdx & bowlerIdx are player embeddings
categorical_features = ['batsmanIdx','bowlerIdx']

# Create a empty dataframe for synthetic data
df1 = pd.DataFrame()

# Loop for 12 iterations. Minimize generator & discriminator loss
for i in range(12):
    print(i)
    train_dataset = df.sample(frac=0.05)
    seed(33)

    ctgan = CTGAN(epochs=20,verbose=True,generator_lr=.001,discriminator_lr=.001,batch_size=1000)
    ctgan.fit(train_dataset, categorical_features)

    # Generate synthetic data
    samples = ctgan.sample(30000)

   # Concatenate the synthetic data after each iteration
    df1 = pd.concat([df1,samples])
    print(samples.head())
    print(df1.shape)

# Output the synthetic data to file
df1.to_csv("output1.csv",index=False)

0
Epoch 1, Loss G:  8.3825,Loss D: -0.6159
Epoch 2, Loss G:  3.5117,Loss D: -0.3016
Epoch 3, Loss G:  2.1619,Loss D: -0.5713
Epoch 4, Loss G:  0.9847,Loss D:  0.1010
Epoch 5, Loss G:  0.6198,Loss D:  0.0789
Epoch 6, Loss G:  0.1710,Loss D:  0.0959
Epoch 7, Loss G:  0.3236,Loss D: -0.1554
Epoch 8, Loss G:  0.2317,Loss D: -0.0765
Epoch 9, Loss G: -0.0127,Loss D:  0.0275
Epoch 10, Loss G:  0.1477,Loss D: -0.0353
Epoch 11, Loss G:  0.0997,Loss D: -0.0129
Epoch 12, Loss G:  0.0066,Loss D: -0.0486
Epoch 13, Loss G:  0.0351,Loss D: -0.0805
Epoch 14, Loss G: -0.1399,Loss D: -0.0021
Epoch 15, Loss G: -0.1503,Loss D: -0.0518
Epoch 16, Loss G: -0.2306,Loss D: -0.0234
Epoch 17, Loss G: -0.2986,Loss D:  0.0469
Epoch 18, Loss G: -0.1941,Loss D: -0.0560
Epoch 19, Loss G: -0.3794,Loss D:  0.0000
Epoch 20, Loss G: -0.2763,Loss D:  0.0368
   batsmanIdx  bowlerIdx  ballNum  ballsRemaining  runs   runRate  numWickets  \
0         906        224        8              75    81  1.955153           4   
1        4159        433       17              31   126  1.799280           9   
2         229        351      192              66    82  1.608527           5   
3        1926        962       63               0   117  1.658105           0   
4         286        431      128               1    36  1.605079           0   

   runsMomentum  perfIndex  isWinner  
0      0.146670   6.937595         1  
1      0.160534  10.904346         1  
2      0.516010  11.698128         1  
3      0.380986  11.914613         0  
4      0.112255   5.392120         0  
(30000, 10)
1
Epoch 1, Loss G:  7.9977,Loss D: -0.3592
Epoch 2, Loss G:  3.7418,Loss D: -0.3371
Epoch 3, Loss G:  1.6685,Loss D: -0.3211
Epoch 4, Loss G:  1.0539,Loss D: -0.3495
Epoch 5, Loss G:  0.4664,Loss D: -0.0907
Epoch 6, Loss G:  0.4004,Loss D: -0.1208
Epoch 7, Loss G:  0.3250,Loss D: -0.1482
Epoch 8, Loss G:  0.1753,Loss D:  0.0169
Epoch 9, Loss G:  0.1382,Loss D:  0.0661
Epoch 10, Loss G:  0.1509,Loss D: -0.1023
Epoch 11, Loss G: -0.0235,Loss D:  0.0210
Epoch 12, Loss G: -0.1636,Loss D: -0.0124
Epoch 13, Loss G: -0.3370,Loss D: -0.0185
Epoch 14, Loss G: -0.3054,Loss D: -0.0085
Epoch 15, Loss G: -0.5142,Loss D:  0.0121
Epoch 16, Loss G: -0.3813,Loss D: -0.0921
Epoch 17, Loss G: -0.5838,Loss D:  0.0210
Epoch 18, Loss G: -0.4033,Loss D: -0.0181
Epoch 19, Loss G: -0.5711,Loss D:  0.0269
Epoch 20, Loss G: -0.4828,Loss D: -0.0830
   batsmanIdx  bowlerIdx  ballNum  ballsRemaining  runs   runRate  numWickets  \
0        2202        265      223              39    13  0.868927           0   
1        3641        856       35              59    26  2.236160           6   
2         676       2903      218              93    16  0.460693           1   
3        3482       3459       44             117   102  0.851471           8   
4        3046       3076       59               5    84  1.016824           2   

   runsMomentum  perfIndex  isWinner  
0      0.138586   4.733462         0  
1      0.124453   5.146831         1  
2      0.273168  10.106869         0  
3      0.129520   5.361127         0  
4      1.083525  25.677574         1  
(60000, 10)
...
...
11
Epoch 1, Loss G:  8.8362,Loss D: -0.7111
Epoch 2, Loss G:  4.1322,Loss D: -0.8468
Epoch 3, Loss G:  1.2782,Loss D:  0.1245
Epoch 4, Loss G:  1.1135,Loss D: -0.3588
Epoch 5, Loss G:  0.6033,Loss D: -0.1255
Epoch 6, Loss G:  0.6912,Loss D: -0.1906
Epoch 7, Loss G:  0.3340,Loss D: -0.1048
Epoch 8, Loss G:  0.3515,Loss D: -0.0730
Epoch 9, Loss G:  0.1702,Loss D:  0.0237
Epoch 10, Loss G:  0.1064,Loss D:  0.0632
Epoch 11, Loss G:  0.0884,Loss D: -0.0005
Epoch 12, Loss G:  0.0556,Loss D: -0.0607
Epoch 13, Loss G: -0.0917,Loss D: -0.0223
Epoch 14, Loss G: -0.1492,Loss D:  0.0258
Epoch 15, Loss G: -0.0986,Loss D: -0.0112
Epoch 16, Loss G: -0.1428,Loss D: -0.0060
Epoch 17, Loss G: -0.2225,Loss D: -0.0263
Epoch 18, Loss G: -0.2255,Loss D: -0.0328
Epoch 19, Loss G: -0.3482,Loss D:  0.0277
Epoch 20, Loss G: -0.2667,Loss D: -0.0721
   batsmanIdx  bowlerIdx  ballNum  ballsRemaining  runs   runRate  numWickets  \
0         367       1447      129              27    30  1.242120           2   
1        2481       1528      221               4    10  1.344024           2   
2        1034       3116      132              87   153  1.142750           3   
3        1201       2868      151              60   136  1.091638           1   
4        4327       3291      108              89    22  0.842775           2   

   runsMomentum  perfIndex  isWinner  
0      1.978739   6.393691         1  
1      0.539650   6.783990         0  
2      0.107156  12.154197         0  
3      3.193574  11.992059         0  
4      0.127507  12.210876         0  
(360000, 10)

E. Sample of the Synthetic data

synthetic_data = ctgan.sample(20000)
print(synthetic_data.head(100))

    batsmanIdx  bowlerIdx  ballNum  ballsRemaining  runs    runRate  \
0         1073       3059       72              72   149   2.230236   
1         3769       1443      106               7   137   0.881409   
2          448       3048      166               6   220   1.092504   
3         2969       1244      103              82   207  12.314862   
4          180       1372      125             111    14   1.310051   
..         ...        ...      ...             ...   ...        ...   
95        1521       1040      153               6   166   1.097363   
96        2366         62       25             114   119   0.910642   
97        3506       1736      100             118   140   1.640921   
98        3343       2347       47              54    50   0.696462   
99        1957       2888      136              27   153   1.315565   

    numWickets  runsMomentum  perfIndex  isWinner  
0            0      0.111707  17.466925         0  
1            1      0.130352  14.274113         0  
2            1      0.173541  11.076731         1  
3            1      0.218977   6.239951         0  
4            4      2.829380   9.183323         1  
..         ...           ...        ...       ...  
95           0      0.223437   7.011180         0  
96           1      0.451371  16.908120         1  
97           5      0.156936   9.217205         0  
98           6      0.124536   6.273091         0  
99           1      0.249329  14.221554         0  

[100 rows x 10 columns]

F. Evaluating the synthetic T20 match data

Here the quality of the synthetic data set is evaluated.

a) Statistical evaluation

Read the real T20 match data
Read the generated T20 synthetic match data

import pandas as pd

# Read the T20 match and synthetic match data
df = pd.read_csv('/kaggle/input/cricket1/t20.csv').  #1.2 million rows
synthetic=pd.read_csv('/kaggle/input/synthetic/synthetic.csv')   #300K

# Randomly sample 1000 rows, and generate stats
df1=df.sample(n=1000)
real=df1.describe()
realData_stats=real.transpose
print(realData_stats)

synthetic1=synthetic.sample(n=1000)
synthetic=synthetic1.describe()
syntheticData_stats=synthetic.transpose
syntheticData_stats

a) Stats of real T20 match data

<bound method DataFrame.transpose of         batsmanIdx    bowlerIdx      ballNum  ballsRemaining         runs  \
count  1000.000000  1000.000000  1000.000000     1000.000000  1000.000000   
mean   2323.940000  1776.481000   118.165000       59.236000    77.649000   
std    1329.703046  1011.470703    70.564291       35.312934    49.098763   
min       8.000000    13.000000     1.000000        1.000000    -2.000000   
25%    1134.750000   850.000000    58.000000       28.750000    39.000000   
50%    2265.000000  1781.500000   117.000000       59.000000    72.000000   
75%    3510.000000  2662.250000   178.000000       89.000000   111.000000   
max    4738.000000  3481.000000   265.000000      127.000000   246.000000   

           runRate   numWickets  runsMomentum    perfIndex     isWinner  
count  1000.000000  1000.000000   1000.000000  1000.000000  1000.000000  
mean      1.734979     2.614000      0.310568     9.580386     0.499000  
std       5.698104     2.267189      0.686171     4.530856     0.500249  
min      -2.000000     0.000000      0.071429     0.000000     0.000000  
25%       1.009063     1.000000      0.105769     6.666667     0.000000  
50%       1.272727     2.000000      0.141026     9.236842     0.000000  
75%       1.546891     4.000000      0.250000    12.146735     1.000000  
max     166.000000    10.000000     10.000000    30.800000     1.000000

b) Stats of Synthetic T20 match data

     
           batsmanIdx    bowlerIdx      ballNum  ballsRemaining         runs  \
count  1000.000000  1000.000000  1000.000000     1000.000000  1000.000000   
mean   2304.135000  1760.776000   116.081000       50.102000    74.357000   
std    1342.348684  1003.496003    72.019228       35.795236    48.103446   
min       2.000000    15.000000    -4.000000       -2.000000    -1.000000   
25%    1093.000000   881.000000    46.000000       18.000000    30.000000   
50%    2219.500000  1763.500000   116.000000       45.000000    75.000000   
75%    3496.500000  2644.750000   180.250000       77.000000   112.000000   
max    4718.000000  3481.000000   253.000000      124.000000   222.000000   

           runRate   numWickets  runsMomentum    perfIndex     isWinner  
count  1000.000000  1000.000000   1000.000000  1000.000000  1000.000000  
mean      1.637225     3.096000      0.336540     9.278073     0.507000  
std       1.691060     2.640408      0.502346     4.727677     0.500201  
min      -4.388339     0.000000      0.083351    -0.902991     0.000000  
25%       1.077789     1.000000      0.115770     5.731931     0.000000  
50%       1.369655     2.000000      0.163085     9.104328     1.000000  
75%       1.660477     5.000000      0.311586    12.619318     1.000000  
max      23.757001    10.000000      4.630908    29.829497     1.000000

c) Plotting the Generator and Discriminator loss

import pandas as pd

# CTGAN prints out a new line for each epoch
epochs_output = str(output).split('\n')

# CTGAN separates the values with commas
raw_values = [line.split(',') for line in epochs_output]
loss_values = pd.DataFrame(raw_values)[:-1] # convert to df and delete last row (empty)

# Rename columns
loss_values.columns = ['Epoch', 'Generator Loss', 'Discriminator Loss']

# Extract the numbers from each column 
loss_values['Epoch'] = loss_values['Epoch'].str.extract('(\d+)').astype(int)
loss_values['Generator Loss'] = loss_values['Generator Loss'].str.extract('([-+]?\d*\.\d+|\d+)').astype(float)
loss_values['Discriminator Loss'] = loss_values['Discriminator Loss'].str.extract('([-+]?\d*\.\d+|\d+)').astype(float)

# the result is a row for each epoch that contains the generator and discriminator loss
loss_values.head()

	Epoch	Generator Loss	Discriminator Loss
0	1	8.0158	-0.3840
1	2	4.6748	-0.9589
2	3	1.1503	-0.0066
3	4	1.5593	-0.8148
4	5	0.6734	-0.1425
5	6	0.5342	-0.2202
6	7	0.4539	-0.1462
7	8	0.2907	-0.0155
8	9	0.2399	0.0172
9	10	0.1520	-0.0236

import plotly.graph_objects as go

# Plot loss function
fig = go.Figure(data=[go.Scatter(x=loss_values['Epoch'], y=loss_values['Generator Loss'], name='Generator Loss'),
                      go.Scatter(x=loss_values['Epoch'], y=loss_values['Discriminator Loss'], name='Discriminator Loss')])


# Update the layout for best viewing
fig.update_layout(template='plotly_white',
                    legend_orientation="h",
                    legend=dict(x=0, y=1.1))

title = 'CTGAN loss function for T20 dataset - ' 
fig.update_layout(title=title, xaxis_title='Epoch', yaxis_title='Loss')
fig.show()

G. Qualitative evaluation of Synthetic data

a) Quality of continuous columns in synthetic data

KSComplement -This metric computes the similarity of a real column vs. a synthetic column in terms of the column shapes.The KSComplement uses the Kolmogorov-Smirnov statistic. Closer to 1.0 is good and 0 is worst

from sdmetrics.single_column import KSComplement
numerical_columns=['ballNum','ballsRemaining','runs','runRate','numWickets','runsMomentum','perfIndex']
total_score = 0
for column_name in numerical_columns:
    column_score = KSComplement.compute(df[column_name], synthetic[column_name])
    total_score += column_score
    print('Column:', column_name, ', Score: ', column_score)

print('\nAverage: ', total_score/len(numerical_columns))

Column: ballNum , Score:  0.9502754283367316
Column: ballsRemaining , Score:  0.8770284103276166
Column: runs , Score:  0.9136464248633367
Column: runRate , Score:  0.9183841670732166
Column: numWickets , Score:  0.9016209114638712
Column: runsMomentum , Score:  0.8773491702213716
Column: perfIndex , Score:  0.9173808852778924

Average:  0.9079550567948624

b) Quality of categorical columns

This statistic measures the quality of generated categorical columns. 1 is best and 0 is worst

categorical_columns=['batsmanIdx','bowlerIdx']
from sdmetrics.single_column import TVComplement

total_score = 0
for column_name in categorical_columns:
    column_score = TVComplement.compute(df[column_name], synthetic[column_name])
    total_score += column_score
    print('Column:', column_name, ', Score: ', column_score)

print('\nAverage: ', total_score/len(categorical_columns))

Column: batsmanIdx , Score:  0.8436263499539245
Column: bowlerIdx , Score:  0.7356177407921669

Average:  0.7896220453730457

The performance is decent but not excellent. I was unable to execute more epochs as it it required larger than the memory allowed

c) Correlation similarity

This metric measures the correlation between a pair of numerical columns and computes the similarity between the real and synthetic data – it compares the trends of 2D distributions. Best 1.0 and 0.0 is worst

import itertools
from sdmetrics.column_pairs import CorrelationSimilarity

total_score = 0
total_pairs = 0
for pair in itertools.combinations(numerical_columns,2):
    col_A, col_B = pair
    score = CorrelationSimilarity.compute(df[[col_A, col_B]], synthetic[[col_A, col_B]])
    print('Columns:', pair, ' Score:', score)
    total_score += score
    total_pairs += 1

print('\nAverage: ', total_score/total_pairs)

Columns: ('ballNum', 'ballsRemaining')  Score: 0.7153942317384889
Columns: ('ballNum', 'runs')  Score: 0.8838043045134777
Columns: ('ballNum', 'runRate')  Score: 0.8710243133637056
Columns: ('ballNum', 'numWickets')  Score: 0.7978515509750435
Columns: ('ballNum', 'runsMomentum')  Score: 0.8956281260834316
Columns: ('ballNum', 'perfIndex')  Score: 0.9275145840528048
Columns: ('ballsRemaining', 'runs')  Score: 0.9566928975064546
Columns: ('ballsRemaining', 'runRate')  Score: 0.9127313819127167
Columns: ('ballsRemaining', 'numWickets')  Score: 0.6770737279315224
Columns: ('ballsRemaining', 'runsMomentum')  Score: 0.7939260278412358
Columns: ('ballsRemaining', 'perfIndex')  Score: 0.8694582252638351
Columns: ('runs', 'runRate')  Score: 0.999593795992159
Columns: ('runs', 'numWickets')  Score: 0.9510731832916608
Columns: ('runs', 'runsMomentum')  Score: 0.9956131422133428
Columns: ('runs', 'perfIndex')  Score: 0.9742931845536701
Columns: ('runRate', 'numWickets')  Score: 0.8859830711832263
Columns: ('runRate', 'runsMomentum')  Score: 0.9174744874779561
Columns: ('runRate', 'perfIndex')  Score: 0.9491100087911353
Columns: ('numWickets', 'runsMomentum')  Score: 0.8989709776329797
Columns: ('numWickets', 'perfIndex')  Score: 0.7178946968801441
Columns: ('runsMomentum', 'perfIndex')  Score: 0.9744441623018661

Average:  0.8840738134048025

d) Category coverage

This metric measures whether a synthetic column covers all the possible categories that are present in a real column. 1.0 is best , 0 is worst

from sdmetrics.single_column import CategoryCoverage

total_score = 0
for column_name in categorical_columns:
    column_score = CategoryCoverage.compute(df[column_name], synthetic[column_name])
    total_score += column_score
    print('Column:', column_name, ', Score: ', column_score)

print('\nAverage: ', total_score/len(categorical_columns))

Column: batsmanIdx , Score:  0.9533951919021509
Column: bowlerIdx , Score:  0.9913966160022942

Average:  0.9723959039522225

H. Augmenting the T20 match data set

In this final part I augment my T20 match data set with the generated synthetic T20 data set.

import pandas as pd
from numpy import savetxt
import tensorflow as tf
from tensorflow import keras
import pandas as pd
import numpy as np


from keras.layers import Input, Embedding, Flatten, Dense, Reshape, Concatenate, Dropout
from keras.models import Model
import matplotlib.pyplot as plt

# Read real and synthetic data
df = pd.read_csv('/kaggle/input/cricket1/t20.csv')
synthetic=pd.read_csv('/kaggle/input/synthetic/synthetic.csv')

# Augment the data. Concatenate real & synthetic data
df1=pd.concat([df,synthetic])

# Create training and test samples
print("Shape of dataframe=",df1.shape)
train_dataset = df1.sample(frac=0.8,random_state=0)
test_dataset = df1.drop(train_dataset.index)
train_dataset1 = train_dataset[['batsmanIdx','bowlerIdx','ballNum','ballsRemaining','runs','runRate','numWickets','runsMomentum','perfIndex']]
test_dataset1 = test_dataset[['batsmanIdx','bowlerIdx','ballNum','ballsRemaining','runs','runRate','numWickets','runsMomentum','perfIndex']]
train_dataset1
train_labels = train_dataset.pop('isWinner')
test_labels = test_dataset.pop('isWinner')
print(train_dataset1.shape)

a=train_dataset1.describe()
stats=a.transpose
print(a)

a) Create A Deep Learning Model in Keras

from numpy.random import seed
seed(33)
tf.random.set_seed(432)
# create input layers for each of the predictors
batsmanIdx_input = Input(shape=(1,), name='batsmanIdx')
bowlerIdx_input = Input(shape=(1,), name='bowlerIdx')
ballNum_input = Input(shape=(1,), name='ballNum')
ballsRemaining_input = Input(shape=(1,), name='ballsRemaining')
runs_input = Input(shape=(1,), name='runs')
runRate_input = Input(shape=(1,), name='runRate')
numWickets_input = Input(shape=(1,), name='numWickets')
runsMomentum_input = Input(shape=(1,), name='runsMomentum')
perfIndex_input = Input(shape=(1,), name='perfIndex')

no_of_unique_batman=len(df1["batsmanIdx"].unique()) 
print(no_of_unique_batman)
no_of_unique_bowler=len(df1["bowlerIdx"].unique()) 
print(no_of_unique_bowler)

embedding_size_bat = no_of_unique_batman ** (1/4)
print(embedding_size_bat)
embedding_size_bwl = no_of_unique_bowler ** (1/4)
print(embedding_size_bwl)

# create embedding layer for the categorical predictor
batsmanIdx_embedding = Embedding(input_dim=no_of_unique_batman+1, output_dim=16,input_length=1)(batsmanIdx_input)
print(batsmanIdx_embedding)
batsmanIdx_flatten = Flatten()(batsmanIdx_embedding)
print(batsmanIdx_flatten)
bowlerIdx_embedding = Embedding(input_dim=no_of_unique_bowler+1, output_dim=16,input_length=1)(bowlerIdx_input)
bowlerIdx_flatten = Flatten()(bowlerIdx_embedding)
print(bowlerIdx_flatten)
# concatenate all the predictors
x = keras.layers.concatenate([batsmanIdx_flatten,bowlerIdx_flatten, ballNum_input, ballsRemaining_input, runs_input, runRate_input, numWickets_input, runsMomentum_input, perfIndex_input])
print(x.shape)
# add hidden layers
x = Dense(96, activation='relu')(x)
x = Dropout(0.1)(x)
x = Dense(64, activation='relu')(x)
x = Dropout(0.1)(x)
x = Dense(32, activation='relu')(x)
x = Dropout(0.1)(x)
x = Dense(16, activation='relu')(x)
x = Dropout(0.1)(x)
x = Dense(8, activation='relu')(x)
x = Dropout(0.1)(x)
# add output layer
output = Dense(1, activation='sigmoid', name='output')(x)
print(output.shape)
# create model
model = Model(inputs=[batsmanIdx_input,bowlerIdx_input, ballNum_input, ballsRemaining_input, runs_input, runRate_input, numWickets_input, runsMomentum_input, perfIndex_input], outputs=output)
model.summary()
# compile model
#optimizer=keras.optimizers.Adam(learning_rate=.01, beta_1=0.1, beta_2=0.999, epsilon=None, decay=0.0, amsgrad=True)
#optimizer=keras.optimizers.RMSprop(learning_rate=0.001, rho=0.2, momentum=0.2, epsilon=1e-07)
#optimizer=keras.optimizers.SGD(learning_rate=.01,momentum=0.1) #- Works without dropout
#optimizer = tf.keras.optimizers.RMSprop(0.01)
#optimizer=keras.optimizers.SGD(learning_rate=.01,momentum=0.1)
#optimizer=keras.optimizers.RMSprop(learning_rate=.005, rho=0.1, momentum=0, epsilon=1e-07)

optimizer=keras.optimizers.Adam(learning_rate=.015, beta_1=0.9, beta_2=0.999, epsilon=1e-07, amsgrad=True)

model.compile(optimizer=optimizer, loss='binary_crossentropy', metrics=['accuracy'])

# train the model
history=model.fit([train_dataset1['batsmanIdx'],train_dataset1['bowlerIdx'],train_dataset1['ballNum'],train_dataset1['ballsRemaining'],train_dataset1['runs'],
           train_dataset1['runRate'],train_dataset1['numWickets'],train_dataset1['runsMomentum'],train_dataset1['perfIndex']], train_labels, epochs=20, batch_size=1024,
          validation_data = ([test_dataset1['batsmanIdx'],test_dataset1['bowlerIdx'],test_dataset1['ballNum'],test_dataset1['ballsRemaining'],test_dataset1['runs'],
           test_dataset1['runRate'],test_dataset1['numWickets'],test_dataset1['runsMomentum'],test_dataset1['perfIndex']],test_labels), verbose=1)

plt.plot(history.history["loss"])
plt.plot(history.history["val_loss"])
plt.title("model loss")
plt.ylabel("loss")
plt.xlabel("epoch")
plt.legend(["train", "test"], loc="upper left")
plt.show()

==================================================================================================
Total params: 144,497
Trainable params: 144,497
Non-trainable params: 0
__________________________________________________________________________________________________
Epoch 1/20
1219/1219 [==============================] - 15s 11ms/step - loss: 0.6285 - accuracy: 0.6372 - val_loss: 0.5164 - val_accuracy: 0.7606
Epoch 2/20
1219/1219 [==============================] - 14s 11ms/step - loss: 0.5594 - accuracy: 0.7121 - val_loss: 0.4920 - val_accuracy: 0.7663
Epoch 3/20
1219/1219 [==============================] - 14s 12ms/step - loss: 0.5338 - accuracy: 0.7244 - val_loss: 0.4541 - val_accuracy: 0.7878
Epoch 4/20
1219/1219 [==============================] - 14s 11ms/step - loss: 0.5176 - accuracy: 0.7317 - val_loss: 0.4226 - val_accuracy: 0.7933
Epoch 5/20
1219/1219 [==============================] - 13s 11ms/step - loss: 0.4966 - accuracy: 0.7420 - val_loss: 0.4547 - val_accuracy: 0.7
...
...
poch 18/20
1219/1219 [==============================] - 14s 11ms/step - loss: 0.4300 - accuracy: 0.7747 - val_loss: 0.3536 - val_accuracy: 0.8288
Epoch 19/20
1219/1219 [==============================] - 14s 12ms/step - loss: 0.4269 - accuracy: 0.7766 - val_loss: 0.3565 - val_accuracy: 0.8302
Epoch 20/20
1219/1219 [==============================] - 14s 11ms/step - loss: 0.4259 - accuracy: 0.7775 - val_loss: 0.3498 - val_accuracy: 0.831

As can be seen the accuracy with augmented dataset is around 0.77, while without it I was getting 0.867 with just the real data. This degradation is probably due to the folllowing reasons

Only a fraction of the dataset was used for training. This was not representative of the data distribution for CTGAN to correctly synthesise data
The number of epochs had to be kept low to prevent Kaggle/Colab from crashing

I. Conclusion

This post shows how we can generate synthetic T20 match data to augment real T20 match data. Assuming we have sufficient processing power we should be able to generate synthetic data for augmenting our data set. This should improve the accuracy of the Win Probabily Deep Learning model.

References

Also see

To see all posts click Index of posts

Player Performance Estimation using AI Collaborative Filtering

1. Introduction

Often times before crucial matches, or in general, we would like to know the performance of a batsman against a bowler or vice-versa, but we may not have the data. We generally have data where different batsmen would have faced different sets of bowlers with certain performance data like ballsFaced, totalRuns, fours, sixes, strike rate and timesOut. Similarly different bowlers would have performance figures(deliveries, runsConceded, economyRate and wicketTaken) against different sets of batsmen. We will never have the data for all batsmen against all bowlers. However, it would be good estimate the performance of batsmen against a bowler, even though we do not have the performance data. This could be done using collaborative filtering which identifies and computes based on the similarity between batsmen vs bowlers & bowlers vs batsmen.

This post shows an approach whereby we can estimate a batsman’s performance against bowlers even though the batsman may not have faced those bowlers, based on his/her performance against other bowlers. It also estimates the performance of bowlers against batsmen using the same approach. This is based on the recommender algorithm which is used to recommend products to customers based on their rating on other products.

This idea came to me while generating the performance of batsmen vs bowlers & vice-versa for 2 IPL teams in this IPL 2022 with my Shiny app GooglyPlusPlus in the optimization tab, I found that there were some batsmen for which there was no data against certain bowlers, probably because they are playing for the first time in their team or because they were new (see picture below)

In the picture above there is no data for Dewald Brevis against Jasprit Bumrah and YS Chahal. Wouldn’t be great to estimate the performance of Brevis against Bumrah or vice-versa? Can we estimate this performance?

While pondering on this problem, I realized that this problem formulation is similar to the problem formulation for the famous Netflix movie recommendation problem, in which user’s ratings for certain movies are known and based on these ratings, the recommender engine can generate ratings for movies not yet seen.

This post estimates a player’s (batsman/bowler) using the recommender engine This post is based on R package recommenderlab

“Michael Hahsler (2021). recommenderlab: Lab for Developing and Testing Recommender Algorithms. R package version 0.2-7. https://github.com/mhahsler/recommenderlab”

Note 1: Thw data for this analysis is taken from Cricsheet after being processed by my R package yorkr.

You can also read this post in RPubs at Player Performance Estimation using AI Collaborative Filtering

A PDF copy of this post is available at Player Performance Estimation using AI Collaborative Filtering.pdf

You can download this R Markdown file and the associated data and perform the analysis yourself using any other recommender engine from Github at playerPerformanceEstimation

Problem statement

In the table below we see a set of bowlers vs a set of batsmen and the number of times the bowlers got these batsmen out.
By knowing the performance of the bowlers against some of the batsmen we can use collaborative filter to determine the missing values. This is done using the recommender engine.

The Recommender Engine works as follows. Let us say that there are feature vectors $x^1$ , $x^2$ and $x^3$ for the 3 bowlers which identify the characteristics of these bowlers (“fast”, “lateral drift through the air”, “movement off the pitch”). Let each batsman be identified by parameter vectors $\theta^1$ , $\theta^2$ and so on

For e.g. consider the following table

Then by assuming an initial estimate for the parameter vector $\theta$ and the feature vector xx we can formulate this as an optimization problem which tries to minimize the error for $\theta^T*x$ This can work very well as the algorithm can determine features which cannot be captured. So for e.g. some particular bowler may have very impressive figures. This could be due to some aspect of the bowling which cannot be captured by the data for e.g. let’s say the bowler uses the ‘scrambled seam’ when he is most effective, with a slightly different arc to the flight. Though the algorithm cannot identify the feature as we know it, but the ML algorithm should pick up intricacies which cannot be captured in data.

Hence the algorithm can be quite effective.

Note: The recommender lab performance is not very good and the Mean Square Error is quite high. Also, the ROC and AUC curves show that not in aLL cases the algorithm is doing a clean job of separating the True positives (TPR) from the False Positives (FPR)

Note: This is similar to the recommendation problem

The collaborative optimization object can be considered as a minimization of both $\theta$ and the features x and can be written as

J( $x^{(1)},x^{(2)},..x^{(n_{u})}$ , $\theta^{(1)},\theta^{(2)},..,\theta^{(n_{m})}$ }= 1/2 $\sum(\theta^{j})^{T}x^{i}- y^{(i,j)})^{2} + \lambda\sum\sum (x_{k}^{i})^{2} + \lambda\sum\sum (_\theta{k}^{j})^{2}$

The collaborative filtering algorithm can be summarized as follows

Initialize $\theta^1$ , $\theta^2$ … $\theta^{n_{u}}$ and the set of features be $x^1$ , $x^2$ , … , $x^{n_{m}}$ to small random values
Minimize J( $\theta^1$ , $\theta^2$ … $\theta^{n_{u}}$ , $x^1$ , $x^2$ , … , $x^{n_{m}}$ ) using gradient descent. For every
j=1,2, … $n_{u}$ , i= 1,2,.., $n_{m}$
$x_{k}^{i}$ := $x_{k}^{i}$ – $\alpha$ ( $\sigma$ $(\theta^j)^T$ ) $x^i$ – $y^(i,j)\theta_{k}^{j} + \lambda x_{k}^i$

&

$\theta_{k}^{i}$ := $\theta_{k}^{i}$ – $\alpha$ ( $\sigma$ $(\theta^j)^T)x^i - y^(i,j)\theta_{k}^{j} + \lambda x_{k}^i$
Hence for a batsman with parameters $\theta$ and a bowler with (learned) features x, predict the “times out” for the player where the value is not known using $\theta^Tx$

The above derivation for the recommender problem is taken from Machine Learning by Prof Andrew Ng at Coursera from the lecture Collaborative filtering

There are 2 main types of Collaborative Filtering(CF) approaches

User based Collaborative Filtering User-based CF is a memory-based algorithm which tries to mimics word-of-mouth by analyzing rating data from many individuals. The assumption is that users with similar preferences will rate items similarly.
Item based Collaborative Filtering Item-based CF is a model-based approach which produces recommendations based on the relationship between items inferred from the rating matrix. The assumption behind this approach is that users will prefer items that are similar to other items they like.

1a. A note on ROC and Precision-Recall curves

A small note on interpreting ROC & Precision-Recall curves in the post below

ROC Curve: The ROC curve plots the True Positive Rate (TPR) against the False Positive Rate (FPR). Ideally the TPR should increase faster than the FPR and the AUC (area under the curve) should be close to 1

Precision-Recall: The precision-recall curve shows the tradeoff between precision and recall for different threshold. A high area under the curve represents both high recall and high precision, where high precision relates to a low false positive rate, and high recall relates to a low false negative rate

library(reshape2)
library(dplyr)
library(ggplot2)
library(recommenderlab)
library(tidyr)
load("recom_data/batsmenVsBowler20_22.rdata")

2. Define recommender lab helper functions

Helper functions for the RMarkdown notebook are created

eval – Gives details of RMSE, MSE and MAE of ML algorithm
evalRecomMethods – Evaluates different recommender methods and plot the ROC and Precision-Recall curves

# This function returns the error for the chosen algorithm and also predicts the estimates
# for the given data
eval <- function(data, train1, k1,given1,goodRating1,recomType1="UBCF"){
  set.seed(2022)
  e<- evaluationScheme(data,
                       method = "split",
                       train = train1,
                       k = k1,
                       given = given1,
                       goodRating = goodRating1)
  
  r1 <- Recommender(getData(e, "train"), recomType1)
  print(r1)
  
  p1 <- predict(r1, getData(e, "known"), type="ratings")
  print(p1)
  
  error = calcPredictionAccuracy(p1, getData(e, "unknown"))
  
  print(error)
  p2 <- predict(r1, data, type="ratingMatrix")
  p2
}
# This function will evaluate the different recommender algorithms and plot the AUC and ROC curves
evalRecomMethods <- function(data,k1,given1,goodRating1){
  set.seed(2022)
  e<- evaluationScheme(data,
                       method = "cross",
                       k = k1,
                       given = given1,
                       goodRating = goodRating1)
  
  models_to_evaluate <- list(
    `IBCF Cosinus` = list(name = "IBCF", 
                          param = list(method = "cosine")),
    `IBCF Pearson` = list(name = "IBCF", 
                          param = list(method = "pearson")),
    `UBCF Cosinus` = list(name = "UBCF",
                          param = list(method = "cosine")),
    `UBCF Pearson` = list(name = "UBCF",
                          param = list(method = "pearson")),
    `Zufälliger Vorschlag` = list(name = "RANDOM", param=NULL)
  )
  
  n_recommendations <- c(1, 5, seq(10, 100, 10))
  list_results <- evaluate(x = e, 
                           method = models_to_evaluate, 
                           n = n_recommendations)
  plot(list_results, annotate=c(1,3), legend="bottomright")
  plot(list_results, "prec/rec", annotate=3, legend="topleft")
}

3. Batsman performance estimation

The section below regenerates the performance for batsmen based on incomplete data for the different fields in the data frame namely balls faced, fours, sixes, strike rate, times out. The recommender lab allows one to test several different algorithms all at once namely

User based – Cosine similarity method, Pearson similarity
Item based – Cosine similarity method, Pearson similarity
Popular
Random
SVD and a few others

3a. Batting dataframe

head(df)

##   batsman1         bowler1 ballsFaced totalRuns fours sixes  SR timesOut
## 1 A Badoni        A Mishra          0         0     0     0 NaN        0
## 2 A Badoni        A Nortje          0         0     0     0 NaN        0
## 3 A Badoni         A Zampa          0         0     0     0 NaN        0
## 4 A Badoni     Abdul Samad          0         0     0     0 NaN        0
## 5 A Badoni Abhishek Sharma          0         0     0     0 NaN        0
## 6 A Badoni      AD Russell          0         0     0     0 NaN        0

3b Data set and data preparation

For this analysis the data from Cricsheet has been processed using my R package yorkr to obtain the following 2 data sets – batsmenVsBowler – This dataset will contain the performance of the batsmen against the bowler and will capture a) ballsFaced b) totalRuns c) Fours d) Sixes e) SR f) timesOut – bowlerVsBatsmen – This data set will contain the performance of the bowler against the difference batsmen and will include a) deliveries b) runsConceded c) EconomyRate d) wicketsTaken

Obviously many rows/columns will be empty

This is a large data set and hence I have filtered for the period > Jan 2020 and < Dec 2022 which gives 2 datasets a) batsmanVsBowler20_22.rdata b) bowlerVsBatsman20_22.rdata

I also have 2 other datasets of all batsmen and bowlers in these 2 dataset in the files c) all-batsmen20_22.rds d) all-bowlers20_22.rds

You can download the data and this RMarkdown notebook from Github at PlayerPerformanceEstimation

Feel free to download and analyze the data and use any recommendation engine you choose

3c. Exploratory analysis

Initially an exploratory analysis is done on the data

df3 <- select(df, batsman1,bowler1,timesOut)
df6 <- xtabs(timesOut ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
print(df8[1:10,1:10])

##                 A Mishra A Nortje A Zampa Abdul Samad Abhishek Sharma
## A Badoni              NA       NA      NA          NA              NA
## A Manohar             NA       NA      NA          NA              NA
## A Nortje              NA       NA      NA          NA              NA
## AB de Villiers        NA        4       3          NA              NA
## Abdul Samad           NA       NA      NA          NA              NA
## Abhishek Sharma       NA       NA      NA          NA              NA
## AD Russell             1       NA      NA          NA              NA
## AF Milne              NA       NA      NA          NA              NA
## AJ Finch              NA       NA      NA          NA               3
## AJ Tye                NA       NA      NA          NA              NA
##                 AD Russell AF Milne AJ Tye AK Markram Akash Deep
## A Badoni                NA       NA     NA         NA         NA
## A Manohar               NA       NA     NA         NA         NA
## A Nortje                NA       NA     NA         NA         NA
## AB de Villiers           3       NA      3         NA         NA
## Abdul Samad             NA       NA     NA         NA         NA
## Abhishek Sharma         NA       NA     NA         NA         NA
## AD Russell              NA       NA      6         NA         NA
## AF Milne                NA       NA     NA         NA         NA
## AJ Finch                NA       NA     NA         NA         NA
## AJ Tye                  NA       NA     NA         NA         NA

The dots below represent data for which there is no performance data. These cells need to be estimated by the algorithm

set.seed(2022)
r <- as(df8,"realRatingMatrix")
getRatingMatrix(r)[1:15,1:15]

## 15 x 15 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 15 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                               
## A Badoni         . . . . . . . . . . . . . . .
## A Manohar        . . . . . . . . . . . . . . .
## A Nortje         . . . . . . . . . . . . . . .
## AB de Villiers   . 4 3 . . 3 . 3 . . . 4 3 . .
## Abdul Samad      . . . . . . . . . . . . . . .
## Abhishek Sharma  . . . . . . . . . . . 1 . . .
## AD Russell       1 . . . . . . 6 . . . 3 3 3 .
## AF Milne         . . . . . . . . . . . . . . .
## AJ Finch         . . . . 3 . . . . . . 1 . . .
## AJ Tye           . . . . . . . . . . . 1 . . .
## AK Markram       . . . 3 . . . . . . . . . . .
## AM Rahane        9 . . . . 3 . 3 . . . 3 3 . .
## Anmolpreet Singh . . . . . . . . . . . . . . .
## Anuj Rawat       . . . . . . . . . . . . . . .
## AR Patel         . . . . . . . 1 . . . . . . .

r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:15,1:15]

## 15 x 15 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 15 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                              
## AB de Villiers  . 4 3 . . 3 . 3 . . . 4 3 . .
## Abdul Samad     . . . . . . . . . . . . . . .
## Abhishek Sharma . . . . . . . . . . . 1 . . .
## AD Russell      1 . . . . . . 6 . . . 3 3 3 .
## AJ Finch        . . . . 3 . . . . . . 1 . . .
## AM Rahane       9 . . . . 3 . 3 . . . 3 3 . .
## AR Patel        . . . . . . . 1 . . . . . . .
## AT Rayudu       2 . . . . . 1 . . . . 3 . . .
## B Kumar         3 . 3 . . . . . . . . . . 3 .
## BA Stokes       . . . . . . 3 4 . . . 3 . . .
## CA Lynn         . . . . . . . 9 . . . 3 . . .
## CH Gayle        . . . . . 6 . 3 . . . 6 . . .
## CH Morris       . 3 . . . . . . . . . 3 . . .
## D Padikkal      . 4 . . . 3 . . . . . . 3 . .
## DA Miller       . . . . . 3 . . . . . 3 . . .

# Get the summary of the data
summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   3.000   3.463   4.000  21.000

# Normalize the data
r0_m <- normalize(r0)
getRatingMatrix(r0_m)[1:15,1:15]

## 15 x 15 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 15 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                                                       
## AB de Villiers   .         -0.7857143 -1.7857143 .  .       -1.7857143
## Abdul Samad      .          .          .         .  .        .        
## Abhishek Sharma  .          .          .         .  .        .        
## AD Russell      -2.6562500  .          .         .  .        .        
## AJ Finch         .          .          .         . -0.03125  .        
## AM Rahane        4.6041667  .          .         .  .       -1.3958333
## AR Patel         .          .          .         .  .        .        
## AT Rayudu       -2.1363636  .          .         .  .        .        
## B Kumar          0.3636364  .          0.3636364 .  .        .        
## BA Stokes        .          .          .         .  .        .        
## CA Lynn          .          .          .         .  .        .        
## CH Gayle         .          .          .         .  .        1.5476190
## CH Morris        .          0.3500000  .         .  .        .        
## D Padikkal       .          0.6250000  .         .  .       -0.3750000
## DA Miller        .          .          .         .  .       -0.7037037
##                                                                              
## AB de Villiers   .         -1.7857143 . . . -0.7857143 -1.785714  .         .
## Abdul Samad      .          .         . . .  .          .         .         .
## Abhishek Sharma  .          .         . . . -1.6000000  .         .         .
## AD Russell       .          2.3437500 . . . -0.6562500 -0.656250 -0.6562500 .
## AJ Finch         .          .         . . . -2.0312500  .         .         .
## AM Rahane        .         -1.3958333 . . . -1.3958333 -1.395833  .         .
## AR Patel         .         -2.3333333 . . .  .          .         .         .
## AT Rayudu       -3.1363636  .         . . . -1.1363636  .         .         .
## B Kumar          .          .         . . .  .          .         0.3636364 .
## BA Stokes       -0.6086957  0.3913043 . . . -0.6086957  .         .         .
## CA Lynn          .          5.3200000 . . . -0.6800000  .         .         .
## CH Gayle         .         -1.4523810 . . .  1.5476190  .         .         .
## CH Morris        .          .         . . .  0.3500000  .         .         .
## D Padikkal       .          .         . . .  .         -0.375000  .         .
## DA Miller        .          .         . . . -0.7037037  .         .         .

4. Create a visual representation of the rating data before and after the normalization

The histograms show the bias in the data is removed after normalization

r0=r[(m=rowCounts(r) > 10),]
getRatingMatrix(r0)[1:15,1:10]

## 15 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                    
## AB de Villiers  . 4 3 . . 3 . 3 . .
## Abdul Samad     . . . . . . . . . .
## Abhishek Sharma . . . . . . . . . .
## AD Russell      1 . . . . . . 6 . .
## AJ Finch        . . . . 3 . . . . .
## AM Rahane       9 . . . . 3 . 3 . .
## AR Patel        . . . . . . . 1 . .
## AT Rayudu       2 . . . . . 1 . . .
## B Kumar         3 . 3 . . . . . . .
## BA Stokes       . . . . . . 3 4 . .
## CA Lynn         . . . . . . . 9 . .
## CH Gayle        . . . . . 6 . 3 . .
## CH Morris       . 3 . . . . . . . .
## D Padikkal      . 4 . . . 3 . . . .
## DA Miller       . . . . . 3 . . . .

#Plot ratings
image(r0, main = "Raw Ratings")

#Plot normalized ratings
r0_m <- normalize(r0)
getRatingMatrix(r0_m)[1:15,1:15]

## 15 x 15 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 15 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                                                       
## AB de Villiers   .         -0.7857143 -1.7857143 .  .       -1.7857143
## Abdul Samad      .          .          .         .  .        .        
## Abhishek Sharma  .          .          .         .  .        .        
## AD Russell      -2.6562500  .          .         .  .        .        
## AJ Finch         .          .          .         . -0.03125  .        
## AM Rahane        4.6041667  .          .         .  .       -1.3958333
## AR Patel         .          .          .         .  .        .        
## AT Rayudu       -2.1363636  .          .         .  .        .        
## B Kumar          0.3636364  .          0.3636364 .  .        .        
## BA Stokes        .          .          .         .  .        .        
## CA Lynn          .          .          .         .  .        .        
## CH Gayle         .          .          .         .  .        1.5476190
## CH Morris        .          0.3500000  .         .  .        .        
## D Padikkal       .          0.6250000  .         .  .       -0.3750000
## DA Miller        .          .          .         .  .       -0.7037037
##                                                                              
## AB de Villiers   .         -1.7857143 . . . -0.7857143 -1.785714  .         .
## Abdul Samad      .          .         . . .  .          .         .         .
## Abhishek Sharma  .          .         . . . -1.6000000  .         .         .
## AD Russell       .          2.3437500 . . . -0.6562500 -0.656250 -0.6562500 .
## AJ Finch         .          .         . . . -2.0312500  .         .         .
## AM Rahane        .         -1.3958333 . . . -1.3958333 -1.395833  .         .
## AR Patel         .         -2.3333333 . . .  .          .         .         .
## AT Rayudu       -3.1363636  .         . . . -1.1363636  .         .         .
## B Kumar          .          .         . . .  .          .         0.3636364 .
## BA Stokes       -0.6086957  0.3913043 . . . -0.6086957  .         .         .
## CA Lynn          .          5.3200000 . . . -0.6800000  .         .         .
## CH Gayle         .         -1.4523810 . . .  1.5476190  .         .         .
## CH Morris        .          .         . . .  0.3500000  .         .         .
## D Padikkal       .          .         . . .  .         -0.375000  .         .
## DA Miller        .          .         . . . -0.7037037  .         .         .

image(r0_m, main = "Normalized Ratings")

set.seed(1234)
hist(getRatings(r0), breaks=25)

hist(getRatings(r0_m), breaks=25)

4a. Data for analysis

The data frame of the batsman vs bowlers from the period 2020 -2022 is read as a dataframe. To remove rows with very low number of ratings(timesOut, SR, Fours, Sixes etc), the rows are filtered so that there are at least more 10 values in the row. For the player estimation the dataframe is converted into a wide-format as a matrix (m x n) of batsman x bowler with each of the columns of the dataframe i.e. timesOut, SR, fours or sixes. These different matrices can be considered as a rating matrix for estimation.

A similar approach is taken for estimating bowler performance. Here a wide form matrix (m x n) of bowler x batsman is created for each of the columns of deliveries, runsConceded, ER, wicketsTaken

5. Batsman’s times Out

The code below estimates the number of times the batsmen would lose his/her wicket to the bowler. As discussed in the algorithm above, the recommendation engine will make an initial estimate features for the bowler and an initial estimate for the parameter vector for the batsmen. Then using gradient descent the recommender engine will determine the feature and parameter values such that the over Mean Squared Error is minimum

From the plot for the different algorithms it can be seen that UBCF performs the best. However the AUC & ROC curves are not optimal and the AUC> 0.5

df3 <- select(df, batsman1,bowler1,timesOut)
df6 <- xtabs(timesOut ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
# Filter only rows where the row count is > 10
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                    
## AB de Villiers  . 4 3 . . 3 . 3 . .
## Abdul Samad     . . . . . . . . . .
## Abhishek Sharma . . . . . . . . . .
## AD Russell      1 . . . . . . 6 . .
## AJ Finch        . . . . 3 . . . . .
## AM Rahane       9 . . . . 3 . 3 . .
## AR Patel        . . . . . . . 1 . .
## AT Rayudu       2 . . . . . 1 . . .
## B Kumar         3 . 3 . . . . . . .
## BA Stokes       . . . . . . 3 4 . .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   3.000   3.463   4.000  21.000

# Evaluate the different plotting methods
evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

#Evaluate the error
a=eval(r0[1:dim(r0)[1]],0.8,k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 70 users.
## 18 x 145 rating matrix of class 'realRatingMatrix' with 1755 ratings.
##     RMSE      MSE      MAE 
## 2.069027 4.280872 1.496388

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
m=as(c,"data.frame")
names(m) =c("batsman","bowler","TimesOut")

6. Batsman’s Strike rate

This section deals with the Strike rate of batsmen versus bowlers and estimates the values for those where the data is incomplete using UBCF method.

Even here all the algorithms do not perform too efficiently. I did try out a few variations but could not lower the error (suggestions welcome!!)

df3 <- select(df, batsman1,bowler1,SR)
df6 <- xtabs(SR ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                                                           
## AB de Villiers   96.8254 171.4286  33.33333  . 66.66667 223.07692   .     
## Abdul Samad       .      228.0000   .        .  .       100.00000   .     
## Abhishek Sharma 150.0000   .        .        .  .        66.66667   .     
## AD Russell      111.4286   .        .        .  .         .         .     
## AJ Finch        250.0000 116.6667   .        . 50.00000  85.71429 112.5000
## AJ Tye            .        .        .        .  .         .       100.0000
## AK Markram        .        .        .       50  .         .         .     
## AM Rahane       121.1111   .        .        .  .       113.82979 117.9487
## AR Patel        183.3333   .      200.00000  .  .       433.33333   .     
## AT Rayudu       126.5432 200.0000 122.22222  .  .       105.55556   .     
##                                
## AB de Villiers  109.52381 .   .
## Abdul Samad       .       .   .
## Abhishek Sharma   .       .   .
## AD Russell      195.45455 .   .
## AJ Finch          .       .   .
## AJ Tye            .       .   .
## AK Markram        .       .   .
## AM Rahane        33.33333 . 200
## AR Patel        171.42857 .   .
## AT Rayudu       204.76190 .   .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   5.882  85.714 116.667 128.529 160.606 600.000

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

a=eval(r0[1:dim(r0)[1]],0.8, k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 105 users.
## 27 x 145 rating matrix of class 'realRatingMatrix' with 3220 ratings.
##       RMSE        MSE        MAE 
##   77.71979 6040.36508   58.58484

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
n=as(c,"data.frame")
names(n) =c("batsman","bowler","SR")

7. Batsman’s Sixes

The snippet of code estimes the sixes of the batsman against bowlers. The ROC and AUC curve for UBCF looks a lot better here, as it significantly greater than 0.5

df3 <- select(df, batsman1,bowler1,sixes)
df6 <- xtabs(sixes ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                      
## AB de Villiers  3 3 . . . 18 .  3 . .
## AD Russell      3 . . . .  . . 12 . .
## AJ Finch        2 . . . .  . .  . . .
## AM Rahane       7 . . . .  3 1  . . .
## AR Patel        4 . 3 . .  6 .  1 . .
## AT Rayudu       5 2 . . .  . .  1 . .
## BA Stokes       . . . . .  . .  . . .
## CA Lynn         . . . . .  . .  9 . .
## CH Gayle       17 . . . . 17 .  . . .
## CH Morris       . . 3 . .  . .  . . .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    3.00    3.00    4.68    6.00   33.00

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

## Timing stopped at: 0.003 0 0.002

a=eval(r0[1:dim(r0)[1]],0.8, k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 52 users.
## 14 x 145 rating matrix of class 'realRatingMatrix' with 1634 ratings.
##      RMSE       MSE       MAE 
##  3.529922 12.460350  2.532122

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
o=as(c,"data.frame")
names(o) =c("batsman","bowler","Sixes")

8. Batsman’s Fours

The code below estimates 4s for the batsmen

df3 <- select(df, batsman1,bowler1,fours)
df6 <- xtabs(fours ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                      
## AB de Villiers   . 1 . . . 24 . 3 . .
## Abhishek Sharma  . . . . .  . . . . .
## AD Russell       1 . . . .  . . 9 . .
## AJ Finch         . 1 . . .  3 2 . . .
## AK Markram       . . . . .  . . . . .
## AM Rahane       11 . . . .  8 7 . . 3
## AR Patel         . . . . .  . . 3 . .
## AT Rayudu       11 2 3 . .  6 . 6 . .
## BA Stokes        1 . . . .  . . . . .
## CA Lynn          . . . . .  . . 6 . .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   4.000   6.339   9.000  55.000

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

## Timing stopped at: 0.008 0 0.008

## Warning in .local(x, method, ...): 
##   Recommender 'UBCF Pearson' has failed and has been removed from the results!

a=eval(r0[1:dim(r0)[1]],0.8, k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 67 users.
## 17 x 145 rating matrix of class 'realRatingMatrix' with 2083 ratings.
##      RMSE       MSE       MAE 
##  5.486661 30.103447  4.060990

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
p=as(c,"data.frame")
names(p) =c("batsman","bowler","Fours")

9. Batsman’s Total Runs

The code below estimates the total runs that would have scored by the batsman against different bowlers

df3 <- select(df, batsman1,bowler1,totalRuns)
df6 <- xtabs(totalRuns ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                          
## A Badoni         .  . . . .   . .   . . .
## A Manohar        .  . . . .   . .   . . .
## A Nortje         .  . . . .   . .   . . .
## AB de Villiers  61 36 3 . 6 261 .  69 . .
## Abdul Samad      . 57 . . .  12 .   . . .
## Abhishek Sharma  3  . . . .   6 .   . . .
## AD Russell      39  . . . .   . . 129 . .
## AF Milne         .  . . . .   . .   . . .
## AJ Finch        15  7 . . 3  18 9   . . .
## AJ Tye           .  . . . .   . 4   . . .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   24.00   41.36   54.00  452.00

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given1=7,goodRating1=median(getRatings(r0)))

a=eval(r0[1:dim(r0)[1]],0.8, k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 105 users.
## 27 x 145 rating matrix of class 'realRatingMatrix' with 3256 ratings.
##       RMSE        MSE        MAE 
##   41.50985 1723.06788   29.52958

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
q=as(c,"data.frame")
names(q) =c("batsman","bowler","TotalRuns")

10. Batsman’s Balls Faced

The snippet estimates the balls faced by batsmen versus bowlers

df3 <- select(df, batsman1,bowler1,ballsFaced)
df6 <- xtabs(ballsFaced ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Mishra', 'A Nortje', 'A Zampa' ... ]]

##                                         
## A Badoni         .  . . . .   . .  . . .
## A Manohar        .  . . . .   . .  . . .
## A Nortje         .  . . . .   . .  . . .
## AB de Villiers  63 21 9 . 9 117 . 63 . .
## Abdul Samad      . 25 . . .  12 .  . . .
## Abhishek Sharma  2  . . . .   9 .  . . .
## AD Russell      35  . . . .   . . 66 . .
## AF Milne         .  . . . .   . .  . . .
## AJ Finch         6  6 . . 6  21 8  . . .
## AJ Tye           .  . . . .   9 4  . . .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   18.00   30.21   39.00  384.00

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

a=eval(r0[1:dim(r0)[1]],0.8, k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 112 users.
## 28 x 145 rating matrix of class 'realRatingMatrix' with 3378 ratings.
##       RMSE        MSE        MAE 
##   33.91251 1150.05835   23.39439

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
r=as(c,"data.frame")
names(r) =c("batsman","bowler","BallsFaced")

11. Generate the Batsmen Performance Estimate

This code generates the estimated dataframe with known and ‘predicted’ values

a1=merge(m,n,by=c("batsman","bowler"))
a2=merge(a1,o,by=c("batsman","bowler"))
a3=merge(a2,p,by=c("batsman","bowler"))
a4=merge(a3,q,by=c("batsman","bowler"))
a5=merge(a4,r,by=c("batsman","bowler"))
a6= select(a5, batsman,bowler,BallsFaced,TotalRuns,Fours, Sixes, SR,TimesOut)
head(a6)

##          batsman          bowler BallsFaced TotalRuns Fours Sixes  SR TimesOut
## 1 AB de Villiers        A Mishra         94       124     7     5 144        5
## 2 AB de Villiers        A Nortje         26        42     4     3 148        3
## 3 AB de Villiers         A Zampa         28        42     5     7 106        4
## 4 AB de Villiers Abhishek Sharma         22        28     0    10 136        5
## 5 AB de Villiers      AD Russell         70       135    14    12 207        4
## 6 AB de Villiers        AF Milne         31        45     6     6 130        3

12. Bowler analysis

Just like the batsman performance estimation we can consider the bowler’s performances also for estimation. Consider the following table

As in the batsman analysis, for every batsman a set of features like (“strong backfoot player”, “360 degree player”,“Power hitter”) can be estimated with a set of initial values. Also every bowler will have an associated parameter vector θθ. Different bowlers will have performance data for different set of batsmen. Based on the initial estimate of the features and the parameters, gradient descent can be used to minimize actual values {for e.g. wicketsTaken(ratings)}.

load("recom_data/bowlerVsBatsman20_22.rdata")

12a. Bowler dataframe

Inspecting the bowler dataframe

head(df2)

##    bowler1        batsman1 balls runsConceded       ER wicketTaken
## 1 A Mishra        A Badoni     0            0 0.000000           0
## 2 A Mishra       A Manohar     0            0 0.000000           0
## 3 A Mishra        A Nortje     0            0 0.000000           0
## 4 A Mishra  AB de Villiers    63           61 5.809524           0
## 5 A Mishra     Abdul Samad     0            0 0.000000           0
## 6 A Mishra Abhishek Sharma     2            3 9.000000           0

names(df2)

## [1] "bowler1"      "batsman1"     "balls"        "runsConceded" "ER"          
## [6] "wicketTaken"

13. Balls bowled by bowler

The below section estimates the balls bowled for each bowler. We can see that UBCF Pearson and UBCF Cosine both perform well

df3 <- select(df2, bowler1,batsman1,balls)
df6 <- xtabs(balls ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Badoni', 'A Manohar', 'A Nortje' ... ]]

##                                          
## A Mishra        . . .  63  .  2 35 .  6 .
## A Nortje        . . .  21 25  .  . .  6 .
## A Zampa         . . .   9  .  .  . .  . .
## Abhishek Sharma . . .   9  .  .  . .  6 .
## AD Russell      . . . 117 12  9  . . 21 9
## AF Milne        . . .   .  .  .  . .  8 4
## AJ Tye          . . .  63  .  . 66 .  . .
## Akash Deep      . . .   .  .  .  . .  . .
## AR Patel        . . . 188  5  1 84 . 29 5
## Arshdeep Singh  . . .   6  6 24 18 . 12 .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   18.00   29.61   36.00  384.00

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

a=eval(r0[1:dim(r0)[1]],0.8,k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 96 users.
## 24 x 195 rating matrix of class 'realRatingMatrix' with 3954 ratings.
##      RMSE       MSE       MAE 
##  30.72284 943.89294  19.89204

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
s=as(c,"data.frame")
names(s) =c("bowler","batsman","BallsBowled")

14. Runs conceded by bowler

This section estimates the runs conceded by the bowler. The UBCF Cosinus algorithm performs the best with TPR increasing fastewr than FPR

df3 <- select(df2, bowler1,batsman1,runsConceded)
df6 <- xtabs(runsConceded ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Badoni', 'A Manohar', 'A Nortje' ... ]]

##                                            
## A Mishra        . . .  61  .  3  41 . 15  .
## A Nortje        . . .  36 57  .   . .  8  .
## A Zampa         . . .   3  .  .   . .  .  .
## Abhishek Sharma . . .   6  .  .   . .  3  .
## AD Russell      . . . 276 12  6   . . 21  .
## AF Milne        . . .   .  .  .   . . 10  4
## AJ Tye          . . .  69  .  . 138 .  .  .
## Akash Deep      . . .   .  .  .   . .  .  .
## AR Patel        . . . 205  5  . 165 . 33 13
## Arshdeep Singh  . . .  18  3 51  51 .  6  .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    9.00   24.00   41.34   54.00  458.00

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

## Timing stopped at: 0.004 0 0.004

## Warning in .local(x, method, ...): 
##   Recommender 'UBCF Pearson' has failed and has been removed from the results!

a=eval(r0[1:dim(r0)[1]],0.8,k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 95 users.
## 24 x 195 rating matrix of class 'realRatingMatrix' with 3820 ratings.
##       RMSE        MSE        MAE 
##   43.16674 1863.36749   30.32709

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
t=as(c,"data.frame")
names(t) =c("bowler","batsman","RunsConceded")

15. Economy Rate of the bowler

This section computes the economy rate of the bowler. The performance is not all that good

df3 <- select(df2, bowler1,batsman1,ER)
df6 <- xtabs(ER ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Badoni', 'A Manohar', 'A Nortje' ... ]]

##                                                                       
## A Mishra        . . .  5.809524  .     9.00  7.028571 . 15.000000  .  
## A Nortje        . . . 10.285714 13.68  .     .        .  8.000000  .  
## A Zampa         . . .  2.000000  .     .     .        .  .         .  
## Abhishek Sharma . . .  4.000000  .     .     .        .  3.000000  .  
## AD Russell      . . . 14.153846  6.00  4.00  .        .  6.000000  .  
## AF Milne        . . .  .         .     .     .        .  7.500000  6.0
## AJ Tye          . . .  6.571429  .     .    12.545455 .  .         .  
## Akash Deep      . . .  .         .     .     .        .  .         .  
## AR Patel        . . .  6.542553  6.00  .    11.785714 .  6.827586 15.6
## Arshdeep Singh  . . . 18.000000  3.00 12.75 17.000000 .  3.000000  .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3529  5.2500  7.1126  7.8139  9.8000 36.0000

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

## Timing stopped at: 0.003 0 0.004

## Warning in .local(x, method, ...): 
##   Recommender 'UBCF Pearson' has failed and has been removed from the results!

a=eval(r0[1:dim(r0)[1]],0.8,k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 95 users.
## 24 x 195 rating matrix of class 'realRatingMatrix' with 3839 ratings.
##      RMSE       MSE       MAE 
##  4.380680 19.190356  3.316556

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
u=as(c,"data.frame")
names(u) =c("bowler","batsman","EconomyRate")

16. Wickets Taken by bowler

The code below computes the wickets taken by the bowler versus different batsmen

df3 <- select(df2, bowler1,batsman1,wicketTaken)
df6 <- xtabs(wicketTaken ~ ., df3)
df7 <- as.data.frame.matrix(df6)
df8 <- data.matrix(df7)
df8[df8 == 0] <- NA
r <- as(df8,"realRatingMatrix")
r0=r[(rowCounts(r) > 10),]
getRatingMatrix(r0)[1:10,1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 10 column names 'A Badoni', 'A Manohar', 'A Nortje' ... ]]

##                                   
## A Mishra       . . . . . . 1 . . .
## A Nortje       . . . 4 . . . . . .
## A Zampa        . . . 3 . . . . . .
## AD Russell     . . . 3 . . . . . .
## AJ Tye         . . . 3 . . 6 . . .
## AR Patel       . . . 4 . 1 3 . 1 1
## Arshdeep Singh . . . 3 . . 3 . . .
## AS Rajpoot     . . . . . . 3 . . .
## Avesh Khan     . . . . . . 1 . 3 .
## B Kumar        . . . 9 . . 3 . 1 .

summary(getRatings(r0))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   3.000   3.423   3.000  21.000

evalRecomMethods(r0[1:dim(r0)[1]],k1=5,given=7,goodRating1=median(getRatings(r0)))

## Timing stopped at: 0.003 0 0.003

## Warning in .local(x, method, ...): 
##   Recommender 'UBCF Pearson' has failed and has been removed from the results!

a=eval(r0[1:dim(r0)[1]],0.8,k1=5,given1=7,goodRating1=median(getRatings(r0)),"UBCF")

## Recommender of type 'UBCF' for 'realRatingMatrix' 
## learned using 64 users.
## 16 x 195 rating matrix of class 'realRatingMatrix' with 1908 ratings.
##     RMSE      MSE      MAE 
## 2.672677 7.143203 1.956934

b=round(as(a,"matrix")[1:10,1:10])
c <- as(b,"realRatingMatrix")
v=as(c,"data.frame")
names(v) =c("bowler","batsman","WicketTaken")

17. Generate the Bowler Performance estmiate

The entire dataframe is regenerated with known and ‘predicted’ values

r1=merge(s,t,by=c("bowler","batsman"))
r2=merge(r1,u,by=c("bowler","batsman"))
r3=merge(r2,v,by=c("bowler","batsman"))
r4= select(r3,bowler, batsman, BallsBowled,RunsConceded,EconomyRate, WicketTaken)
head(r4)

##     bowler         batsman BallsBowled RunsConceded EconomyRate WicketTaken
## 1 A Mishra  AB de Villiers         102          144           8           4
## 2 A Mishra     Abdul Samad          13           20           7           4
## 3 A Mishra Abhishek Sharma          14           26           8           2
## 4 A Mishra      AD Russell          47           85           9           3
## 5 A Mishra        AJ Finch          45           61          11           4
## 6 A Mishra          AJ Tye          14           20           5           4

18. Conclusion

This post showed an approach for performing the Batsmen Performance Estimate & Bowler Performance Estimate. The performance of the recommender engine could have been better. In any case, I think this approach will work for player estimation provided the recommender algorithm is able to achieve a high degree of accuracy. This will be a good way to estimate as the algorithm will be able to determine features and nuances of batsmen and bowlers which cannot be captured by data.

References

Also see

To see all posts click Index of posts

Understanding Neural Style Transfer with Tensorflow and Keras

Neural Style Transfer (NST) is a fascinating area of Deep Learning and Convolutional Neural Networks. NST is an interesting technique, in which the style from an image, known as the ‘style image’ is transferred to another image ‘content image’ and we get a third a image which is a generated image which has the content of the original image and the style of another image.

NST can be used to reimagine how famous painters like Van Gogh, Claude Monet or a Picasso would have visualised a scenery or architecture. NST uses Convolutional Neural Networks (CNNs) to achieve this artistic style transfer from one image to another. NST was originally implemented by Gati et al., in their paper Neural Algorithm of Artistic Style. Convolutional Neural Networks have been very successful in image classification image recognition et cetera. CNN networks have also been have also generated very interesting pictures using Neural Style Transfer which will be shown in this post. An interesting aspect of CNN’s is that the first couple of layers in the CNN capture basic features of the image like edges and pixel values. But as we go deeper into the CNN, the network captures higher level features of the input image.

To get started with Neural Style transfer we will be using the VGG19 pre-trained network. The VGG19 CNN is a compact pre-trained your network which can be used for performing the NST. However, we could have also used Resnet or InceptionV3 networks for this purpose but these are very large networks. The idea of using a network trained on a different task and applying it to a new task is called transfer learning.

What needs to be done to transfer the style from one of the image to another image. This brings us to the question – What is ‘style’? What is it that distinguishes Van Gogh’s painting or Picasso’s cubist art. Convolutional Neural Networks capture basic features in the lower layers and much more complex features in the deeper layers. Style can be computed by taking the correlation of the feature maps in a layer L. This is my interpretation of how style is captured. Since style is intrinsic to the image, it implies that the style feature would exist across all the filters in a layer. Hence, to pick up this style we would need to get the correlation of the filters across channels of a lawyer. This is computed mathematically, using the Gram matrix which calculates the correlation of the activation of a the filter by the style image and generated image

To transfer the style from one image to the content image we need to do two parallel operations while doing forward propagation
– Compute the content loss between the source image and the generated image
– Compute the style loss between the style image and the generated image
– Finally we need to compute the total loss

In order to get transfer the style from the ‘style’ image to the ‘content ‘image resulting in a ‘generated’ image the total loss has to be minimised. Therefore backward propagation with gradient descent is done to minimise the total loss comprising of the content and style loss.

Initially we make the Generated Image ‘G’ the same as the source image ‘S’

The content loss at layer ‘l’

$L_{content} = 1/2 \sum_{i}^{j} ( F^{l}_{i,j} - P^{l}_{i,j})^{2}$

where $F^{l}_{i,j}$ and $P^{l}_{i,j}$ represent the activations at layer ‘l’ in a filter i, at position ‘j’. The intuition is that the activations will be same for similar source and generated image. We need to minimise the content loss so that the generated stylized image is as close to the original image as possible. An intermediate layer of VGG19 block5_conv2 is used

The Style layers that are are used are

style_layers = [‘block1_conv1’,

‘block2_conv1’,

‘block3_conv1’,

‘block4_conv1’,

‘block5_conv1’]

To compute the Style Loss the Gram matrix needs to be computed. The Gram Matrix is computed by unrolling the filters as shown below (source: Convolutional Neural Networks by Prof Andrew Ng, Coursera). The result is a matrix of size $n_{c}$ x $n_{c}$ where $n_{c}$ is the number of channels

The above diagram shows the filters of height $n_{H}$ and width $n_{W}$ with $n_{C}$ channels

The contribution of layer ‘l’ to style loss is given by

$L^{'}_{style} = \frac{\sum_{i}^{j} (G^{2}_{i,j} - A^l{i,j})^2}{4N^{2}_{l}M^{2}_{l}}$

where $G_{i,j}$ and $A_{i,j}$ are the Gram matrices of the style and generated images respectively. By minimising the distance in the gram matrices of the style and generated image we can ensure that generated image is a stylized version of the original image similar to the style image

The total loss is given by

$L_{total} = \alpha L_{content} + \beta L_{style}$

Back propagation with gradient descent works to minimise the content loss between the source and generated image, while the style loss tries to minimise the discrepancies in the style of the style image and generated image. Running through forward and backpropagation through several epochs successfully transfers the style from the style image to the source image.

You can check the Notebook at Neural Style Transfer

Note: The code in this notebook is largely based on the Neural Style Transfer tutorial from Tensorflow, though I may have taken some changes from other blogs. I also made a few changes to the code in this tutorial, like removing the scaling factor, or the class definition (Personally, I belong to the old school (C language) and am not much in love with the ‘self.”..All references are included below

Note: Here is a interesting thought. Could we do a Neural Style Transfer in music? Imagine Carlos Santana playing ‘Hotel California’ or Brian May style in ‘Another brick in the wall’. While our first reaction would be that it may not sound good as we are used to style of these songs, we may be surprised by a possible style transfer. This is definitely music to the ears!

Here are few runs from this

A) Run 1

1. Neural Style Transfer – a) Content Image – My portrait. b) Style Image – Wassily Kadinsky Oil on canvas, 1913, Vassily Kadinsky’s composition

2. Result of Neural Style Transfer

2) Run 2

a) Content Image – Portrait of my parents b) Style Image – Vincent Van Gogh’s ,Starry Night Oil on canvas 1889

2. Result of Neural Style Transfer

Run 3

1. Content Image – Caesar 2 (Masai Mara- 20 Jun 2018). Style Image – The Great Wave at Kanagawa – Katsushika Hokosai, 1826-1833

Screenshot 2020-04-12 at 12.40.44 PM

2. Result of Neural Style Transfer

lkg

Run 4

1. Content Image – Junagarh Fort , Rajasthan Sep 2016 b) Style Image – Le Pont Japonais by Claude Monet, Oil on canvas, 1920

2. Result of Neural Style Transfer

Neural Style Transfer is a very ingenious idea which shows that we can segregate the style of a painting and transfer to another image.

References

1. A Neural Algorithm of Artistic Style, Leon A. Gatys, Alexander S. Ecker, Matthias Bethge
2. Neural style transfer
3. Neural Style Transfer: Creating Art with Deep Learning using tf.keras and eager execution
4. Convolutional Neural Network, DeepLearning.AI Specialization, Prof Andrew Ng
5. Intuitive Guide to Neural Style Transfer

Using Reinforcement Learning to solve Gridworld

“Take up one idea. Make that one idea your life — think of it, dream of it, live on that idea. Let the brain, muscles, nerves, every part of your body, be full of that idea, and just leave every other idea alone. This is the way to success.”

– Swami Vivekananda

“Be the change you want to see in the world”

– Mahatma Gandhi

“If you want to shine like the sun, first burn like the sun”

-Shri A.P.J Abdul Kalam

Reinforcement Learning

Reinforcement Learning (RL) involves decision making under uncertainty which tries to maximize return over successive states.There are four main elements of a Reinforcement Learning system: a policy, a reward signal, a value function. The policy is a mapping from the states to actions or a probability distribution of actions. Every action the agent takes results in a numerical reward. The agent’s sole purpose is to maximize the reward in the long run.

Reinforcement Learning is very different from Supervised, Unsupervised and Semi-supervised learning where the data is either labeled, unlabeled or partially labeled and the learning algorithm tries to learn the target values from the input features which is then used either for inference or prediction. In unsupervised the intention is to extract patterns from the data. In Reinforcement Learning the agent/robot takes action in each state based on the reward it would get for a particular action in a specific state with the goal of maximizing the reward. In many ways Reinforcement Learning is similar to how human beings and animals learn. Every action we take is with the goal of increasing our overall happiness, contentment, money,fame, power over the opposite!

RL has been used very effectively in many situations, the most famous is AlphaGo from Deep Mind, the first computer program to defeat a professional Go player in the Go game, which is supposed to be extremely complex. Also AlphaZero, from DeepMind has a higher ELO rating than that of Stockfish and was able to beat Stockfish 800+ times in 1000 chess matches. Take a look at DeepMind

In this post, I use some of the basic concepts of Reinforcment Learning to solve Grids (mazes). With this we can solve mazes, with arbitrary size, shape and complexity fairly easily. The RL algorithm can find the optimal path through the maze. Incidentally, I recollect recursive algorithms in Data Structures book which take a much more complex route with a lot of back tracking to solve maze problems

Reinforcement Learning involves decision making under uncertainty which tries to maximize return over successive states.There are four main elements of a Reinforcement Learning system: a policy, a reward signal, a value function. The policy is a mapping from the states to actions or a probability distribution of actions. Every action the agent takes results in a numerical reward. The agent’s sole purpose is to maximize the reward in the long run.

The reward indicates the immediate return, a value function specifies the return in the long run. Value of a state is the expected reward that an agent can accrue.

The agent/robot takes an action in At in state St and moves to state S’t anf gets a reward Rt+1 as shown

An agent will seek to maximize the overall return as it transition across states
The expected return can be expressed as
$G_{t} = R_{t+1} + \gamma G_{t+1}$ where $G_{t}$ is the expected return in time t and the discounted expected return $G_{t+1}$ in time t+1

A policy is a mapping from states to probabilities of selecting each possible action. If the agent is following policy $\pi$ at time t, then $\pi(a|s)$ is the probability that $A_{t}$ = a if $S_{t}$ = s.

The value function of a state s under a policy $\pi$ , denoted $v_{\pi}(s)$ , is the expected return when starting in s and following $\pi$ thereafter

This can be written as

$v_{\pi}(s) = E_{\pi}[G_{t} |S_{t}=s] = E_{\pi}[\sum_{k=0}^{k=Inf} \gamma^{k}R_{t+k+1}|S_{t}=s]$

= $E_{\pi}[R_{t+1} + \gamma G_{t+1} |S_{t}=s]$

$v_{\pi}(s)=\sum_{a} \pi(a|s) \sum_{s',r} p(s',r|s,a)[r+\gamma*v_{\pi}(s')]$

Similarly the action value function gives the expected return when taking an action ‘a’ in state ‘s’
$q_{\pi}(s,a)= \sum_{s',r} p(s',r|s,a)[r+\gamma*\pi(a|s)q_{\pi}(s',a')]$

These are Bellman’s equation for the state value function

The Bellman equations give the equation for each of the state

The Bellman optimality equations give the optimal policy of choosing specific actions in specific states to achieve the maximum reward and reach the goal efficiently. They are given as

$v_{*}(s)=max_{a}\sum_{s',r} p(s',r|s,a)[r+\gamma*v_{*}(s')]$

$q_{*}(s,a)=\sum_{s',r} p(s',r|s,a)[r+\gamma*max_{a}q_{*}(s',a')]$

The Bellman equations cannot be used directly in goal directed problems and dynamic programming is used instead where the value functions are computed iteratively

n this post I solve Grids using Reinforcement Learning. In the problem below the Maze has 2 end states as shown in the corner. There are four possible actions in each state up, down, right and left. If an action in a state takes it out of the grid then the agent remains in the same state. All actions have a reward of -1 while the end states have a reward of 0

This is shown as

where the reward for any transition is Rt=−1 except the transition to the end states at the corner which have a reward of 0. The policy is a uniform policy with all actions being equi-probable with a probability of 1/4 or 0.25

You can fork/clone the code from my Github repository – Gridworld

Note: This post shows 3 different grids each with slightly more complexity and uses 3 methods

a) Bellman Update

b) Greedification

c) Bellman Optimality Update

with dynamic programming to solve the Grids

1. Gridworld-1

In [1]:

import numpy as np
import random

In [2]:

gamma = 1 # discounting rate
gridSize = 4
rewardValue = -1
terminationStates = [[0,0], [gridSize-1, gridSize-1]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

The action value provides the next state for a given action in a state and the accrued reward

In [3]:

def actionValue(initialPosition,action):
    if initialPosition in terminationStates:
        finalPosition = initialPosition
        reward=0
    else:
        #Compute final position
        finalPosition = np.array(initialPosition) + np.array(action)
        reward= rewardValue
    # If the action moves the finalPosition out of the grid, stay in same cell
    if -1 in finalPosition or gridSize in finalPosition:
        finalPosition = initialPosition
        reward= rewardValue
    
    #print(finalPosition)
    return finalPosition, reward

1a. Bellman Update

In [4]:

# Initialize valueMap and valueMap1
valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]

In [5]:

def policy_evaluation(numIterations,gamma,theta,valueMap):
    for i in range(numIterations):
        delta=0
        for state in states:
            weightedRewards=0
            for action in actions:
                finalPosition,reward = actionValue(state,action)
                weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
            valueMap1[state[0],state[1]]=weightedRewards
            delta =max(delta,abs(weightedRewards-valueMap[state[0],state[1]]))
        valueMap = np.copy(valueMap1)
        if(delta < 0.01):                                                
            print(valueMap)
            break

In [6]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
policy_evaluation(1000,1,0.001,valueMap)

[[  0.         -13.89528403 -19.84482978 -21.82635535]
 [-13.89528403 -17.86330422 -19.84586777 -19.84482978]
 [-19.84482978 -19.84586777 -17.86330422 -13.89528403]
 [-21.82635535 -19.84482978 -13.89528403   0.        ]]

Findings

The valueMap is the result of several sweeps through all the states. It can be seen that the cells in the corner state have a higher value. We can start on any cell in the grid and move in the direction which is greater than the current state and we will reach the end state

1b. Greedify

The previous alogirthm while it works is somewhat inefficient as we have to sweep over the states to compute the state value function. The approach below works on the same problem but after each computation of the value function, a greedifications takes place to ensure that the action with the highest return is selected after which the policy π is followed

To make the transitions clearer I also create another grid which shows the path from any cell to the end states as

‘u’ – up

‘d’ – down

‘r’ – right

‘l’ – left

Important note: If there are several alternative actions with equal value then the algorithm will break the tie randomly

In [7]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [8]:

# Compute the value state function for the Grid
def policy_evaluate(states,actions,gamma,valueMap):
    #print("iterations=",i)
    for state in states:
        weightedRewards=0
        for action in actions:
            finalPosition,reward = actionValue(state,action)
            weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Set the computed weighted rewards to valueMap1
        valueMap1[state[0],state[1]]=weightedRewards
    # Copy to original valueMap
    valueMap = np.copy(valueMap1)
    return(valueMap)

In [9]:

def argmax(q_values):
    idx=np.argmax(q_values)
    return(np.random.choice(np.where(a==a[idx])[0].tolist()))


# Compute the best action in each state
def greedify_policy(state,pi,pi1,gamma,valueMap):  
        q_values=np.zeros(len(actions))
        for idx,action in enumerate(actions):
            finalPosition,reward = actionValue(state,action)
            q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Find the index of the action for which the q_value is 
        idx=q_values.argmax()
        pi[state[0],state[1]]=idx 
        if(idx == 0):
            pi1[state[0],state[1]]='u'
        elif(idx == 1):
            pi1[state[0],state[1]]='d'
        elif(idx == 2):
            pi1[state[0],state[1]]='r'
        elif(idx == 3):
            pi1[state[0],state[1]]='l'

In [10]:

def improve_policy(pi, pi1,gamma,valueMap):
    policy_stable = True
    for state in states:
        old = pi[state].copy()
        # Greedify policy for state
        greedify_policy(state,pi,pi1,gamma,valueMap)
        if not np.array_equal(pi[state], old):
            policy_stable = False
    print(pi)
    print(pi1)
    return pi, pi1, policy_stable

In [11]:

def policy_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    pi = np.ones((gridSize,gridSize))/4
    pi1 = np.chararray((gridSize, gridSize))
    pi1[:] = 'a'
    policy_stable = False
    print("here")
    while not policy_stable:
        valueMap = policy_evaluate(states,actions,gamma,valueMap)
        pi,pi1, policy_stable = improve_policy(pi,pi1,  gamma,valueMap)
    return valueMap, pi,pi1

In [12]:

theta=0.1
valueMap, pi,pi1 = policy_iteration(gamma, theta)

[[0. 3. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 1.]
 [0. 0. 2. 0.]]
[[b'u' b'l' b'u' b'u']
 [b'u' b'u' b'u' b'u']
 [b'u' b'u' b'u' b'd']
 [b'u' b'u' b'r' b'u']]
[[0. 3. 3. 0.]
 [0. 0. 0. 1.]
 [0. 0. 1. 1.]
 [0. 2. 2. 0.]]
[[b'u' b'l' b'l' b'u']
 [b'u' b'u' b'u' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'r' b'r' b'u']]
[[0. 3. 3. 1.]
 [0. 0. 1. 1.]
 [0. 0. 1. 1.]
 [0. 2. 2. 0.]]
[[b'u' b'l' b'l' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'r' b'r' b'u']]
[[0. 3. 3. 1.]
 [0. 0. 1. 1.]
 [0. 0. 1. 1.]
 [0. 2. 2. 0.]]
[[b'u' b'l' b'l' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'r' b'r' b'u']]

Findings

From the above valueMap we can see that greedification solves this much faster as below

1c. Bellman Optimality update

The Bellman optimality update directly updates the value state function for the action that results in the maximum return in a state

In [13]:

gamma = 1 # discounting rate
rewardValue = -1
gridSize = 4
terminationStates = [[0,0], [gridSize-1, gridSize-1]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

In [14]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [15]:

def bellman_optimality_update(valueMap, state, gamma):

    q_values=np.zeros(len(actions))
    
    for idx,action in enumerate(actions):
        finalPosition,reward = actionValue(state,action)
        q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
    # Find the index of the action for which the q_value is 
    idx=q_values.argmax()
            
    max = np.argmax(q_values)
    valueMap[state[0],state[1]] = q_values[max]    
    #print(q_values[max])

In [16]:

def value_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    while True:
        delta = 0
        for state in states:
            v_old=valueMap[state[0],state[1]]
            bellman_optimality_update(valueMap, state, gamma)
            delta = max(delta, abs(v_old - valueMap[state[0],state[1]]))
        if delta < theta:
            break
    pi = np.ones((gridSize,gridSize))/4
    for state in states:
        greedify_policy(state,pi,pi1,gamma,valueMap)
    print(pi)
    print(pi1)
    return valueMap, pi,pi1

In [17]:

gamma = 1
theta = 0.01
valueMap,pi,pi1=value_iteration(gamma, theta)
pi
pi1

[[0. 3. 3. 1.]
 [0. 0. 0. 1.]
 [0. 0. 1. 1.]
 [0. 2. 2. 0.]]
[[b'u' b'l' b'l' b'd']
 [b'u' b'u' b'u' b'd']
 [b'u' b'u' b'd' b'd']
 [b'u' b'r' b'r' b'u']]

Out[17]:

chararray([[b'u', b'l', b'l', b'd'],
           [b'u', b'u', b'u', b'd'],
           [b'u', b'u', b'd', b'd'],
           [b'u', b'r', b'r', b'u']], dtype='|S1')

Findings

The above valueMap shows the optimal path from any state

2.Gridworld 2

To make the problem more interesting, I created a 2nd grid which has more interesting structure as shown below <img src=”fig5.png”

The end state is the grey cell. Transitions to the black cells have a negative reward of -10. All other transitions have a reward of -1, while the end state has a reward of 0

In [2]:

##2a. Bellman Update

In [3]:

gamma = 1 # discounting rate
gridSize = 4

terminationStates = [[0,0]]
#terminationStates = [[0,0]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

In [4]:

rewardValue = np.zeros((gridSize,gridSize)) -1
rewardValue[0]=np.array([-1,-10,-10,-10])
rewardValue[2]=np.array([-10,-10,-10,-1])
rewardValue

Out[4]:

array([[ -1., -10., -10., -10.],
       [ -1.,  -1.,  -1.,  -1.],
       [-10., -10., -10.,  -1.],
       [ -1.,  -1.,  -1.,  -1.]])

In [5]:

def actionValue(initialPosition,action):
    if initialPosition in terminationStates:
        finalPosition = initialPosition
        reward=0
    else:
        #Compute final position
        finalPosition = np.array(initialPosition) + np.array(action)
        
        # If the action moves the finalPosition out of the grid, stay in same cell
        if -1 in finalPosition or gridSize in finalPosition:
                finalPosition = initialPosition
                reward= rewardValue[finalPosition[0],finalPosition[1]]
        else:
                reward= rewardValue[finalPosition[0],finalPosition[1]]
    
    #print(finalPosition)
    return finalPosition, reward

In [6]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]

In [7]:

def policy_evaluation(numIterations,gamma,theta,valueMap):
    for i in range(numIterations):
        delta=0
        #print("iterations=",i)
        for state in states:
            weightedRewards=0
            for action in actions:
                finalPosition,reward = actionValue(state,action)
                #print("reward=",reward,"valueMap=",valueMap[finalPosition[0],finalPosition][1])
                weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
            #print(weightedRewards)
            valueMap1[state[0],state[1]]=weightedRewards
            #print("wr=",weightedRewards,"va=",valueMap[state[0],state[1]]) 
            delta =max(delta,abs(weightedRewards-valueMap[state[0],state[1]]))
        valueMap = np.copy(valueMap1)
        #print(valueMap1)
        if(delta < 0.01):
            print(delta)                                                   
            print(valueMap)
            break

In [8]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
policy_evaluation(1000,1,0.0001,valueMap)

0.009862596190146178
[[   0.         -137.28514189 -209.19560831 -239.01378395]
 [-129.2494276  -180.67825796 -220.31626448 -237.86482779]
 [-194.08846546 -213.88769305 -231.5579035  -241.29920147]
 [-217.15664109 -227.25768494 -237.76348718 -241.51200989]]

2b. Greedify

In [9]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [10]:

# Compute the value state function for the Grid
def policy_evaluate(states,actions,gamma,valueMap):
    #print("iterations=",i)
    for state in states:
        weightedRewards=0
        for action in actions:
            finalPosition,reward = actionValue(state,action)
            weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Set the computed weighted rewards to valueMap1
        valueMap1[state[0],state[1]]=weightedRewards
    # Copy to original valueMap
    valueMap = np.copy(valueMap1)
    return(valueMap)

In [11]:

def argmax(q_values):
    idx=np.argmax(q_values)
    return(np.random.choice(np.where(a==a[idx])[0].tolist()))


# Compute the best action in each state
def greedify_policy(state,pi,pi1,gamma,valueMap):  
        q_values=np.zeros(len(actions))
        for idx,action in enumerate(actions):
            finalPosition,reward = actionValue(state,action)
            q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Find the index of the action for which the q_value is 
        idx=q_values.argmax()
        pi[state[0],state[1]]=idx 
        if(idx == 0):
            pi1[state[0],state[1]]='u'
        elif(idx == 1):
            pi1[state[0],state[1]]='d'
        elif(idx == 2):
            pi1[state[0],state[1]]='r'
        elif(idx == 3):
            pi1[state[0],state[1]]='l'

In [12]:

def improve_policy(pi, pi1,gamma,valueMap):
    policy_stable = True
    for state in states:
        old = pi[state].copy()
        # Greedify policy for state
        greedify_policy(state,pi,pi1,gamma,valueMap)
        if not np.array_equal(pi[state], old):
            policy_stable = False
    print(pi)
    print(pi1)
    return pi, pi1, policy_stable

In [13]:

def policy_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    pi = np.ones((gridSize,gridSize))/4
    pi1 = np.chararray((gridSize, gridSize))
    pi1[:] = 'a'
    policy_stable = False
    print("here")
    while not policy_stable:
        valueMap = policy_evaluate(states,actions,gamma,valueMap)
        pi,pi1, policy_stable = improve_policy(pi,pi1,  gamma,valueMap)
    return valueMap, pi,pi1

In [14]:

theta=0.1
valueMap, pi,pi1 = policy_iteration(gamma, theta)

here
[[0. 3. 1. 1.]
 [0. 3. 2. 1.]
 [0. 1. 1. 1.]
 [1. 1. 2. 1.]]
[[b'u' b'l' b'd' b'd']
 [b'u' b'l' b'r' b'd']
 [b'u' b'd' b'd' b'd']
 [b'd' b'd' b'r' b'd']]
[[0. 3. 1. 1.]
 [0. 3. 2. 1.]
 [0. 1. 1. 1.]
 [1. 2. 2. 1.]]
[[b'u' b'l' b'd' b'd']
 [b'u' b'l' b'r' b'd']
 [b'u' b'd' b'd' b'd']
 [b'd' b'r' b'r' b'd']]
[[0. 3. 1. 1.]
 [0. 3. 2. 1.]
 [0. 1. 1. 1.]
 [2. 2. 2. 1.]]
[[b'u' b'l' b'd' b'd']
 [b'u' b'l' b'r' b'd']
 [b'u' b'd' b'd' b'd']
 [b'r' b'r' b'r' b'd']]
[[0. 3. 1. 1.]
 [0. 3. 2. 1.]
 [0. 1. 1. 1.]
 [2. 2. 2. 1.]]
[[b'u' b'l' b'd' b'd']
 [b'u' b'l' b'r' b'd']
 [b'u' b'd' b'd' b'd']
 [b'r' b'r' b'r' b'd']]

In [15]:

## 2c. Bellman Optimality update

In [16]:

gamma = 1 # discounting rate
rewardValue = np.zeros((gridSize,gridSize)) -1
rewardValue[0]=np.array([-1,-10,-10,-10])
rewardValue[2]=np.array([-10,-10,-10,-1])
rewardValue
gridSize = 4
terminationStates = [[0,0]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

In [17]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [18]:

2c. Bellman Optimality Update

def bellman_optimality_update(valueMap, state, gamma):

    q_values=np.zeros(len(actions))
    
    for idx,action in enumerate(actions):
        finalPosition,reward = actionValue(state,action)
        q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
    # Find the index of the action for which the q_value is 
    idx=q_values.argmax()
            
    max = np.argmax(q_values)
    valueMap[state[0],state[1]] = q_values[max]    
    #print(q_values[max])

In [19]:

def value_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    while True:
        delta = 0
        for state in states:
            v_old=valueMap[state[0],state[1]]
            bellman_optimality_update(valueMap, state, gamma)
            delta = max(delta, abs(v_old - valueMap[state[0],state[1]]))
        if delta < theta:
            break
    pi = np.ones((gridSize,gridSize))/4
    for state in states:
        greedify_policy(state,pi,pi1,gamma,valueMap)
    print(pi)
    print(pi1)
    return valueMap, pi,pi1

In [20]:

gamma = 1
theta = 0.000001
valueMap,pi,pi1=value_iteration(gamma, theta)
pi
pi1

[[0. 3. 1. 1.]
 [0. 3. 3. 3.]
 [0. 0. 0. 0.]
 [2. 2. 2. 0.]]
[[b'u' b'l' b'd' b'd']
 [b'u' b'l' b'l' b'l']
 [b'u' b'u' b'u' b'u']
 [b'r' b'r' b'r' b'u']]

Out[20]:

chararray([[b'u', b'l', b'd', b'd'],
           [b'u', b'l', b'l', b'l'],
           [b'u', b'u', b'u', b'u'],
           [b'r', b'r', b'r', b'u']], dtype='|S1')

Findings

The above shows the path from any cell to the stop cell as

3. Another maze

This is the third grid world which I create where the green cell is the end state and has a reward of 0. Transitions to the black cell will receive a reward of -10 and all other transitions will receive a reward of -1

In [2]:

gamma = 1 # discounting rate
gridSize = 5
terminationStates = [[2,2]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

In [3]:

rewardValue = np.zeros((gridSize,gridSize)) -1
rewardValue[1]=np.array([-1,-10,-1,-10,-1])
rewardValue[3]=np.array([-1,-10,-1,-10,-1])
rewardValue

Out[3]:

array([[ -1.,  -1.,  -1.,  -1.,  -1.],
       [ -1., -10.,  -1., -10.,  -1.],
       [ -1.,  -1.,  -1.,  -1.,  -1.],
       [ -1., -10.,  -1., -10.,  -1.],
       [ -1.,  -1.,  -1.,  -1.,  -1.]])

In [4]:

3a. Bellman Update

def actionValue(initialPosition,action):
    if initialPosition in terminationStates:
        finalPosition = initialPosition
        reward=0
    else:
        #Compute final position
        finalPosition = np.array(initialPosition) + np.array(action)
        
        # If the action moves the finalPosition out of the grid, stay in same cell
        if -1 in finalPosition or gridSize in finalPosition:
                finalPosition = initialPosition
                reward= rewardValue[finalPosition[0],finalPosition[1]]
        else:
                reward= rewardValue[finalPosition[0],finalPosition[1]]
    
    #print(finalPosition)
    return finalPosition, reward

In [5]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]

In [6]:

def policy_evaluation(numIterations,gamma,theta,valueMap):
    for i in range(numIterations):
        delta=0
        #print("iterations=",i)
        for state in states:
            weightedRewards=0
            for action in actions:
                finalPosition,reward = actionValue(state,action)
                #print("reward=",reward,"valueMap=",valueMap[finalPosition[0],finalPosition][1])
                weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
            #print(weightedRewards)
            valueMap1[state[0],state[1]]=weightedRewards
            #print("wr=",weightedRewards,"va=",valueMap[state[0],state[1]]) 
            delta =max(delta,abs(weightedRewards-valueMap[state[0],state[1]]))
        valueMap = np.copy(valueMap1)
        #print(valueMap1)
        if(delta < 0.01):
            print(delta)                                                   
            print(valueMap)
            break

In [7]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
policy_evaluation(1000,1,0.0001,valueMap)

0.009697101372182715
[[-82.49768079 -80.51647225 -74.9345659  -80.51647225 -82.49768079]
 [-80.51647225 -71.15241689 -59.80375072 -71.15241689 -80.51647225]
 [-74.9345659  -59.80375072   0.         -59.80375072 -74.9345659 ]
 [-80.51647225 -71.15241689 -59.80375072 -71.15241689 -80.51647225]
 [-82.49768079 -80.51647225 -74.9345659  -80.51647225 -82.49768079]]

3b. Greedify

In [8]:

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [9]:

# Compute the value state function for the Grid
def policy_evaluate(states,actions,gamma,valueMap):
    #print("iterations=",i)
    for state in states:
        weightedRewards=0
        for action in actions:
            finalPosition,reward = actionValue(state,action)
            weightedRewards += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Set the computed weighted rewards to valueMap1
        valueMap1[state[0],state[1]]=weightedRewards
    # Copy to original valueMap
    valueMap = np.copy(valueMap1)
    return(valueMap)

In [10]:

def argmax(q_values):
    idx=np.argmax(q_values)
    return(np.random.choice(np.where(a==a[idx])[0].tolist()))


# Compute the best action in each state
def greedify_policy(state,pi,pi1,gamma,valueMap):  
        q_values=np.zeros(len(actions))
        for idx,action in enumerate(actions):
            finalPosition,reward = actionValue(state,action)
            q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
        # Find the index of the action for which the q_value is 
        idx=q_values.argmax()
        pi[state[0],state[1]]=idx 
        if(idx == 0):
            pi1[state[0],state[1]]='u'
        elif(idx == 1):
            pi1[state[0],state[1]]='d'
        elif(idx == 2):
            pi1[state[0],state[1]]='r'
        elif(idx == 3):
            pi1[state[0],state[1]]='l'

In [11]:

def improve_policy(pi, pi1,gamma,valueMap):
    policy_stable = True
    for state in states:
        old = pi[state].copy()
        # Greedify policy for state
        greedify_policy(state,pi,pi1,gamma,valueMap)
        if not np.array_equal(pi[state], old):
            policy_stable = False
    print(pi)
    print(pi1)
    return pi, pi1, policy_stable

In [12]:

def policy_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    pi = np.ones((gridSize,gridSize))/4
    pi1 = np.chararray((gridSize, gridSize))
    pi1[:] = 'a'
    policy_stable = False
    print("here")
    while not policy_stable:
        valueMap = policy_evaluate(states,actions,gamma,valueMap)
        pi,pi1, policy_stable = improve_policy(pi,pi1,  gamma,valueMap)
    return valueMap, pi,pi1

In [13]:

theta=0.1
valueMap, pi,pi1 = policy_iteration(gamma, theta)

here
[[0. 2. 0. 2. 0.]
 [0. 0. 1. 0. 0.]
 [3. 2. 0. 3. 2.]
 [0. 1. 0. 1. 0.]
 [1. 2. 1. 2. 1.]]
[[b'u' b'r' b'u' b'r' b'u']
 [b'u' b'u' b'd' b'u' b'u']
 [b'l' b'r' b'u' b'l' b'r']
 [b'u' b'd' b'u' b'd' b'u']
 [b'd' b'r' b'd' b'r' b'd']]
[[0. 3. 0. 2. 0.]
 [0. 0. 1. 0. 0.]
 [3. 2. 0. 3. 2.]
 [1. 1. 0. 1. 1.]
 [1. 3. 1. 2. 1.]]
[[b'u' b'l' b'u' b'r' b'u']
 [b'u' b'u' b'd' b'u' b'u']
 [b'l' b'r' b'u' b'l' b'r']
 [b'd' b'd' b'u' b'd' b'd']
 [b'd' b'l' b'd' b'r' b'd']]
[[0. 3. 0. 2. 0.]
 [0. 0. 1. 0. 0.]
 [3. 2. 0. 3. 2.]
 [1. 1. 0. 1. 1.]
 [1. 3. 1. 2. 1.]]
[[b'u' b'l' b'u' b'r' b'u']
 [b'u' b'u' b'd' b'u' b'u']
 [b'l' b'r' b'u' b'l' b'r']
 [b'd' b'd' b'u' b'd' b'd']
 [b'd' b'l' b'd' b'r' b'd']]

In [14]:

gamma = 1 # discounting rate
gridSize=5
rewardValue = np.zeros((gridSize,gridSize)) -1
rewardValue = np.zeros((gridSize,gridSize)) -1
rewardValue[1]=np.array([-1,-10,-1,-10,-1])
rewardValue[3]=np.array([-1,-10,-1,-10,-1])
print(rewardValue)


terminationStates = [[2,2]]
actions = [[-1, 0], [1, 0], [0, 1], [0, -1]]
numIterations = 1000

[[ -1.  -1.  -1.  -1.  -1.]
 [ -1. -10.  -1. -10.  -1.]
 [ -1.  -1.  -1.  -1.  -1.]
 [ -1. -10.  -1. -10.  -1.]
 [ -1.  -1.  -1.  -1.  -1.]]

In [15]:

3c. Bellman Optimality Update

valueMap = np.zeros((gridSize, gridSize))
valueMap1 = np.zeros((gridSize, gridSize))
states = [[i, j] for i in range(gridSize) for j in range(gridSize)]
pi = np.ones((gridSize,gridSize))/4
pi1 = np.chararray((gridSize, gridSize))
pi1[:] = 'a'

In [16]:

def bellman_optimality_update(valueMap, state, gamma):

    q_values=np.zeros(len(actions))
    
    for idx,action in enumerate(actions):
        finalPosition,reward = actionValue(state,action)
        q_values[idx] += 1/4* (reward + gamma * valueMap[finalPosition[0],finalPosition][1])
    # Find the index of the action for which the q_value is 
    idx=q_values.argmax()
            
    max = np.argmax(q_values)
    valueMap[state[0],state[1]] = q_values[max]    
    #print(q_values[max])

In [17]:

def value_iteration(gamma, theta):
    valueMap = np.zeros((gridSize, gridSize))
    while True:
        delta = 0
        for state in states:
            v_old=valueMap[state[0],state[1]]
            bellman_optimality_update(valueMap, state, gamma)
            delta = max(delta, abs(v_old - valueMap[state[0],state[1]]))
        if delta < theta:
            break
    pi = np.ones((gridSize,gridSize))/4
    for state in states:
        greedify_policy(state,pi,pi1,gamma,valueMap)
    print(pi)
    print(pi1)
    return valueMap, pi,pi1

In [18]:

gamma = 1
theta = 0.000001
valueMap,pi,pi1=value_iteration(gamma, theta)
pi
pi1

[[1. 2. 1. 3. 1.]
 [1. 1. 1. 1. 1.]
 [2. 2. 0. 3. 3.]
 [0. 0. 0. 0. 0.]
 [0. 2. 0. 3. 0.]]
[[b'd' b'r' b'd' b'l' b'd']
 [b'd' b'd' b'd' b'd' b'd']
 [b'r' b'r' b'u' b'l' b'l']
 [b'u' b'u' b'u' b'u' b'u']
 [b'u' b'r' b'u' b'l' b'u']]

Out[18]:

chararray([[b'd', b'r', b'd', b'l', b'd'],
           [b'd', b'd', b'd', b'd', b'd'],
           [b'r', b'r', b'u', b'l', b'l'],
           [b'u', b'u', b'u', b'u', b'u'],
           [b'u', b'r', b'u', b'l', b'u']], dtype='|S1')

Findings

We can see that the Bellman Optimality Update correctly finds the path the to end node which we can see from the valueMap1 above which is

Conclusion:

We can see how with the Bellman equations implemented iteratively with dynamic programming we can solve mazes of arbitrary shapes and complexities as long as we correctly choose the reward for the transitions

References
1. Reinforcement learning – An introduction by Richard S. Sutton and Andrew G Barto
2. Reinforcement learning (RL) 101 with Python Blog by Gerard Martinez
3. Reinforcement Learning Demystified: Solving MDPs with Dynamic Programming Blog by Mohammed Ashraf

You may also like

1. My book ‘Deep Learning from first principles:Second Edition’ now on Amazon
2. Big Data-4: Webserver log analysis with RDDs, Pyspark, SparkR and SparklyR
3. Practical Machine Learning with R and Python – Part 3
3. Pitching yorkpy…on the middle and outside off-stump to IPL – Part 2
4. Sixer – R package cricketr’s new Shiny avatar
5. Natural language processing: What would Shakespeare say?
6. Getting started with Tensorflow, Keras in Python and R

To see all posts click Index of posts

Getting started with Tensorflow, Keras in Python and R

The Pale Blue Dot

“From this distant vantage point, the Earth might not seem of any particular interest. But for us, it’s different. Consider again that dot. That’s here, that’s home, that’s us. On it everyone you love, everyone you know, everyone you ever heard of, every human being who ever was, lived out their lives. The aggregate of our joy and suffering, thousands of confident religions, ideologies, and economic doctrines, every hunter and forager, every hero and coward, every creator and destroyer of civilization, every king and peasant, every young couple in love, every mother and father, hopeful child, inventor and explorer, every teacher of morals, every corrupt politician, every “superstar,” every “supreme leader,” every saint and sinner in the history of our species lived there—on the mote of dust suspended in a sunbeam.”

Carl Sagan

Tensorflow and Keras are Deep Learning frameworks that really simplify a lot of things to the user. If you are familiar with Machine Learning and Deep Learning concepts then Tensorflow and Keras are really a playground to realize your ideas. In this post I show how you can get started with Tensorflow in both Python and R

Tensorflow in Python

For tensorflow in Python, I found Google’s Colab an ideal environment for running your Deep Learning code. This is an Google’s research project where you can execute your code on GPUs, TPUs etc

Tensorflow in R (RStudio)

To execute tensorflow in R (RStudio) you need to install tensorflow and keras as shown below
In this post I show how to get started with Tensorflow and Keras in R.

# Install Tensorflow in RStudio
#install_tensorflow()
# Install Keras
#install_packages("keras")
library(tensorflow)
libary(keras)

This post takes 3 different Machine Learning problems and uses the
Tensorflow/Keras framework to solve it

Note:
You can view the Google Colab notebook at Tensorflow in Python
The RMarkdown file has been published at RPubs and can be accessed
at Getting started with Tensorflow in R

Checkout my book ‘Deep Learning from first principles: Second Edition – In vectorized Python, R and Octave’. My book starts with the implementation of a simple 2-layer Neural Network and works its way to a generic L-Layer Deep Learning Network, with all the bells and whistles. The derivations have been discussed in detail. The code has been extensively commented and included in its entirety in the Appendix sections. My book is available on Amazon as paperback ($14.99) and in kindle version($9.99/Rs449).

1. Multivariate regression with Tensorflow – Python

This code performs multivariate regression using Tensorflow and keras on the advent of Parkinson disease through sound recordings see Parkinson Speech Dataset with Multiple Types of Sound Recordings Data Set . The clinician’s motorUPDRS score has to be predicted from the set of features

In [0]:

# Import tensorflow
import tensorflow as tf
from tensorflow import keras

In [2]:

#Get the data rom the UCI Machine Learning repository
dataset = keras.utils.get_file("parkinsons_updrs.data", "https://archive.ics.uci.edu/ml/machine-learning-databases/parkinsons/telemonitoring/parkinsons_updrs.data")

Downloading data from https://archive.ics.uci.edu/ml/machine-learning-databases/parkinsons/telemonitoring/parkinsons_updrs.data
917504/911261 [==============================] - 0s 0us/step

In [3]:

# Read the CSV file 
import pandas as pd
parkinsons = pd.read_csv(dataset, na_values = "?", comment='\t',
                      sep=",", skipinitialspace=True)
print(parkinsons.shape)
print(parkinsons.columns)
#Check if there are any NAs in the rows
parkinsons.isna().sum()

(5875, 22)
Index(['subject#', 'age', 'sex', 'test_time', 'motor_UPDRS', 'total_UPDRS',
       'Jitter(%)', 'Jitter(Abs)', 'Jitter:RAP', 'Jitter:PPQ5', 'Jitter:DDP',
       'Shimmer', 'Shimmer(dB)', 'Shimmer:APQ3', 'Shimmer:APQ5',
       'Shimmer:APQ11', 'Shimmer:DDA', 'NHR', 'HNR', 'RPDE', 'DFA', 'PPE'],
      dtype='object')

Out[3]:

subject#         0
age              0
sex              0
test_time        0
motor_UPDRS      0
total_UPDRS      0
Jitter(%)        0
Jitter(Abs)      0
Jitter:RAP       0
Jitter:PPQ5      0
Jitter:DDP       0
Shimmer          0
Shimmer(dB)      0
Shimmer:APQ3     0
Shimmer:APQ5     0
Shimmer:APQ11    0
Shimmer:DDA      0
NHR              0
HNR              0
RPDE             0
DFA              0
PPE              0
dtype: int64

Note: To see how to create dummy variables see my post Practical Machine Learning with R and Python – Part 2

In [4]:

# Drop the columns subject number as it is not relevant
parkinsons1=parkinsons.drop(['subject#'],axis=1)

# Create dummy variables for sex (M/F)
parkinsons2=pd.get_dummies(parkinsons1,columns=['sex'])
parkinsons2.head()

Out[4]
age test_time motor_UPDRS total_UPDRS Jitter(%) Jitter(Abs) Jitter:RAP Jitter:PPQ5 Jitter:DDP Shimmer Shimmer(dB) Shimmer:APQ3 Shimmer:APQ5 Shimmer:APQ11 Shimmer:DDA NHR HNR RPDE DFA PPE sex_0 sex_1
0 72 5.6431 28.199 34.398 0.00662 0.000034 0.00401 0.00317 0.01204 0.02565 0.230 0.01438 0.01309 0.01662 0.04314 0.014290 21.640 0.41888 0.54842 0.16006 1 0
1 72 12.6660 28.447 34.894 0.00300 0.000017 0.00132 0.00150 0.00395 0.02024 0.179 0.00994 0.01072 0.01689 0.02982 0.011112 27.183 0.43493 0.56477 0.10810 1 0
2 72 19.6810 28.695 35.389 0.00481 0.000025 0.00205 0.00208 0.00616 0.01675 0.181 0.00734 0.00844 0.01458 0.02202 0.020220 23.047 0.46222 0.54405 0.21014 1 0
3 72 25.6470 28.905 35.810 0.00528 0.000027 0.00191 0.00264 0.00573 0.02309 0.327 0.01106 0.01265 0.01963 0.03317 0.027837 24.445 0.48730 0.57794 0.33277 1 0
4 72 33.6420 29.187 36.375 0.00335 0.000020 0.00093 0.00130 0.00278 0.01703 0.176 0.00679 0.00929 0.01819 0.02036 0.011625 26.126 0.47188 0.56122 0.19361 1 0

# Create a training and test data set with 80%/20%
train_dataset = parkinsons2.sample(frac=0.8,random_state=0)
test_dataset = parkinsons2.drop(train_dataset.index)

# Select columns
train_dataset1= train_dataset[['age', 'test_time', 'Jitter(%)', 'Jitter(Abs)',
       'Jitter:RAP', 'Jitter:PPQ5', 'Jitter:DDP', 'Shimmer', 'Shimmer(dB)',
       'Shimmer:APQ3', 'Shimmer:APQ5', 'Shimmer:APQ11', 'Shimmer:DDA', 'NHR',
       'HNR', 'RPDE', 'DFA', 'PPE', 'sex_0', 'sex_1']]
test_dataset1= test_dataset[['age','test_time', 'Jitter(%)', 'Jitter(Abs)',
       'Jitter:RAP', 'Jitter:PPQ5', 'Jitter:DDP', 'Shimmer', 'Shimmer(dB)',
       'Shimmer:APQ3', 'Shimmer:APQ5', 'Shimmer:APQ11', 'Shimmer:DDA', 'NHR',
       'HNR', 'RPDE', 'DFA', 'PPE', 'sex_0', 'sex_1']]

In [7]:

# Generate the statistics of the columns for use in normalization of the data
train_stats = train_dataset1.describe()
train_stats = train_stats.transpose()
train_stats

Out[7]:

count mean std min 25% 50% 75% max
age 4700.0 64.792766 8.870401 36.000000 58.000000 65.000000 72.000000 85.000000
test_time 4700.0 93.399490 53.630411 -4.262500 46.852250 93.405000 139.367500 215.490000
Jitter(%) 4700.0 0.006136 0.005612 0.000830 0.003560 0.004900 0.006770 0.099990
Jitter(Abs) 4700.0 0.000044 0.000036 0.000002 0.000022 0.000034 0.000053 0.000396
Jitter:RAP 4700.0 0.002969 0.003089 0.000330 0.001570 0.002235 0.003260 0.057540
Jitter:PPQ5 4700.0 0.003271 0.003760 0.000430 0.001810 0.002480 0.003460 0.069560
Jitter:DDP 4700.0 0.008908 0.009267 0.000980 0.004710 0.006705 0.009790 0.172630
Shimmer 4700.0 0.033992 0.025922 0.003060 0.019020 0.027385 0.039810 0.268630
Shimmer(dB) 4700.0 0.310487 0.231016 0.026000 0.175000 0.251000 0.363250 2.107000
Shimmer:APQ3 4700.0 0.017125 0.013275 0.001610 0.009190 0.013615 0.020562 0.162670
Shimmer:APQ5 4700.0 0.020151 0.016848 0.001940 0.010750 0.015785 0.023733 0.167020
Shimmer:APQ11 4700.0 0.027508 0.020270 0.002490 0.015630 0.022685 0.032713 0.275460
Shimmer:DDA 4700.0 0.051375 0.039826 0.004840 0.027567 0.040845 0.061683 0.488020
NHR 4700.0 0.032116 0.060206 0.000304 0.010827 0.018403 0.031452 0.748260
HNR 4700.0 21.704631 4.288853 1.659000 19.447750 21.973000 24.445250 37.187000
RPDE 4700.0 0.542549 0.100212 0.151020 0.471235 0.543490 0.614335 0.966080
DFA 4700.0 0.653015 0.070446 0.514040 0.596470 0.643285 0.710618 0.865600
PPE 4700.0 0.219559 0.091506 0.021983 0.156470 0.205340 0.264017 0.731730
sex_0 4700.0 0.681489 0.465948 0.000000 0.000000 1.000000 1.000000 1.000000
sex_1 4700.0 0.318511 0.465948 0.000000 0.000000 0.000000 1.000000 1.000000

In [0]:

# Create the target variable
train_labels = train_dataset.pop('motor_UPDRS')
test_labels = test_dataset.pop('motor_UPDRS')

In [0]:

# Normalize the data by subtracting the mean and dividing by the standard deviation
def normalize(x):
  return (x - train_stats['mean']) / train_stats['std']

# Create normalized training and test data
normalized_train_data = normalize(train_dataset1)
normalized_test_data = normalize(test_dataset1)

In [0]:

# Create a Deep Learning model with keras
model = tf.keras.Sequential([
    keras.layers.Dense(6, activation=tf.nn.relu, input_shape=[len(train_dataset1.keys())]),
    keras.layers.Dense(9, activation=tf.nn.relu),
    keras.layers.Dense(6,activation=tf.nn.relu),
    keras.layers.Dense(1)
  ])

# Use the Adam optimizer with a learning rate of 0.01
optimizer=keras.optimizers.Adam(lr=.01, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.0, amsgrad=False)

# Set the metrics required to be Mean Absolute Error and Mean Squared Error.For regression, the loss is mean_squared_error
model.compile(loss='mean_squared_error',
                optimizer=optimizer,
                metrics=['mean_absolute_error', 'mean_squared_error'])

In [0]:

# Create a model
history=model.fit(
  normalized_train_data, train_labels,
  epochs=1000, validation_data = (normalized_test_data,test_labels), verbose=0)

In [26]:

hist = pd.DataFrame(history.history)
hist['epoch'] = history.epoch
hist.tail()

Out[26]:

	loss	mean_absolute_error	mean_squared_error	val_loss	val_mean_absolute_error	val_mean_squared_error	epoch
995	15.773989	2.936990	15.773988	16.980803	3.028168	16.980803	995
996	15.238623	2.873420	15.238622	17.458752	3.101033	17.458752	996
997	15.437594	2.895500	15.437593	16.926016	2.971508	16.926018	997
998	15.867891	2.943521	15.867892	16.950249	2.985036	16.950249	998
999	15.846878	2.938914	15.846880	17.095623	3.014504	17.095625	999

In [30]:

def plot_history(history):
  hist = pd.DataFrame(history.history)
  hist['epoch'] = history.epoch

  plt.figure()
  plt.xlabel('Epoch')
  plt.ylabel('Mean Abs Error')
  plt.plot(hist['epoch'], hist['mean_absolute_error'],
           label='Train Error')
  plt.plot(hist['epoch'], hist['val_mean_absolute_error'],
           label = 'Val Error')
  plt.ylim([2,5])
  plt.legend()

  plt.figure()
  plt.xlabel('Epoch')
  plt.ylabel('Mean Square Error ')
  plt.plot(hist['epoch'], hist['mean_squared_error'],
           label='Train Error')
  plt.plot(hist['epoch'], hist['val_mean_squared_error'],
           label = 'Val Error')
  plt.ylim([10,40])
  plt.legend()
  plt.show()


plot_history(history)

Observation

It can be seen that the mean absolute error is on an average about +/- 4.0. The validation error also is about the same. This can be reduced by playing around with the hyperparamaters and increasing the number of iterations

1a. Multivariate Regression in Tensorflow – R

# Install Tensorflow in RStudio
#install_tensorflow()
# Install Keras
#install_packages("keras")
library(tensorflow)

library(keras)

library(dplyr)

library(dummies)

## dummies-1.5.6 provided by Decision Patterns

library(tensorflow)
library(keras)

Multivariate regression

This code performs multivariate regression using Tensorflow and keras on the advent of Parkinson disease through sound recordings see Parkinson Speech Dataset with Multiple Types of Sound Recordings Data Set. The clinician’s motorUPDRS score has to be predicted from the set of features.

Read the data

# Download the Parkinson's data from UCI Machine Learning repository
dataset <- read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/parkinsons/telemonitoring/parkinsons_updrs.data")

# Set the column names
names(dataset) <- c("subject","age", "sex", "test_time","motor_UPDRS","total_UPDRS","Jitter","Jitter.Abs",
                 "Jitter.RAP","Jitter.PPQ5","Jitter.DDP","Shimmer", "Shimmer.dB", "Shimmer.APQ3",
                 "Shimmer.APQ5","Shimmer.APQ11","Shimmer.DDA", "NHR","HNR", "RPDE", "DFA","PPE")

# Remove the column 'subject' as it is not relevant to analysis
dataset1 <- subset(dataset, select = -c(subject))

# Make the column 'sex' as a factor for using dummies
dataset1$sex=as.factor(dataset1$sex)
# Add dummy variables for categorical cariable 'sex'
dataset2 <- dummy.data.frame(dataset1, sep = ".")

## Warning in model.matrix.default(~x - 1, model.frame(~x - 1), contrasts =
## FALSE): non-list contrasts argument ignored

dataset3 <- na.omit(dataset2)

Split the data as training and test in 80/20

## Split data 80% training and 20% test
sample_size <- floor(0.8 * nrow(dataset3))

## set the seed to make your partition reproducible
set.seed(12)
train_index <- sample(seq_len(nrow(dataset3)), size = sample_size)

train_dataset <- dataset3[train_index, ]
test_dataset <- dataset3[-train_index, ]

train_data <- train_dataset %>% select(sex.0,sex.1,age, test_time,Jitter,Jitter.Abs,Jitter.PPQ5,Jitter.DDP,
                              Shimmer, Shimmer.dB,Shimmer.APQ3,Shimmer.APQ11,
                              Shimmer.DDA,NHR,HNR,RPDE,DFA,PPE)

train_labels <- select(train_dataset,motor_UPDRS)
test_data <- test_dataset %>% select(sex.0,sex.1,age, test_time,Jitter,Jitter.Abs,Jitter.PPQ5,Jitter.DDP,
                              Shimmer, Shimmer.dB,Shimmer.APQ3,Shimmer.APQ11,
                              Shimmer.DDA,NHR,HNR,RPDE,DFA,PPE)
test_labels <- select(test_dataset,motor_UPDRS)

Normalize the data

 # Normalize the data by subtracting the mean and dividing by the standard deviation
normalize<-function(x) {
  y<-(x - mean(x)) / sd(x)
  return(y)
}

normalized_train_data <-apply(train_data,2,normalize)
# Convert to matrix
train_labels <- as.matrix(train_labels)
normalized_test_data <- apply(test_data,2,normalize)
test_labels <- as.matrix(test_labels)

Create the Deep Learning Model

model <- keras_model_sequential()
model %>% 
  layer_dense(units = 6, activation = 'relu', input_shape = dim(normalized_train_data)[2]) %>% 
  layer_dense(units = 9, activation = 'relu') %>%
  layer_dense(units = 6, activation = 'relu') %>%
  layer_dense(units = 1)

# Set the metrics required to be Mean Absolute Error and Mean Squared Error.For regression, the loss is 
# mean_squared_error
model %>% compile(
  loss = 'mean_squared_error',
  optimizer = optimizer_rmsprop(),
  metrics = c('mean_absolute_error','mean_squared_error')
)

# Fit the model
# Use the test data for validation
history <- model %>% fit(
  normalized_train_data, train_labels, 
  epochs = 30, batch_size = 128, 
  validation_data = list(normalized_test_data,test_labels)
)

Plot mean squared error, mean absolute error and loss for training data and test data

plot(history)

Fig1

2. Binary classification in Tensorflow – Python

This is a simple binary classification problem from UCI Machine Learning repository and deals with data on Breast cancer from the Univ. of Wisconsin Breast Cancer Wisconsin (Diagnostic) Data Set bold text

In [31]:

import tensorflow as tf
from tensorflow import keras
import pandas as pd
# Read the data set from UCI ML site
dataset_path = keras.utils.get_file("breast-cancer-wisconsin.data", "https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data")
raw_dataset = pd.read_csv(dataset_path, sep=",", na_values = "?", skipinitialspace=True,)
dataset = raw_dataset.copy()

#Check for Null and drop
dataset.isna().sum()
dataset = dataset.dropna()
dataset.isna().sum()

# Set the column names
dataset.columns = ["id","thickness",	"cellsize",	"cellshape","adhesion","epicellsize",
                    "barenuclei","chromatin","normalnucleoli","mitoses","class"]
dataset.head()

Downloading data from https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data
24576/19889 [=====================================] - 0s 1us/step
id	thickness	cellsize	cellshape	adhesion	epicellsize	barenuclei	chromatin	normalnucleoli	mitoses	class
0	1002945	5	4	4	5	7	10.0	3	2	1	2
1	1015425	3	1	1	1	2	2.0	3	1	1	2
2	1016277	6	8	8	1	3	4.0	3	7	1	2
3	1017023	4	1	1	3	2	1.0	3	1	1	2
4	1017122	8	10	10	8	7	10.0	9	7	1	4

# Create a training/test set in the ratio 80/20
train_dataset = dataset.sample(frac=0.8,random_state=0)
test_dataset = dataset.drop(train_dataset.index)

# Set the training and test set
train_dataset1= train_dataset[['thickness','cellsize','cellshape','adhesion',
                'epicellsize', 'barenuclei', 'chromatin', 'normalnucleoli','mitoses']]
test_dataset1=test_dataset[['thickness','cellsize','cellshape','adhesion',
                'epicellsize', 'barenuclei', 'chromatin', 'normalnucleoli','mitoses']]

In [34]:

# Generate the stats for each column to be used for normalization
train_stats = train_dataset1.describe()
train_stats = train_stats.transpose()
train_stats

Out[34]:

	count	mean	std	min	25%	50%	75%	max
thickness	546.0	4.430403	2.812768	1.0	2.0	4.0	6.0	10.0
cellsize	546.0	3.179487	3.083668	1.0	1.0	1.0	5.0	10.0
cellshape	546.0	3.225275	3.005588	1.0	1.0	1.0	5.0	10.0
adhesion	546.0	2.921245	2.937144	1.0	1.0	1.0	4.0	10.0
epicellsize	546.0	3.261905	2.252643	1.0	2.0	2.0	4.0	10.0
barenuclei	546.0	3.560440	3.651946	1.0	1.0	1.0	7.0	10.0
chromatin	546.0	3.483516	2.492687	1.0	2.0	3.0	5.0	10.0
normalnucleoli	546.0	2.875458	3.064305	1.0	1.0	1.0	4.0	10.0
mitoses	546.0	1.609890	1.736762	1.0	1.0	1.0	1.0	10.0

In [0]:

# Create target variables
train_labels = train_dataset.pop('class')
test_labels = test_dataset.pop('class')

In [0]:

# Set the target variables as 0 or 1
train_labels[train_labels==2] =0 # benign
train_labels[train_labels==4] =1 # malignant

test_labels[test_labels==2] =0 # benign
test_labels[test_labels==4] =1 # malignant

In [0]:

# Normalize by subtracting mean and dividing by standard deviation
def normalize(x):
  return (x - train_stats['mean']) / train_stats['std']

# Convert columns to numeric
train_dataset1 = train_dataset1.apply(pd.to_numeric)
test_dataset1 = test_dataset1.apply(pd.to_numeric)

# Normalize
normalized_train_data = normalize(train_dataset1)
normalized_test_data = normalize(test_dataset1)

In [0]:

# Create a model
model = tf.keras.Sequential([
    keras.layers.Dense(6, activation=tf.nn.relu, input_shape=[len(train_dataset1.keys())]),
    keras.layers.Dense(9, activation=tf.nn.relu),
    keras.layers.Dense(6,activation=tf.nn.relu),
    keras.layers.Dense(1)
  ])

# Use the RMSProp optimizer
optimizer = tf.keras.optimizers.RMSprop(0.01)

# Since this is binary classification use binary_crossentropy
model.compile(loss='binary_crossentropy',
                optimizer=optimizer,
                metrics=['acc'])


# Fit a model
history=model.fit(
  normalized_train_data, train_labels,
  epochs=1000, validation_data=(normalized_test_data,test_labels), verbose=0)

In [55]:

hist = pd.DataFrame(history.history)
hist['epoch'] = history.epoch
hist.tail()

	loss	acc	val_loss	val_acc	epoch
995	0.112499	0.992674	0.454739	0.970588	995
996	0.112499	0.992674	0.454739	0.970588	996
997	0.112499	0.992674	0.454739	0.970588	997
998	0.112499	0.992674	0.454739	0.970588	998
999	0.112499	0.992674	0.454739	0.970588	999

In [58]:

# Plot training and test accuracy 
plt.plot(history.history['acc'])
plt.plot(history.history['val_acc'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.ylim([0.9,1])
plt.show()












# Plot training and test loss
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.ylim([0,0.5])
plt.show()

2a. Binary classification in Tensorflow -R

This is a simple binary classification problem from UCI Machine Learning repository and deals with data on Breast cancer from the Univ. of Wisconsin Breast Cancer Wisconsin (Diagnostic) Data Set

# Read the data for Breast cancer (Wisconsin)
dataset <- read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data")

# Rename the columns
names(dataset) <- c("id","thickness",   "cellsize", "cellshape","adhesion","epicellsize",
                    "barenuclei","chromatin","normalnucleoli","mitoses","class")

# Remove the columns id and class
dataset1 <- subset(dataset, select = -c(id, class))
dataset2 <- na.omit(dataset1)

# Convert the column to numeric
dataset2$barenuclei <- as.numeric(dataset2$barenuclei)

Normalize the data

train_data <-apply(dataset2,2,normalize)
train_labels <- as.matrix(select(dataset,class))

# Set the target variables as 0 or 1 as it binary classification
train_labels[train_labels==2,]=0
train_labels[train_labels==4,]=1

Create the Deep Learning model

model <- keras_model_sequential()
model %>% 
  layer_dense(units = 6, activation = 'relu', input_shape = dim(train_data)[2]) %>% 
  layer_dense(units = 9, activation = 'relu') %>%
  layer_dense(units = 6, activation = 'relu') %>%
  layer_dense(units = 1)

# Since this is a binary classification we use binary cross entropy
model %>% compile(
  loss = 'binary_crossentropy',
  optimizer = optimizer_rmsprop(),
  metrics = c('accuracy')  # Metrics is accuracy
)

Fit the model. Use 20% of data for validation

history <- model %>% fit(
  train_data, train_labels, 
  epochs = 30, batch_size = 128, 
  validation_split = 0.2
)

Plot the accuracy and loss for training and validation data

plot(history)

3. MNIST in Tensorflow – Python

This takes the famous MNIST handwritten digits . It ca be seen that Tensorflow and Keras make short work of this famous problem of the late 1980s

# Download MNIST data
mnist=tf.keras.datasets.mnist
# Set training and test data and labels
(training_images,training_labels),(test_images,test_labels)=mnist.load_data()

print(training_images.shape)
print(test_images.shape)

(60000, 28, 28)
(10000, 28, 28)

In [61]:

# Plot a sample image from MNIST and show contents
import matplotlib.pyplot as plt
plt.imshow(training_images[1])
print(training_images[1])
[[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 51 159 253
159 50 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 48 238 252 252
252 237 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 54 227 253 252 239
233 252 57 6 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 10 60 224 252 253 252 202
84 252 253 122 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 163 252 252 252 253 252 252
96 189 253 167 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 51 238 253 253 190 114 253 228
47 79 255 168 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 48 238 252 252 179 12 75 121 21
0 0 253 243 50 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 38 165 253 233 208 84 0 0 0 0
0 0 253 252 165 0 0 0 0 0]
[ 0 0 0 0 0 0 0 7 178 252 240 71 19 28 0 0 0 0
0 0 253 252 195 0 0 0 0 0]
[ 0 0 0 0 0 0 0 57 252 252 63 0 0 0 0 0 0 0
0 0 253 252 195 0 0 0 0 0]
[ 0 0 0 0 0 0 0 198 253 190 0 0 0 0 0 0 0 0
0 0 255 253 196 0 0 0 0 0]
[ 0 0 0 0 0 0 76 246 252 112 0 0 0 0 0 0 0 0
0 0 253 252 148 0 0 0 0 0]
[ 0 0 0 0 0 0 85 252 230 25 0 0 0 0 0 0 0 0
7 135 253 186 12 0 0 0 0 0]
[ 0 0 0 0 0 0 85 252 223 0 0 0 0 0 0 0 0 7
131 252 225 71 0 0 0 0 0 0]
[ 0 0 0 0 0 0 85 252 145 0 0 0 0 0 0 0 48 165
252 173 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 86 253 225 0 0 0 0 0 0 114 238 253
162 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 85 252 249 146 48 29 85 178 225 253 223 167
56 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 85 252 252 252 229 215 252 252 252 196 130 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 28 199 252 252 253 252 252 233 145 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 25 128 252 253 252 141 37 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0]]

# Normalize the images by dividing by 255.0
training_images = training_images/255.0
test_images = test_images/255.0

# Create a Sequential Keras model
model = tf.keras.models.Sequential([tf.keras.layers.Flatten(),
                                   tf.keras.layers.Dense(1024,activation=tf.nn.relu),
                                   tf.keras.layers.Dense(10,activation=tf.nn.softmax)])
model.compile(optimizer='adam',loss='sparse_categorical_crossentropy',metrics=['accuracy'])

In [68]:

history=model.fit(training_images,training_labels,validation_data=(test_images, test_labels), epochs=5, verbose=1)

Train on 60000 samples, validate on 10000 samples
Epoch 1/5
60000/60000 [==============================] - 17s 291us/sample - loss: 0.0020 - acc: 0.9999 - val_loss: 0.0719 - val_acc: 0.9810
Epoch 2/5
60000/60000 [==============================] - 17s 284us/sample - loss: 0.0021 - acc: 0.9998 - val_loss: 0.0705 - val_acc: 0.9821
Epoch 3/5
60000/60000 [==============================] - 17s 286us/sample - loss: 0.0017 - acc: 0.9999 - val_loss: 0.0729 - val_acc: 0.9805
Epoch 4/5
60000/60000 [==============================] - 17s 284us/sample - loss: 0.0014 - acc: 0.9999 - val_loss: 0.0762 - val_acc: 0.9804
Epoch 5/5
60000/60000 [==============================] - 17s 280us/sample - loss: 0.0015 - acc: 0.9999 - val_loss: 0.0735 - val_acc: 0.9812

Fig 1

Fig 2

MNIST in Tensorflow – R

The following code uses Tensorflow to learn MNIST’s handwritten digits ### Load MNIST data

mnist <- dataset_mnist()
x_train <- mnist$train$x
y_train <- mnist$train$y
x_test <- mnist$test$x
y_test <- mnist$test$y

Reshape and rescale

# Reshape the array
x_train <- array_reshape(x_train, c(nrow(x_train), 784))
x_test <- array_reshape(x_test, c(nrow(x_test), 784))
# Rescale
x_train <- x_train / 255
x_test <- x_test / 255

Convert out put to One Hot encoded format

y_train <- to_categorical(y_train, 10)
y_test <- to_categorical(y_test, 10)

Fit the model

Use the softmax activation for recognizing 10 digits and categorical cross entropy for loss

model <- keras_model_sequential() 
model %>% 
  layer_dense(units = 256, activation = 'relu', input_shape = c(784)) %>% 
  layer_dense(units = 128, activation = 'relu') %>%
  layer_dense(units = 10, activation = 'softmax') # Use softmax

model %>% compile(
  loss = 'categorical_crossentropy',
  optimizer = optimizer_rmsprop(),
  metrics = c('accuracy')
)

Fit the model

Note: A smaller number of epochs has been used. For better performance increase number of epochs

history <- model %>% fit(
  x_train, y_train, 
  epochs = 5, batch_size = 128, 
  validation_data = list(x_test,y_test)
)

Plot the accuracy and loss for training and test data

plot(history)

Conclusion
This post shows how to use Tensorflow and Keras in both Python & R
Hope you have fun with Tensorflow!!

To see all posts click Index of posts

My presentations on ‘Elements of Neural Networks & Deep Learning’ -Parts 6,7,8

This is the final set of presentations in my series ‘Elements of Neural Networks and Deep Learning’. This set follows the earlier 2 sets of presentations namely
1. My presentations on ‘Elements of Neural Networks & Deep Learning’ -Part1,2,3
2. My presentations on ‘Elements of Neural Networks & Deep Learning’ -Parts 4,5

In this final set of presentations I discuss initialization methods, regularization techniques including dropout. Next I also discuss gradient descent optimization methods like momentum, rmsprop, adam etc. Lastly, I briefly also touch on hyper-parameter tuning approaches. The corresponding implementations are available in vectorized R, Python and Octave are available in my book ‘Deep Learning from first principles:Second edition- In vectorized Python, R and Octave‘

1. Elements of Neural Networks and Deep Learning – Part 6
This part discusses initialization methods specifically like He and Xavier. The presentation also focuses on how to prevent over-fitting using regularization. Lastly the dropout method of regularization is also discusses

The corresponding implementations in vectorized R, Python and Octave of the above discussed methods are available in my post Deep Learning from first principles in Python, R and Octave – Part 6

2. Elements of Neural Networks and Deep Learning – Part 7
This presentation introduces exponentially weighted moving average and shows how this is used in different approaches to gradient descent optimization. The key techniques discussed are learning rate decay, momentum method, rmsprop and adam.

The equivalent implementations of the gradient descent optimization techniques in R, Python and Octave can be seen in my post Deep Learning from first principles in Python, R and Octave – Part 7

3. Elements of Neural Networks and Deep Learning – Part 8
This last part touches upon hyper-parameter tuning in Deep Learning networks

This concludes this series of presentations on “Elements of Neural Networks and Deep Learning’

Important note: Do check out my later version of these videos at Take 4+: Presentations on ‘Elements of Neural Networks and Deep Learning’ – Parts 1-8 . These have more content and also include some corrections. Check it out!

To see all posts click Index of posts

Deep Learning from first principles in Python, R and Octave – Part 4

In this 4th post of my series on Deep Learning from first principles in Python, R and Octave – Part 4, I explore the details of creating a multi-class classifier using the Softmax activation unit in a neural network. The earlier posts in this series were

1. Deep Learning from first principles in Python, R and Octave – Part 1. In this post I implemented logistic regression as a simple Neural Network in vectorized Python, R and Octave
2. Deep Learning from first principles in Python, R and Octave – Part 2. This 2nd part implemented the most elementary neural network with 1 hidden layer and any number of activation units in the hidden layer with sigmoid activation at the output layer
3. Deep Learning from first principles in Python, R and Octave – Part 3. The 3rd implemented a multi-layer Deep Learning network with an arbitrary number if hidden layers and activation units per hidden layer. The output layer was for binary classification which was based on the sigmoid unit. This multi-layer deep network was implemented in vectorized Python, R and Octave.

You can download the PDF version of this book from Github at https://github.com/tvganesh/DeepLearningBook-2ndEd

This 4th part takes a swing at multi-class classification and uses the Softmax as the activation unit in the output layer. Inclusion of the Softmax activation unit in the activation layer requires us to compute the derivative of Softmax, or rather the “Jacobian” of the Softmax function, besides also computing the log loss for this Softmax activation during back propagation. Since the derivation of the Jacobian of a Softmax and the computation of the Cross Entropy/log loss is very involved, I have implemented a basic neural network with just 1 hidden layer with the Softmax activation at the output layer. I also perform multi-class classification based on the ‘spiral’ data set from CS231n Convolutional Neural Networks Stanford course, to test the performance and correctness of the implementations in Python, R and Octave. You can clone download the code for the Python, R and Octave implementations from Github at Deep Learning – Part 4

Note: A detailed discussion of the derivation below can also be seen in my video presentation Neural Networks 5

The Softmax function takes an N dimensional vector as input and generates a N dimensional vector as output.
The Softmax function is given by
$S_{j}= \frac{e_{j}}{\sum_{i}^{N}e_{k}}$
There is a probabilistic interpretation of the Softmax, since the sum of the Softmax values of a set of vectors will always add up to 1, given that each Softmax value is divided by the total of all values.

As mentioned earlier, the Softmax takes a vector input and returns a vector of outputs. For e.g. the Softmax of a vector a=[1, 3, 6] is another vector S=[0.0063,0.0471,0.9464]. Notice that vector output is proportional to the input vector. Also, taking the derivative of a vector by another vector, is known as the Jacobian. By the way, The Matrix Calculus You Need For Deep Learning by Terence Parr and Jeremy Howard, is very good paper that distills all the main mathematical concepts for Deep Learning in one place.

Let us take a simple 2 layered neural network with just 2 activation units in the hidden layer is shown below

$Z_{1}^{1} =W_{11}^{1}x_{1} + W_{21}^{1}x_{2} + b_{1}^{1}$
$Z_{2}^{1} =W_{12}^{1}x_{1} + W_{22}^{1}x_{2} + b_{2}^{1}$
and
$A_{1}^{1} = g'(Z_{1}^{1})$
$A_{2}^{1} = g'(Z_{2}^{1})$
where g'() is the activation unit in the hidden layer which can be a relu, sigmoid or a
tanh function

Note: The superscript denotes the layer. The above denotes the equation for layer 1
of the neural network. For layer 2 with the Softmax activation, the equations are
$Z_{1}^{2} =W_{11}^{2}x_{1} + W_{21}^{2}x_{2} + b_{1}^{2}$
$Z_{2}^{2} =W_{12}^{2}x_{1} + W_{22}^{2}x_{2} + b_{2}^{2}$
and
$A_{1}^{2} = S(Z_{1}^{2})$
$A_{2}^{2} = S(Z_{2}^{2})$
where S() is the Softmax activation function
$S=\begin{pmatrix} S(Z_{1}^{2})\\ S(Z_{2}^{2}) \end{pmatrix}$
$S=\begin{pmatrix} \frac{e^{Z1}}{e^{Z1}+e^{Z2}}\\ \frac{e^{Z2}}{e^{Z1}+e^{Z2}} \end{pmatrix}$

The Jacobian of the softmax ‘S’ is given by
$\begin{pmatrix} \frac {\partial S_{1}}{\partial Z_{1}} & \frac {\partial S_{1}}{\partial Z_{2}}\\ \frac {\partial S_{2}}{\partial Z_{1}} & \frac {\partial S_{2}}{\partial Z_{2}} \end{pmatrix}$
$\begin{pmatrix} \frac{\partial}{\partial Z_{1}} \frac {e^{Z1}}{e^{Z1}+ e^{Z2}} & \frac{\partial}{\partial Z_{2}} \frac {e^{Z1}}{e^{Z1}+ e^{Z2}}\\ \frac{\partial}{\partial Z_{1}} \frac {e^{Z2}}{e^{Z1}+ e^{Z2}} & \frac{\partial}{\partial Z_{2}} \frac {e^{Z2}}{e^{Z1}+ e^{Z2}} \end{pmatrix}$ – (A)

Now the ‘division-rule’ of derivatives is as follows. If u and v are functions of x, then
$\frac{d}{dx} \frac {u}{v} =\frac {vdu -udv}{v^{2}}$
Using this to compute each element of the above Jacobian matrix, we see that
when i=j we have
$\frac {\partial}{\partial Z1}\frac{e^{Z1}}{e^{Z1}+e^{Z2}} = \frac {\sum e^{Z1} - e^{Z1^{2}}}{\sum ^{2}}$
and when $i \neq j$
$\frac {\partial}{\partial Z1}\frac{e^{Z2}}{e^{Z1}+e^{Z2}} = \frac {0 - e^{z1}e^{Z2}}{\sum ^{2}}$
This is of the general form
$\frac {\partial S_{j}}{\partial z_{i}} = S_{i}( 1-S_{j})$ when i=j
and
$\frac {\partial S_{j}}{\partial z_{i}} = -S_{i}S_{j}$ when $i \neq j$
Note: Since the Softmax essentially gives the probability the following
notation is also used
$\frac {\partial p_{j}}{\partial z_{i}} = p_{i}( 1-p_{j})$ when i=j
and
$\frac {\partial p_{j}}{\partial z_{i}} = -p_{i}p_{j} when i \neq j$
If you throw the “Kronecker delta” into the equation, then the above equations can be expressed even more concisely as
$\frac {\partial p_{j}}{\partial z_{i}} = p_{i} (\delta_{ij} - p_{j})$
where $\delta_{ij} = 1$ when i=j and 0 when $i \neq j$

This reduces the Jacobian of the simple 2 output softmax vectors equation (A) as
$\begin{pmatrix} p_{1}(1-p_{1}) & -p_{1}p_{2} \\ -p_{2}p_{1} & p_{2}(1-p_{2}) \end{pmatrix}$
The loss of Softmax is given by
$L = -\sum y_{i} log(p_{i})$
For the 2 valued Softmax output this is
$\frac {dL}{dp1} = -\frac {y_{1}}{p_{1}}$
$\frac {dL}{dp2} = -\frac {y_{2}}{p_{2}}$
Using the chain rule we can write
$\frac {\partial L}{\partial w_{pq}} = \sum _{i}\frac {\partial L}{\partial p_{i}} \frac {\partial p_{i}}{\partial w_{pq}}$ (1)
and
$\frac {\partial p_{i}}{\partial w_{pq}} = \sum _{k}\frac {\partial p_{i}}{\partial z_{k}} \frac {\partial z_{k}}{\partial w_{pq}}$ (2)
In expanded form this is
$\frac {\partial L}{\partial w_{pq}} = \sum _{i}\frac {\partial L}{\partial p_{i}} \sum _{k}\frac {\partial p_{i}}{\partial z_{k}} \frac {\partial z_{k}}{\partial w_{pq}}$
Also
$\frac {\partial L}{\partial Z_{i}} =\sum _{i} \frac {\partial L}{\partial p} \frac {\partial p}{\partial Z_{i}}$
Therefore
$\frac {\partial L}{\partial Z_{1}} =\frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial Z_{1}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial Z_{1}}$
$\frac {\partial L}{\partial z_{1}}=-\frac {y1}{p1} p1(1-p1) - \frac {y2}{p2}*(-p_{2}p_{1})$
Since
$\frac {\partial p_{j}}{\partial z_{i}} = p_{i}( 1-p_{j})$ when i=j
and
$\frac {\partial p_{j}}{\partial z_{i}} = -p_{i}p_{j}$ when $i \neq j$
which simplifies to
$\frac {\partial L}{\partial Z_{1}} = -y_{1} + y_{1}p_{1} + y_{2}p_{1} =$
$p_{1}\sum (y_{1} + y_2) - y_{1}$
$\frac {\partial L}{\partial Z_{1}}= p_{1} - y_{1}$
Since
$\sum_{i} y_{i} =1$
Similarly
$\frac {\partial L}{\partial Z_{2}} =\frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial Z_{2}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial Z_{2}}$
$\frac {\partial L}{\partial z_{2}}=-\frac {y1}{p1}*(p_{1}p_{2}) - \frac {y2}{p2}*p_{2}(1-p_{2})$
$y_{1}p_{2} + y_{2}p_{2} - y_{2}$
$\frac {\partial L}{\partial Z_{2}} =p_{2}\sum (y_{1} + y_2) - y_{2}\\ = p_{2} - y_{2}$
In general this is of the form
$\frac {\partial L}{\partial z_{i}} = p_{i} -y_{i}$
For e.g if the probabilities computed were p=[0.1, 0.7, 0.2] then this implies that the class with probability 0.7 is the likely class. This would imply that the ‘One hot encoding’ for yi would be yi=[0,1,0] therefore the gradient pi-yi = [0.1,-0.3,0.2]

<strong>Note: Further, we could extend this derivation for a Softmax activation output that outputs 3 classes
$S=\begin{pmatrix} \frac{e^{z1}}{e^{z1}+e^{z2}+e^{z3}}\\ \frac{e^{z2}}{e^{z1}+e^{z2}+e^{z3}} \\ \frac{e^{z3}}{e^{z1}+e^{z2}+e^{z3}} \end{pmatrix}$

We could derive
$\frac {\partial L}{\partial z1}= \frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial z_{1}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial z_{1}} +\frac {\partial L}{\partial p_{3}} \frac {\partial p_{3}}{\partial z_{1}}$ which similarly reduces to
$\frac {\partial L}{\partial z_{1}}=-\frac {y1}{p1} p1(1-p1) - \frac {y2}{p2}*(-p_{2}p_{1}) - \frac {y3}{p3}*(-p_{3}p_{1})$
$-y_{1}+ y_{1}p_{1} + y_{2}p_{1} + y_{3}p1 = p_{1}\sum (y_{1} + y_2 + y_3) - y_{1} = p_{1} - y_{1}$
Interestingly, despite the lengthy derivations the final result is simple and intuitive!

As seen in my post ‘Deep Learning from first principles with Python, R and Octave – Part 3 the key equations for forward and backward propagation are

Forward propagation equations layer 1
$Z_{1} = W_{1}X +b_{1}$ and $A_{1} = g(Z_{1})$
Forward propagation equations layer 1
$Z_{2} = W_{2}A_{1} +b_{2}$ and $A_{2} = S(Z_{2})$

Using the result (A) in the back propagation equations below we have
Backward propagation equations layer 2
$\partial L/\partial W_{2} =\partial L/\partial Z_{2}*A_{1}=(p_{2}-y_{2})*A_{1}$
$\partial L/\partial b_{2} =\partial L/\partial Z_{2}=p_{2}-y_{2}$
$\partial L/\partial A_{1} = \partial L/\partial Z_{2} * W_{2}=(p_{2}-y_{2})*W_{2}$
Backward propagation equations layer 1
$\partial L/\partial W_{1} =\partial L/\partial Z_{1} *A_{0}=(p_{1}-y_{1})*A_{0}$
$\partial L/\partial b_{1} =\partial L/\partial Z_{1}=(p_{1}-y_{1})$

2.0 Spiral data set

As I mentioned earlier, I will be using the ‘spiral’ data from CS231n Convolutional Neural Networks to ensure that my vectorized implementations in Python, R and Octave are correct. Here is the ‘spiral’ data set.

import numpy as np
import matplotlib.pyplot as plt
import os
os.chdir("C:/junk/dl-4/dl-4")
exec(open("././DLfunctions41.py").read())

# Create an input data set - Taken from CS231n Convolutional Neural networks
# http://cs231n.github.io/neural-networks-case-study/
N = 100 # number of points per class
D = 2 # dimensionality
K = 3 # number of classes
X = np.zeros((N*K,D)) # data matrix (each row = single example)
y = np.zeros(N*K, dtype='uint8') # class labels
for j in range(K):
  ix = range(N*j,N*(j+1))
  r = np.linspace(0.0,1,N) # radius
  t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta
  X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
  y[ix] = j
# Plot the data
plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral)
plt.savefig("fig1.png", bbox_inches='tight')

The implementations of the vectorized Python, R and Octave code are shown diagrammatically below

2.1 Multi-class classification with Softmax – Python code

A simple 2 layer Neural network with a single hidden layer , with 100 Relu activation units in the hidden layer and the Softmax activation unit in the output layer is used for multi-class classification. This Deep Learning Network, plots the non-linear boundary of the 3 classes as shown below

import numpy as np
import matplotlib.pyplot as plt
import os
os.chdir("C:/junk/dl-4/dl-4")
exec(open("././DLfunctions41.py").read())

# Read the input data
N = 100 # number of points per class
D = 2 # dimensionality
K = 3 # number of classes
X = np.zeros((N*K,D)) # data matrix (each row = single example)
y = np.zeros(N*K, dtype='uint8') # class labels
for j in range(K):
  ix = range(N*j,N*(j+1))
  r = np.linspace(0.0,1,N) # radius
  t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta
  X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
  y[ix] = j
  
# Set the number of features, hidden units in hidden layer and number of classess
numHidden=100 # No of hidden units in hidden layer
numFeats= 2 # dimensionality
numOutput = 3 # number of classes

# Initialize the model
parameters=initializeModel(numFeats,numHidden,numOutput)
W1= parameters['W1']
b1= parameters['b1']
W2= parameters['W2']
b2= parameters['b2']

# Set the learning rate
learningRate=0.6 

# Initialize losses
losses=[]
# Perform Gradient descent
for i in range(10000):
    # Forward propagation through hidden layer with Relu units
    A1,cache1= layerActivationForward(X.T,W1,b1,'relu')
    
    # Forward propagation through output layer with Softmax
    A2,cache2 = layerActivationForward(A1,W2,b2,'softmax')
    
    # No of training examples
    numTraining = X.shape[0]
    # Compute log probs. Take the log prob of correct class based on output y
    correct_logprobs = -np.log(A2[range(numTraining),y])
    # Conpute loss
    loss = np.sum(correct_logprobs)/numTraining
    
    # Print the loss
    if i % 1000 == 0:
        print("iteration %d: loss %f" % (i, loss))
        losses.append(loss)

    dA=0

    # Backward  propagation through output layer with Softmax
    dA1,dW2,db2 = layerActivationBackward(dA, cache2, y, activationFunc='softmax')
    # Backward  propagation through hidden layer with Relu unit
    dA0,dW1,db1 = layerActivationBackward(dA1.T, cache1, y, activationFunc='relu')
    
    #Update paramaters with the learning rate
    W1 += -learningRate * dW1
    b1 += -learningRate * db1
    W2 += -learningRate * dW2.T
    b2 += -learningRate * db2.T

#Plot losses vs iterations  
i=np.arange(0,10000,1000)
plt.plot(i,losses)

plt.xlabel('Iterations')
plt.ylabel('Loss')
plt.title('Losses vs Iterations')
plt.savefig("fig2.png", bbox="tight")

#Compute the multi-class Confusion Matrix
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score

# We need to determine the predicted values from the learnt data
# Forward propagation through hidden layer with Relu units
A1,cache1= layerActivationForward(X.T,W1,b1,'relu')
    
# Forward propagation through output layer with Softmax
A2,cache2 = layerActivationForward(A1,W2,b2,'softmax')
#Compute predicted values from weights and biases
yhat=np.argmax(A2, axis=1)

a=confusion_matrix(y.T,yhat.T)
print("Multi-class Confusion Matrix")
print(a)

## iteration 0: loss 1.098507
## iteration 1000: loss 0.214611
## iteration 2000: loss 0.043622
## iteration 3000: loss 0.032525
## iteration 4000: loss 0.025108
## iteration 5000: loss 0.021365
## iteration 6000: loss 0.019046
## iteration 7000: loss 0.017475
## iteration 8000: loss 0.016359
## iteration 9000: loss 0.015703
## Multi-class Confusion Matrix
## [[ 99   1   0]
##  [  0 100   0]
##  [  0   1  99]]

Check out my compact and minimal book “Practical Machine Learning with R and Python:Second edition- Machine Learning in stereo” available in Amazon in paperback($10.99) and kindle($7.99) versions. My book includes implementations of key ML algorithms and associated measures and metrics. The book is ideal for anybody who is familiar with the concepts and would like a quick reference to the different ML algorithms that can be applied to problems and how to select the best model. Pick your copy today!!

2.2 Multi-class classification with Softmax – R code

The spiral data set created with Python was saved, and is used as the input with R code. The R Neural Network seems to perform much,much slower than both Python and Octave. Not sure why! Incidentally the computation of loss and the softmax derivative are identical for both R and Octave. yet R is much slower. To compute the softmax derivative I create matrices for the One Hot Encoded yi and then stack them before subtracting pi-yi. I am sure there is a more elegant and more efficient way to do this, much like Python. Any suggestions?

library(ggplot2)
library(dplyr)
library(RColorBrewer)

source("DLfunctions41.R")
# Read the spiral dataset
Z <- as.matrix(read.csv("spiral.csv",header=FALSE)) 
Z1=data.frame(Z)
#Plot the dataset
ggplot(Z1,aes(x=V1,y=V2,col=V3)) +geom_point() + 
  scale_colour_gradientn(colours = brewer.pal(10, "Spectral"))

# Setup the data
X <- Z[,1:2]
y <- Z[,3]
X1 <- t(X)
Y1 <- t(y)

# Initialize number of features, number of hidden units in hidden layer and
# number of classes
numFeats<-2 # No features
numHidden<-100 # No of hidden units
numOutput<-3 # No of classes

# Initialize model
parameters <-initializeModel(numFeats, numHidden,numOutput)

W1 <-parameters[['W1']]
b1 <-parameters[['b1']]
W2 <-parameters[['W2']]
b2 <-parameters[['b2']]

# Set the learning rate
learningRate <- 0.5
# Initialize losses
losses <- NULL
# Perform gradient descent
for(i in 0:9000){

# Forward propagation through hidden layer with Relu units
retvals <- layerActivationForward(X1,W1,b1,'relu')
A1 <- retvals[['A']]
cache1 <- retvals[['cache']]
forward_cache1 <- cache1[['forward_cache1']]
activation_cache <- cache1[['activation_cache']]

# Forward propagation through output layer with Softmax units
retvals = layerActivationForward(A1,W2,b2,'softmax')
A2 <- retvals[['A']]
cache2 <- retvals[['cache']]
forward_cache2 <- cache2[['forward_cache1']]
activation_cache2 <- cache2[['activation_cache']]

# No oftraining examples
numTraining <- dim(X)[1]
dA <-0

# Select the elements where the y values are 0, 1 or 2 and make a vector
a=c(A2[y==0,1],A2[y==1,2],A2[y==2,3])
# Take log
correct_probs = -log(a)
# Compute loss
loss= sum(correct_probs)/numTraining

if(i %% 1000 == 0){
sprintf("iteration %d: loss %f",i, loss)
print(loss)
}
# Backward propagation through output layer with Softmax units
retvals = layerActivationBackward(dA, cache2, y, activationFunc='softmax')
dA1 = retvals[['dA_prev']]
dW2= retvals[['dW']]
db2= retvals[['db']]
# Backward propagation through hidden layer with Relu units
retvals = layerActivationBackward(t(dA1), cache1, y, activationFunc='relu')
dA0 = retvals[['dA_prev']]
dW1= retvals[['dW']]
db1= retvals[['db']]

# Update parameters
W1 <- W1 - learningRate * dW1
b1 <- b1 - learningRate * db1
W2 <- W2 - learningRate * t(dW2)
b2 <- b2 - learningRate * t(db2)
}

## [1] 1.212487
## [1] 0.5740867
## [1] 0.4048824
## [1] 0.3561941
## [1] 0.2509576
## [1] 0.7351063
## [1] 0.2066114
## [1] 0.2065875
## [1] 0.2151943
## [1] 0.1318807

#Create iterations
iterations <- seq(0,10)
#df=data.frame(iterations,losses)
ggplot(df,aes(x=iterations,y=losses)) + geom_point() + geom_line(color="blue") +
    ggtitle("Losses vs iterations") + xlab("Iterations") + ylab("Loss")

plotDecisionBoundary(Z,W1,b1,W2,b2)

Multi-class Confusion Matrix

library(caret)
library(e1071)

# Forward propagation through hidden layer with Relu units
retvals <- layerActivationForward(X1,W1,b1,'relu')
A1 <- retvals[['A']]

# Forward propagation through output layer with Softmax units
retvals = layerActivationForward(A1,W2,b2,'softmax')
A2 <- retvals[['A']]
yhat <- apply(A2, 1,which.max) -1
Confusion Matrix and Statistics
          Reference
Prediction  0  1  2
         0 97  0  1
         1  2 96  4
         2  1  4 95

Overall Statistics                                        
               Accuracy : 0.96            
                 95% CI : (0.9312, 0.9792)
    No Information Rate : 0.3333          
    P-Value [Acc > NIR] : <2e-16          
                                          
                  Kappa : 0.94            
 Mcnemar's Test P-Value : 0.5724          
Statistics by Class:

                     Class: 0 Class: 1 Class: 2
Sensitivity            0.9700   0.9600   0.9500
Specificity            0.9950   0.9700   0.9750
Pos Pred Value         0.9898   0.9412   0.9500
Neg Pred Value         0.9851   0.9798   0.9750
Prevalence             0.3333   0.3333   0.3333
Detection Rate         0.3233   0.3200   0.3167
Detection Prevalence   0.3267   0.3400   0.3333
Balanced Accuracy      0.9825   0.9650   0.9625

My book “Practical Machine Learning with R and Python” includes the implementation for many Machine Learning algorithms and associated metrics. Pick up your copy today!

2.3 Multi-class classification with Softmax – Octave code

A 2 layer Neural network with the Softmax activation unit in the output layer is constructed in Octave. The same spiral data set is used for Octave also
source("DL41functions.m") # Read the spiral data data=csvread("spiral.csv"); # Setup the data X=data(:,1:2); Y=data(:,3); # Set the number of features, number of hidden units in hidden layer and number of classes numFeats=2; #No features numHidden=100; # No of hidden units numOutput=3; # No of classes # Initialize model [W1 b1 W2 b2] = initializeModel(numFeats,numHidden,numOutput); # Initialize losses losses=[] #Initialize learningRate learningRate=0.5; for k =1:10000 # Forward propagation through hidden layer with Relu units [A1,cache1 activation_cache1]= layerActivationForward(X',W1,b1,activationFunc ='relu'); # Forward propagation through output layer with Softmax units [A2,cache2 activation_cache2] = layerActivationForward(A1,W2,b2,activationFunc='softmax'); # No of training examples numTraining = size(X)(1); # Select rows where Y=0,1,and 2 and concatenate to a long vector a=[A2(Y==0,1) ;A2(Y==1,2) ;A2(Y==2,3)]; #Select the correct column for log prob correct_probs = -log(a); #Compute log loss loss= sum(correct_probs)/numTraining; if(mod(k,1000) == 0) disp(loss); losses=[losses loss]; endif dA=0; # Backward propagation through output layer with Softmax units [dA1 dW2 db2] = layerActivationBackward(dA, cache2, activation_cache2,Y,activationFunc='softmax'); # Backward propagation through hidden layer with Relu units [dA0,dW1,db1] = layerActivationBackward(dA1', cache1, activation_cache1, Y, activationFunc='relu'); #Update parameters W1 += -learningRate * dW1; b1 += -learningRate * db1; W2 += -learningRate * dW2'; b2 += -learningRate * db2'; endfor # Plot Losses vs Iterations iterations=0:1000:9000 plotCostVsIterations(iterations,losses) # Plot the decision boundary plotDecisionBoundary( X,Y,W1,b1,W2,b2)

The code for the Python, R and Octave implementations can be downloaded from Github at Deep Learning – Part 4

Conclusion

In this post I have implemented a 2 layer Neural Network with the Softmax classifier. In Part 3, I implemented a multi-layer Deep Learning Network. I intend to include the Softmax activation unit into the generalized multi-layer Deep Network along with the other activation units of sigmoid,tanh and relu.

Stick around, I’ll be back!!
Watch this space!

References
1. Deep Learning Specialization
2. Neural Networks for Machine Learning
3. CS231 Convolutional Neural Networks for Visual Recognition
4. Eli Bendersky’s Website – The Softmax function and its derivative
5. Cross Validated – Backpropagation with Softmax / Cross Entropy
6. Stackoverflow – CS231n: How to calculate gradient for Softmax loss function?
7. Math Stack Exchange – Derivative of Softmax
8. The Matrix Calculus for Deep Learning

To see all post click Index of posts

The setup

How it works

Part 1: Bursts, or finding the moment a story breaks

1. Counting hashtags finds bots

2. Names, not hashtags

3. One story, seven events

4. Need for embeddings

5. The burst formula

7. The burst is short

Part 2: From a burst to a story

1. Following a story for a day

2. Twelve checkpoints, facts and claims

3. Many names for one thing

4. Captions that say only what’s new

5. When following goes wrong

A detour to Chennai: a story told in two languages

Play the stories

Share:

Like this:

Share:

Like this:

Share:

Like this:

Share:

Like this:

1. Introduction

Problem statement

1a. A note on ROC and Precision-Recall curves

2. Define recommender lab helper functions

3. Batsman performance estimation

3a. Batting dataframe

3b Data set and data preparation

3c. Exploratory analysis

4. Create a visual representation of the rating data before and after the normalization

4a. Data for analysis

5. Batsman’s times Out

6. Batsman’s Strike rate

7. Batsman’s Sixes

8. Batsman’s Fours

9. Batsman’s Total Runs

10. Batsman’s Balls Faced

11. Generate the Batsmen Performance Estimate

12. Bowler analysis

12a. Bowler dataframe

13. Balls bowled by bowler

14. Runs conceded by bowler

15. Economy Rate of the bowler

16. Wickets Taken by bowler

17. Generate the Bowler Performance estmiate

18. Conclusion

References

Also see

Share:

Like this:

A) Run 1

Run 3

Run 4

References

See also

Share:

Like this:

Reinforcement Learning

1. Gridworld-1

1a. Bellman Update

Findings

1b. Greedify

Findings

1c. Bellman Optimality update

Findings

2.Gridworld 2

##2a. Bellman Update

2b. Greedify

2c. Bellman Optimality Update

3. Another maze

3a. Bellman Update

3b. Greedify

3c. Bellman Optimality Update

Findings

Conclusion:

Share: