# Predicting Molecular Properties EDA

## Special purpose of this EDA
I am not a physicist or chemist, so do not have domain specific knowledge relating to this competition. 

Despite this, I have come up with a cunning plan;

1. Spend some time thinking about what intuitively might affect the scalar coupling constant described in the competition rules.
2. For each of the features identified, explore how it is related to the outcome in the training dataset. 
3. Where features show some systematic relationship to the outcome variable, explore that further by 
using the features in a selection of basic machine learning models.
4. Use my best model as the basis for making predictions using the test data.

In this kernel, I will complete parts 1 and 2, and prepare an output dataset with the relevant features that
I will use in a second kernel that will complete steps 3  and 4. 

## Identifying some really simple features to investigate
The data is described on the competition page, here I will ensure that I can read in the training dataset and
verify the names and types of the feature and outcome variables are what is expected;
```{r get the data}
library(tidyverse)
train <- read_csv('../input/champs-scalar-coupling/train.csv')
names(train)
str(train)
dim(train)
```
We are asked to predict the value of the scalar coupling constant for some pairs of atoms in some molecules. 

So for a given pair, we have information about the molecule they are part of, and information about the two
atoms in the pair. It seems likely that no single piece of information about the pair and their host molecule
will have a great deal of predictive power, but maybe a combination of simple pieces of information could be 
useful. The features I have chosen to look at are;

1. Number of atoms in the host molecule.
2. The atomic numbers of each atom in the pair.
3. The distance between the atoms in the pair.  

## Exploring how the really simple features are related to the scalar coupling constant

### Number of atoms in the host molecule
First read in the structures.csv file;
```{r get structure data}
structures <- read_csv('../input/champs-scalar-coupling/structures.csv')
names(structures)
```

The following code should give the number of atoms in each molecule, by summarising the structures file;
```{r summarise structures file to get number of atoms per molecule}
atom_counts <- structures %>% group_by(molecule_name) %>% summarise(no_of_atoms = n())
```
Now to combine the atom_counts data with the data for the scalar coupling constant data for each atom pair. 
```{r combine the atom_counts data with the training data}
train <- train %>% inner_join(atom_counts)
head(train)
```
Finally, a scatterplot of scalar coupling contant against number of atoms per molecule can be plotted;
```{r scatterplot scalar coupling constant against number of atoms per molecule}
library(ggplot2)
ggplot(train, aes(x=no_of_atoms, y=scalar_coupling_constant)) + geom_point() + geom_smooth()
```

There appears to be a slight trend towards decreasing values of the coupling constant as the number
of atoms in a host molecule decreases. To improve confidence that this is a real effect, it would
be useful to investiagte how the number of atoms per molecule is distributed across all of the molecules. 
A histogram would be usful here.
```{r histogram of number of atoms per molecule}
ggplot(atom_counts, aes(x=no_of_atoms)) + geom_histogram()
```

The distribution of number of atoms per molecule looks vaguely normal apart from a couple of strange
looking gaps where the are no molecules counted at all. Since the overall, but slight, trend towards
decreasing value of the coupling constant seems to be independent of the number density of atoms
across the distribution, I will take this as a real, but small, effect. 

### Atomic number of the atoms in the pair
Atomic number is number of protons in the nucleus of an atom and hence identifies the specific element.
Larger atomic number indicates more protons, more electrons, more complicated charge distributions
around an atom, and might intuitively result in a larger affect on the value of the scalar coupling
constant. Lets see how scalar coupling constant is related to the sum of the atomic number of the atoms
in a pair. 

First we have to get the atomic number data for each element. 
```{r get periodic table of elements}
#install.packages("PeriodicTable")
#library(PeriodicTable)
periodicTable <- read_csv('../input/periodic-table/PeriodicTable.csv')
head(periodicTable)
```

Now to get the atomic number for each atom into the structures data frame. 
```{r add atomic number column to structures file}
# extract atomic number from periodic table
atomic_no_table <- select(periodicTable,Symbol,`Atomic Number`)
atomic_no_table <- rename(atomic_no_table, atomic_no = `Atomic Number`)
atomic_no_table <- rename(atomic_no_table, atom = Symbol)

# join the atomic_no_table to the structures table on "atom" column
structures_extended <- structures %>% inner_join(atomic_no_table)

# Check that the atomic number for each atom is now in the table
head(structures_extended)

# get the atomic numbers for each atom in a pair into the train dataset, add these to get the required feature
train <- train %>% inner_join(structures_extended, by=c("molecule_name"="molecule_name", "atom_index_0"="atom_index"))
train <- rename(train,x_0 = x, y_0 = y, z_0 = z, atomic_no_0 = atomic_no)

train <- train %>% inner_join(structures_extended, by=c("molecule_name"="molecule_name", "atom_index_1"="atom_index"))
train <- rename(train,x_1 = x, y_1 = y, z_1 = z, atomic_no_1 = atomic_no)

train <- mutate(train, atomic_no_sum = atomic_no_0 + atomic_no_1)

# plot the sum of atomic number in each pair against the scalar coupling constant
ggplot(train, aes(x=atomic_no_sum, y=scalar_coupling_constant)) + geom_point() + geom_smooth()
```

The sum of atomic numbers in a pair does not appear to be partiulcarly helpful, but perhaps the 
individual atomic numbers indicating the different types of atoms, might be useful in combination
with other features. 

### distance between the atoms in a pair
In this section, the distance between the atoms in a pair will be calculated, and a scatterplot of
distance against scalar coupling constant made. 

```{r calculate and investigate distance between atoms in a pair}
train <- mutate(train,distance=sqrt((x_1-x_0)**2 + (y_1-y_0)**2 + (z_1-z_0)**2))
                
# plot distance against scalar coupling constant
ggplot(train, aes(x=distance, y=scalar_coupling_constant)) + geom_point() + geom_smooth()
```
The plot has some structure to it, but not a linear one. Looking ahead, a linear model using this feature
would not seem to be applicable. When it comes to exploring machine learning models later, linear
models will not be explored. 

Finally, the new version of the training dataset with the added features will be written out, to be 
processed by part 2 of this EDA. 

```{r write out dataset for EDA part 2}
names(train)
dim(train)
train <- select(train, id, atom_index_0, atom_index_1, atomic_no_0, atomic_no_1, molecule_name,
               distance, no_of_atoms, scalar_coupling_constant)
write_csv(train, path="train_expanded.csv")
```