text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 012d94dd-8736-46d9-b49d-aae96b9a8320
using PlutoTest

# ╔═╡ d02132f0-b0e8-11eb-12cb-353922cb782a
md"""
# [Day 16: Aunt Sue](https://adventofcode.com/2015/day/16)

Your Aunt Sue has given you a wonderful gift, and you'd like to send her a thank you card. However, there's a small problem: she signed it "From, Aunt Sue".

You have 500 Aunts named "Sue".

So, to avoid sending the card to the wrong person, you need to figure out which Aunt Sue (which you conveniently number 1 to 500, for sanity) gave you the gift. You open the present and, as luck would have it, good ol' Aunt Sue got you a My First Crime Scene Analysis Machine! Just what you wanted. Or needed, as the case may be.

The My First Crime Scene Analysis Machine (MFCSAM for short) can detect a few specific compounds in a given sample, as well as how many distinct kinds of those compounds there are. According to the instructions, these are what the MFCSAM can detect:

-   `children`, by human DNA age analysis.
-   `cats`. It doesn't differentiate individual breeds.
-   Several seemingly random breeds of dog: [samoyeds](https://en.wikipedia.org/wiki/Samoyed_%28dog%29), [pomeranians](https://en.wikipedia.org/wiki/Pomeranian_%28dog%29), [akitas](https://en.wikipedia.org/wiki/Akita_%28dog%29), and [vizslas](https://en.wikipedia.org/wiki/Vizsla).
-   `goldfish`. No other kinds of fish.
-   `trees`, all in one group.
-   `cars`, presumably by exhaust or gasoline or something.
-   `perfumes`, which is handy, since many of your Aunts Sue wear a few kinds.

In fact, many of your Aunts Sue have many of these. You put the wrapping from the gift into the MFCSAM. It beeps inquisitively at you a few times and then prints out a message on [ticker tape](https://en.wikipedia.org/wiki/Ticker_tape):

```
children: 3
cats: 7
samoyeds: 2
pomeranians: 3
akitas: 0
vizslas: 0
goldfish: 5
trees: 3
cars: 2
perfumes: 1

```

You make a list of the things you can remember about each Aunt Sue. Things missing from your list aren't zero - you simply don't remember the value.

**What is the number of the Sue that got you the gift?**
"""

# ╔═╡ 56d67502-d140-4066-95e5-44dd2641ed7f
puzzle_input = "Sue 1: children: 1, cars: 8, vizslas: 7
Sue 2: akitas: 10, perfumes: 10, children: 5
Sue 3: cars: 5, pomeranians: 4, vizslas: 1
Sue 4: goldfish: 5, children: 8, perfumes: 3
Sue 5: vizslas: 2, akitas: 7, perfumes: 6
Sue 6: vizslas: 0, akitas: 1, perfumes: 2
Sue 7: perfumes: 8, cars: 4, goldfish: 10
Sue 8: perfumes: 7, children: 2, cats: 1
Sue 9: pomeranians: 3, goldfish: 10, trees: 10
Sue 10: akitas: 7, trees: 8, pomeranians: 4
Sue 11: goldfish: 1, perfumes: 4, cars: 6
Sue 12: samoyeds: 6, trees: 6, perfumes: 2
Sue 13: akitas: 10, pomeranians: 0, vizslas: 2
Sue 14: cars: 2, perfumes: 3, children: 4
Sue 15: goldfish: 2, children: 8, cars: 5
Sue 16: goldfish: 9, cars: 0, vizslas: 5
Sue 17: cats: 5, trees: 6, perfumes: 6
Sue 18: cars: 0, perfumes: 8, pomeranians: 7
Sue 19: trees: 2, goldfish: 5, perfumes: 4
Sue 20: akitas: 4, vizslas: 4, trees: 0
Sue 21: pomeranians: 7, trees: 0, goldfish: 10
Sue 22: cars: 4, vizslas: 0, perfumes: 3
Sue 23: vizslas: 8, trees: 1, akitas: 2
Sue 24: children: 7, trees: 0, akitas: 1
Sue 25: goldfish: 3, akitas: 2, trees: 2
Sue 26: pomeranians: 4, vizslas: 4, samoyeds: 2
Sue 27: cars: 0, trees: 8, akitas: 5
Sue 28: perfumes: 6, cats: 0, cars: 2
Sue 29: trees: 7, akitas: 1, vizslas: 1
Sue 30: perfumes: 9, cars: 9, trees: 10
Sue 31: pomeranians: 5, akitas: 9, samoyeds: 1
Sue 32: pomeranians: 10, vizslas: 5, goldfish: 5
Sue 33: vizslas: 2, akitas: 3, trees: 7
Sue 34: goldfish: 10, perfumes: 0, samoyeds: 7
Sue 35: akitas: 6, cats: 7, perfumes: 10
Sue 36: pomeranians: 8, vizslas: 7, akitas: 6
Sue 37: goldfish: 2, cars: 10, children: 7
Sue 38: goldfish: 2, perfumes: 3, cars: 7
Sue 39: trees: 9, vizslas: 10, cars: 5
Sue 40: goldfish: 1, pomeranians: 0, trees: 2
Sue 41: trees: 2, goldfish: 6, vizslas: 3
Sue 42: akitas: 1, cars: 3, vizslas: 3
Sue 43: akitas: 1, pomeranians: 1, vizslas: 3
Sue 44: goldfish: 7, akitas: 3, vizslas: 10
Sue 45: akitas: 8, samoyeds: 8, goldfish: 2
Sue 46: trees: 0, vizslas: 4, cars: 9
Sue 47: cars: 9, trees: 10, perfumes: 4
Sue 48: akitas: 0, vizslas: 5, perfumes: 4
Sue 49: goldfish: 9, trees: 1, cars: 4
Sue 50: goldfish: 2, perfumes: 5, cars: 2
Sue 51: samoyeds: 1, goldfish: 2, perfumes: 7
Sue 52: cars: 0, perfumes: 4, goldfish: 8
Sue 53: goldfish: 9, vizslas: 2, akitas: 9
Sue 54: trees: 1, goldfish: 9, children: 5
Sue 55: cars: 0, akitas: 5, trees: 4
Sue 56: trees: 4, samoyeds: 5, children: 9
Sue 57: children: 0, vizslas: 8, cars: 3
Sue 58: trees: 4, pomeranians: 5, akitas: 5
Sue 59: vizslas: 10, cats: 3, children: 2
Sue 60: cats: 6, vizslas: 2, cars: 2
Sue 61: akitas: 1, vizslas: 0, children: 4
Sue 62: akitas: 4, trees: 9, children: 10
Sue 63: pomeranians: 6, vizslas: 6, cars: 4
Sue 64: perfumes: 8, pomeranians: 1, children: 8
Sue 65: perfumes: 3, goldfish: 6, trees: 5
Sue 66: goldfish: 10, akitas: 8, vizslas: 4
Sue 67: vizslas: 10, samoyeds: 3, trees: 2
Sue 68: samoyeds: 4, cars: 7, perfumes: 3
Sue 69: perfumes: 2, goldfish: 0, trees: 2
Sue 70: trees: 8, vizslas: 7, akitas: 6
Sue 71: cars: 2, children: 7, perfumes: 3
Sue 72: cars: 1, akitas: 9, perfumes: 0
Sue 73: vizslas: 4, akitas: 7, cars: 5
Sue 74: samoyeds: 3, cars: 3, akitas: 2
Sue 75: trees: 2, cars: 1, vizslas: 7
Sue 76: samoyeds: 9, perfumes: 1, trees: 6
Sue 77: trees: 6, perfumes: 10, cars: 7
Sue 78: trees: 0, children: 8, vizslas: 5
Sue 79: vizslas: 0, trees: 0, samoyeds: 1
Sue 80: trees: 6, goldfish: 8, perfumes: 0
Sue 81: samoyeds: 8, pomeranians: 6, akitas: 5
Sue 82: vizslas: 6, perfumes: 9, akitas: 4
Sue 83: cats: 0, vizslas: 3, pomeranians: 10
Sue 84: cars: 4, perfumes: 6, samoyeds: 5
Sue 85: vizslas: 7, trees: 5, goldfish: 7
Sue 86: goldfish: 2, trees: 2, vizslas: 1
Sue 87: trees: 6, goldfish: 10, pomeranians: 4
Sue 88: vizslas: 1, akitas: 0, perfumes: 8
Sue 89: goldfish: 8, akitas: 3, vizslas: 7
Sue 90: vizslas: 9, akitas: 7, perfumes: 9
Sue 91: children: 7, cars: 7, trees: 9
Sue 92: vizslas: 10, akitas: 8, goldfish: 1
Sue 93: goldfish: 7, vizslas: 2, pomeranians: 0
Sue 94: cats: 2, samoyeds: 6, pomeranians: 3
Sue 95: samoyeds: 4, children: 4, pomeranians: 10
Sue 96: pomeranians: 9, cats: 1, goldfish: 3
Sue 97: trees: 1, akitas: 6, goldfish: 1
Sue 98: vizslas: 7, akitas: 2, perfumes: 7
Sue 99: pomeranians: 6, perfumes: 2, trees: 1
Sue 100: cars: 3, children: 9, trees: 10
Sue 101: children: 0, perfumes: 0, vizslas: 3
Sue 102: cars: 4, goldfish: 5, children: 2
Sue 103: pomeranians: 3, perfumes: 7, cats: 8
Sue 104: akitas: 0, perfumes: 5, vizslas: 5
Sue 105: akitas: 7, vizslas: 2, samoyeds: 8
Sue 106: goldfish: 7, perfumes: 0, cats: 8
Sue 107: cats: 6, pomeranians: 9, cars: 6
Sue 108: akitas: 3, vizslas: 10, cats: 5
Sue 109: akitas: 10, perfumes: 2, cars: 7
Sue 110: goldfish: 7, pomeranians: 1, trees: 1
Sue 111: akitas: 10, samoyeds: 6, vizslas: 6
Sue 112: cats: 6, akitas: 7, trees: 9
Sue 113: akitas: 1, trees: 9, vizslas: 8
Sue 114: vizslas: 2, cats: 1, cars: 4
Sue 115: akitas: 0, trees: 5, goldfish: 7
Sue 116: goldfish: 2, trees: 10, akitas: 2
Sue 117: cars: 4, goldfish: 10, perfumes: 5
Sue 118: cars: 5, perfumes: 6, trees: 0
Sue 119: perfumes: 5, vizslas: 1, cats: 0
Sue 120: perfumes: 8, akitas: 9, vizslas: 4
Sue 121: samoyeds: 2, vizslas: 7, perfumes: 6
Sue 122: children: 6, trees: 9, perfumes: 2
Sue 123: cars: 7, akitas: 0, pomeranians: 0
Sue 124: akitas: 7, cats: 8, vizslas: 5
Sue 125: goldfish: 3, trees: 1, cars: 4
Sue 126: cars: 4, perfumes: 3, akitas: 0
Sue 127: children: 10, vizslas: 5, akitas: 9
Sue 128: akitas: 3, samoyeds: 2, cats: 8
Sue 129: cats: 8, akitas: 1, vizslas: 8
Sue 130: trees: 4, cars: 6, perfumes: 6
Sue 131: akitas: 7, perfumes: 6, goldfish: 9
Sue 132: akitas: 6, vizslas: 7, trees: 1
Sue 133: akitas: 5, vizslas: 7, children: 9
Sue 134: cars: 8, goldfish: 4, pomeranians: 4
Sue 135: samoyeds: 1, pomeranians: 6, akitas: 4
Sue 136: perfumes: 10, goldfish: 1, cars: 3
Sue 137: cars: 3, samoyeds: 6, vizslas: 7
Sue 138: samoyeds: 10, akitas: 3, perfumes: 4
Sue 139: perfumes: 10, vizslas: 2, goldfish: 7
Sue 140: samoyeds: 7, cars: 1, trees: 2
Sue 141: children: 6, cats: 5, cars: 9
Sue 142: cats: 0, trees: 1, akitas: 10
Sue 143: samoyeds: 4, cars: 0, children: 7
Sue 144: trees: 0, cars: 4, perfumes: 8
Sue 145: goldfish: 7, cars: 5, trees: 1
Sue 146: perfumes: 7, cars: 7, goldfish: 0
Sue 147: trees: 2, goldfish: 7, vizslas: 5
Sue 148: samoyeds: 8, perfumes: 1, trees: 0
Sue 149: vizslas: 2, samoyeds: 5, trees: 0
Sue 150: akitas: 4, perfumes: 4, pomeranians: 2
Sue 151: trees: 2, cars: 0, goldfish: 10
Sue 152: goldfish: 7, vizslas: 0, trees: 0
Sue 153: children: 9, cats: 0, pomeranians: 10
Sue 154: cars: 6, goldfish: 10, akitas: 5
Sue 155: perfumes: 9, trees: 2, akitas: 3
Sue 156: pomeranians: 9, perfumes: 5, cars: 9
Sue 157: akitas: 0, trees: 2, cars: 7
Sue 158: goldfish: 10, trees: 8, akitas: 7
Sue 159: akitas: 5, trees: 10, cars: 10
Sue 160: akitas: 3, trees: 5, cars: 8
Sue 161: samoyeds: 2, cars: 7, perfumes: 4
Sue 162: cars: 6, vizslas: 10, pomeranians: 5
Sue 163: cars: 10, perfumes: 6, vizslas: 9
Sue 164: pomeranians: 7, cars: 4, vizslas: 2
Sue 165: goldfish: 9, vizslas: 3, trees: 1
Sue 166: goldfish: 1, samoyeds: 3, trees: 1
Sue 167: vizslas: 4, goldfish: 7, cats: 5
Sue 168: children: 1, cars: 5, samoyeds: 7
Sue 169: trees: 1, samoyeds: 3, goldfish: 6
Sue 170: goldfish: 2, cars: 3, perfumes: 9
Sue 171: cars: 4, goldfish: 0, trees: 6
Sue 172: cats: 8, perfumes: 6, trees: 1
Sue 173: akitas: 9, goldfish: 7, cars: 10
Sue 174: vizslas: 2, trees: 0, akitas: 1
Sue 175: perfumes: 3, vizslas: 8, akitas: 4
Sue 176: perfumes: 0, akitas: 6, goldfish: 3
Sue 177: perfumes: 6, children: 1, goldfish: 10
Sue 178: cars: 5, vizslas: 3, children: 10
Sue 179: perfumes: 3, trees: 8, cats: 9
Sue 180: perfumes: 8, vizslas: 4, trees: 7
Sue 181: perfumes: 7, vizslas: 9, samoyeds: 4
Sue 182: vizslas: 9, trees: 4, pomeranians: 4
Sue 183: trees: 9, cars: 3, goldfish: 5
Sue 184: perfumes: 2, cars: 4, vizslas: 3
Sue 185: children: 10, akitas: 10, cats: 9
Sue 186: cars: 5, samoyeds: 0, trees: 0
Sue 187: trees: 2, goldfish: 3, cars: 4
Sue 188: goldfish: 3, vizslas: 1, cats: 6
Sue 189: trees: 2, pomeranians: 10, cars: 7
Sue 190: perfumes: 10, akitas: 3, samoyeds: 0
Sue 191: cats: 5, vizslas: 6, akitas: 6
Sue 192: samoyeds: 5, trees: 1, perfumes: 8
Sue 193: pomeranians: 0, akitas: 9, cats: 0
Sue 194: trees: 1, goldfish: 0, perfumes: 10
Sue 195: perfumes: 2, akitas: 7, cars: 5
Sue 196: perfumes: 5, samoyeds: 8, cars: 1
Sue 197: vizslas: 2, pomeranians: 9, trees: 1
Sue 198: trees: 8, vizslas: 6, children: 8
Sue 199: pomeranians: 4, cars: 7, vizslas: 5
Sue 200: trees: 0, perfumes: 10, akitas: 10
Sue 201: cats: 9, akitas: 4, vizslas: 0
Sue 202: goldfish: 9, pomeranians: 9, cats: 6
Sue 203: cars: 5, perfumes: 5, trees: 2
Sue 204: pomeranians: 7, children: 2, akitas: 6
Sue 205: samoyeds: 7, pomeranians: 7, children: 6
Sue 206: trees: 1, cars: 1, pomeranians: 4
Sue 207: goldfish: 2, perfumes: 5, trees: 0
Sue 208: perfumes: 2, samoyeds: 4, trees: 1
Sue 209: cars: 8, perfumes: 6, goldfish: 9
Sue 210: perfumes: 4, cars: 8, samoyeds: 3
Sue 211: perfumes: 2, cars: 8, trees: 9
Sue 212: trees: 7, perfumes: 2, akitas: 5
Sue 213: children: 3, goldfish: 5, vizslas: 0
Sue 214: akitas: 6, goldfish: 0, children: 0
Sue 215: trees: 8, akitas: 3, goldfish: 1
Sue 216: goldfish: 6, perfumes: 8, akitas: 3
Sue 217: children: 7, trees: 2, vizslas: 6
Sue 218: goldfish: 8, samoyeds: 4, pomeranians: 6
Sue 219: goldfish: 8, samoyeds: 0, children: 9
Sue 220: perfumes: 1, cars: 8, vizslas: 6
Sue 221: perfumes: 9, cars: 10, children: 10
Sue 222: perfumes: 9, vizslas: 1, trees: 0
Sue 223: goldfish: 1, akitas: 2, vizslas: 8
Sue 224: samoyeds: 8, akitas: 7, vizslas: 4
Sue 225: goldfish: 1, cars: 4, perfumes: 10
Sue 226: goldfish: 9, trees: 4, perfumes: 5
Sue 227: vizslas: 5, trees: 4, goldfish: 7
Sue 228: cars: 1, cats: 10, perfumes: 4
Sue 229: vizslas: 8, cars: 10, akitas: 4
Sue 230: cats: 1, children: 8, vizslas: 3
Sue 231: perfumes: 7, cats: 6, samoyeds: 7
Sue 232: cars: 3, children: 9, perfumes: 7
Sue 233: vizslas: 1, samoyeds: 2, children: 2
Sue 234: trees: 1, samoyeds: 8, children: 2
Sue 235: trees: 6, akitas: 9, goldfish: 7
Sue 236: children: 10, trees: 0, samoyeds: 8
Sue 237: pomeranians: 4, trees: 1, goldfish: 2
Sue 238: vizslas: 4, akitas: 2, cars: 0
Sue 239: goldfish: 9, cars: 10, perfumes: 4
Sue 240: perfumes: 3, vizslas: 6, trees: 6
Sue 241: pomeranians: 6, akitas: 4, trees: 2
Sue 242: cars: 8, perfumes: 5, children: 7
Sue 243: trees: 4, perfumes: 7, cars: 3
Sue 244: perfumes: 6, akitas: 1, vizslas: 7
Sue 245: akitas: 3, perfumes: 9, samoyeds: 0
Sue 246: pomeranians: 3, vizslas: 9, samoyeds: 1
Sue 247: cars: 0, goldfish: 7, cats: 2
Sue 248: trees: 5, goldfish: 6, perfumes: 3
Sue 249: trees: 0, pomeranians: 7, perfumes: 9
Sue 250: cars: 9, trees: 1, goldfish: 10
Sue 251: perfumes: 3, cars: 8, trees: 7
Sue 252: cars: 5, akitas: 7, trees: 8
Sue 253: perfumes: 7, akitas: 3, trees: 8
Sue 254: goldfish: 8, samoyeds: 1, vizslas: 7
Sue 255: perfumes: 3, cars: 4, children: 6
Sue 256: perfumes: 9, trees: 8, children: 7
Sue 257: trees: 8, children: 6, cars: 4
Sue 258: vizslas: 1, trees: 10, goldfish: 9
Sue 259: vizslas: 5, trees: 6, goldfish: 9
Sue 260: trees: 0, goldfish: 6, cars: 7
Sue 261: cars: 1, perfumes: 4, goldfish: 9
Sue 262: cars: 7, goldfish: 9, cats: 9
Sue 263: cars: 0, children: 5, goldfish: 8
Sue 264: cars: 2, akitas: 8, trees: 0
Sue 265: perfumes: 9, children: 8, samoyeds: 7
Sue 266: cats: 1, children: 1, vizslas: 10
Sue 267: vizslas: 8, children: 2, trees: 6
Sue 268: akitas: 10, vizslas: 3, cats: 2
Sue 269: children: 4, goldfish: 1, cats: 6
Sue 270: vizslas: 5, cars: 9, akitas: 9
Sue 271: vizslas: 5, children: 4, akitas: 3
Sue 272: cars: 1, goldfish: 0, vizslas: 0
Sue 273: goldfish: 10, samoyeds: 1, akitas: 2
Sue 274: goldfish: 10, children: 2, pomeranians: 0
Sue 275: children: 0, vizslas: 1, samoyeds: 6
Sue 276: children: 1, vizslas: 3, samoyeds: 1
Sue 277: perfumes: 4, cats: 6, children: 10
Sue 278: pomeranians: 7, goldfish: 3, cars: 4
Sue 279: perfumes: 5, goldfish: 9, trees: 7
Sue 280: goldfish: 6, trees: 5, perfumes: 8
Sue 281: cars: 2, akitas: 1, vizslas: 7
Sue 282: vizslas: 4, akitas: 3, children: 8
Sue 283: pomeranians: 8, akitas: 9, vizslas: 4
Sue 284: samoyeds: 10, trees: 10, pomeranians: 2
Sue 285: akitas: 9, perfumes: 7, goldfish: 6
Sue 286: akitas: 2, vizslas: 7, goldfish: 10
Sue 287: pomeranians: 8, cars: 6, samoyeds: 5
Sue 288: pomeranians: 1, trees: 0, goldfish: 0
Sue 289: trees: 10, samoyeds: 1, children: 0
Sue 290: cats: 10, samoyeds: 6, trees: 0
Sue 291: vizslas: 9, trees: 6, goldfish: 5
Sue 292: cats: 4, perfumes: 8, cars: 3
Sue 293: goldfish: 10, perfumes: 10, cats: 0
Sue 294: cats: 7, trees: 6, akitas: 4
Sue 295: vizslas: 8, cars: 1, akitas: 6
Sue 296: vizslas: 5, akitas: 10, trees: 1
Sue 297: pomeranians: 8, samoyeds: 5, vizslas: 4
Sue 298: perfumes: 10, children: 5, vizslas: 2
Sue 299: cars: 10, akitas: 7, cats: 5
Sue 300: trees: 1, perfumes: 7, cars: 7
Sue 301: cars: 9, vizslas: 1, perfumes: 3
Sue 302: perfumes: 9, vizslas: 1, akitas: 5
Sue 303: akitas: 9, trees: 1, goldfish: 10
Sue 304: children: 10, vizslas: 6, pomeranians: 8
Sue 305: trees: 3, goldfish: 6, cats: 9
Sue 306: cars: 5, perfumes: 9, vizslas: 5
Sue 307: children: 0, goldfish: 7, trees: 2
Sue 308: trees: 9, samoyeds: 4, cars: 0
Sue 309: cats: 8, vizslas: 2, perfumes: 3
Sue 310: cars: 6, pomeranians: 6, vizslas: 6
Sue 311: vizslas: 6, akitas: 7, cats: 10
Sue 312: trees: 0, goldfish: 7, cars: 0
Sue 313: perfumes: 5, akitas: 5, cars: 2
Sue 314: akitas: 10, vizslas: 3, samoyeds: 8
Sue 315: cars: 3, perfumes: 1, goldfish: 8
Sue 316: pomeranians: 6, goldfish: 9, perfumes: 1
Sue 317: goldfish: 4, akitas: 6, cars: 2
Sue 318: perfumes: 8, vizslas: 8, akitas: 0
Sue 319: akitas: 10, cars: 5, vizslas: 6
Sue 320: vizslas: 4, akitas: 3, cats: 4
Sue 321: goldfish: 4, akitas: 8, cars: 8
Sue 322: pomeranians: 5, vizslas: 7, cats: 1
Sue 323: perfumes: 1, trees: 6, goldfish: 0
Sue 324: goldfish: 6, trees: 10, cars: 10
Sue 325: akitas: 2, samoyeds: 6, trees: 9
Sue 326: vizslas: 4, akitas: 7, cars: 9
Sue 327: children: 3, perfumes: 4, cars: 1
Sue 328: akitas: 9, perfumes: 6, cars: 10
Sue 329: perfumes: 2, goldfish: 0, trees: 1
Sue 330: vizslas: 10, pomeranians: 7, goldfish: 6
Sue 331: trees: 3, vizslas: 8, cars: 3
Sue 332: akitas: 2, cats: 1, goldfish: 8
Sue 333: cars: 6, trees: 2, vizslas: 0
Sue 334: samoyeds: 7, cars: 7, trees: 3
Sue 335: cats: 7, children: 1, perfumes: 8
Sue 336: akitas: 5, goldfish: 10, vizslas: 5
Sue 337: cats: 3, vizslas: 0, akitas: 10
Sue 338: perfumes: 8, cars: 1, trees: 8
Sue 339: cars: 4, samoyeds: 8, children: 2
Sue 340: goldfish: 9, pomeranians: 1, samoyeds: 1
Sue 341: akitas: 3, trees: 0, goldfish: 2
Sue 342: perfumes: 4, vizslas: 8, pomeranians: 9
Sue 343: akitas: 4, cars: 5, goldfish: 4
Sue 344: samoyeds: 5, cats: 4, trees: 0
Sue 345: samoyeds: 4, cars: 8, akitas: 2
Sue 346: akitas: 3, vizslas: 10, perfumes: 10
Sue 347: goldfish: 10, akitas: 4, cars: 1
Sue 348: perfumes: 10, cats: 4, vizslas: 5
Sue 349: akitas: 2, vizslas: 4, cars: 7
Sue 350: akitas: 5, vizslas: 5, cars: 6
Sue 351: vizslas: 8, perfumes: 6, cars: 3
Sue 352: cars: 10, vizslas: 0, goldfish: 10
Sue 353: cars: 10, perfumes: 5, children: 7
Sue 354: vizslas: 6, akitas: 3, samoyeds: 9
Sue 355: akitas: 2, perfumes: 7, cars: 10
Sue 356: cars: 10, perfumes: 7, children: 6
Sue 357: akitas: 4, cars: 8, trees: 1
Sue 358: trees: 2, cars: 1, goldfish: 2
Sue 359: vizslas: 5, cars: 9, trees: 4
Sue 360: perfumes: 4, akitas: 3, cars: 3
Sue 361: children: 3, akitas: 2, cats: 5
Sue 362: cars: 8, cats: 4, akitas: 10
Sue 363: cats: 2, trees: 1, vizslas: 4
Sue 364: vizslas: 2, pomeranians: 5, samoyeds: 9
Sue 365: samoyeds: 2, akitas: 7, goldfish: 9
Sue 366: goldfish: 8, trees: 7, cats: 2
Sue 367: perfumes: 2, vizslas: 6, trees: 5
Sue 368: cars: 5, samoyeds: 0, perfumes: 6
Sue 369: samoyeds: 10, trees: 10, vizslas: 1
Sue 370: trees: 2, vizslas: 3, cars: 4
Sue 371: akitas: 6, pomeranians: 2, cats: 4
Sue 372: trees: 2, perfumes: 3, goldfish: 9
Sue 373: vizslas: 5, children: 0, pomeranians: 6
Sue 374: trees: 1, vizslas: 8, perfumes: 10
Sue 375: cars: 0, akitas: 6, children: 0
Sue 376: akitas: 1, vizslas: 0, trees: 0
Sue 377: samoyeds: 10, cats: 5, pomeranians: 0
Sue 378: goldfish: 3, pomeranians: 7, cats: 7
Sue 379: perfumes: 0, cats: 0, trees: 8
Sue 380: perfumes: 4, samoyeds: 1, akitas: 7
Sue 381: akitas: 4, pomeranians: 2, children: 4
Sue 382: vizslas: 9, akitas: 4, trees: 10
Sue 383: trees: 1, vizslas: 10, akitas: 6
Sue 384: trees: 3, akitas: 8, goldfish: 3
Sue 385: goldfish: 6, perfumes: 2, children: 9
Sue 386: children: 10, akitas: 7, goldfish: 7
Sue 387: goldfish: 3, vizslas: 10, perfumes: 5
Sue 388: children: 4, trees: 0, cars: 2
Sue 389: trees: 0, cats: 3, goldfish: 10
Sue 390: samoyeds: 9, pomeranians: 0, cats: 6
Sue 391: samoyeds: 10, trees: 3, akitas: 4
Sue 392: akitas: 9, goldfish: 10, perfumes: 7
Sue 393: goldfish: 6, cars: 2, akitas: 9
Sue 394: trees: 4, goldfish: 9, vizslas: 7
Sue 395: vizslas: 4, samoyeds: 1, goldfish: 6
Sue 396: vizslas: 5, cats: 0, samoyeds: 1
Sue 397: goldfish: 7, cats: 0, trees: 7
Sue 398: cars: 10, akitas: 1, vizslas: 7
Sue 399: samoyeds: 10, cats: 6, goldfish: 6
Sue 400: cats: 6, samoyeds: 0, trees: 2
Sue 401: trees: 1, children: 4, goldfish: 2
Sue 402: cats: 8, vizslas: 4, children: 3
Sue 403: cars: 9, perfumes: 8, pomeranians: 2
Sue 404: goldfish: 8, trees: 2, cars: 5
Sue 405: perfumes: 1, pomeranians: 5, vizslas: 5
Sue 406: perfumes: 6, trees: 2, pomeranians: 6
Sue 407: trees: 0, goldfish: 6, cars: 6
Sue 408: trees: 0, samoyeds: 7, goldfish: 9
Sue 409: samoyeds: 10, goldfish: 6, pomeranians: 0
Sue 410: perfumes: 5, vizslas: 6, trees: 0
Sue 411: goldfish: 2, trees: 2, pomeranians: 0
Sue 412: pomeranians: 4, perfumes: 8, cats: 8
Sue 413: vizslas: 4, cars: 5, akitas: 1
Sue 414: perfumes: 2, trees: 8, goldfish: 7
Sue 415: akitas: 3, trees: 1, perfumes: 3
Sue 416: cars: 7, trees: 1, perfumes: 8
Sue 417: cars: 5, goldfish: 5, trees: 1
Sue 418: cars: 9, goldfish: 4, samoyeds: 2
Sue 419: pomeranians: 8, akitas: 1, goldfish: 6
Sue 420: cars: 0, cats: 0, children: 8
Sue 421: akitas: 10, goldfish: 1, vizslas: 8
Sue 422: children: 8, vizslas: 6, samoyeds: 10
Sue 423: samoyeds: 3, goldfish: 10, vizslas: 8
Sue 424: cars: 3, children: 7, goldfish: 4
Sue 425: cars: 9, perfumes: 9, goldfish: 8
Sue 426: akitas: 5, trees: 10, vizslas: 10
Sue 427: vizslas: 10, cars: 3, akitas: 7
Sue 428: cats: 6, perfumes: 5, goldfish: 10
Sue 429: goldfish: 7, trees: 5, vizslas: 10
Sue 430: perfumes: 3, trees: 7, cars: 3
Sue 431: cars: 2, vizslas: 1, akitas: 6
Sue 432: pomeranians: 8, perfumes: 5, cars: 3
Sue 433: children: 8, cars: 0, perfumes: 7
Sue 434: samoyeds: 0, vizslas: 9, akitas: 10
Sue 435: akitas: 3, vizslas: 8, cats: 4
Sue 436: goldfish: 5, trees: 8, samoyeds: 8
Sue 437: cars: 10, samoyeds: 9, goldfish: 7
Sue 438: samoyeds: 5, akitas: 7, perfumes: 9
Sue 439: goldfish: 10, perfumes: 5, cars: 0
Sue 440: pomeranians: 1, samoyeds: 9, children: 4
Sue 441: vizslas: 4, perfumes: 2, cats: 5
Sue 442: trees: 0, pomeranians: 3, cars: 7
Sue 443: akitas: 0, cars: 2, vizslas: 10
Sue 444: children: 1, akitas: 9, trees: 0
Sue 445: cars: 5, perfumes: 7, goldfish: 9
Sue 446: akitas: 0, perfumes: 1, vizslas: 2
Sue 447: vizslas: 7, perfumes: 0, cars: 5
Sue 448: vizslas: 6, goldfish: 10, trees: 0
Sue 449: cars: 7, vizslas: 7, trees: 3
Sue 450: pomeranians: 4, akitas: 4, vizslas: 8
Sue 451: cats: 4, perfumes: 8, children: 3
Sue 452: samoyeds: 8, akitas: 9, cars: 1
Sue 453: cars: 8, akitas: 5, vizslas: 2
Sue 454: vizslas: 9, perfumes: 4, akitas: 4
Sue 455: akitas: 3, goldfish: 2, vizslas: 6
Sue 456: cars: 4, perfumes: 5, goldfish: 10
Sue 457: trees: 9, pomeranians: 4, goldfish: 10
Sue 458: pomeranians: 1, perfumes: 9, children: 6
Sue 459: samoyeds: 0, goldfish: 8, vizslas: 6
Sue 460: cars: 10, goldfish: 8, samoyeds: 8
Sue 461: akitas: 8, goldfish: 9, vizslas: 2
Sue 462: cars: 1, vizslas: 2, akitas: 8
Sue 463: goldfish: 2, akitas: 4, samoyeds: 10
Sue 464: children: 5, perfumes: 5, cars: 5
Sue 465: perfumes: 9, trees: 0, samoyeds: 6
Sue 466: akitas: 5, goldfish: 3, cats: 6
Sue 467: perfumes: 3, goldfish: 0, trees: 4
Sue 468: goldfish: 2, children: 4, trees: 1
Sue 469: cars: 0, perfumes: 8, children: 7
Sue 470: vizslas: 8, cats: 5, samoyeds: 9
Sue 471: pomeranians: 7, trees: 2, goldfish: 3
Sue 472: goldfish: 8, akitas: 4, perfumes: 5
Sue 473: perfumes: 2, pomeranians: 3, cars: 8
Sue 474: samoyeds: 0, akitas: 7, pomeranians: 6
Sue 475: vizslas: 7, perfumes: 1, trees: 6
Sue 476: vizslas: 3, samoyeds: 1, perfumes: 10
Sue 477: cars: 6, perfumes: 5, vizslas: 2
Sue 478: pomeranians: 1, goldfish: 3, akitas: 7
Sue 479: goldfish: 10, trees: 0, cars: 3
Sue 480: cats: 3, akitas: 5, vizslas: 8
Sue 481: pomeranians: 5, vizslas: 2, trees: 3
Sue 482: cars: 8, samoyeds: 10, goldfish: 10
Sue 483: pomeranians: 3, vizslas: 6, goldfish: 5
Sue 484: perfumes: 7, vizslas: 4, akitas: 7
Sue 485: goldfish: 1, trees: 0, perfumes: 10
Sue 486: goldfish: 6, perfumes: 0, akitas: 10
Sue 487: cats: 2, akitas: 10, trees: 1
Sue 488: akitas: 1, goldfish: 3, cars: 7
Sue 489: goldfish: 3, akitas: 6, vizslas: 6
Sue 490: goldfish: 8, perfumes: 2, akitas: 2
Sue 491: trees: 4, vizslas: 8, perfumes: 6
Sue 492: cars: 9, perfumes: 3, cats: 0
Sue 493: trees: 3, vizslas: 6, goldfish: 7
Sue 494: trees: 8, samoyeds: 1, perfumes: 5
Sue 495: children: 9, akitas: 8, vizslas: 4
Sue 496: vizslas: 2, pomeranians: 1, perfumes: 7
Sue 497: trees: 2, akitas: 4, vizslas: 6
Sue 498: akitas: 8, pomeranians: 7, trees: 0
Sue 499: perfumes: 6, goldfish: 3, vizslas: 7
Sue 500: cars: 1, perfumes: 6, vizslas: 1"

# ╔═╡ 3278c05d-c412-457e-b6c4-37d76d576cd4
function parsesues(input)
	sues = Dict()

	for line in split(input, "\n")
		colonpos = findfirst(":", line)[1]
		id = parse(Int, line[5: colonpos-1])

		data = Dict()
		for info in split(line[colonpos+2:end], ", ")
			parts = split(info, ": ")

			data[parts[1]] = parse(Int, parts[2])
		end

		sues[id] = data
	end
	sues
end

# ╔═╡ 255547d8-d696-4f62-95a1-56b2a0bfdb9e
firstindex

# ╔═╡ 18ce336d-a5cb-42d2-8c2f-5e3eb26d9321
puzzle_sues = parsesues(puzzle_input)

# ╔═╡ 1a49f050-69ff-4404-a494-2b8e35fdaebb
puzzle_facts = Dict(
	"children" => 3,
	"cats" => 7,
	"samoyeds" => 2,
	"pomeranians" => 3,
	"akitas" => 0,
	"vizslas" => 0,
	"goldfish" => 5,
	"trees" => 3,
	"cars" => 2,
	"perfumes" => 1
)

# ╔═╡ b46d289c-e2a1-45b1-9f2e-cac5e8093510
function findsue1(sues, facts)
	valid_sues = filter(id -> begin
		for key in keys(sues[id])
			if facts[key] != sues[id][key]
				return false
			end
		end

		true
	end, keys(sues))

	@assert length(valid_sues) == 1
	return first(valid_sues)
end

# ╔═╡ 3c7002f4-f218-4329-b672-01fe36023d28
part1 = findsue1(puzzle_sues, puzzle_facts)

# ╔═╡ d9dc6c87-6b3d-44e0-8ab5-45358a23c23d
@test part1 == 213

# ╔═╡ 39160fe9-a017-44f0-b298-7d7397f4efd7
md"Your puzzle answer was `213`."

# ╔═╡ 7ae737f3-dfe7-4509-b98c-c6b827e774a9
md"""
# Part Two

As you're about to send the thank you note, something in the MFCSAM's instructions catches your eye. Apparently, it has an outdated retroencabulator, and so the output from the machine isn't exact values - some of them indicate ranges.

In particular, the cats and trees readings indicates that there are greater than that many (due to the unpredictable nuclear decay of cat dander and tree pollen), while the pomeranians and goldfish readings indicate that there are fewer than that many (due to the modial interaction of magnetoreluctance).

**What is the number of the real Aunt Sue?**
"""

# ╔═╡ e6ab6c98-3053-440d-9913-c044a7e97853
function findsue2(sues, facts)
	valid_sues = filter(id -> begin
		for key in keys(sues[id])
			fact = facts[key]
			value = sues[id][key]

			if key == "cats" || key == "trees"
				if value < fact
					return false
				end
			elseif key == "pomeranians" || key == "goldfish"
				if value > fact
					return false
				end
			else
				if fact != value
					return false
				end
			end
		end

		true
	end, keys(sues))

	# @assert length(valid_sues) == 1
	return collect(valid_sues)
end

# ╔═╡ afcbab33-d3fa-4928-9171-4dafa83d8536
part2 = findsue2(puzzle_sues, puzzle_facts)[2] # why two results?

# ╔═╡ 6129d2e9-331e-4760-bd12-a037369ccab9
@test part2 == 323

# ╔═╡ e2edeffa-752d-4242-b5cf-1798221a0038
md"Your puzzle answer was `323`."

# ╔═╡ Cell order:
# ╠═012d94dd-8736-46d9-b49d-aae96b9a8320
# ╟─d02132f0-b0e8-11eb-12cb-353922cb782a
# ╟─56d67502-d140-4066-95e5-44dd2641ed7f
# ╠═3278c05d-c412-457e-b6c4-37d76d576cd4
# ╠═255547d8-d696-4f62-95a1-56b2a0bfdb9e
# ╠═18ce336d-a5cb-42d2-8c2f-5e3eb26d9321
# ╠═1a49f050-69ff-4404-a494-2b8e35fdaebb
# ╠═b46d289c-e2a1-45b1-9f2e-cac5e8093510
# ╠═3c7002f4-f218-4329-b672-01fe36023d28
# ╠═d9dc6c87-6b3d-44e0-8ab5-45358a23c23d
# ╟─39160fe9-a017-44f0-b298-7d7397f4efd7
# ╟─7ae737f3-dfe7-4509-b98c-c6b827e774a9
# ╠═e6ab6c98-3053-440d-9913-c044a7e97853
# ╠═afcbab33-d3fa-4928-9171-4dafa83d8536
# ╠═6129d2e9-331e-4760-bd12-a037369ccab9
# ╟─e2edeffa-752d-4242-b5cf-1798221a0038


text:
"""
Implements tabulated MACs from various literature sources including
`Henke1974`, `Henke1982`, `Bastin1988`, `Bastin1989`, `Henke1993`,
`Bastin1997`, `Ruste1979`, `Kohlhaas1970`, `Weisweiler1975`,
`Bastin1990`, `Poml2020`, `Ruste1975`, `Farthing1990` and `Sabbatucci2016`.

The custom MACs are stored in the file `data/specialmacs.csv`.
"""
struct CustomMAC <: NeXLAlgorithm end

# The file '../data/specialmacs.csv' is the place to add custom MACs from any source.

let custommacs = nothing
    global function getcustommacs()
        if isnothing(custommacs)
            @info "Loading custom MACs."
            custommacs = Dict{Tuple{Symbol,Int,CharXRay},Float64}()
            for row in CSV.File(
                joinpath(dirname(pathof(@__MODULE__)), "..", "data", "specialmacs.csv"),
                header = 1,
            )
                custommacs[(Symbol(row[3]), row[1], parse(CharXRay, row[2]))] = row[4]
            end
        end
        return custommacs
    end
end

mac(elm::Element, cxr::CharXRay, model::Symbol, def = missing) =
    get(getcustommacs(), (model, z(elm), cxr), def)

function getcustommacs(
    model::Symbol,
    withsimilar = true,
)::Vector{Tuple{Element,CharXRay,Float64}}
    macs = getcustommacs()
    tmp = [
        (elements[key[2]], key[3], macs[key]) for
        key in filter(k -> k[1] == model, keys(macs))
    ]
    if withsimilar
        res = []
        for macs in tmp
            cxrs = characteristic(element(macs[2]), transitionsbyshell[shell(macs[2])])
            append!(res, [(macs[1], cxr, macs[3]) for cxr in cxrs])
        end
    else
        res = tmp
    end
    return res
end

function NeXLUncertainties.asa(::Type{DataFrame}, ::Type{CustomMAC}, model::Symbol)
    macs = getcustommacs()
    df = DataFrame(Element=String[], CharXRay=CharXRay[], MAC=Float64[])
    for key in filter(k -> k[1] == model, keys(macs))
        push!(df, (symbol(elements[key[2]]), key[3], macs[key]))
    end
    return df
end

ext: jl


text:
# tools for displaying latin squares

# alphabets

abc = [ "$ch" for ch in 'a':'z']
ABC = uppercase.(abc)
greek = [ "$(Char(t))" for t in 945:969 if t != 962 ]
GREEK = uppercase.(greek)

export print_latin

"""
`print_latin(A)` prints out the Latin square `A` using uppercase
Latin letters.

`print_latin(A,B)` prints out the pair of (presumably orthogonal)
Latin squares using upper case Latin letters for `A` and lower case Greek
letters for `B`.

Both versions permit alternate character sets:
* `print_latin(A,charlist)`
* `print_latin(A,B,charlist1,charlist2)`

The following pre-built character sets are available:
* `LatinSquares.abc`: lower case Latin letters
* `LatinSquares.ABC`: upper case Latin letters
* `LatinSquares.greek`: lower case Greek letters
* `LatinSquares.GREEK`: upper case Greek letters
"""

function print_latin(A::Matrix{Int}, chars::Array{String,1}=ABC)
    r,c = size(A)
    for i=1:r
        for j=1:c
            print(chars[A[i,j]])
            if j < c
                print(" ")
            else
                println()
            end
        end
    end
    nothing
end

function print_latin(A::Matrix{Int}, B::Matrix{Int}, chars1::Array{String}=ABC, chars2::Array{String}=greek)
    r,c = size(A)
    for i=1:r
        for j=1:c
            print(chars1[A[i,j]]*chars2[B[i,j]])
            if j<c
                print(" ")
            else
                println()
            end
        end
    end
    nothing
end


text:
mutable struct TimerXMLElement <: MyXMLElement
    el::XMLOrNothing
    mcmc_el::MCMCXMLElement
    filename::String

    TimerXMLElement(mcmc_el::MCMCXMLElement) = new(nothing, mcmc_el, "")
    TimerXMLElement(mcmc_el::MCMCXMLElement, filename::String) = new(nothing, mcmc_el, filename)
end

function make_xml(t_el::TimerXMLElement)
    el = new_element(bn.REPORT)

    if length(t_el.filename) > 0
        set_attribute(el, bn.FILENAME, t_el.filename)
    end

    prop_el = new_child(el, bn.PROPERTY)
    set_attribute(prop_el, bn.NAME, bn.TIMER)
    add_ref_el(prop_el, t_el.mcmc_el.el)

    t_el.el = el
    return el
end

function set_filename(t::TimerXMLElement, fn::String)
    t.filename = fn
end


text:
function calcBounds(nbit::Int,numlist::Array{Float64,1};overflow::Bool=false,underflow::Bool=true)

    numlist[1] == 0 || throw(error("The first element of x should be 0."))
    numlist[end] == Inf || throw(error("The last element of x should be Inf."))
    numlist == sort(numlist) || throw(error("Elements of x are expected to be in an ascending order."))

    if nbit == 8
        bounds = bounds8
    elseif nbit == 16
        bounds = bounds16
    else
        throw(error("Only 8/16bit supported."))
    end

    bounds[1] = 0.0
    bounds[end] = floatmax(Float64)     # no overflow rounding mode

    # use the arithmetic mean
    for i in 2:length(bounds)-1
        bounds[i] = 0.5*(numlist[i-1] + numlist[i])
    end

    # tie to even / round half to even by shifting every second bound a tiny bit
    for i in 2:2:length(bounds)-1
        bounds[i] = nextfloat(bounds[i])
    end

    # UNDERFLOW ROUNDING MODE
    if ~underflow
        bounds[2] = floatmin(Float64)
    end

    #TODO this has no effect, always no overflow somehow
    # OVERFLOW ROUNDING MODE
    if overflow
        # overflow for everything > 2 times the largest representable number
        bounds[end] = 2*numlist[end-1]
    else
        bounds[end] = floatmax(Float64)
    end

    return nothing
end

function setOverflow8(overflow::Bool)
    calcBounds(8,sonum8,overflow=overflow)
end

function setOverflow16(overflow::Bool)
    calcBounds(16,sonum16,overflow=overflow)
end

function setUnderflow8(underflow::Bool)
    calcBounds(8,sonum8,underflow=underflow)
end

function setUnderflow16(underflow::Bool)
    calcBounds(16,sonum16,underflow=underflow)
end


text:
using Requests
import Requests: get, post, put, delete, options

function safe(name)
	result = replace(name, ":", " ")
	result = replace(result, "\\", " ")
	result = replace(result, "/", " ")
	result = replace(result, "*", " ")
	result = replace(result, "?", " ")
	result = replace(result, "|", " ")
	result = replace(result, "<", " ")
	result = replace(result, ">", " ")
	result = replace(result, "\"", " ")
	result = replace(result, r"[^\x20-\x7e]", "")
	result = strip(result)
	return result
end

function urlsafe(title)
	result = replace(title, r"[^\x20-\x7e]", "")
	result = strip(result)
	result = replace(result, " ", "+")
	result = replace(result, ",", "")
	result = replace(result, "!", "")
	result = replace(result, "'", "")
	result = replace(result, "-", "")
	return result
end

function main(dataPath)
	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	movies = movies[1]

	@printf "Finished loading, going into loop\n"
	@printf "%s/%s" 0 length(movies[:,1])

	for i = 1:length(movies[:, 1])
		if i % 1_000 == 0
			@printf "\n%s/%s" i length(movies[:,1])
		elseif i % 100 == 0
			@printf "."
		end

		title = movies[i, 2]

		if !isfile("movies/$(safe(title))/3.jpg")
			data = []
			num = 1
			try
				data = Requests.json(get("https://www.googleapis.com/youtube/v3/search?part=snippet&maxResults=50&q=$(urlsafe(title))+trailer&key=AIzaSyDRmOXB5UYs51X_1qNO_5iXn_gsc8OBINw"))

				while data["items"][num]["id"]["kind"] != "youtube#video"
					num += 1
				end

				if !isdir("movies/$(safe(title))")
					mkdir("movies/$(safe(title))")
				end

				id = data["items"][num]["id"]["videoId"]
				for j = 1:3
					image = get("https://i.ytimg.com/vi/$(id)/$(j).jpg")
					save(image, "movies/$(safe(title))/$(j).jpg")
				end
			catch ex
				if data != [] && haskey(data, "items") && length(data["items"]) < num
					if isdir("movies/$(safe(title))")
						@show safe(title)
					end
					continue
				end

				@show safe(title)
				#throw(ex)
				i -= 1
				continue
			end
		end
	end

	@printf "Done :)\n"
end

main("../ml-latest/fixed/")


text:
using Revise
using ADCME
using PyCall
using LinearAlgebra
using AdFem
using PyPlot
using Random
Random.seed!(233)

# TODO: specify your input parameters
σ1 = rand(10, 3)
σ2 = zeros(10, 3)
ε1 = rand(10, 3)
ε2 = rand(10, 3)
μ = rand(10)
η = rand(10)
λ = rand(10)
Δt = rand()

u1 = update_stress_viscosity(ε2, ε1, σ1, η, μ, λ, Δt)

u2 = update_stress_viscosity(constant(ε2), ε1, σ1, η, μ, λ, Δt)

sess = Session(); init(sess)
@show run(sess, u2)-u1

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(update_stress_viscosity(ε2, ε1, σ1, η, μ, x, Δt)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(10))
v_ = rand(10)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")
savefig("gradtest.png")


text:
function filtered(filters::Vector{Function}, b...)
    np = length(b)
    [(filters[idx])(b[idx]) for idx in 1:np]
end
function filtered_fn_call(m::Function, filters::Vector{Function}, b...)
    filt = filtered(filters, b...)
    m(filt...)
end

function map(m::Base.Callable, bf::Block...)
    blks = [blocks(x) for x in bf]
    filters = [x.filter for x in bf]
    fc = (b...)->filtered_fn_call(m, filters, b...)
    map(fc, blks...)
end

function pmap(m, bf::Block...; kwargs...)
    affs = [x.affinity for x in bf]
    filters = [x.filter for x in bf]
    fc = (b...)->filtered_fn_call(m, filters, b...)
    # without any affinities, fall back to default pmap
    (sum([((x == no_affinity) ? 0 : 1) for x in affs]) == 0) && (return block_pmap(fc, [blocks(b) for b in bf]...; kwargs...))

    blks = [x.block for x in bf]
    prepare = [x.prepare for x in bf]

    n1 = length(affs[1])
    n = length(affs)

    # check if all parts are of same number of blocks
    for idx in 2:n
        (length(affs[idx]) != n1) && error("all parts must be of the same number of blocks")
    end

    afflist = [intersect([affs[x][j] for x in 1:n]...) for j in 1:n1]

    # check if at least one common processor for each block
    for idx in 1:n1
        isempty(afflist[idx]) && error("no common processor for block $(idx)")
    end

    donelist = falses(n1)
    np = nworkers()
    results = cell(n1)

    # function to produce the next work item from the queue.
    # in this case it's just an index.
    nextidx(p) = begin
        for (idx,a) in enumerate(afflist)
            if !donelist[idx] && (p in a)
                donelist[idx] = true
                return idx
            end
        end
        0
    end

    fetch_results = true
    for (n,v) in kwargs
        (n == :fetch_results) && (fetch_results = v)
    end
    @sync begin
        wrkrs = workers()
        for p=1:np
            wpid = wrkrs[p]
            @async begin
                while true
                    idx = nextidx(wpid)
                    (idx == 0) && break
                    #println("processor $(wpid) got idx $idx ... preparing data...")
                    d = [(prepare[bid])((blks[bid])[idx]) for bid in 1:length(blks)]
                    #println("processor $(wpid) preapred idx $idx")
                    results[idx] = fetch_results ? remotecall_fetch(fc, wpid, d...) : remotecall_wait(fc, wpid, d...)
                end
            end
        end
    end
    results
end

pmapreduce(m, r, bf::Block...) = reduce(r, pmap(m, bf...))
function pmapreduce(m, r, v0, bf::Block...)
    results = pmap(m, bf...)
    ret = v0
    for res in results
        ret = r(ret, res)
    end
    ret
    #reduce(r, v0, pmap(m, bf...))
end


text:
using StaticArrays
Base.@kwdef struct DenseSheeps{L} <: Sheeps{L}
    # # TODO: gridsize needed?
    capacity::Int = 1
    # available_sheeps::Base.RefValue{Int} = Ref(0)
    current_sheep::Base.RefValue{Int} = Ref(0)
    diffusion_candidates::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[]
    grid::SizedArray{Tuple{L,L},Int64,2,2} = SizedMatrix{L,L}( zeros(Int,L,L) )
end

function DenseSheeps{L}( total_sheep; kwargs...  ) where L
    sheep = DenseSheeps{L}(;kwargs...)
    for _ in 1:total_sheep
        i = rand(1:L)
        j = rand(1:L)
        n_sheep = getNSheep( sheep, i, j )
        setNSheep!( sheep, i, j, n_sheep + 1 )
    end # for
    sheep.current_sheep[] = total_sheep
    kickSheep!( sheep )
    return sheep
end # function

function DenseSheeps(gridsize; n_sheeps, kwargs...)
    return DenseSheeps{gridsize}(n_sheeps; kwargs...)
end # function

##< iterator interface
function Base.iterate( sheep::DenseSheeps{L}, state = firstindex(sheep.grid) ) where L
    for i in range(state,lastindex(sheep.grid), step = 1)
        n_sheep = sheep.grid[i]
        if n_sheep != 0
            return (n_sheep, CartesianIndices(sheep.grid)[i]), i+1
        end # if
    end # for
    return nothing
end # function

Base.eltype( sheep::DenseSheeps ) = Tuple{ Base.eltype( sheep.grid ), CartesianIndex{2} }
Base.length( sheep::DenseSheeps ) = Base.length( findall(!iszero, sheep.grid) )
##>

function _setNSheep!( sheep::DenseSheeps, i, j, n )
    old = sheep.grid[i,j]
    sheep.grid[i,j] = n
    if n > sheep.capacity
        push!( sheep.diffusion_candidates, (i,j) )
    end # if
    return old
end # function

function getNSheep( sheep::DenseSheeps, i, j )
    return sheep.grid[i,j]
end

"""
Returns random index of occupied grid-field.
"""
function getRandomSheep( sheep::DenseSheeps )
    rand( findall(!iszero, sheep.grid) )
end # function


text:
"""
    ncdim_def(ds, name, val, attrib = Dict())
    ncdim_def(ds, dim::NcDim)
    ncdim_def(ds, dims::Vector{NcDim})

# Examples
```julia
## 1. by existing nc file
# f = "/mnt/j/Researches/Obs/Obs_ChinaHW_cluster/ncell_16/clusterId_HItasmax_movTRS_Obs.nc"
dims = nc_dims(f)
ds = nc_open("f2.nc", "c")
ncdim_def(ds, dims);
println(ds)
close(ds)

## 2. by assigned variables
# dates = nc_date(ds)
ds = nc_open("f3.nc", "c")
ncdim_def(ds, "time",
    CFTime.timeencode.(dates, "days since 1970-01-01", eltype(dates)),
    Dict("units" => "days since 1970-01-01",
        "calendar" => "proleptic_gregorian", "long_name" => "time"))
close(ds)
```
"""
function ncdim_def(ds, name, val, attrib = Dict(); verbose = false)
    # x = NcDim(name, val, attrib)
    # val = val |> collect
    if name in keys(ds.dim)
        verbose && @warn "Dimension `$name`: exist!"
        return
    end
    defDim(ds, name, length(val))
    defVar(ds, name, val, (name,); attrib = attrib)
end

# NcDim
function ncdim_def(ds, dim::NcDim)
    ncdim_def(ds, dim.name, dim.vals, dim.atts)
end

function ncdim_def(ds, dims::Vector{NcDim})
    for i = 1:length(dims)
        dim = dims[i]
        ncdim_def(ds, dim.name, dim.vals, dim.atts)
    end
end


text:
module StatsFunctions

export logsumexp

function logsumexp(a, b)
    x = max(a,b)
    y = min(a,b)

    x + log1p(exp(y-x))
end

function logsumexp(xs::Array{T,N}) where {T <: Number, N}
    xm = maximum(xs)
    xm + log(sum(exp.(xs.-xm)))
end

end


text:
############################################################################################
"""
$(SIGNATURES)
Compute cross-chain diagnostics for parameter.

# Examples
```julia
```

"""
function chainparamdiagnostics(arr::Array{T,3}, computingtime::F) where {T<:Real,F<:Real}
    ess, rhat = MCMCDiagnosticTools.ess_rhat(arr)
    gelmandiagnostic = MCMCDiagnosticTools.gelmandiag(arr)
    vals = map(
        iter -> (
            ESS=ess[iter],
            ESSperSec=ess[iter] / computingtime,
            Rhat=rhat[iter],
            GelmanPSRF=gelmandiagnostic.psrf[iter],
        ),
        eachindex(ess),
    )
    return vals
end

function singlechaindiagnostics(arr::Array{T,3}, Nparams::Integer) where {T<:Real}
    vals = map(
        iter -> (
            ESS=NaN,
            ESSperSec=NaN,
            Rhat=NaN,
            GelmanPSRF=NaN),
        Base.OneTo(Nparams)
    )
    return vals
end

"""
$(SIGNATURES)
Compute in-chain diagnostics for parameter.

# Examples
```julia
```

"""
function paramdiagnostics(vec::AbstractVector{T}) where {T<:Real}
    vals = (
        Mean=Statistics.mean(vec),
        MCSE=MCMCDiagnosticTools.mcse(vec),
        StdDev=Statistics.std(vec),
    )
    return vals
end
function paramdiagnostics(arr::Array{T,3}) where {T<:Real}
    return map(
        iter -> paramdiagnostics(vec(view(arr, :, iter, :))), Base.OneTo(size(arr, 2))
    )
end

############################################################################################
"""
$(SIGNATURES)
Compute quantiles provided by `printdefault` for parameter.

# Examples
```julia
```

"""
function paramquantiles(vec::AbstractVector{T}, printdefault::PrintDefault) where {T<:Real}
    @unpack Ndigits, quantiles = printdefault
    vals = (Quantiles=Statistics.quantile(vec, quantiles),)
    return vals
end
function paramquantiles(arr::Array{T,3}, printdefault::PrintDefault) where {T<:Real}
    return map(
        iter -> paramquantiles(vec(view(arr, :, iter, :)), printdefault),
        Base.OneTo(size(arr, 2)),
    )
end

############################################################################################
"""
$(SIGNATURES)
Merge all statistics for all parameter.

# Examples
```julia
```

"""
function mergediagnostics(paramdiagnostic, chainparamdiagnostic, paramquantiles)
    Nparam = size(paramdiagnostic, 1)
    return map(
        param -> merge(
            paramdiagnostic[param], chainparamdiagnostic[param], paramquantiles[param]
        ),
        Base.OneTo(Nparam),
    )
end

############################################################################################
"""
$(SIGNATURES)
Return summary for trace parameter chains. `Model` defines flattening type of parameter,
`sym` defines parameter to be flattened, `backend` may be Val(:text), or Val(:latex).

# Examples
```julia
```

"""
function chainsummary(
    trace::Trace,
    model::ModelWrapper,
    sym::S,
    backend, #i.e., Val(:text), or Val(:latex)
    burnin::Integer,
    thinning::Integer,
    printdefault::PrintDefault=PrintDefault();
    kwargs...,
) where {S<:Union{Symbol,NTuple{k,Symbol} where k}}
    ## Assign utility values
    @unpack Ndigits, quantiles = printdefault
    @unpack default = model.info.reconstruct
    @unpack progress = trace.info
    tagged = Tagged(model, sym)
    Nparams = length(tagged)
    paramnames = ModelWrappers.paramnames(
        default, tagged.info.constraint, subset(model.val, tagged.parameter)
    )
    computingtime = progress.enabled ? (progress.tlast - progress.tinit) : NaN
    arr3D = trace_to_3DArray(trace, model, tagged, burnin, thinning)
    ## Compute summary statistics
    #!NOTE If more than 1 chain used, can use cross-chain diagnostics
    if trace.info.sampling.Nchains > 1
        #NOTE: Most of the computing time of the function happens here
        chainparamdiagnostic = chainparamdiagnostics(arr3D, computingtime)
    else
        chainparamdiagnostic = singlechaindiagnostics(arr3D, Nparams)
    end
    paramdiagnostic = paramdiagnostics(arr3D)
    paramquantile = paramquantiles(arr3D, printdefault)
    diag = mergediagnostics(paramdiagnostic, chainparamdiagnostic, paramquantile)
    ## Assign row and header variables for table
    statsnames = union(keys(paramdiagnostic[begin]), keys(chainparamdiagnostic[begin]))
    quantilenames = string.("Q", round.(quantiles .* 100; digits=1))
    tablenames = union(string.(statsnames), quantilenames)
    Nstats = length(tablenames)
    ## Create table
    _reconstruct = ModelWrappers.ReConstructor(diag)
    diag_flattened = flatten(_reconstruct, diag)
    table = round.(reshape(diag_flattened, Nstats, Nparams)'; digits=Ndigits)
    ## Print table
    PrettyTables.pretty_table(
        table, backend=backend, header=tablenames, row_names=paramnames, kwargs...
    )
    ## Return table arguments
    return table, tablenames, paramnames
end

############################################################################################
#export
export chainsummary


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

include("linearcode.jl")

abstract type AbstractHammingCode <: AbstractLinearCode end

function construct_ham_matrix(r::Int, q::Int)
    ncols = Int(floor((q^r - 1) / (q - 1)))
    M = Matrix{Int}(undef, r, ncols)

    for i in 1:ncols
        M[:, i] = reverse(digits(parse(Int, string(i, base = q)), pad = r), dims = 1)
    end

    return M
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
function test_factorization_machines()
    println("-- Testing Factorization Machines-based recommender")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    num_factors = 2
    learning_rate = 15e-4
    max_iter = 100

    user = 1
    topk = 4
    items = [i for i in 1:8]

    mf = MatrixFactorization(data, num_factors)
    build!(mf, learning_rate=learning_rate, max_iter=max_iter)
    mf_rec = recommend(mf, user, topk, items)

    fm = FactorizationMachines(data, num_factors)
    build!(fm, learning_rate=learning_rate, max_iter=max_iter)
    fm_rec = recommend(fm, user, topk, items)

    # if no contextual features are considered in FMs, the recommendation
    # result must be same as matrix factorization.
    @test [item for (item, score) in mf_rec] == [item for (item, score) in fm_rec]
end

function test_factorization_machines_with_attributes()
    println("-- Testing Factorization Machines-based recommender with contextual user/item attributes")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    set_user_attribute(data, 1, [1, 0, 20])
    set_user_attribute(data, 2, [0, 1, 23])
    set_user_attribute(data, 3, [0, 1, 12])

    set_item_attribute(data, 1, [1, 1, 0])
    set_item_attribute(data, 2, [0, 1, 1])
    set_item_attribute(data, 3, [1, 1, 1])
    set_item_attribute(data, 4, [1, 0, 0])
    set_item_attribute(data, 5, [1, 0, 1])
    set_item_attribute(data, 6, [0, 0, 1])
    set_item_attribute(data, 7, [0, 1, 0])
    set_item_attribute(data, 8, [1, 1, 1])

    recommender = FactorizationMachines(data, 2)
    build!(recommender, reg_w0=0.3, reg_w=0.3, reg_V=0.3, learning_rate=3e-10, random_init=true)

    rec = recommend(recommender, 1, 4, [i for i in 1:8])
    @test size(rec, 1) == 4  # top-4 recos
end

test_factorization_machines()
test_factorization_machines_with_attributes()


text:
#Thin QR implementation

using LinearAlgebra

"""
Allocate matrices and vectors for the QR factorization of A

# Arguments
- `A::Array{Float64,2}`: matrix of size m x n such that m < n

# Return

- `QR::Array{Float64,2}`: matrix of size m+1 x n
- `v::Array{Float64,1}`: vector of size m
- `u::Array{Float64,1}`: vector of size m

"""
function allocate_matrices(A)
    (m,n) = size(A)

    if m < n
        error("m must be greater or equal than n")
    end

    QR = zeros(m+1,n)
    v = zeros(m)
    u = zeros(m)
    return QR,v,u
end

"""
Performs QR factorization of A

# Arguments
- `A::Array{Float64,2}`: matrix of size m x n such that m < n.
- `QR::Array{Float64,2}`: matrix of size m+1 x n. The result of the factorization is stored in QR such that:
    - the upper triangular part contains R. Therefore, R = triu(QR)[1:n,:]
    - the sub-diagonal part contains the reflectors v. Hence, V = tril(QR,-1)[2:n+1,:]
- `v::Array{Float64,1}`: vector of size m
- `u::Array{Float64,1}`: vector of size m
"""
function qr_factorization!(A, QR::Array{Float64,2}, v::Array{Float64,1}, u::Array{Float64,1} )
    (m,n) = size(A)

    #copying elements of A in R
    for j = 1 : n
        for i = 1 : m
            QR[i,j] = A[i,j]
        end
    end

    #Total Complexity: O(n(3m + 2mn)) = O(3nm + 2mn^2) = O(2n^2m) + O(nm)
    @inbounds @views for j = 1 : min(m-1,n)

        #copying j-th column of R into v
        #Complexity: O(m)
        for i = 1 : m
            v[i] = i < j ? 0 : QR[i,j]
        end

        #calculating householders
        s = norm(v)

        if v[j] >= 0
            s = -s
        end

        v[j] = v[j] - s

        norm_v = norm(v)

        #Complexity: O(m)
        @. v = v/norm_v

        #calculating R
        # R = R - 2v*(v'*R)

        t = j+1

        QR[j,j] = s
        @. QR[t:m, j] = 0

        #Complexity: O(mn)
        u[t:n]' .= v[j:m]' * QR[j:m, t:n]

        #Complexity: O(mn)
        @. QR[j:m, t:n] = QR[j:m, t:n] - (2*v[j:m])*u[t:n]'

        #copying householder reflector into QR matrix
        @. QR[j+1:m+1, j] = v[j : m]
    end

    #if m== n then QR[m+1,n]=0
    if m == n
        QR[m+1,n] = 0.0
    end

end

"""
Performs Q^t * A and stores the result in W (Q^t is Q transposed)

# Arguments
- `QR::Array{Float64,2}`: matrix of size m+1 x n. The result of the factorization is stored in QR such that:
    - the upper triangular part contains R. Therefore, R = triu(QR)[1:n,:]
    - the sub-diagonal part contains the reflectors v. Hence, V = tril(QR,-1)[2:n+1,:]
- `A::Array{Float64,2}`: matrix of size m x n.
- `W::Array{Float64,2}`:matrix of size m x n in which the result is stored.
"""
function Q_t_times_A!(QR,A,W)
    #Taking householder vectors
    @views V = tril(QR, -1)[2:end,:]
    (m,n) = size(V)

    W .= 0
    W .= A .- 2 .* V[:,1] .* (V[:,1]'*A)

    #Complexity: O(mn^2)
    num_iteration = n
    if V[m,n] == 0.0
        num_iteration -= 1
    end

    for j = 2:num_iteration
        W .= W .- 2 .* V[:,j] .* (V[:,j]'*W)
    end
end

"""
Given the QR factorization returns the Q factor

# Arguments
- `QR::Array{Float64,2}`: matrix of size m+1 x n. The result of the factorization is stored in QR such that:
    - the upper triangular part contains R. Therefore, R = triu(QR)[1:n,:]
    - the sub-diagonal part contains the reflectors v. Hence, V = tril(QR,-1)[2:n+1,:]
- `Q::Array{Float64,2}`:matrix of size m x n in which the result is stored.
"""
function get_Q!(QR, Q)
    @views V = tril(QR, -1)[2:end,:]
    (m,n) = size(V)

    #initialize Q
    @views @inbounds Q .= (I - 2 .* V[:,n] .* V[:,n]')[:,1:n]

    @views @inbounds for j in n-1:-1:1
        Q[j:end,j:end] .= (I - 2 .* V[:,j] .* V[:,j]')[j:end,j:end] * Q[j:end,j:end]
    end
end


text:
# Because there exists the `FiniteDiff.jl` package with quite different purposes,
# this module is not expected to be reexported.
module FiniteDiff

using ImageCore
using ImageCore.MappedArrays: of_eltype

"""
Although stored as an array, image can also be viewed as a function from discrete grid space
Zᴺ to continuous space R if it is gray image, to C if it is complex-valued image
(MRI rawdata), to Rᴺ if it is colorant image, etc.
This module provides the discrete
version of gradient-related operators by viewing image arrays as functions.

This module provides:

- forward/backward difference [`fdiff`](@ref) are the Images-flavor of `Base.diff`
- gradient operator [`fgradient`](@ref) and its adjoint via keyword `adjoint=true`.
- divergence operator [`fdiv`](@ref) computes the sum of discrete derivatives of vector
  fields.
- laplacian operator [`flaplacian`](@ref) is the divergence of the gradient fields.

Every function in this module has its in-place version.
"""
FiniteDiff

export
    fdiff, fdiff!,
    fdiv, fdiv!,
    fgradient, fgradient!,
    flaplacian, flaplacian!

"""
    fdiff(A::AbstractArray; dims::Int, rev=false, boundary=:periodic)

A one-dimension finite difference operator on array `A`. Unlike `Base.diff`, this function doesn't
shrink the array size.

Take vector as an example, it computes `(A[2]-A[1], A[3]-A[2], ..., A[1]-A[end])`.

# Keywords

- `rev::Bool`
  If `rev==true`, then it computes the backward difference
  `(A[end]-A[1], A[1]-A[2], ..., A[end-1]-A[end])`.
- `boundary`
  By default it computes periodically in the boundary, i.e., `:periodic`.
  In some cases, one can fill zero values with `boundary=:zero`.

# Examples

```jldoctest; setup=:(using ImageBase.FiniteDiff: fdiff)
julia> A = [2 4 8; 3 9 27; 4 16 64]
3×3 $(Matrix{Int}):
 2   4   8
 3   9  27
 4  16  64

julia> diff(A, dims=2) # this function exists in Base
3×2 $(Matrix{Int}):
  2   4
  6  18
 12  48

julia> fdiff(A, dims=2)
3×3 $(Matrix{Int}):
  2   4   -6
  6  18  -24
 12  48  -60

julia> fdiff(A, dims=2, rev=true) # reverse diff
3×3 $(Matrix{Int}):
  -6   2   4
 -24   6  18
 -60  12  48

julia> fdiff(A, dims=2, boundary=:zero) # fill boundary with zeros
3×3 $(Matrix{Int}):
  2   4  0
  6  18  0
 12  48  0
```

See also [`fdiff!`](@ref) for the in-place version.
"""
fdiff(A::AbstractArray; kwargs...) = fdiff!(similar(A, maybe_floattype(eltype(A))), A; kwargs...)

"""
    fdiff!(dst::AbstractArray, src::AbstractArray; dims::Int, rev=false, boundary=:periodic)

The in-place version of [`fdiff`](@ref)
"""
function fdiff!(dst::AbstractArray, src::AbstractArray;
        dims=_fdiff_default_dims(src),
        rev=false,
        boundary::Symbol=:periodic)
    isnothing(dims) && throw(UndefKeywordError(:dims))
    axes(dst) == axes(src) || throw(ArgumentError("axes of all input arrays should be equal. Instead they are $(axes(dst)) and $(axes(src))."))
    N = ndims(src)
    1 <= dims <= N || throw(ArgumentError("dimension $dims out of range (1:$N)"))

    src = of_eltype(maybe_floattype(eltype(dst)), src)
    r = axes(src)
    r0 = ntuple(i -> i == dims ? UnitRange(first(r[i]), last(r[i]) - 1) : UnitRange(r[i]), N)
    r1 = ntuple(i -> i == dims ? UnitRange(first(r[i])+1, last(r[i])) : UnitRange(r[i]), N)

    d0 = ntuple(i -> i == dims ? UnitRange(last(r[i]), last(r[i])) : UnitRange(r[i]), N)
    d1 = ntuple(i -> i == dims ? UnitRange(first(r[i]), first(r[i])) : UnitRange(r[i]), N)

    if rev
        dst[r1...] .= view(src, r1...) .- view(src, r0...)
        if boundary == :periodic
            dst[d1...] .= view(src, d1...) .- view(src, d0...)
        elseif boundary == :zero
            dst[d1...] .= zero(eltype(dst))
        else
            throw(ArgumentError("Wrong boundary condition $boundary"))
        end
    else
        dst[r0...] .= view(src, r1...) .- view(src, r0...)
        if boundary == :periodic
            dst[d0...] .= view(src, d1...) .- view(src, d0...)
        elseif boundary == :zero
            dst[d0...] .= zero(eltype(dst))
        else
            throw(ArgumentError("Wrong boundary condition $boundary"))
        end
    end

    return dst
end

_fdiff_default_dims(A) = nothing
_fdiff_default_dims(A::AbstractVector) = 1

maybe_floattype(::Type{T}) where T = T
maybe_floattype(::Type{T}) where T<:FixedPoint = floattype(T)
maybe_floattype(::Type{CT}) where CT<:Color = base_color_type(CT){maybe_floattype(eltype(CT))}

"""
    fdiv(Vs...)

Compute the divergence of vector fields `Vs`.

See also [`fdiv!`](@ref) for the in-place version.
"""
function fdiv(V₁::AbstractArray{T}, Vs::AbstractArray{T}...) where T
    fdiv!(similar(V₁, maybe_floattype(T)), V₁, Vs...)
end
fdiv(Vs::Tuple) = fdiv(Vs...)

"""
    fdiv!(out, Vs...)

The in-place version of divergence operator [`fdiv`](@ref).
"""
function fdiv!(out::AbstractArray, V₁::AbstractArray, Vs::AbstractArray...)
    # This is the optimized version of `sum(v->fgradient(v; ajoint=true), (V₁, Vs...))`
    # by removing unnecessary memory allocations.
    all(v->axes(v) == axes(out), (V₁, Vs...)) || throw(ArgumentError("All axes of vector fields Vs and X should be the same."))

    # TODO(johnnychen94): for better performance, we can eliminate this `tmp` allocation by fusing multiple `fdiff` in the inner loop.
    out .= fdiff(V₁; dims=1, rev=true, boundary=:periodic)
    tmp = similar(out)
    for i in 1:length(Vs)
        out .+= fdiff!(tmp, Vs[i]; dims=i+1, rev=true, boundary=:periodic)
    end
    return out
end
fdiv!(out::AbstractArray, Vs::Tuple) = fdiv!(out, Vs...)

"""
    flaplacian(X::AbstractArray)

The discrete laplacian operator, i.e., the divergence of the gradient fields of `X`.

See also [`flaplacian!`](@ref) for the in-place version.
"""
flaplacian(X::AbstractArray) = flaplacian!(similar(X, maybe_floattype(eltype(X))), X)

"""
    flaplacian!(out, X)
    flaplacian!(out, ∇X::Tuple, X)

The in-place version of the laplacian operator [`flaplacian`](@ref).

!!! tip Avoiding allocations
    The two-argument method will allocate memory to store the intermediate
    gradient fields `∇X`. If you call this repeatedly with images of consistent size and type,
    consider using the three-argument form with pre-allocated memory for `∇X`,
    which will eliminate allocation by this function.
"""
flaplacian!(out, X::AbstractArray) = fdiv!(out, fgradient(X))
flaplacian!(out, ∇X::Tuple, X::AbstractArray) = fdiv!(out, fgradient!(∇X, X))

"""
    fgradient(X::AbstractArray; adjoint=false, boundary=:periodic) -> (∂₁X, ∂₂X, ..., ∂ₙX)

Computes the gradient fields of `X`. If `adjoint==true` then it computes the adjoint gradient
fields.

Each gradient vector is computed as forward difference along specific dimension, e.g.,
[`∂ᵢX = fdiff(X, dims=i)`](@ref fdiff). The `boundary` keyword is passed to `fdiff` to
specify the behavior on boundary.

Mathematically, the adjoint operator ∂ᵢ' of ∂ᵢ is defined as `<∂ᵢu, v> := <u, ∂ᵢ'v>`, where
`<a, b>` is the inner product of vector `a` and `b`.

See also the in-place version [`fgradient!(X)`](@ref) to reuse the allocated memory.
"""
function fgradient(X::AbstractArray{T,N}; adjoint::Bool=false, boundary=:periodic) where {T,N}
    fgradient!(ntuple(i->similar(X, maybe_floattype(T)), N), X; adjoint=adjoint, boundary=boundary)
end

"""
    fgradient!(∇X::Tuple, X::AbstractArray; adjoint=false, boundary=:periodic)

The in-place version of (adjoint) gradient operator [`fgradient`](@ref).

The input `∇X = (∂₁X, ∂₂X, ..., ∂ₙX)` is a tuple of arrays that are similar to `X`, i.e.,
`eltype(∂ᵢX) == eltype(X)` and `axes(∂ᵢX)  == axes(X)` for all `i`.
"""
function fgradient!(∇X::NTuple{N, <:AbstractArray}, X; adjoint::Bool=false, boundary=:periodic) where N
    all(v->axes(v) == axes(X), ∇X) || throw(ArgumentError("All axes of vector fields ∇X and X should be the same."))
    for i in 1:N
        if adjoint
            # the negative adjoint of gradient operator for forward difference is the backward difference
            # see also
            # Getreuer, Pascal. "Rudin-Osher-Fatemi total variation denoising using split Bregman." _Image Processing On Line_ 2 (2012): 74-95.
            fdiff!(∇X[i], X, dims=i, rev=true, boundary=boundary)
            # TODO(johnnychen94): ideally we can get avoid flipping the signs for better performance.
            @. ∇X[i] = -∇X[i]
        else
            fdiff!(∇X[i], X, dims=i, boundary=boundary)
        end
    end
    return ∇X
end

if VERSION < v"1.1"
    isnothing(x) = x === nothing
end

end # module


text:
module test_julia_hashing

using Test

# hash
# == isequal
@test hash(2) == hash(2)
@test !((==) isa Core.Builtin)

# objectid
# ===
@test objectid(2) == objectid(2) == ccall(:jl_object_id, UInt, (Any,), 2)
@test (===) isa Core.Builtin

end # module test_julia_hashing


text:
### Packages used
using Statistics
using NaNStatistics
using Optim
using DelimitedFiles

### Structs used
struct hyperParameters
    sampleRate              ::Int64          # [Hz] Sample rate.
    unloadingFitRange       ::Int64          # [-] Vector samples lengths (measured from start of unloading) to be included in unload fit.
    unloadingFitFunction    ::String         # [string] Function to use when fitting.
    compensateCreep         ::Bool           # [Bool] Compensate creep using Feng's method, y/n
    constrainHead           ::Int            # Experimental, not implemented
    constrainTail           ::Int            # Experimental, not implemented
    machineCompliance       ::Float64        # Machine compliance
end

struct control
    plotMode                ::Bool           # Activates plotting of intermediate results
    verboseMode             ::Bool           # Verbose output
end

struct metaInfoExperimentalSeries
    designatedName          ::String         # Designated name
    relativeHumidity        ::Float64        # Relative humidity
    indenterType            ::String
    indentationNormal       ::String
    springConstant          ::Float64
    areaFile                ::String
    targetDir               ::String
    thermalHoldTime         ::Int64          # Should be changed to something of unit "second".
end

"""
ReadIBW(filename)

Reads a .IBW binary file, extracting the position and force signal of the AFM indenter.

filename is the string to a file which either needs to be absolute or on the path.
"""
function readIBW(filename)

    A = open( filename, "r")
    firstByte = read(A, Int8)
    close(A)

    if firstByte == 0
        machineFormat = 'b'
    else
        machineFormat = 'l'
    end

    A = open( filename, "r")
    version         = read(A, Int16)
    checksum        = ltoh(read(A,UInt16)) # Does not work
    wfmSize         = ltoh(read(A, Int32))
    formulaSize     = ltoh(read(A, Int32))
    noteSize        = ltoh(read(A, Int32)) # Does not work
    dataEUnitsSize  = ltoh(read(A, Int32))
    dimEUnitsSize   = Array{Int32, 1}(undef, 4); read!(ltoh(A),dimEUnitsSize)
    dimLabelsSize   = Array{Int32, 1}(undef, 4); read!(ltoh(A),dimLabelsSize)
    sIndicesSize    = read(A, Int32)
    optionSize1     = read(A, Int32)
    optionSize2     = read(A, Int32)
    ignore          = read(A, Int32)
    CreationDate    = read(A, UInt32)
    modData         = read(A, UInt32)
    npnts           = read(A, Int32)                  # WORKS!

    type = read(A, Int16)
    if type == 2
        datatype = "single"
    elseif type == 4
        datatype = "double"
    else
        println("ERROR")
    end

    read(A, Int16)
    for aLoop = 1:6
        read(A, Char)
    end
    read(A, Int16)

    for aLoop = 1:32
        read(A, Char)
    end
    for aLoop = 1:6
        read(A, Int32)
    end

    sfA = Array{Float64, 1}(undef, 4)
    read!(A, sfA)
    sfB = Array{Float64, 1}(undef, 4)
    read!(A, sfB)

    dUnits = Array{Char, 1}(undef,4)
    for aLoop = 1:4
        dUnits[aLoop] = read(A, Char)
    end
    dUnits = join(dUnits)

    xUnits  = Array{Char, 1}(undef, 16)
    for aLoop = 1:16
        xUnits[aLoop] = read(A, Char)
    end
    xUnits = join(xUnits)

    fsValid = read(A, Int16)
    whpad3 = read(A, Int16)
    read(A, Float64)
    read(A, Float64)
    for aLoop = 1:26
        read(A, Int32)
    end
    for aLoop = 1:3
        read(A, Int16)
    end
    read(A, Char)
    read(A, Char)
    read(A, Int32)
    read(A, Int16)
    read(A, Int16)
    read(A, Int32)
    modDate = read(A, Int32)
    creationDate = read(A, Int32)
    wdata = Array{Float32, 1}(undef, npnts); read!(A,wdata)

    close(A)

    x0 = Float64(sfB[1]*fsValid)
    dx = Float64(sfA[1]*fsValid)
    return wdata#, npnts# , dUnits, x0, dx, xUnits

end

#data, dUnits, x0, dx, xUnits = ReadIBW(filename)
"""
IBWtoTXT(filename)

Takes a filename, calls ReadIBW and formats the returned data into x and y signals where
x is the position of the sensor and y is the deflection of the sensor.
"""
function IBWtoTXT(filename::String)
    data = readIBW(filename)
    lengthOfData = length(data) ÷ 3

    #t = data[1:lengthOfData]
    y = data[lengthOfData+1:2*lengthOfData]
    x = data[2*lengthOfData+1:3*lengthOfData]
    return [x y]
end

function dataPreProcessing(filename::String)
    xy0 = IBWtoTXT(filename)
    # Import signal
    xy0 .*= 1e9
    # Convert to nano-meters

    return xy0
end

function offsetAndDriftCompensation(xy::Matrix{Float32})
    endRangeBaseFit = 25;
    breakForContact = true;
    aLoop = 0;
    basefit = 0

    while breakForContact
        aLoop += 1

        sIdx = min(10+endRangeBaseFit*(aLoop-1), size(xy,1)-endRangeBaseFit);
        eIdx = min(10+endRangeBaseFit*(aLoop), size(xy,1)-1);

        dPdt = 100*(xy[eIdx,2] - xy[sIdx,2]) / endRangeBaseFit;

        if dPdt > 3.3  && aLoop > 3
            basefit = endRangeBaseFit*(aLoop-1);
            breakForContact = false;
            #println("if")
        elseif endRangeBaseFit*aLoop > 1000
            basefit = endRangeBaseFit*7;
            breakForContact = false;
            #println("elseif")
        end
    end

    zeroLineM = ones(length(xy[10:basefit,1])) \ xy[10:basefit,2]
    xy[:,2] .-= zeroLineM

    # Locate the point where the ramp begins.
    noise = max(0.1, std(xy[10:basefit,2]) );
    rampStartIdx = max(1, findfirst(x -> x > 4.0*noise, xy[:,2]) - 1);
    xZero = xy[rampStartIdx,1];
    xy[:,1] .-= xZero;                                       # Shift the curve by the value at the start of the ramp.
    xy[:,2] .-= xy[rampStartIdx,2];

    zeroLineD = xy[10:basefit,1] \ xy[10:basefit,2];

    zeroLine = [zeroLineD zeroLineM];

    return xy , basefit , zeroLine , rampStartIdx, xZero
end

function subdirImport(placeToLook::String,stringToMatch::String)
    return filter(x-> contains(x, stringToMatch) , readdir(placeToLook))
end

function findStartOfHold(xy::Matrix{Float32}, directionOfSearch::String)
    # 1. Determine range of deflection values.
    # 2. Bin the values (heuristic bin size at the moment)
    # 3. Determine the most common deflection value at the highest load levels under the assumption that this bin will contain
    #    the hold sequence.
    # 4. Determine the mean value of all values larger than this bin value.
    # 5. Find the first time the vector exceeds this value.
    # 6. This is taken as the first value in the hold sequence.

    sensorRange = maximum(xy[:,2]) - minimum(xy[:,2])
    vecLengthTemp = Int64(round(4.0*sensorRange))
    edgesOfHist = range(minimum(xy[:,2]), maximum(xy[:,2]), length = vecLengthTemp)
    histTemp = histcounts(xy[:,2] , edgesOfHist)

    tailVecStart = Int64(round(0.9*vecLengthTemp))

    peakIdx = argmax(histTemp[tailVecStart:end])
    peakIdx += tailVecStart-1

    idxTemp = map(x -> x > edgesOfHist[peakIdx], xy[:,2])

    meanOfPlateau = mean(xy[idxTemp,2])

    if cmp(directionOfSearch,"first") == 0
        returnIdx = findfirst(x -> x ≥ meanOfPlateau, xy[:,2])
    elseif cmp(directionOfSearch,"last") == 0
        returnIdx = findlast(x -> x ≥ meanOfPlateau, xy[:,2])
    end

    return returnIdx

end

function determineThermalCreep(xy::Matrix{Float32},sampleRate::Int64,thermalHoldTime::Int64,ctrl::control)
    sensorRange = maximum(xy[:,2]) - minimum(xy[:,2])
    vecLengthTemp = Int64(round(0.1*sensorRange));
    edgesOfHist = range(minimum(xy[:,2]), maximum(xy[:,2]), length = vecLengthTemp)
    peakIdx = argmax(histcounts(xy[:,2] , edgesOfHist))

    idxTemp = (xy[:,2] .< edgesOfHist[peakIdx+1]) .& (xy[:,2] .> edgesOfHist[max(1,peakIdx-1)])
    meanOfPlateau = mean(xy[idxTemp,2])
    stdOfPlateau = std(xy[idxTemp,2])

    noiseMultiplier = 5.0; # 15 # OBS HARD CODED AND SHOULD BE PULLED OUT!
    thermalHoldStartIdx = findlast(x -> x > meanOfPlateau+noiseMultiplier*stdOfPlateau, xy[:,2])
    thermalHoldEndIdx = findlast(x -> x > meanOfPlateau-noiseMultiplier*stdOfPlateau, xy[:,2])

    thermalHoldStartIdx += sampleRate
    thermalHoldEndIdx -= sampleRate

    # Ensure that the thermal hold sequence contains at least 25 seconds (out of the 30 secounds
    # specified).
    if thermalHoldStartIdx > thermalHoldEndIdx
        ctrl.verboseMode && println("Missed thermal hold. Increasing search range.")

        while thermalHoldTime+thermalHoldStartIdx > thermalHoldEndIdx
            noiseMultiplier += 1.0
            thermalHoldStartIdx = findlast(x -> x > meanOfPlateau+noiseMultiplier*stdOfPlateau, xy[:,2])
            thermalHoldEndIdx = findlast(x -> x > meanOfPlateau-noiseMultiplier*stdOfPlateau, xy[:,2])

            try
            thermalHoldStartIdx += sampleRate
            thermalHoldEndIdx -= sampleRate
            catch
                return 0.0
                println("Failure in the termal hold calculation")
            end
        end
        ctrl.verboseMode && println("Thermal hold found using multiplier $noiseMultiplier")
    end

    # Fit a function of displacement (due to thermal fluctuation)
    # h_thermal(time) = A1 + A2*time^A3
    thermalHoldDisplacement = xy[thermalHoldStartIdx:thermalHoldEndIdx,1];
    thermalHoldTime = collect(1:length(thermalHoldDisplacement))./sampleRate;

    deltaDisp = (xy[thermalHoldEndIdx,1] - xy[thermalHoldStartIdx,1])

    thermalCreepFun(x) = xy[thermalHoldStartIdx,1] .+ x[1].*(thermalHoldTime)#.^x[2]
    thermalHoldMinFun(x) = sqrt(sum((   (thermalCreepFun(x) .- thermalHoldDisplacement)./thermalHoldDisplacement ).^2))
    result = optimize(thermalHoldMinFun, [deltaDisp], BFGS())
    thermal_p = result.minimizer

    # Estimate the thermal drift rate by taking the median of the differentiated h_thermal
    # dhtdt = d(h_thermal(time))/d(time)
    # The functional form accounts for any viscous effects lingering from the unloading at the start
    # of the thermal hold, while the median provides a roboust average of the thermal drift rate.
    #dhtdt = median(thermal_p[1] .* thermal_p[2].*thermalHoldTime.^(thermal_p[2] - 1));
    dhtdt = thermal_p[1];
end

function determineCreepDuringHold(xy_hold,sampleRate::Int64)
    # Fit the creep during the hold time at maximum load to a linear spring-dashpot model
    #
    # Equation (17c) - not labeled in [1]
    #
    # h(t) = h_i + \beta * t.^(1/3)
    #
    #   h(t)    - Displacement as a function of time during hold sequence.
    #   h_i     - Fitting constant
    #   \beta   - Fitting constant
    #   t       - Time

    holdTimeVals = collect(1:length(xy_hold[:,1])) ./ sampleRate
    # Generate time signal

    deltaDisp = (xy_hold[end,1] - xy_hold[1,1])

    # Define fitting functions
    hOftFun(x) = xy_hold[1,1] .+ x[1].*holdTimeVals.^x[2]
    minFcn(x) = sqrt( sum( ( (hOftFun(x) - xy_hold[:,1])./ xy_hold[:,1] ).^2 ) )
    holdDrift = optimize(minFcn, [deltaDisp 0.33], BFGS())
    crp_p = holdDrift.minimizer

    h_dot_tot = crp_p[1] * crp_p[2]*holdTimeVals[end]^(crp_p[2] - 1.0)

    return max(h_dot_tot, 0.0)
end

function modulusfitter(indentationSet::metaInfoExperimentalSeries,hyperParameters,ctrl::control,resultFile::String)
    xy = dataPreProcessing(indentationSet.targetDir*resultFile)
    # Load raw displacements

    #ctrl.plotMode && display(plot([xy[1:100:end,1]],[xy[1:100:end,2]]))

    xy , basefit , zeroLine , rampStartIdx, xZero  = offsetAndDriftCompensation(xy)
    # Find initial contact
    #ctrl.plotMode && display(plot!(xy[1:100:end,1],xy[1:100:end,2]))

    xy[:,1] .-= xy[:,2]
    xy[:,2] .*= indentationSet.springConstant
    xy[:,1] .-= hyperParameters.machineCompliance.*xy[:,2];
    # Convert displacement-deflection matrix to indentation-force matrix

    #################################################################
    # Determine the start of the hold time at circa max force.
    #
    # 1. Determine range of deflection values.
    # 2. Bin the values (heuristic bin size at the moment)
    # 3. Determine the most common deflection value at the highest load levels under the assumption that this bin will contain
    #    the hold sequence.
    # 4. Determine the mean value of all values larger than this bin value.
    # 5. Find the first time the vector exceeds this value.
    # 6. This is taken as the first value in the hold sequence.
    holdStartIdx = findStartOfHold(xy,"first")
    ctrl.plotMode && display(plot(xy[:,1], xy[:,2], xlims = (0.0, maximum(xy[:,1])), xlab = "Indentation [nm]", ylab = "Force [uN]", legend = false))
    # ctrl.plotMode && display(plot!([xy[holdStartIdx,1]], [xy[holdStartIdx,2]],
    #                          seriestype = :scatter, lab = "Start of hold", legend = :topleft))

    # Split into loading and unloading.
    #xy_load = xy[1:holdStartIdx,:];
    xy_unld1 = xy[holdStartIdx:end,:];

    #Determine the end of the hold time.
    unloadStartIdx = findStartOfHold(xy_unld1,"last")
    #ctrl.plotMode && display(plot!([xy_unld1[unloadStartIdx,1]], [xy_unld1[unloadStartIdx,2]],seriestype = :scatter))

    # Split into two new pieces
    xy_hold = xy_unld1[1:unloadStartIdx-1,:];
    xy_unld = xy_unld1[unloadStartIdx:end,:];

    # Accept only indentations that had positive creep. "Negative" creep (indenter moves outwards
    # during hold sequence) can occur if the thermal drift is substantial, but this typically
    # indicates that the system was not in equilibrium (since the thermal drift dominates the creep)
    # and furthermore it messes up the mathematical framework if you accept such indentations (see
    # Cheng & Cheng articles.)
    condition1 = xy[holdStartIdx,1] < xy_unld1[unloadStartIdx,1]

    # Accept only monotonously increasing load-displacement curves. A curve may show weird behaviour
    # and our solution is to simply drop the curve in that case.
    condition2 = minimum(xy[rampStartIdx:holdStartIdx,1]) > (0.0-eps())
    condition3 = maximum(xy_unld[:,1]) < (xy_unld[1,1]+0.5)

    ctrl.plotMode && display(title!(string(condition3)))
    sleep(1.0)

    if condition1 && condition2 && condition3
        xy_unld5 = xy_unld[1:Int64(round(2000*0.95)),:];  # OBS 2000 is hard coded!
        # Make sure that the thermal hold sequence is not included in the unloading curve.

        dhtdt = determineThermalCreep(xy,hyperParameters.sampleRate,indentationSet.thermalHoldTime,ctrl)
        #dhtdt == 0.0 && return 0.0

        # Fitting of the unloading curve.
        stiffness_fit = Array{Float64}(undef,1)
        dispVals = xy_unld5[1:hyperParameters.unloadingFitRange,1]
        forceVals = xy_unld5[1:hyperParameters.unloadingFitRange,2]
        Fmax = xy_unld5[1,2]               # Maximum force during unloading

        if cmp(hyperParameters.unloadingFitFunction,"Oliver-Pharr") == 0
            Dmax = xy_unld5[1,1]
            # Maximum indentation depth during unloading

            function unloadFitFun(fitCoefs)
                return fitCoefs[1].*(dispVals .- fitCoefs[2]).^fitCoefs[3] .- forceVals
            end
            function unloadFitMinFun(fitCoefs)
                sqrt( sum( (unloadFitFun(fitCoefs) ./ forceVals).^2 ) )
            end

            lx = [0.0, 0.0 , 0.0]; ux = [Inf, minimum(dispVals)-1e-2 , Inf];
            dfc = TwiceDifferentiableConstraints(lx, ux)
            resultFit = optimize(unloadFitMinFun, dfc, [1.0, 1.0, 1.0], IPNewton())
            uld_p = resultFit.minimizer
            stiffness_fit = uld_p[1]*uld_p[3]*(Dmax - uld_p[2]).^(uld_p[3] - 1)

        elseif cmp(hyperParameters.unloadingFitFunction, "Feng") == 0

            unloadFitFun2(fitCoefs) = fitCoefs[1] .+ fitCoefs[2].*forceVals.^0.5 + fitCoefs[3].*forceVals.^fitCoefs[4] .- dispVals
            unloadFitMinFun2(fitCoefs) = sqrt(sum( (unloadFitFun2(fitCoefs) ./ dispVals).^2) )
            resultFit = optimize(unloadFitMinFun2, [1.0 1.0 1.0 1.0], BFGS())
            uld_p = resultFit.minimizer
            stiffness_fit = inv(( 0.5*uld_p[2].*Fmax.^-0.5 + uld_p[4]*uld_p[3].*Fmax.^(uld_p[4] - 1.0) ))

        end

        h_dot_tot = determineCreepDuringHold(xy_hold,hyperParameters.sampleRate)
        dPdt = [1/hyperParameters.sampleRate .* collect(0:(length(xy_unld5[:,1])-1)) ones(length(xy_unld5[:,1]))] \ xy_unld5[:,2]

        if hyperParameters.compensateCreep
            stiffness = inv(1/stiffness_fit + h_dot_tot/(abs(dPdt[1])));
        else
            stiffness = stiffness_fit;
        end

        # Equation (2) in [1]
        maxIndentation = median(xy_unld5[1,1]) - dhtdt*(length(xy[rampStartIdx:holdStartIdx,1])+length(xy_hold[:,1]))/hyperParameters.sampleRate;  #%OBS OBS OBS

        if cmp(indentationSet.indenterType,"pyramid") == 0
            x0 = maxIndentation - 0.72*Fmax/stiffness;
        elseif cmp(indentationSet.indenterType,"hemisphere") == 0
            x0 = maxIndentation - 0.75*Fmax/stiffness;
        end
        x0 < 0.0 && return 0.0

        area_xy = readdlm(indentationSet.areaFile, ' ', Float64, '\n')
        # % Determine the area by loading the calibration data and fitting a polynom to the data.

        if (x0 > 100.0)
            area_fit_end = length(area_xy[:,1])
        elseif (x0 < 100.0)
            area_fit_end = findfirst( x -> x > 100,area_xy[:,1])
        end

        tempVec = area_xy[1:area_fit_end,1]
        p_area = [tempVec.^2 tempVec tempVec.^0.5 tempVec.^0.25 tempVec.^0.125] \ area_xy[1:area_fit_end,2]

        unloadArea = [x0^2 x0 x0^0.5 x0^0.25 x0^0.125] * p_area
        unloadArea = unloadArea[1]
        unloadArea < 0.0 && return 0.0

        # % Equation (1) in [1]
        Er = sqrt(pi)/(2.0)/sqrt(unloadArea) / (1.0/stiffness )
        if cmp(indentationSet.indenterType,"pyramid") == 0
            Er = Er/1.05;
        end
        return Er
    else
        return 0.0
    end

    println(Er)
    #return Er
end

function calculateMachineCompliance(indentationSet::metaInfoExperimentalSeries,hyperParameters,ctrl::control)
    resultNames = subdirImport(indentationSet.targetDir,".ibw");     # Find the .ibw files in targetDir
    resultNames = resultNames[6:end]
    ap = []
    bp = []
    for file in resultNames
        println(file)
        tempS , tempA = extractSingleComplianceExperiment(indentationSet,hyperParameters,ctrl,file)
        push!(ap,tempS)
        push!(bp,tempA)
        #println([tempS tempA])
    end

    squaredInverseArea = 1.0 ./sqrt.(bp)
    #return squaredInverseArea
    #return ones(size(ap))
    return [squaredInverseArea[:] ap]

    effectiveCompliance = [squaredInverseArea[:] ones(size(ap))] \ ap;
    return effectiveCompliance

    return effectiveCompliance[2]
end

function extractSingleComplianceExperiment(indentationSet::metaInfoExperimentalSeries,hyperParameters,ctrl::control,resultFile::String)

    xy = dataPreProcessing(indentationSet.targetDir*resultFile)
    # Load raw displacements

    #ctrl.plotMode && display(plot([xy[1:100:end,1]],[xy[1:100:end,2]]))

    xy , basefit , zeroLine , rampStartIdx, xZero  = offsetAndDriftCompensation(xy)
    # Find initial contact
    #ctrl.plotMode && display(plot!(xy[1:100:end,1],xy[1:100:end,2]))

    xy[:,1] .-= xy[:,2]
    xy[:,2] .*= indentationSet.springConstant
    # Convert displacement-deflection matrix to indentation-force matrix

    holdStartIdx = findStartOfHold(xy,"first")

    # Split into loading and unloading.
    xy_unld1 = xy[holdStartIdx:end,:];

    #Determine the end of the hold time.
    unloadStartIdx = findStartOfHold(xy_unld1,"last")

    # Split into two new pieces
    xy_hold = xy_unld1[1:unloadStartIdx-1,:];
    xy_unld = xy_unld1[unloadStartIdx:end,:];

    xy_unld5 = xy_unld[1:min(Int64(round(2000*0.95)),size(xy_unld,1)),:];  # OBS 2000 is hard coded!

    # Fitting of the unloading curve.
    stiffness_fit = Array{Float64}(undef,1)
    dispVals = xy_unld5[1:min(hyperParameters.unloadingFitRange,size(xy_unld5,1)),1]
    forceVals = xy_unld5[1:min(hyperParameters.unloadingFitRange,size(xy_unld5,1)),2]
    Fmax = xy_unld5[1,2]               # Maximum force during unloading

    if cmp(hyperParameters.unloadingFitFunction,"Oliver-Pharr") == 0
        Dmax = xy_unld5[1,1]
        # Maximum indentation depth during unloading

        function unloadFitFun(fitCoefs)
            return fitCoefs[1].*(dispVals .- fitCoefs[2]).^fitCoefs[3] .- forceVals
        end
        function unloadFitMinFun(fitCoefs)
            sqrt( sum( (unloadFitFun(fitCoefs) ./ forceVals).^2 ) )
        end

        lx = [-Inf, -Inf , 0.0]; ux = [Inf, minimum(dispVals)-1e-2 , Inf];
        dfc = TwiceDifferentiableConstraints(lx, ux)
        resultFit = optimize(unloadFitMinFun, dfc, [1.0, 1.0, 1.0], IPNewton())
        uld_p = resultFit.minimizer
        println(uld_p)
        stiffness_fit = uld_p[1]*uld_p[3]*(Dmax - uld_p[2]).^(uld_p[3] - 1)

    elseif cmp(hyperParameters.unloadingFitFunction, "Feng") == 0

        unloadFitFun2(fitCoefs) = fitCoefs[1] .+ fitCoefs[2].*forceVals.^0.5 + fitCoefs[3].*forceVals.^fitCoefs[4] .- dispVals
        unloadFitMinFun2(fitCoefs) = sqrt(sum( (unloadFitFun2(fitCoefs) ./ dispVals).^2) )
        resultFit = optimize(unloadFitMinFun2, [1.0 1.0 1.0 1.0], BFGS())
        uld_p = resultFit.minimizer
        println(uld_p)
        stiffness_fit = inv(( 0.5*uld_p[2].*Fmax.^-0.5 + uld_p[4]*uld_p[3].*Fmax.^(uld_p[4] - 1.0) ))

    end
    stiffness = stiffness_fit
    dhtdt = 0;

    maxIndentation = median(xy_unld5[1,1]) - dhtdt*(length(xy[rampStartIdx:holdStartIdx,1])+length(xy_hold[:,1]))/hyperParameters.sampleRate;  #%OBS OBS OBS

    if cmp(indentationSet.indenterType,"pyramid") == 0
        x0 = maxIndentation - 0.72*Fmax/stiffness;
    elseif cmp(indentationSet.indenterType,"hemisphere") == 0
        x0 = maxIndentation - 0.75*Fmax/stiffness;
    end

    area_xy = readdlm(indentationSet.areaFile, ' ', Float64, '\n')
    # % Determine the area by loading the calibration data and fitting a polynom to the data.

    if (x0 > 100.0)
        area_fit_end = length(area_xy[:,1])
    elseif (x0 < 100.0)
        area_fit_end = findfirst( x -> x > 100,area_xy[:,1])
    end

    tempVec = area_xy[1:area_fit_end,1]
    p_area = [tempVec.^2 tempVec tempVec.^0.5 tempVec.^0.25 tempVec.^0.125] \ area_xy[1:area_fit_end,2]

    unloadArea = [x0^2 x0 x0^0.5 x0^0.25 x0^0.125] * p_area
    unloadArea = unloadArea[1]

    return 1/stiffness_fit , unloadArea
    # Assign outputs
end

max_stars_repo_name: abrandberg/afmIndentation

max_issues_repo_name: abrandberg/afmIndentation

max_forks_repo_name: abrandberg/afmIndentation


text:
mutable struct FuncReferType
    ref::Ptr{UInt8}
    doi::Ptr{UInt8}
    bibtex::Ptr{UInt8}
end

xc_func_reference_get_ref(frt::Ref{FuncReferType}) =
    ccall((:xc_func_reference_get_ref, libxc), Ptr{UInt8}, (Ref{FuncReferType},), frt)
xc_func_reference_get_doi(frt::Ref{FuncReferType}) =
    ccall((:xc_func_reference_get_doi, libxc), Ptr{UInt8}, (Ref{FuncReferType},), frt)
xc_func_reference_get_bibtex(frt::Ref{FuncReferType}) =
    ccall((:xc_func_reference_get_bibtex, libxc), Ptr{UInt8}, (Ref{FuncReferType},), frt)

mutable struct FuncParamsType
    value::Float64
    description::Ptr{UInt8}
end

const MAX_REFERENCES = 5

mutable struct XCFuncInfoType
    number::Cint
    kind::Cint
    name::Ptr{UInt8}
    family::Cint
    refs::NTuple{MAX_REFERENCES, Ref{FuncReferType}}

    flags::Cint

    dens_threshold::Float64

    n_ext_params::Cint
    ext_params::Ref{FuncParamsType}
    set_ext_params::Ptr{Cvoid}

    init::Ptr{Cvoid}
    end_::Ptr{Cvoid}
    lda::Ptr{Cvoid}
    gga::Ptr{Cvoid}
    mgga::Ptr{Cvoid}
end

xc_func_info_get_number(info::Ref{XCFuncInfoType}) =
    ccall( (:xc_func_info_get_number, libxc), Cint, (Ref{XCFuncInfoType},), info)
xc_func_info_get_kind(info::Ref{XCFuncInfoType}) =
    ccall( (:xc_func_info_get_kind, libxc), Cint, (Ref{XCFuncInfoType},), info)
xc_func_info_get_name(info::Ref{XCFuncInfoType}) = unsafe_string(
    ccall( (:xc_func_info_get_name, libxc), Ptr{UInt8}, (Ref{XCFuncInfoType},), info) )
xc_func_info_get_family(info::Ref{XCFuncInfoType}) =
    ccall( (:xc_func_info_get_family, libxc), Cint, (Ref{XCFuncInfoType},), info)
xc_func_info_get_flags(info::Ref{XCFuncInfoType}) =
    ccall( (:xc_func_info_get_flags, libxc), Cint, (Ref{XCFuncInfoType},), info)
xc_func_info_get_references(info::Ref{XCFuncInfoType}, number::Int64) =
    ccall( (:xc_func_info_get_references, libxc), Ref{FuncReferType}, (Ref{XCFuncInfoType},Cint), info, Cint(number))
xc_func_info_get_n_ext_params(info::Ref{XCFuncInfoType}) =
    ccall( (:xc_func_info_get_n_ext_params, libxc), Cint, (Ref{XCFuncInfoType},), info)
xc_func_info_get_ext_params_description(info::Ref{XCFuncInfoType}, number::Int64) = unsafe_string(
    ccall( (:xc_func_info_get_ext_params_description, libxc), Ref{FuncReferType}, (Ref{XCFuncInfoType},Cint), info, Cint(number)) )
xc_func_info_get_ext_params_default_value(info::Ref{XCFuncInfoType}, number::Int64) =
    ccall( (:xc_func_info_get_ext_params_default_value, libxc), Float64, (Ref{XCFuncInfoType},Cint), info, Cint(number))

mutable struct _XCFuncType
    # info::Ref{XCFuncInfoType}
    nspin::Cint
    n_func_aux::Cint
    # func_aux::Ptr{Ref{XCFuncType}}
    mix_coef::Ptr{Float64}

    cam_omega::Float64
    cam_alpha::Float64
    cam_beta::Float64

    nlc_b::Float64
    nlc_C::Float64

    n_rho::Cint
    n_sigma::Cint
    n_tau::Cint
    n_lapl::Cint

    n_zk::Cint

    n_vrho::Cint
    n_vsigma::Cint
    n_vtau::Cint
    n_vlapl::Cint

    n_v2rho2::Cint
    n_v2sigma2::Cint
    n_v2tau2::Cint
    n_v2lapl2::Cint
    n_v2rhosigma::Cint
    n_v2rhotau::Cint
    n_v2rholapl::Cint
    n_v2sigmatau::Cint
    n_v2sigmalapl::Cint
    n_v2lapltau::Cint

    n_v3rho3::Cint
    n_v3rho2sigma::Cint
    n_v3rhosigma2::Cint
    n_v3sigma3::Cint

    params::Ptr{Cvoid}
    dens_threshold::Float64
end

mutable struct XCFuncType
end

xc_func_alloc() = ccall( (:xc_func_alloc, libxc), Ptr{XCFuncType}, () )
xc_func_init(p::Ptr{XCFuncType}, functional::Integer, nspin::Integer) =
    ccall( (:xc_func_init, libxc), Cint, (Ptr{XCFuncType}, Cint, Cint), p, functional, nspin)
xc_func_end(p::Ptr{XCFuncType}) = ccall( (:xc_func_end, libxc), Cvoid, (Ptr{XCFuncType},), p)
xc_func_free(p::Ptr{XCFuncType}) = ccall( (:xc_func_free, libxc), Cvoid, (Ptr{XCFuncType},), p)
xc_func_get_info(p::Ptr{XCFuncType}) = ccall( (:xc_func_get_info, libxc), Ptr{XCFuncInfoType}, (Ref{XCFuncType},), p )
xc_func_set_dens_threshold(p::Ptr{XCFuncType}, dens_threshold::Float64) = ccall( (:xc_func_set_dens_threshold, libxc), Cvoid, (Ptr{XCFuncType}, Float64), p, dens_threshold)
xc_func_set_ext_params(p::Ptr{XCFuncType}, ext_params::Ptr{Float64}) = ccall( (:xc_func_set_ext_params, libxc), Cvoid, (Ptr{XCFuncType}, Ptr{Float64}), p, ext_params)

######################################################################
##  LDA-xc
######################################################################

xc_lda!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, zk::Array{Float64, 1}, vrho::Array{Float64, 1}, v2rho2::Array{Float64, 1}, v3rho3::Array{Float64, 1}) =
    ccall( (:xc_lda, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), p, np, rho, zk, vrho, v2rho2, v3rho3)
xc_lda_exc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, zk::Array{Float64, 1}) =
    ccall( (:xc_lda_exc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}), p, np, rho, zk)
xc_lda_exc_vxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, zk::Array{Float64, 1}, vrho::Array{Float64, 1}) =
    ccall( (:xc_lda_exc_vxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), p, np, rho, zk, vrho)
xc_lda_vxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, vrho::Array{Float64, 1}) =
    ccall( (:xc_lda_vxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}), p, np, rho, vrho)
xc_lda_fxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, v2rho2::Array{Float64, 1}) =
    ccall( (:xc_lda_fxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}), p, np, rho, v2rho2)
xc_lda_kxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, v3rho3::Array{Float64, 1}) =
    ccall( (:xc_lda_kxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64}), p, np, rho, v3rho3)

######################################################################
##  GGA-xc
######################################################################

xc_gga!(p::Ptr{XCFuncType}, np::Int,
        rho::Array{Float64, 1}, sigma::Array{Float64, 1},
        zk::Array{Float64, 1}, vrho::Array{Float64, 1}, vsigma::Array{Float64, 1},
        v2rho2::Array{Float64, 1}, v2rhosigma::Array{Float64, 1}, v2sigma2::Array{Float64, 1},
        v3rho3::Array{Float64, 1}, v3rho2sigma::Array{Float64, 1}, v3rhosigma2::Array{Float64, 1}, v3sigma3::Array{Float64, 1}) =
    ccall( (:xc_gga, libxc), Cvoid, (Ptr{XCFuncType}, Cint,
                                     Ptr{Float64}, Ptr{Float64},
                                     Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                                     Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                                     Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                     p, np,
                                     rho, sigma,
                                     zk, vrho, vsigma,
                                     v2rho2, v2rhosigma, v2sigma2,
                                     v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)

xc_gga_exc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, sigma::Array{Float64, 1},
            zk::Array{Float64, 1}) =
    ccall( (:xc_gga_exc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                         Ptr{Float64}),
                                         p, np, rho, sigma,
                                         zk)

xc_gga_exc_vxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, sigma::Array{Float64, 1},
                zk::Array{Float64, 1}, vrho::Array{Float64, 1}, vsigma::Array{Float64, 1}) =
    ccall( (:xc_gga_exc_vxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                             Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                             p, np, rho, sigma,
                                             zk, vrho, vsigma)

xc_gga_vxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, sigma::Array{Float64, 1},
            vrho::Array{Float64, 1}, vsigma::Array{Float64, 1}) =
    ccall( (:xc_gga_vxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                         Ptr{Float64}, Ptr{Float64}),
                                         p, np, rho, sigma,
                                         vrho, vsigma)

xc_gga_fxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, sigma::Array{Float64, 1},
            v2rho2::Array{Float64, 1}, v2rhosigma::Array{Float64, 1}, v2sigma2::Array{Float64, 1}) =
    ccall( (:xc_gga_fxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                         Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                         p, np, rho, sigma,
                                         v2rho2, v2rhosigma, v2sigma2)

xc_gga_kxc!(p::Ptr{XCFuncType}, np::Int, rho::Array{Float64, 1}, sigma::Array{Float64, 1},
            v3rho3::Array{Float64, 1}, v3rho2sigma::Array{Float64, 1}, v3rhosigma2::Array{Float64, 1}, v3sigma3::Array{Float64, 1}) =
    ccall( (:xc_gga_kxc, libxc), Cvoid, (Ptr{XCFuncType}, Cint, Ptr{Float64}, Ptr{Float64},
                                         Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
                                         p, np, rho, sigma,
                                         v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3)

# TODO
#
# void xc_gga_lb_modified  (const xc_func_type *p, int np, const double *rho, const double *sigma,
# 			   double r, double *vrho);
#
# double xc_gga_ak13_get_asymptotic (double homo);
#
# double xc_hyb_exx_coef(const xc_func_type *p);
# void  xc_hyb_cam_coef(const xc_func_type *p, double *omega, double *alpha, double *beta);
# void  xc_nlc_coef(const xc_func_type *p, double *nlc_b, double *nlc_C);

######################################################################
##  meta GGA-xc
######################################################################
# void xc_mgga        (const xc_func_type *p, int np,
# 		      const double *rho, const double *sigma, const double *lapl, const double *tau,
# 		      double *zk, double *vrho, double *vsigma, double *vlapl, double *vtau,
# 		      double *v2rho2, double *v2sigma2, double *v2lapl2, double *v2tau2,
# 		      double *v2rhosigma, double *v2rholapl, double *v2rhotau,
# 		      double *v2sigmalapl, double *v2sigmatau, double *v2lapltau);
# void xc_mgga_exc    (const xc_func_type *p, int np,
# 		      const double *rho, const double *sigma, const double *lapl, const double *tau,
# 		      double *zk);
# void xc_mgga_exc_vxc(const xc_func_type *p, int np,
# 		      const double *rho, const double *sigma, const double *lapl, const double *tau,
# 		      double *zk, double *vrho, double *vsigma, double *vlapl, double *vtau);
# void xc_mgga_vxc    (const xc_func_type *p, int np,
# 		      const double *rho, const double *sigma, const double *lapl, const double *tau,
# 		      double *vrho, double *vsigma, double *vlapl, double *vtau);
# void xc_mgga_fxc    (const xc_func_type *p, int np,
# 		      const double *rho, const double *sigma, const double *lapl, const double *tau,
# 		      double *v2rho2, double *v2sigma2, double *v2lapl2, double *v2tau2,
# 		      double *v2rhosigma, double *v2rholapl, double *v2rhotau,
# 		      double *v2sigmalapl, double *v2sigmatau, double *v2lapltau);

######################################################################
# /* Functionals that are defined as mixtures of others */
######################################################################

# void xc_mix_func
#   (const xc_func_type *func, int np,
#    const double *rho, const double *sigma, const double *lapl, const double *tau,
#    double *zk, double *vrho, double *vsigma, double *vlapl, double *vtau,
#    double *v2rho2, double *v2sigma2, double *v2lapl2, double *v2tau2,
#    double *v2rhosigma, double *v2rholapl, double *v2rhotau,
#    double *v2sigmalapl, double *v2sigmatau, double *v2lapltau);


text:
# DESCRIPTION
# Model order reduction based on residualization. Keep rigid-body states and
# long and lat flapping angles.
#
# ------------------------------------------------------------------------------

function ModRed10(A,B)
    idxss=[collect(1:8); collect(15:16)]
    idxfs=[collect(13:14); collect(17:31); collect(33:23)]
    Ared=A[idxss,idxss]-A[idxss,idxfs]/A[idxfs,idxfs]*A[idxfs,idxss]
    Bred=B[idxss,1:4]-A[idxss,idxfs]/A[idxfs,idxfs]*B[idxfs,1:4]
    return Ared, Bred
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text: latexify_default_dt(::Type{Symbolics.Num}) = (Symbolics.@variables(h); h)


text:
print("\ttst_basic.jl")

# generate random instance
srand(0)
m = 2^4
n = 2^4
A = hcat(sprand(m, n, 0.1), speye(m))
u = 10.0 * sprand(n, 1.0)
b = 1.0 + rand(m)
c = vcat(1.0 - 2.0 * rand(n), zeros(m))

model = Tulip.Model(A, b, c, u.nzind, u.nzval)
status = Tulip.solve!(model, verbose=0, tol=10.0^-8)

@test status == :Optimal

objval = dot(model.sol.x, model.c)
const OBJVAL = -15.989485968940144
@test abs(objval - OBJVAL) <= 10.0^-8

println("\tPassed.")