text:
##
## function get_transcripts_from_gff(record::FASTA.Record,
##    transDict::Dict)::Array{FASTA.Record,1}
##

dcAll = Dict("Seq" => [String("1-5,8-12,+"), String("1-5,+"), String("1-20,+")])
append!(dcAll["Seq"],[String("1-5,8-12,-"), String("1-5,-"), String("1-20,-")])
append!(dcAll["Seq"],[String("1-20,5-10,+"), String("1-20,5-10,-")])
outRecordsAll=Array{FASTA.Record,1}()

# + strand testing
record = FASTA.Record("Seq", dna"GTGGTGGCGCTCACCGACTG")
dc = Dict("Seq" => [String("1-5,8-12,+")])
recordOut = FASTA.Record("Seq", dna"GTGGTCGCTC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-5,+")])
recordOut = FASTA.Record("Seq", dna"GTGGT")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-20,+")])
recordOut = FASTA.Record("Seq", dna"GTGGTGGCGCTCACCGACTG")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

# - strand testing
dc = Dict("Seq" => [String("1-5,8-12,-")])
recordOut = FASTA.Record("Seq", dna"GAGCGACCAC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-5,-")])
recordOut = FASTA.Record("Seq", dna"ACCAC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-20,-")])
recordOut = FASTA.Record("Seq", dna"CAGTCGGTGAGCGCCACCAC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-20,5-10,+")])
recordOut = FASTA.Record("Seq", dna"GTGGTGGCGCTCACCGACTGTGGCGC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

dc = Dict("Seq" => [String("1-20,5-10,-")])
recordOut = FASTA.Record("Seq", dna"GCGCCACAGTCGGTGAGCGCCACCAC")
outRecords = Array{FASTA.Record,1}()
push!(outRecords,recordOut)
push!(outRecordsAll,recordOut)
@test outRecords == get_transcripts_from_dict(record,dc)

# Test all records in list
@test outRecordsAll == get_transcripts_from_dict(record,dcAll)

##
##  function get_transcripts_from_gff(fa::String, gff::String)::Array{FASTA.Record,1}
##

# Creating Temp Files (gff and fasta) for testing.
f1,file_stream = mktemp(tempdir())

writer = GFF3.Writer(file_stream)
write(writer,GFF3.Record(b"chr1	HAVANA	gene	1	100	.	+	.	ID=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	1	100	.	+	.	ID=ENST00000456328.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	1	20	.	+	.	ID=exon:ENST00000456328.2:1;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	30	60	.	+	.	ID=exon:ENST00000456328.2:2;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	65	80	.	+	.	ID=exon:ENST00000456328.2:3;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	2	98	.	+	.	ID=ENST00000450305.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	2	15	.	+	.	ID=exon:ENST00000450305.2:1;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	20	35	.	+	.	ID=exon:ENST00000450305.2:2;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	37	45	.	+	.	ID=exon:ENST00000450305.2:3;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	50	65	.	+	.	ID=exon:ENST00000450305.2:4;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	70	80	.	+	.	ID=exon:ENST00000450305.2:5;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	85	95	.	+	.	ID=exon:ENST00000450305.2:6;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	gene	1	100	.	+	.	ID=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	1	100	.	-	.	ID=ENST00000456328.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	1	20	.	-	.	ID=exon:ENST00000456328.2:1;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	30	60	.	-	.	ID=exon:ENST00000456328.2:2;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	65	80	.	-	.	ID=exon:ENST00000456328.2:3;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450305.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	2	15	.	-	.	ID=exon:ENST00000450305.2:1;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	20	35	.	-	.	ID=exon:ENST00000450305.2:2;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	37	45	.	-	.	ID=exon:ENST00000450305.2:3;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	50	65	.	-	.	ID=exon:ENST00000450305.2:4;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	70	80	.	-	.	ID=exon:ENST00000450305.2:5;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	85	95	.	-	.	ID=exon:ENST00000450305.2:6;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450304.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	2	15	.	-	.	ID=exon:ENST00000450304.2:1;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	20	35	.	-	.	ID=exon:ENST00000450304.2:2;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	37	45	.	-	.	ID=exon:ENST00000450304.2:3;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	50	65	.	-	.	ID=exon:ENST00000450304.2:4;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	70	80	.	-	.	ID=exon:ENST00000450304.2:5;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	85	95	.	-	.	ID=exon:ENST00000450304.2:6;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	2	98	.	+	.	ID=ENST00000450303.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	2	98	.	+	.	ID=exon:ENST00000450302.2:1;Parent=ENST00000450303.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450300.2;Parent=ENST00000450300.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr1	HAVANA	exon	2	98	.	-	.	ID=exon:ENST00000450300.2:1;Parent=ENST00000450300.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	gene	1	100	.	+	.	ID=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	1	100	.	+	.	ID=ENST00000456328.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	1	20	.	+	.	ID=exon:ENST00000456328.2:1;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	30	60	.	+	.	ID=exon:ENST00000456328.2:2;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	65	80	.	+	.	ID=exon:ENST00000456328.2:3;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	2	98	.	+	.	ID=ENST00000450305.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	2	15	.	+	.	ID=exon:ENST00000450305.2:1;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	20	35	.	+	.	ID=exon:ENST00000450305.2:2;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	37	45	.	+	.	ID=exon:ENST00000450305.2:3;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	50	65	.	+	.	ID=exon:ENST00000450305.2:4;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	70	80	.	+	.	ID=exon:ENST00000450305.2:5;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	85	95	.	+	.	ID=exon:ENST00000450305.2:6;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	gene	1	100	.	+	.	ID=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	1	100	.	-	.	ID=ENST00000456328.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	1	20	.	-	.	ID=exon:ENST00000456328.2:1;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	30	60	.	-	.	ID=exon:ENST00000456328.2:2;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	65	80	.	-	.	ID=exon:ENST00000456328.2:3;Parent=ENST00000456328.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450305.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	2	15	.	-	.	ID=exon:ENST00000450305.2:1;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	20	35	.	-	.	ID=exon:ENST00000450305.2:2;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	37	45	.	-	.	ID=exon:ENST00000450305.2:3;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	50	65	.	-	.	ID=exon:ENST00000450305.2:4;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	70	80	.	-	.	ID=exon:ENST00000450305.2:5;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	85	95	.	-	.	ID=exon:ENST00000450305.2:6;Parent=ENST00000450305.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450304.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	2	15	.	-	.	ID=exon:ENST00000450304.2:1;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	20	35	.	-	.	ID=exon:ENST00000450304.2:2;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	37	45	.	-	.	ID=exon:ENST00000450304.2:3;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	50	65	.	-	.	ID=exon:ENST00000450304.2:4;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	70	80	.	-	.	ID=exon:ENST00000450304.2:5;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	85	95	.	-	.	ID=exon:ENST00000450304.2:6;Parent=ENST00000450304.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	2	98	.	+	.	ID=ENST00000450303.2;Parent=ENSG00000223972.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	2	98	.	+	.	ID=exon:ENST00000450302.2:1;Parent=ENST00000450303.2;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	transcript	2	98	.	-	.	ID=ENST00000450300.2;Parent=ENST00000450300.5;gene_id=ENSG00000223972.5"))
write(writer,GFF3.Record(b"chr2	HAVANA	exon	2	98	.	-	.	ID=exon:ENST00000450300.2:1;Parent=ENST00000450300.2;gene_id=ENSG00000223972.5"))
close(writer)

f2,file_stream = mktemp(tempdir())
writer = FASTA.Writer(file_stream)
write(writer,FASTA.Record("chr1", dna"TACTGCCTGCATAAGGAGAACGGAGTTGCCAAGGACGAAAGCGACTCTAGGTTCTAACCGTCGACTTTGGCGGAAAGGTTTCACTCAGGAAGCAGACACT"))
write(writer,FASTA.Record("chr2", dna"TACTGCCTGCATAAGGAGAACGGAGTTGCCAAGGACGAAAGCGACTCTAGGTTCTAACCGTCGACTTTGGCGGAAAGGTTTCACTCAGGAAGCAGACACT"))
close(writer)

outRecords = Array{FASTA.Record,1}()
push!(outRecords,FASTA.Record("Seq", dna"TACTGCCTGCATAAGGAGAACAAGGACGAAAGCGACTCTAGGTTCTAACCGCTTTGGCGGAAAGGTT"))
push!(outRecords,FASTA.Record("Seq", dna"ACTGCCTGCATAAGACGGAGTTGCCAAGGAGAAAGCGACGGTTCTAACCGTCGACGCGGAAAGGTTTCAGGAAGCAG"))
push!(outRecords,FASTA.Record("Seq", dna"AACCTTTCCGCCAAAGCGGTTAGAACCTAGAGTCGCTTTCGTCCTTGTTCTCCTTATGCAGGCAGTA"))
push!(outRecords,FASTA.Record("Seq", dna"CTGCTTCCTGAAACCTTTCCGCGTCGACGGTTAGAACCGTCGCTTTCTCCTTGGCAACTCCGTCTTATGCAGGCAGT"))
push!(outRecords,FASTA.Record("Seq", dna"ACTGCCTGCATAAGGAGAACGGAGTTGCCAAGGACGAAAGCGACTCTAGGTTCTAACCGTCGACTTTGGCGGAAAGGTTTCACTCAGGAAGCAGACA"))
push!(outRecords,FASTA.Record("Seq", dna"TGTCTGCTTCCTGAGTGAAACCTTTCCGCCAAAGTCGACGGTTAGAACCTAGAGTCGCTTTCGTCCTTGGCAACTCCGTTCTCCTTATGCAGGCAGT"))
push!(outRecords,FASTA.Record("Seq", dna"TACTGCCTGCATAAGGAGAACAAGGACGAAAGCGACTCTAGGTTCTAACCGCTTTGGCGGAAAGGTT"))
push!(outRecords,FASTA.Record("Seq", dna"ACTGCCTGCATAAGACGGAGTTGCCAAGGAGAAAGCGACGGTTCTAACCGTCGACGCGGAAAGGTTTCAGGAAGCAG"))
push!(outRecords,FASTA.Record("Seq", dna"AACCTTTCCGCCAAAGCGGTTAGAACCTAGAGTCGCTTTCGTCCTTGTTCTCCTTATGCAGGCAGTA"))
push!(outRecords,FASTA.Record("Seq", dna"CTGCTTCCTGAAACCTTTCCGCGTCGACGGTTAGAACCGTCGCTTTCTCCTTGGCAACTCCGTCTTATGCAGGCAGT"))
push!(outRecords,FASTA.Record("Seq", dna"ACTGCCTGCATAAGGAGAACGGAGTTGCCAAGGACGAAAGCGACTCTAGGTTCTAACCGTCGACTTTGGCGGAAAGGTTTCACTCAGGAAGCAGACA"))
push!(outRecords,FASTA.Record("Seq", dna"TGTCTGCTTCCTGAGTGAAACCTTTCCGCCAAAGTCGACGGTTAGAACCTAGAGTCGCTTTCGTCCTTGGCAACTCCGTTCTCCTTATGCAGGCAGT"))

# Testing
@test outRecords == get_transcripts_from_gff(f2, f1)

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
#=
This is Julia implementation of Go Board by Tejan Karmali.
Reference: Python version of the implementation of Go by
(https://github.com/tensorflow/minigo)
=#

import Base: ==, deepcopy, show
using IterTools: chain

# Cell occupants
#@enum CELL_OCC WHITE = -1 EMPTY BLACK  FILL KO UNKNOWN
WHITE, EMPTY, BLACK, FILL, KO, UNKNOWN = collect(-1:4)

# Represents "group not found" in the LibertyTracker object
MISSING_GROUP_ID = -1

struct PlayerMove
  color
  move
end

function place_stones!(board, color, stones)
  for s in stones
    board[s...] = color
  end
end

function find_reached(board, c, env::GoEnv)
  color = board[c...]
  chain = Set([c])
  reached = Set()
  frontier = [c]
  while length(frontier) != 0
    current = pop!(frontier)
    push!(chain, current)
    for n in env.NEIGHBORS[current]
      if board[n...] == color && n ∉ chain
        push!(frontier, n)
      elseif board[n...] != color
        push!(reached, n)
      end
    end
  end
  return chain, reached
end

function is_koish(board, c, env::GoEnv)
  # Check if c is surrounded on all sides by 1 color, and return that color
  if board[c...] != EMPTY return nothing end
  neighs = Set(board[n...] for n in env.NEIGHBORS[c...])
  if length(neighs) == 1 && EMPTY ∉ neighs
    return collect(neighs)[1]
  else
    return nothing
  end
end

function is_eyeish(board, c, env)
  # Check if c is an eye, for the purpose of restricting MC rollouts.
  color = is_koish(board, c, env)
  if color == nothing
    return nothing
  end

  diag_faults = 0
  diag = env.DIAGONALS[c...]
  if length(diag) < 4
    diag_faults += 1
  end

  for d in diag
    if board[d...] ∉ (color, EMPTY)
      diag_faults += 1
    end
  end
  if diag_faults > 1
    return nothing
  else
    return color
  end
end

# making a (x,y) coordinate for board from flattened index x
coordify(x; N = 19) = 1 + (x - 1) % N, 1 + (x - 1) ÷ N

get_first_n(x::Array{Int8, 3}, n) = size(x, 3) < n ? x : x[:, :, 1:n]

mutable struct Group
  id::Int
  stones::Set{NTuple{2, Int}}
  liberties::Set{NTuple{2, Int}}
  color::Int
end

function ==(g1::Group, g2::Group)
  return g1.stones == g2.stones && g1.liberties == g2.liberties && g1.color == g2.color
end

mutable struct LibertyTracker
  group_index::Array{Int16, 2}
  groups::Dict{Int, Group}
  liberty_cache::Array{UInt8, 2}
  max_group_id::Int16

  function LibertyTracker(env::GoEnv; group_index = nothing, groups = nothing,
    liberty_cache = nothing, max_group_id = 1)
    # group_index: a NxN array of group_ids. -1 means no group
    # groups: a dict of group_id to groups
    # liberty_cache: a NxN array of liberty counts
    N = env.N
    grp_idx = group_index != nothing ? group_index : -ones(Int16, N, N)
    grps = groups != nothing ? groups : Dict{Int, Group}()
    lib_cache = liberty_cache != nothing ? liberty_cache : zeros(UInt8, N, N)
    new(grp_idx, grps, lib_cache, max_group_id)
  end
end

function deepcopy(lib_trac::LibertyTracker)
  new_group_index = deepcopy(lib_trac.group_index)
  new_lib_cache = deepcopy(lib_trac.liberty_cache)
  new_groups = Dict{Int, Group}(
      group.id => Group(group.id, Set(group.stones), Set(group.liberties), group.color)
      for group in values(lib_trac.groups)
  )
  N = size(new_lib_cache, 1)
  return LibertyTracker(GoEnv(N), group_index = new_group_index,
                  groups = new_groups,
                  liberty_cache = new_lib_cache,
                  max_group_id = lib_trac.max_group_id)
end

function from_board(board, env::GoEnv)
  board = deepcopy(board)
  curr_group_id = 0
  lib_tracker = LibertyTracker(env)
  for color ∈ (WHITE, BLACK)
    while color ∈ board
      curr_group_id += 1
      coord_flat = findfirst(x -> x == color, board[:])
      coord = coordify(coord_flat, N = env.N)
      chain, reached = find_reached(board, coord, env)
      liberties = Set(r for r in reached if board[r...] == EMPTY)
      new_group = Group(curr_group_id, chain, liberties, color)
      lib_tracker.groups[curr_group_id] = new_group
      for s in chain
        lib_tracker.group_index[s...] = curr_group_id
      end
      place_stones!(board, FILL, chain)
    end
  end

  lib_tracker.max_group_id = curr_group_id

  liberty_counts = zeros(UInt8, env.N, env.N)
  for group in values(lib_tracker.groups)
    num_libs = length(group.liberties)
    for s in group.stones
      liberty_counts[s...] = num_libs
    end
  end
  lib_tracker.liberty_cache = liberty_counts

  return lib_tracker
end

function _merge_from_played!(lib_trac::LibertyTracker, color, played, libs, other_group_ids)
  stones = Set([played])
  liberties = Set{NTuple{2, Int}}(libs)
  for group_id in other_group_ids
    other = pop!(lib_trac.groups, group_id)
    union!(stones, other.stones)
    union!(liberties, other.liberties)
  end

  if !isempty(other_group_ids)
    setdiff!(liberties, Set([played]))
  end

  @assert setdiff(stones, liberties) == stones

  lib_trac.max_group_id += 1
  result = Group(lib_trac.max_group_id, stones, liberties, color)
  lib_trac.groups[result.id] = result

  for s in result.stones
    lib_trac.group_index[s...] = result.id
    lib_trac.liberty_cache[s...] = length(result.liberties)
  end
  return result
end

function _update_liberties!(lib_trac::LibertyTracker, group_id::Integer;
  add = Set{NTuple{2, Int}}(), remove = Set{NTuple{2, Int}}())

  group = lib_trac.groups[group_id]
  new_libs = setdiff(group.liberties ∪ add, remove)
  lib_trac.groups[group_id] = Group(group_id, group.stones, new_libs, group.color)

  new_lib_count = length(new_libs)
  for s in lib_trac.groups[group_id].stones
    lib_trac.liberty_cache[s...] = new_lib_count
  end
end

function _capture_group!(lib_trac::LibertyTracker, group_id::Integer)
  dead_group = lib_trac.groups[group_id]
  delete!(lib_trac.groups, group_id)
  for s in dead_group.stones
    lib_trac.group_index[s...] = MISSING_GROUP_ID
    lib_trac.liberty_cache[s...] = 0
  end
  return dead_group.stones
end

function _handle_captures!(lib_trac::LibertyTracker, captured_stones, env::GoEnv)
  for s in captured_stones
    for n in env.NEIGHBORS[s...]
      group_id = lib_trac.group_index[n...]
      if group_id != MISSING_GROUP_ID
        _update_liberties!(lib_trac, group_id; add = Set([s]))
      end
    end
  end
end

function add_stone!(lib_trac::LibertyTracker, color, c, env::GoEnv)
  @assert lib_trac.group_index[c...] == MISSING_GROUP_ID
  captured_stones = Set()
  opponent_neighboring_group_ids = Set{Int}()
  friendly_neighboring_group_ids = Set{Int}()
  empty_neighbors = Set{NTuple{2, Int}}()

  for n in env.NEIGHBORS[c...]
    neighbor_group_id = lib_trac.group_index[n...]
    if neighbor_group_id != MISSING_GROUP_ID
      neighbor_group = lib_trac.groups[neighbor_group_id...]
      if neighbor_group.color == color
        push!(friendly_neighboring_group_ids, neighbor_group_id)
      else
        push!(opponent_neighboring_group_ids, neighbor_group_id)
      end
    else
      push!(empty_neighbors, n)
    end
  end

  new_group = _merge_from_played!(lib_trac, color, c, empty_neighbors, friendly_neighboring_group_ids)

  # new_group becomes stale as _update_liberties and
  # _handle_captures are called; must refetch with lib_trac.groups[new_group.id]
  for group_id in opponent_neighboring_group_ids
    neighbor_group = lib_trac.groups[group_id]
    if length(neighbor_group.liberties) == 1
      captured = _capture_group!(lib_trac, group_id)
      union!(captured_stones, captured)
    else
      _update_liberties!(lib_trac, group_id; remove = Set([c]))
    end
  end
  _handle_captures!(lib_trac, captured_stones, env)

  # suicide is illegal
  if length(lib_trac.groups[new_group.id].liberties) == 0
    throw(IllegalMove())
  end

  return captured_stones
end

mutable struct GoPosition <: Position
  #=
  board: an array
  n: an int representing moves played so far
  komi: a float, representing points given to the second player.
  caps: a (int, int) tuple of captures for B, W.
  lib_tracker: a LibertyTracker object
  ko: a Move
  recent: a tuple of PlayerMoves, such that recent[end] is the last move.
  board_deltas: an array of shape (N, N, n) representing changes
      made to the board at each move (played move and captures).
      Should satisfy next_pos.board - next_pos.board_deltas[:, :, 0] == pos.board
  to_play: BLACK or WHITE
  =#
  env::GoEnv
  board::Array{Int8, 2}
  n::Int
  komi::Float32
  caps::NTuple{2, Int}
  lib_tracker::LibertyTracker
  ko
  recent::Vector{PlayerMove}
  board_deltas::Array{Int8, 3}
  to_play::Int
  done::Bool

  function GoPosition(env::GoEnv; board = nothing, n = 0, komi = 7.5,
    caps = (0, 0), lib_tracker = nothing, ko = nothing,
    recent = Vector{PlayerMove}(), board_deltas = nothing, to_play = BLACK)

    b = board != nothing ? board : deepcopy(env.EMPTY_BOARD)
    lib_trac = lib_tracker != nothing ? lib_tracker : from_board(b, env)
    bd = board_deltas != nothing ? board_deltas : zeros(Int8, env.N, env.N, 0)
    new(env, b, n, komi, caps, lib_trac, ko, recent, bd, to_play, false)
  end
end

function deepcopy(pos::GoPosition)
  new_board = deepcopy(pos.board)
  new_lib_tracker = deepcopy(pos.lib_tracker)
  new_recent = deepcopy(pos.recent)
  return GoPosition(pos.env; board = new_board, n = pos.n, komi = pos.komi, caps = pos.caps,
                  lib_tracker = new_lib_tracker, ko = pos.ko, recent = new_recent,
                  board_deltas = pos.board_deltas, to_play = pos.to_play)
end

function show(io::IO, pos::GoPosition)
  pretty_print_map = Dict{Int, String}([
            WHITE => "O",
            EMPTY => ".",
            BLACK => "X",
            FILL => "#",
            KO => "*"
        ])
  board = deepcopy(pos.board)
  N = size(board, 1)

  captures = pos.caps
  if pos.ko != nothing
    place_stones!(board, KO, [pos.ko])
  end

  raw_board_contents = []
  for i = 1:N
    row = []
    for j = 1:N
      appended = length(pos.recent) != 0 && (i, j) == pos.recent[end].move ? "<" : " "
      push!(row, pretty_print_map[board[i,j]] * appended)
    end
    push!(raw_board_contents, join(row))
  end

  row_labels_fmt = [@sprintf("%2d ", i) for i = N:-1:1]
  row_labels = [string(i) for i = N:-1:1]
  annotated_board_contents = [join(r) for r in zip(row_labels_fmt, raw_board_contents, row_labels)]
  header_footer_rows = ["   " * join("ABCDEFGHJKLMNOPQRST"[1:N], " ") * "   "]
  annotated_board = join(collect(chain(header_footer_rows, annotated_board_contents,
                                        header_footer_rows)), "\n")
  details = "\nMove: $(pos.n). Captures X: $(captures[1]) O: $(captures[2])\n"
  turn = pos.done ? "GAME OVER\n" : "To Play: " * (pos.to_play == BLACK ? "X(BLACK)\n" : "O(WHITE)\n")
  println(annotated_board * details * turn)
end

function is_move_suicidal(pos::GoPosition, move)
  potential_libs = Set()
  for n in pos.env.NEIGHBORS[move...]
    neighbor_group_id = pos.lib_tracker.group_index[n...]
    if neighbor_group_id == MISSING_GROUP_ID
    # at least one liberty after playing here, so not a suicide
      return false
    end
    neighbor_group = pos.lib_tracker.groups[neighbor_group_id...]
    if neighbor_group.color == pos.to_play
      union!(potential_libs, neighbor_group.liberties)
    elseif length(neighbor_group.liberties) == 1
      # would capture an opponent group if they only had one lib.
      return false
    end
  end
  # it's possible to suicide by connecting several friendly groups
  # each of which had one liberty.
  setdiff!(potential_libs, Set([move]))
  return isempty(potential_libs)
end

function is_move_legal(pos::GoPosition, move)
  # Checks that a move is on an empty space, not on ko, and not suicide
  if move == nothing
    return true
  end
  if pos.board[move...] != EMPTY
    return false
  end
  if move == pos.ko
    return false
  end
  if is_move_suicidal(pos, move)
    return false
  end
  return true
end

function all_legal_moves(pos::GoPosition)
  # Returns an array of size N² + 1, with 1 = legal, 0 = illegal
  # by default, every move is legal
  N = size(pos.board, 1)
  legal_moves = ones(Int8, N, N)
  # ...unless there is already a stone there
  legal_moves[pos.board .!= EMPTY] .= 0
  # calculate which spots have 4 stones next to them
  # padding is because the edge always counts as a lost liberty.
  adjacent = ones(Int8, N + 2, N + 2)
  adjacent[2:end - 1, 2:end - 1] = abs.(pos.board)
  num_adjacent_stones = adjacent[1:end - 2, 2:end - 1] + adjacent[2:end - 1, 1:end - 2] +
                        adjacent[3:end, 2:end - 1] + adjacent[2:end - 1, 3:end]
  # Surrounded spots are those that are empty and have 4 adjacent stones.
  surrounded_spots = (pos.board .== EMPTY) .* (num_adjacent_stones .== 4)
  # Such spots are possibly illegal, unless they are capturing something.
  # Iterate over and manually check each spot.
  for c in findall(x -> x != 0, surrounded_spots[:])
    coord = coordify(c; N = N)
    if is_move_suicidal(pos, coord)
      legal_moves[coord...] = 0
    end
  end

  # ...and retaking ko is always illegal
  if pos.ko != nothing
    legal_moves[pos.ko...] = 0
  end

  # and pass is always legal
  return vcat(legal_moves[:], [1])
end

function pass_move!(pos::GoPosition; mutate = false)
  new_pos = mutate ? pos : deepcopy(pos)
  new_pos.n += 1

  N = size(pos.board, 1)

  push!(new_pos.recent, PlayerMove(new_pos.to_play, nothing))
  new_pos.board_deltas = cat(dims=3, zeros(Int8, N, N, 1), get_first_n(new_pos.board_deltas, new_pos.env.planes-2))
  new_pos.to_play *= -1
  new_pos.ko = nothing
  if length(new_pos.recent) > 1 && new_pos.recent[end - 1].move == nothing
    new_pos.done = true
  end
  return new_pos
end

function flip_playerturn!(pos::GoPosition; mutate = false)
  new_pos = mutate ? pos : deepcopy(pos)
  new_pos.ko = nothing
  new_pos.to_play *= -1
  return new_pos
end

get_liberties(pos::GoPosition) = pos.lib_tracker.liberty_cache

function play_move!(pos::GoPosition, c; color = nothing, mutate = false)
  # Obeys CGOS Rules of Play. defIn short:
  # No suicides
  # Chinese/area scoring
  # Positional superko (this is very crudely approximate at the moment.)
  if color == nothing
    color = pos.to_play
  end

  new_pos = mutate ? pos : deepcopy(pos)

  @assert !new_pos.done

  if c == nothing
    new_pos = pass_move!(new_pos; mutate = mutate)
    return new_pos
  end

  if !is_move_legal(pos, c)
    throw(IllegalMove())
  end

  potential_ko = is_koish(new_pos.board, c, pos.env)

  place_stones!(new_pos.board, color, [c])
  captured_stones = add_stone!(new_pos.lib_tracker, color, c, pos.env)
  place_stones!(new_pos.board, EMPTY, captured_stones)

  opp_color = -color

  N = size(pos.board, 1)

  new_board_delta = zeros(Int8, N, N)
  new_board_delta[c...] = color
  place_stones!(new_board_delta, color, captured_stones)

  if length(captured_stones) == 1 && potential_ko == opp_color
    new_ko = collect(captured_stones)[1]
  else
    new_ko = nothing
  end

  if new_pos.to_play == BLACK
    new_caps = (new_pos.caps[1] + length(captured_stones), new_pos.caps[2])
  else
    new_caps = (new_pos.caps[1], new_pos.caps[2] + length(captured_stones))
  end
  new_pos.n += 1
  new_pos.caps = new_caps
  new_pos.ko = new_ko
  push!(new_pos.recent, PlayerMove(color, c))
  # keep a rolling history of last 7 deltas - that's all we'll need to
  # extract the last 8 board states.

  new_pos.board_deltas = cat(dims = 3, reshape(new_board_delta, N, N, 1),
			get_first_n(new_pos.board_deltas, new_pos.env.planes-2))
  new_pos.to_play *= -1
  return new_pos
end

function score(pos::GoPosition)
  # Return score from B perspective. If W is winning, score is negative.
  working_board = deepcopy(pos.board)
  N = size(pos.board, 1)
  while EMPTY ∈ working_board
    unassigned_spaces = findfirst(x -> x == EMPTY, working_board[:])
    c = coordify(unassigned_spaces, N =  N)
    territory, borders = find_reached(working_board, c, pos.env)
    border_colors = Set(working_board[b...] for b in borders)
    X_border = BLACK ∈ border_colors
    O_border = WHITE ∈ border_colors
    if X_border && ! O_border
      territory_color = BLACK
    elseif O_border && ! X_border
      territory_color = WHITE
    else
      territory_color = UNKNOWN # dame, or seki
    end
    place_stones!(working_board, territory_color, territory)
  end
  return count(x -> x == BLACK, working_board) -
  count(x -> x == WHITE, working_board) - pos.komi
end

function result(pos::GoPosition)
  points = score(pos)
  if points > 0
    return 1
  elseif points < 0
    return -1
  else
    return 0
  end
end

function result_string(pos::GoPosition)
  points = score(pos)
  if points > 0
    return "B+" * @sprintf("%.1f", points)
  elseif points < 0
    return "W+" *  @sprintf("%.1f", abs(points))
  else
    return "DRAW"
  end
end

function replay_position(pos::GoPosition, result)
  #=
  Wrapper for a GoPosition which replays its history.
  Assumes an empty start position! (i.e. no handicap, and history must be exhaustive.)

  Result must be passed in, since a resign cannot be inferred from position
  history alone.

  for position_w_context in replay_position(position):
    print(position_w_context.position)
  =#
  pos.n == length(pos.recent) ? nothing : throw(AssertionError("GoPosition history is incomplete"))
  replay_buffer = Vector{PositionWithContext}()
  dummy_pos = GoPosition(pos.env; komi = pos.komi)

  for player_move in pos.recent
    color, next_move = player_move.color, player_move.move
    push!(replay_buffer, PositionWithContext(dummy_pos, next_move, result))
    dummy_pos = play_move!(dummy_pos, next_move, color = color)
  end
  return replay_buffer
end

ext: jl


text:
# Xpress model

#################################################
#
#  model type & constructors
#
#################################################

type Model
    ptr_model::Ptr{Void}
    callback::Any
    finalize_env::Bool

    function Model(p::Ptr{Void}; finalize_env::Bool=false)
        model = new(p, nothing, finalize_env)
        finalizer(model, m -> (free_model(m)) )
        model
    end
end

# If environment is tied to this model, the, we can safely finalize
# both at the same time, working around the Julia GC.
#function Model(env::Env; finalize_env::Bool=false)
function Model(; finalize_env::Bool=false)

    #@assert is_valid(env)

    a = Array(Ptr{Void}, 1)
    ret = @xprs_ccall(createprob, Cint, ( Ptr{Ptr{Void}},), a )
    #println(a)
    if ret != 0
        error("It was not possible to create a model, try running Env() and then create the model again.")
        #throw(XpressError(env, ret))
    end

    m = Model(a[1]; finalize_env = finalize_env)

    load_empty(m)

    return m
end

function Model(name::Compat.ASCIIString, sense::Symbol)
    model = Model()
    if sense != :minimize
        set_sense!(model, sense)
    end
    model
end
function Model(name::Compat.ASCIIString)
    model = Model()
    set_sense!(model, :minimize)
    model
end

#################################################
#
#  model error handling
#
#################################################
function get_error_msg(m::Model)
    #@assert env.ptr_env == 1
    out = Array(Cchar, 512)
    ret = @xprs_ccall(getlasterror, Cint, (Ptr{Void},Ptr{Cchar}),
        m.ptr_model, out)

    error( unsafe_string(pointer(out))  )
end
function get_error_msg(m::Model,ret::Int)
    #@assert env.ptr_env == 1
    out = Array(Cint,1)
    out2 = @xprs_ccall(getintattrib, Cint,(Ptr{Void}, Cint, Ptr{Cint}),
        m.ptr_model, ret , out)
end
##############################################
# Error handling
##############################################
type XpressError
    code::Int
    msg::Compat.ASCIIString
end
#function XpressError(ret::Int,m::Model)#, code::Integer)
#    XpressError( ret, get_error_msg(m) )#convert(Int, code), get_error_msg(env))
#end
function XpressError(m::Model)#, code::Integer)
    XpressError( 0, get_error_msg(m) )#convert(Int, code), get_error_msg(env))
end
function XpressError(m::Model,ret::Int)
    XpressError( get_error_msg(m,ret), "get message from optimizer manual" )
end

#################################################
#
#  model manipulation
#
#################################################

Base.unsafe_convert(ty::Type{Ptr{Void}}, model::Model) = model.ptr_model::Ptr{Void}

function free_model(model::Model)

    if model.ptr_model != C_NULL
        ret = @xprs_ccall(destroyprob, Cint, (Ptr{Void},), model.ptr_model)
        if ret != 0
            throw(XpressError(model))
        end
        model.ptr_model = C_NULL
    end
end

function copy(model_src::Model)
    # only copies problem not callbacks and controls
    if model_src.ptr_model != C_NULL
        model_dest = Model()#env)
        ret = @xprs_ccall(copyprob, Cint, (Ptr{Void},Ptr{Void},Ptr{UInt8}),
            model_dest.ptr_model, model_src.ptr_model, "")
        if ret != 0
            throw(XpressError(model_src))
        end
    end
    model_dest
end

# read / write file
#=
XPRSprob a-> Ptr{Void}
const char *a -> Ptr{UInt8}  -> "" : ASCIIString
int -> Cint
const char a[] -> Ptr{Cchar} -> Cchar[]
=#
function load_empty(model::Model)
#=int XPRS_CC XPRSloadlp(XPRSprob prob, const char *probname, int ncol, int
      nrow, const char qrtype[], const double rhs[], const double range[],
      const double obj[], const int mstart[], const int mnel[], const int
      mrwind[], const double dmatval[], const double dlb[], const double
      dub[]);
=#
    ret = @xprs_ccall(loadlp,Cint,
        (Ptr{Void},
            Ptr{UInt8},
            Cint,
            Cint,
            Ptr{Cchar},#type
            Ptr{Float64},#rhs
            Ptr{Float64},#range
            Ptr{Float64},#obj
            Ptr{Cint},#mastart
            Ptr{Cint},#mnel
            Ptr{Cint},#mrwind
            Ptr{Float64},#dmat
            Ptr{Float64},#lb
            Ptr{Float64}#ub
            ),model.ptr_model,"",0,0,Cchar[],Float64[],Float64[],Float64[],
                Cint[],Cint[],Cint[],Float64[],Float64[],Float64[]
    )
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end
function read_model(model::Model, filename::Compat.ASCIIString)
    #@assert is_valid(model.env)
    flags = ""
    ret = @xprs_ccall(readprob, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function write_model(model::Model, filename::Compat.ASCIIString)
    flags = ""
    ret = @xprs_ccall(writeprob, Cint, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

# function presolve_model(model::Model)
#     ret = @grb_ccall(presolvemodel, Ptr{Void}, (Ptr{Void},), model.ptr_model)
#     if ret == C_NULL
#         # Presumably failed to return a model
#         error("presolve_model did not return a model")
#     end
#     return Model(model.env, ret)
# end


text:
include("HW7.jl")
l_0 = 0.3
l_1 = 0.6
l_2 = 0.4
## Answer for 3.b
twist_1 = [0 0 1 0 0 0]';
twist_2 = [0 0 1 l_1 0 0]';
twist_3 = [0 0 1 l_1+l_2 0 0]';
twist_4 = [0 0 0 0 0 1]';
twists = [twist_1 twist_2 twist_3 twist_4];
R = Matrix(I,3,3)
t = [0 l_1+l_2 l_0];
T=RpToTrans(R,t)
angles=[pi/3 -pi/2 pi/4 -0.35]
## Answer for 3.b
g = FwdKin(T,twists,angles)
## Answer for 3.c
Js = JacobianSpace(twists,angles)
## Answer for 3.d
l_rod = 0.15;
P_rod = [0 l_1+l_2+l_rod l_0 1]';
Q_dot = [-0.2 1 0.5 -2]';

V_s = JacobianSpace(twists,angles)*Q_dot;
V_rod = VecTose3( V_s ) * P_rod;

max_stars_repo_name: veysiadn/AdvancedRoboticsJulia

max_issues_repo_name: veysiadn/AdvancedRoboticsJulia

max_forks_repo_name: veysiadn/AdvancedRoboticsJulia


text:
#!/usr/bin/env julia

# Draw lines with various dash styles.

using Compose
using Color

function draw_lines(dash_patterns)
    if length(dash_patterns) == 0
        canvas()
    else
        compose(compose(canvas(), lines((0, 0), (1, 0)), strokedash(dash_patterns[1])),
                compose(canvas(0, 0.1), draw_lines(dash_patterns[2:end])))
    end
end

patterns = Array[
    [5mm, 5mm],
    [5mm, 10mm],
    [10mm, 5mm],
    [5mm, 1mm],
    [1mm, 5mm],
    [0.9mm],
    [15mm, 10mm, 5mm],
    [15mm, 10mm, 5mm, 10mm],
    [15mm, 10mm, 5mm, 10mm, 15mm],
    [5mm, 5mm, 1mm, 5mm]]

c = draw_lines(patterns)

img = PDF("dash.pdf", 4inch, 4(sqrt(3)/2)inch)
draw(img, compose(c, stroke(color("black")), linewidth(1mm)))


text:
export molarmass
# function to calculate the molar mass of an ionic liquid
function molarmass(anion, cation, pdb_dir)
  atcation = readPDB("$pdb_dir/$(cation)_VSIL.pdb")
  atanion  = readPDB("$pdb_dir/$(anion)_VSIL.pdb")
  MMIL = mass(atcation) + mass(atanion)
  return MMIL
end


text:
function example_blank(nwin)
    for i = 1:nwin
        Pages.launch("http://localhost:$(Pages.port)/examples/blanks")
        wait(Pages.conditions["connected"])
    end

    Pages.broadcast("/examples/blank","say","Hello everyone!")
    for (sid,s) in Pages.connections["/examples/blank"]
        Pages.message("/examples/comms",sid,"say","You are connection #$(sid).")
    end
end


text:
module batch
srcDir  = dirname(@__FILE__)*"\\";

include(srcDir*"Model\\R_sep_Model.jl")
include(srcDir*"Model\\C_sep_Model.jl")
include(srcDir*"Model\\Utils_Model.jl")
include(srcDir*"sep_Compute.jl")
using .R_sep_Model
using .C_sep_Model
using .Utils_Model
using .sep_Compute

using JuMP

export batch_model
default_sdir = dirname(dirname(@__FILE__))*"\\assets\\bounds\\"
""" Produces models for a set of constraints."""
function batch_model(t,ρdic,df;sdir = default_sdir,con_list= ["S∞sG","S₂sG","S₂₁sG"])
    for ex in keys(ρdic)
        ρ = ρdic[ex]
        d = filter(:Example => ==(ex),df).Size[1]
        for con ∈  con_list, isR ∈ ['R', 'C']
            (!filter(:Example => ==(ex),df).isReal[1] && (isR == 'R')) ? continue : 0
            println("Currently $isR-modeling example $ex with constraints $con")
            sep_mod = (isR == 'R') ? R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list=con) : C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list=con,noBlock = true)

            Utils_Model.export_model(sep_mod, sdir*"$(ex)_$(con)_$(isR)_$t.dat-s")
        end
    end
end

# """ Reads all .dat-s files and computes bounds storing result summary in a .csv file"""
function batch_Computeξₜˢᵉᵖ(boundsDir = default_sdir)
    datsFiles = [file for file in readdir(boundsDir,join = true) if contains(file,".dat-s")]
    file_loc = boundsDir*"Summary.csv"
    touch(file_loc)
    open(file_loc,"a") do io
        write(io, "|Example|Model|Constraint|Primal|Dual|obj_val|\n")
    end
    for file in datsFiles
        name = split(basename(file),".")[1]
        ex   = split(name,"_")[1]
        con  = split(name,"_")[2]
        num  = split(name,"_")[3]
        # try
        mod  = Utils_Model.read_model(file)
        omod = sep_Compute.Computeξₜˢᵉᵖ(mod)
        # catch
        #     open(file_loc,"a") do io; write(io, "|$ex|$con|err.|err.|err.| \n"); end
        #     continue
        # end
        pstat   = JuMP.primal_status(omod)
        dstat   = JuMP.dual_status(omod)
        ov_temp = round(JuMP.objective_value(omod),digits=3)
        ov = (string(pstat) == string(dstat) == "FEASIBLE_POINT") ?
             ov_temp :
             (contains(string(pstat)*string(dstat),"CER") ? "*" : NaN )

        open(file_loc,"a") do io
            write(io, "|$ex|$num|$con|$pstat|$dstat|$ov|\n")
        end
    end
end

using CSV
using DataFrames

function unstack_constraints(df,bDir = default_sdir)
    temp_df  = CSV.read("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\assets\\bounds\\"*"Summary.csv", DataFrame)
    temp2_df = select(innerjoin(df, temp_df, on = :Example),[:Example,:isReal,:isSeparable,:Model,:Size,:Bi_rank,:Constraint,:obj_val])
    temp4_df = unstack(temp2_df, :Model, :obj_val)

    temp4R_df = select(temp4_df,[:Example,:isReal,:isSeparable,:Size,:Bi_rank,:Constraint,:R])
    temp4Rus_df = unstack(temp4R_df, :Constraint, :R)
    rename!(temp4Rus_df,Dict(:S₂sG => :RS₂sG, :S₂₁sG => :RS₂₁sG, :S∞sG=>:RS∞sG));

    temp4C_df = select(temp4_df,[:Example,:isReal,:isSeparable,:Size,:Bi_rank,:Constraint,:C])
    temp4Cus_df = unstack(temp4C_df, :Constraint, :C)
    rename!(temp4Cus_df,Dict(:S₂sG => :CS₂sG, :S₂₁sG => :CS₂₁sG, :S∞sG=>:CS∞sG));

    temp5_df = innerjoin(temp4Cus_df,select(temp4Rus_df,[:Example,:RS₂sG, :RS₂₁sG, :RS∞sG]), on = :Example)

    CSV.write(bDir*"SummaryUnstacked.csv", temp5_df, delim="&")
    return temp5_df
end

end  # module
## Selecting all real examples
# df_real = filter(:isReal => ==(true),df)
# examples_real = Dict()
# for exa in df_real.ex
#     examples_real[exa] = examples_all[exa]
# end
# ##
# save_dir = sep_rank_proj_path*"\\PreRunModels\\R-models\\t=$(t[1])\\"
# Utils_Model.batch_model(t,examples_real,df_real,save_dir,["S₁sG","S₂sG","S₃sG"],true)
# sep_Compute.batch_Computeξₜˢᵉᵖ(save_dir)
# unstack_constraints(save_dir,df_real)
# ##
# temp_df  = CSV.read("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\PreRunModels\\C-models\\t=3\\SummaryUnstacked.csv", DataFrame, delim = "&")
# temp_dfr = CSV.read("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\PreRunModels\\R-models\\t=3\\SummaryUnstacked.csv", DataFrame, delim = "&")
#
# poes_df = select(innerjoin(temp_dfr, temp_df, on = :ex,makeunique=true),[:ex,:siz,:birank,:S₁sG,:S₂sG,:S₃sG,:S₁sG_1,:S₂sG_1,:S₃sG_1,])
# rename!(poes_df ,Dict(:S₁sG_1 => :RS₁sG, :S₂sG_1 => :RS₂sG, :S₃sG_1 => :RS₃sG))
# CSV.write("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\PreRunModels\\t=2.csv", poes_df, delim="&")

# function quick_run_spec_example(ex,t = (2,2),con_list = "S₁sG",isReal = true)
#     examples     = Examples.get_examples()
#     examples_all = merge(examples["ent"],examples["sep"])
#     df = Examples.get_example_overview(false)
#
#     ρ = examples_all[ex]
#     ρ_meta_data = filter(:ex => ==(ex),df)
#     d =  ρ_meta_data.siz[1]
#     if isReal
#         Lx, mod = R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t,con_list)
#     else
#         Lx, mod = C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t,con_list)
#     end
#     mod_opt = sep_Compute.Computeξₜˢᵉᵖ(mod)
#     return JuMP.values(Lx), mod_opt
# end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
import Pkg
Pkg.instantiate()

Pkg.add("UrlDownload")
Pkg.add("CSV")
Pkg.add("DataFrames")

using UrlDownload, DataFrames, CSV

# url = "https://covidtracking.com/api/v1/states/current.json"
# df = download(url,  "data/covid-current.json")

url = "https://covidtracking.com/api/v1/states/current.json"

df = urldownload(url) |> DataFrame

# url = "https://www.stats.govt.nz/assets/Uploads/Business-price-indexes/Business-price-indexes-December-2019-quarter/Download-data/business-price-indexes-december-2019-quarter-csv.csv"
# df = urldownload(url)
# df = urldownload(url, true)


text:
using Devito, BenchmarkTools

const SUITE = BenchmarkGroup()

nz,ny,nx = 102,101,100
grd = Grid(shape=(nz,ny,nx))
f = Devito.Function(name="f", grid=grd, space_order=8)
d = data(f)
dhalo = data_with_halo(f)
dinhalo = data_with_inhalo(f)
dalloc = data_allocated(f)

_d = zeros(eltype(d), size(d))
_dhalo = zeros(eltype(dhalo),  size(dhalo))
_dinhalo = zeros(eltype(dinhalo), size(dinhalo))
_dalloc = zeros(eltype(dalloc), size(dalloc))

SUITE["function"] = BenchmarkGroup()
SUITE["function"]["data"] = @benchmarkable data($f)
SUITE["function"]["data with halo"] = @benchmarkable data_with_halo($f)
SUITE["function"]["data with inhalo"] = @benchmarkable data_with_inhalo($f)
SUITE["function"]["data allocated"] = @benchmarkable data_allocated($f)
SUITE["function"]["data copy!"] = @benchmarkable copy!($_d, $d)
SUITE["function"]["data with halo copy!"] = @benchmarkable copy!($_dhalo, $dhalo)
SUITE["function"]["data with inhalo copy!"] = @benchmarkable copy!($_dinhalo, $dinhalo)
SUITE["function"]["data allocated copy!"] = @benchmarkable copy!($_dalloc, $dalloc)

stf = SparseTimeFunction(name="stf", grid=grd, npoint=1, nt=1000, coordinates=[50 50 50])
tf = TimeFunction(name="tf", grid=grd, time_order=8, space_order=8)

SUITE["sparse time function"] = BenchmarkGroup()
SUITE["sparse time function"]["data"] = @benchmarkable data($stf)
SUITE["sparse time function"]["data with halo"] = @benchmarkable data_with_halo($stf)
SUITE["sparse time function"]["data with inhalo"] = @benchmarkable data_with_inhalo($stf)
SUITE["sparse time function"]["data allocated"] = @benchmarkable data_allocated($stf)
SUITE["sparse time function"]["inject"] = @benchmarkable inject($stf; field=$(forward(tf)), expr=$stf)
SUITE["sparse time function"]["interpolate"] = @benchmarkable interpolate($stf; expr=tf)

SUITE


text:
# graphstructure.jl
# Authors : Heitor Guimaraes and Luiz Ciafrino
# @brief: Module with some optimization algorithms to solve graph coloring

#=using JuMP, Gurobi

function direct_coloring(G::Array{Array{Int64},1})
	l = length(G)
	colors = zeros(Int,l)
	vertex_color = zeros(Int,l)

	colors[1] = 1
	vertex_color[1] = 1

	for i in 2:l
		unavailable = 0
		for j in G[i]
			if vertex_color[j] >= 1 && vertex_color[j] > unavailable
				unavailable = vertex_color[j]
			end
		end
		vertex_color[i] = unavailable+1
		colors[unavailable+1] = 1

	end
	return colors
end

function gurobi_coloring(G::Array{Array{Int64},1}, E, tmax, r)

	n = length(G)
	h = convert(Int, length(Set(dsatur(G))))
	println("DSATUR greedy coloring $h")
	if h-r > 0
		h = h-r
	end
	println("Solving for DSATUR - $r colors = $h")

	m = Model(solver=GurobiSolver(TimeLimit = tmax, Method=-1, MIPFocus=1, Cuts=3, MIPGap=0.1))

	@defVar(m, x[1:n, 1:h], Bin)
	@defVar(m, y[1:h], Bin)

	@addConstraint(m, constr[i=1:n], sum{x[i,k], k=1:h} == 1 )
	println("Added constrain 1")

	@addConstraint(m, xyconstr[i=1:n, k=1:h], x[i,k] - y[k] <= 0)
	println("Added constrain 2")

	for e in E
		for k in 1:h
			a= e[1]
			b= e[2]
			@addConstraint(m, x[a,k] + x[b,k] <= 1)
		end
	end
	println("Added constrain 3")

	@setObjective(m, Min, sum{y[i], i=1:h})

	@time status = solve(m)

	return getObjectiveValue(m), getValue(y), getValue(x), status
end

function edges(G::Array{Array{Int64},1})
	E = Tuple{Int,Int}[]
	n = length(G)
	for i in 1:n
		for j in G[i]
			push!(E, (i,j))
		end
	end
	return E
end

function load(name)
	println("Loading graph $name")
	g = graphtheoryufrj
	G = g.simpleGraph(name)
	G = G.Graph
	E=edges(G)
	return G,E
end

function run(name, tmax, r)
	println("Coloring graph $name")
	g = graphtheoryufrj
	G = g.simpleGraph(name)
	G = G.Graph
	E=edges(G)

	obj,y,x,status = gurobi_coloring(G, E, tmax, r)
	println("Objective value of $obj")
	return obj,y,x,status
end

function colors(x::Array{Float64,2})
	ret = Array{Int}(size(x)[1])
	for i in 1:size(x)[1]
		for j in 1:size(x)[2]
			if x[i,j] == 1
				ret[i]=j
			end
		end
	end
	ret2 = Array{Int}(size(x)[1])
	c = 0
	for i in Set(ret)
		for j in find(ret .== i)
			ret2[j] = c
		end
		c+=1
	end
	return ret, ret2
end

function coloring_correct(E, color)

	err = 0
	for i in E
		if color[i[1]] == color[i[2]]
			err+=1
		end
	end

	if err == 0
		return true
	else
		return err
	end

end

function color_report(i, tmax, r)
	name = "G$(i).txt"

	t = @elapsed obj,y,x,status = run(name, tmax, r)
	_,sim = colors(x)

	f = open("$(time())color_$(name)", "w")
	write(f, "optimization status = $(status)\n")
	write(f, "chromatic number = $(obj)\n")
	write(f, "runtime = $(t)\n")
	write(f, "$(sim)\n")
	close(f)
end

# l = ["dsjc1000_9.txt", "dsjc1000_5.txt", "dsjc1000_1.txt", "dsjc500_9.txt", "dsjc500_5.txt", "dsjc500_1.txt"]
# for i in l
# 	name = i
# 	g = graphtheoryufrj
# 	println("Graph $(name)")
# 	@time G = g.simpleGraph(name)
# 	G = G.Graph
# 	E=edges(G)
# 	t = @elapsed d = length(Set(dsatur(G)))
# 	println("K = $(d), time = $(t)")
# 	println("=========================================")
# end
# l = []

# for i in 1:20
# 	name = "G$(i).txt"
# 	g = graphtheoryufrj
# 	G = g.simpleGraph(name)
# 	G = G.Graph
# 	E=edges(G)
# 	d = length(Set(dsatur(G)))

# 	push!(l, (d*length(E), name))
# end

# # DRIVER PROGRAM
println("Compile run")
run("g.txt", 10, 0)
println("==================================================")
cd("t3")
color_report(4, 60*10, 1)
color_report(6, 60*10, 1)

## RUN AND SAVE JUST GREEDY REPORT
# cd("t3")
# name = "G18.txt"
# g = graphtheoryufrj
# G = g.simpleGraph(name)
# G = G.Graph
# E=edges(G)
# t = @elapsed d = dsatur(G)
# obj = length(Set(d))
# f = open("color_$(name)", "w")
# write(f, "optimization status = greedy\n")
# write(f, "chromatic number = $(obj)\n")
# write(f, "runtime = $(t)\n")
# write(f, "$(d)\n")
# close(f)
=#


text:
using Revise
using Stheno, BenchmarkTools, Bench, Flux, Zygote
using Zygote: gradient
using BenchmarkTools: median

# Various bits of utility for benchmarking.
include("util.jl")

println("Constructing")
Bench.clear_suites()
benchmarks = @benchset "Stheno" begin

    # @benchset "mean_and_kernel" begin
    #     include("mean_and_kernel/mean.jl")
        include("mean_and_kernel/kernel.jl")
        # include("mean_and_kernel/compose.jl")
        # include("mean_and_kernel/finite.jl")
    # end

    # @benchset "gp" begin
    #     include("gp/abstract_gp.jl")
    # end
end

println("Tuning")
tune!(benchmarks; verbose=true)
println("Running")
results = run(benchmarks; verbose=true)

println("Time:")
pretty_print(median(results))

println()
println("Memory:")
pretty_print(memory(results))

println()
println("Allocs:")
pretty_print(allocs(results))


text:
using EchelleCCFs
using Test

    @testset "Linelists/mask I/O" begin
        mask_fn = joinpath(pkgdir(EchelleCCFs),"data","masks","G2.espresso.mas")
        @test isfile(mask_fn)
        @test_nowarn EchelleCCFs.read_linelist_espresso(mask_fn)
        mask_fn = joinpath(pkgdir(EchelleCCFs),"data","line_lists","VALD_Fe1_DP_rejectTelluricSlope0.0_badLineFilterESPRESSO-strict-NEID-BIS_overlapcutoff6e-05_depthcutoff0.05_allowBlends0_wavesReiners_depthssolar_nbin1depth0.mas")
        @test isfile(mask_fn)
        @test_nowarn EchelleCCFs.read_mask_vald(mask_fn)
        mask_fn = joinpath(pkgdir(EchelleCCFs),"data","masks","HD101501_f75s75v3_mask.csv")
        @test isfile(mask_fn)
        @test_nowarn EchelleCCFs.read_linelist_rvspectml(mask_fn)
    end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
load("Distributions")
using Distributions
load("MCMC")
using MCMC

load("benchmarks/julia/densities.jl")

function mcmc(log_density::Function, x::Vector, sampler::Function, n_chains::Int64, n_samples::Int64)
    xs = zeros(n_samples, length(x))
    gx = log_density(x)
    for iter = 1:n_samples
        x, gx = sampler(log_density, x, gx)
        xs[iter,:] = x
    end
    return xs
end

function multivariate_slice(f::Function, x::Vector, sampler::Function)
    P = length(x)
    fx = f(x)
    for p = 1:length(x)
        # Keep all but the p'th value fixed
        function fp(w)
            x[p] = w
            f(x)
        end

        # Pass this to a univariate slice sampler
        x[p], fx = slice_sampler(x[p], fp, fx)
    end
    x, fx
end
mv_slice(f::Function, x, fx) = multivariate_slice(f, x, slice_sampler)
mv_slice(f::Function, x, fx) = multivariate_slice(f, x, slice_sampler)

function run_australian()
    d = csvread("benchmarks/csv/australian.csv")
    y = d[:,1]
    x = d[:,2:end]
    theta = rand(MultivariateNormal([0. for i in 1:14], .000001 *  eye(14)))
    theta[13] = 0
    samples = mcmc(logistic_regression, theta, mv_slice, 1, 10)
end


text:
function rules_crossover!(Γ₁::ClassificationRules, Γ₂::ClassificationRules, Γ₁′::ClassificationRules, Γ₂′::ClassificationRules)
    i = rand(1:length(Γ₁))
    j = rand(1:length(Γ₂))
    Γ₁′ = deepcopy(Γ₁)
    Γ₂′ = deepcopy(Γ₂)
    _rules(Γ₁′)[i] = _rules(Γ₂)[j]
    _rules(Γ₂′)[j] = _rules(Γ₁)[i]
    return
end

# levels: 50% 30% 10% 10%
#   1) rules level:
#       a) add rule, and
#       b) remove rule
#   2) rule level that replaces a subformula with another subformula
#   3) class level that changes the class
#   4) horizon level
function rules_mutation!(Γ::ClassificationRules)
    level = rand(1:100)

    # rules level
    if level ≤ 50
        # add rule
        if rand(1:2) == 1
            if length(Γ) < Γ.maxnumrules
                push!(_rules(Γ), _rand_rule(Γ.init_relations,Γ.relations,Γ.algebra,Γ.ds,Γ.orders;minmd=Γ.minmd,maxmd=Γ.maxmd,maxdepth=Γ.maxdepth,isroot=true))
            end
            return
        # remove rule
        else
            if length(Γ) > Γ.minnumrules
                deleteat!(_rules(Γ), rand(1:length(Γ)))
                return
            end
        end
    # rule level
    elseif 50 < level && level ≤ 80

        i = rand(1:length(Γ))
        rule = _rules(Γ)[i]

        nodes = _collect_nodes(rule)
        s = reduce(+,[_size(ν) for ν in nodes])
        w = [1 - _size(ν)/s for ν in nodes]./(length(nodes)-1)
        sf = StatsBase.sample(nodes, Weights(w))
        #sf = nodes[rand(1:length(nodes))]

        # upwards modal depth of the selected subformula
        umd = _upwards_md(sf)
        # # upwards modal depth is the previously computed one if the subformula is a leaf;
        # # otherwise, it is minus 1 because we need to replace an internal node which has a modality by definition
        # umd = (_is_leaf(sf) ? umd : umd-1)

        udepth = _upwards_depth(sf)

        # new subtree
        subtree = isdefined(sf, :parent) ?
            _rand_antecedent(Γ.init_relations,Γ.relations,Γ.algebra,Γ.ds,Γ.orders;minmd=Γ.minmd-umd,maxmd=Γ.maxmd-umd,maxdepth=Γ.maxdepth-udepth,isroot=false) :
            _rand_antecedent(Γ.init_relations,Γ.relations,Γ.algebra,Γ.ds,Γ.orders;minmd=Γ.minmd-umd,maxmd=Γ.maxmd-umd,maxdepth=Γ.maxdepth-udepth,isroot=true)
        # set parent node, if exists
        if isdefined(sf, :parent)
            if isdefined(sf.parent, :left) && sf.parent.left == sf
                sf.parent.left = subtree
            elseif isdefined(sf.parent, :right) && sf.parent.right == sf
                sf.parent.right = subtree
            end
            _parent!(subtree, sf.parent)
        else
            _antecedent!(Γ[i],subtree)
        end
    # class level
    elseif 80 < level && level ≤ 90
        i = rand(1:length(Γ))
        C = Γ.ds.unique_classes
        c = C[rand(1:length(C))]

        # TODO: farebbe loop infinito..
        while _consequent(_rules(Γ)[i]) == c
            c = C[rand(1:length(C))]
        end

        _rules(Γ)[i].consequent = c
        return
    # horizon level
    else
        # reduce horizon value by 1, if possible
        if rand(1:2) == 1
            if _horizon(Γ) > Γ.minh
                _horizon!(Γ, _horizon(Γ) - 1)
            end
        # increase horizon value by 1, if possible
        else
            if _horizon(Γ) < Γ.maxh
                _horizon!(Γ, _horizon(Γ) + 1)
            end
        end
    end
end


text:
function bit_sint(bit)
    2^(bit - 1)
end
function bit_uint(bit)
    2^bit
end
function main(arg)
    bit_num = big(arg)
    ans_csv = "Bit,Max_sInt,max_uInt"
    for j in 1:bit_num
        s = 2^(j - 1)
        ans_csv *=  "\n" * string(s) * "," * string(bit_sint(s)) * "," * string(bit_sint(s))
    end
    w = open("max_integer_of_bit" * bit_num * ".csv", "w")
    write(w, ans_csv)
    close(w)
end
map(x -> main(parse(Int, x)), ARGS)


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

using Markdown
using InteractiveUtils

# ╔═╡ b0602873-cda1-4d30-8ec8-1b9f7fdcdb07
using Latexify, PlutoTest

# ╔═╡ eac55bba-cc75-4632-832a-dbb06983be45
md"""
# [Day 2: I Was Told There Would Be No Math](https://adventofcode.com/2015/day/2)
The elves are running low on wrapping paper, and so they need to submit an order for more. They have a list of the dimensions (length `l`, width `w`, and height `h`) of each present, and only want to order exactly as much as they need.

Fortunately, every present is a box (a perfect [right rectangular prism](https://en.wikipedia.org/wiki/Cuboid#Rectangular_cuboid)), which makes calculating the required wrapping paper for each gift a little easier: find the surface area of the box, which is `2*l*w + 2*w*h + 2*h*l`. The elves also need a little extra paper for each present: the area of the smallest side.

For example:

-   A present with dimensions `2x3x4` requires `2*6 + 2*12 + 2*8 = 52` square feet of wrapping paper plus `6` square feet of slack, for a total of `58` square feet.
-   A present with dimensions `1x1x10` requires `2*1 + 2*10 + 2*10 = 42` square feet of wrapping paper plus `1` square foot of slack, for a total of `43` square feet.

**All numbers in the elves' list are in feet. How many total square feet of wrapping paper should they order?**
"""

# ╔═╡ 0568b194-e70a-42b5-a9e8-89e710e61648
puzzle_input = """29x13x26
11x11x14
27x2x5
6x10x13
15x19x10
26x29x15
8x23x6
17x8x26
20x28x3
23x12x24
11x17x3
19x23x28
25x2x25
1x15x3
25x14x4
23x10x23
29x19x7
17x10x13
26x30x4
16x7x16
7x5x27
8x23x6
2x20x2
18x4x24
30x2x26
6x14x23
10x23x9
29x29x22
1x21x14
22x10x13
10x12x10
20x13x11
12x2x14
2x16x29
27x18x26
6x12x20
18x17x8
14x25x1
30x15x22
17x18x7
28x23x24
15x12x25
14x7x20
29x23x8
24x5x22
6x22x8
1x15x26
14x5x1
24x28x28
17x23x23
4x15x7
23x8x11
6x15x1
23x18x13
17x1x26
23x13x17
2x18x8
22x22x1
10x22x6
28x29x20
22x21x25
14x8x23
12x30x14
8x7x5
3x30x15
4x3x29
25x18x3
16x7x16
4x3x8
9x16x30
20x28x3
28x24x6
4x18x2
23x18x5
22x4x30
15x30x9
7x12x12
3x22x29
12x1x9
9x2x25
17x11x10
25x24x7
7x27x26
26x4x12
29x2x26
19x24x12
23x23x3
26x28x16
18x4x16
25x30x18
29x19x19
16x3x27
29x25x29
18x19x5
14x21x30
19x13x26
19x10x15
9x4x7
18x6x6
24x25x29
9x12x27
15x3x22
30x17x21
18x19x28
9x11x12
8x28x22
11x3x4
28x17x20
24x18x15
11x12x13
6x19x24
28x4x5
28x22x23
13x29x2
9x16x15
29x28x1
10x18x30
19x11x12
26x28x25
23x17x13
25x1x21
17x1x27
17x27x28
28x13x15
14x13x25
11x29x7
22x29x5
13x6x14
23x18x13
25x7x17
18x9x20
21x11x2
28x11x13
13x25x1
19x29x25
16x29x4
10x21x10
7x25x17
5x9x3
1x15x6
8x27x29
23x6x30
22x22x29
6x20x30
26x25x29
10x19x19
20x30x9
5x30x24
17x10x27
30x14x30
8x17x4
7x18x6
3x5x4
24x17x15
14x20x17
22x27x15
18x14x15
23x9x11
21x16x29
7x18x21
9x3x29
10x13x4
2x30x4
23x20x4
8x22x21
29x28x4
13x16x25
21x9x11
7x26x26
13x23x30
19x7x10
9x23x21
21x9x17
9x21x15
20x29x22
23x13x15
19x25x2
12x11x30
20x21x6
21x6x17
24x26x9
29x21x29
29x26x16
6x16x1
2x12x6
6x7x20
7x2x22
6x22x4
13x11x27
25x27x14
11x8x6
26x11x14
30x3x29
27x21x20
15x16x26
6x22x10
11x9x25
23x13x6
13x9x3
30x22x13
29x23x14
25x19x6
7x29x11
19x18x5
29x25x13
25x24x27
1x9x12
22x9x17
14x12x28
19x21x17
13x25x17
14x25x12
4x14x30
7x15x28
3x6x25
6x2x16
15x19x11
17x30x20
20x23x7
26x21x6
26x29x24
2x4x30
4x22x18
13x3x28
27x6x21
5x3x27
12x7x11
28x11x9
12x9x2
1x22x20
15x13x28
14x19x16
28x20x3
20x4x9
26x7x26
18x19x25
7x1x13
20x23x29
27x26x8
11x15x15
10x21x23
29x2x11
21x28x20
3x18x23
26x17x17
14x26x17
20x7x17
18x12x8
4x8x8
8x15x23
24x29x5
1x25x8
1x28x17
16x18x13
29x24x22
13x16x10
14x7x16
15x11x29
12x15x19
17x6x28
4x3x9
15x16x8
29x27x11
2x24x20
4x21x3
29x24x27
18x22x22
7x8x18
20x7x8
19x9x2
20x17x2
2x29x10
19x25x1
28x9x3
29x27x20
7x21x7
10x4x22
26x8x5
26x14x1
5x27x9
2x18x3
3x27x17
30x17x23
30x11x20
4x6x7
6x29x27
30x16x20
24x30x28
19x20x26
18x1x25
26x12x12
19x15x29
16x21x24
23x13x26
25x16x10
8x9x18
24x14x1
24x15x21
19x9x14
8x23x11
22x2x16
29x9x26
3x16x25
15x20x30
3x11x12
15x2x3
13x7x4
2x7x27
9x26x11
30x24x19
28x17x21
10x8x2
11x15x26
10x12x20
24x24x27
25x26x16
13x4x20
25x13x11
12x22x3
20x7x1
12x18x6
26x8x20
14x2x7
23x12x1
26x24x24
27x26x23
26x17x5
17x24x2
26x5x6
23x5x1
5x18x30
24x21x19
5x28x11
21x20x14
25x4x22
26x24x11
7x5x8
13x1x30
5x1x6
14x5x2
8x11x7
13x20x1
17x30x14
29x22x10
12x26x3
27x17x3
26x27x4
5x26x17
22x11x19
8x26x3
24x19x22
7x1x4
6x27x30
4x28x14
16x14x18
4x5x20
19x25x4
15x15x1
10x14x14
16x18x24
21x27x15
5x5x10
1x7x13
16x2x8
13x15x11
3x25x10
20x29x8
12x3x2
10x13x12
25x27x1
11x30x19
7x19x13
27x6x18
16x21x19
21x29x5
16x23x12
29x19x15
5x5x10
27x15x1
13x16x22
29x19x5
8x12x9
3x18x5
13x25x3
5x9x21
10x20x16
9x9x11
23x21x1
22x2x15
27x8x13
23x7x3
26x30x15
29x15x16
16x27x13
2x18x9
10x27x8
20x9x25
10x2x17
16x13x13
21x26x1
27x26x24
9x30x16
19x17x28
25x15x1
10x26x6
10x11x11
5x26x25
30x4x15
9x8x23
14x25x7
8x28x8
28x18x24
4x4x25
16x25x11
17x27x8
15x16x9
24x13x21
17x3x27
27x5x26
8x27x12
29x2x8
24x23x30
1x30x21
6x18x20
13x14x12
25x30x23
24x6x24
12x7x21
11x6x8
8x30x30
26x3x12
28x6x5
18x7x1
7x6x20
14x16x18
11x22x15
4x20x10
19x24x19
8x24x11
4x9x10
6x6x22
10x9x29
1x5x28
19x25x29
20x30x3
15x13x13
9x9x24
20x14x29
26x24x13
2x25x8
10x26x2
12x19x12
18x6x20
4x5x14
26x27x10
16x26x20
3x21x15
2x26x18
14x11x17
26x26x25
10x1x11
17x19x19
27x28x26
9x2x10
19x30x15
23x30x14
15x3x20
2x14x22
21x18x8
22x4x29
19x6x29
9x26x29
16x10x9
22x12x22
13x28x14
25x14x28
28x3x30
10x17x1
10x27x22
10x23x19
14x25x9
11x24x8
30x25x10
22x13x28
2x7x6
11x20x8
9x22x14
19x16x9
11x24x4
11x17x2
6x4x10
26x10x10
12x14x5
27x10x3
15x3x6
11x7x19
22x10x12
21x26x10
13x20x3
27x8x8
1x24x23
24x9x22
23x17x23
3x28x19
2x20x28
23x17x24
26x1x4
4x1x12
5x6x16
13x22x13
25x21x21
20x21x12
9x24x25
17x16x12
12x28x9
18x16x27
29x12x2
30x12x15
24x11x10
4x9x22
4x24x5
19x11x5
6x25x6
1x20x17
22x8x21
11x26x4
16x19x3
8x12x8
13x2x18
10x5x11
8x12x17
21x2x5
26x17x26
23x18x17
28x11x14
1x4x27
29x5x28
5x9x10
5x7x25
20x15x27
15x11x17
12x14x1
29x14x4
18x14x18
14x25x24
26x14x18
13x8x11
30x1x23
3x4x12
12x24x9
8x6x16
14x15x30
12x30x8
22x11x18
16x30x28
17x18x4
13x14x23
2x28x8
3x28x30
29x30x8
4x6x26
6x30x17
11x30x30
19x4x3
12x15x20
22x28x4
26x30x2
6x12x7
1x10x5
25x29x7
17x9x18
16x21x29
21x14x7
15x16x11
26x6x15
8x24x7
2x20x4
2x9x3
19x8x13
18x7x22
27x14x17
2x13x8
18x15x26
15x27x27
18x11x15
1x29x20
21x12x11
20x2x15
28x23x9
1x1x17
7x23x9
30x9x27
9x16x18
15x24x28
30x11x18
29x26x10
9x5x25
2x1x19
14x3x14
6x3x6
30x15x20
20x17x27
28x10x9
14x24x28
17x11x6
12x3x6
8x8x15
23x14x21
11x21x7
5x13x30
4x29x25
30x28x24
18x4x9
3x15x6
13x9x19
30x14x7
7x9x9
17x11x26
24x26x13
16x21x16
27x17x25
2x21x11
9x11x27
3x3x7
13x8x14
20x20x26
13x29x22
30x11x1
7x10x19
27x5x9
23x17x15
21x6x13
24x15x16
18x4x14
18x16x6
22x11x18
14x2x5
15x3x7
10x20x29
16x1x10
30x23x1
10x15x11
17x14x5
22x8x13
7x11x28
26x17x3
2x23x2
28x13x19
18x12x28
22x23x16
14x12x1
20x8x19
17x19x13
29x2x12
2x26x27
29x16x4
13x8x18
16x15x30
23x16x2
28x8x27
21x8x23
13x20x26
19x6x17
17x30x15
7x4x30
2x13x30
18x7x19
4x13x27
8x6x5
18x20x25
2x3x30
23x27x13
22x30x4
23x25x25
23x16x19
25x3x1
5x6x15
11x29x12
25x24x7
16x7x20
20x3x2
12x27x15
16x10x12
1x3x14
22x1x26
2x24x18
11x29x16
15x2x9
10x1x24
21x8x11
30x11x23
6x30x21
13x27x29
14x6x5
18x29x19
12x4x28
29x3x14
10x30x28
5x7x15
14x1x10
9x25x14
7x24x18
28x17x21
18x13x25
26x15x1
21x1x19
12x16x21
4x6x13
7x15x26
17x19x5
12x28x2
1x20x19
27x7x5
17x26x8
12x15x19
5x23x10
8x2x8
16x13x12
14x27x1
26x29x3
24x16x14
14x13x13
7x22x23
2x9x30
4x27x8
26x27x15
23x1x6
25x29x18
5x18x1
20x8x20
5x10x25
30x25x15
7x22x25
28x26x17
29x4x1
21x11x27
20x9x8
25x22x12
2x11x11
23x2x16
23x27x20
2x13x28
27x2x24
11x1x17
12x4x27
16x20x22
30x12x10
5x15x4
5x2x27
12x4x25
1x16x4
27x4x4
21x16x3
27x26x3
24x6x6
24x12x12
20x20x25
8x29x2
21x4x5
2x4x8
4x13x19
3x20x10
12x15x16
6x5x4
12x16x20
22x19x17
8x17x22
25x16x15
7x1x19
10x1x7
23x23x5
28x6x12
2x25x12
10x27x12
24x27x19
14x14x20
4x1x5
16x27x29
20x20x24
28x24x30
6x15x15
9x15x30
23x26x3
17x24x21
22x25x25
18x29x10
20x25x1
24x11x16
20x7x21
20x7x9
7x26x2
5x18x1
16x26x28
4x10x18
27x30x21
26x9x9
8x16x14
6x27x8
28x9x20
13x13x4
9x18x16
18x15x18
22x19x14
14x10x17
25x29x11
1x18x19
8x11x26
18x6x14
30x24x13
27x1x27
15x9x3
2x29x17
2x26x21
22x9x9
20x20x20
22x28x2
26x5x16
11x3x14
21x16x16
18x26x7
18x30x6
7x11x12
15x10x2
27x2x16
27x30x24
28x14x24
7x4x8
6x28x15
13x19x1
22x26x30
7x30x24
2x17x21
19x26x2
19x24x15
14x23x2
21x27x15
30x15x14
21x29x5
23x30x2
4x1x2
15x5x13
21x2x30
20x7x16
1x21x25
2x25x1
12x29x5
28x13x16
26x3x12
29x20x23
28x12x20
4x30x8
16x15x16
6x16x29
2x28x13
24x25x2
26x15x22
17x20x11
18x12x7
19x1x18
8x27x13
22x16x8
19x26x17
13x11x10
22x12x3
13x12x14
29x17x9
6x14x10
14x20x10
8x26x9
25x13x22
3x30x25
14x28x1
30x29x12
3x17x15
3x24x14
28x24x22
16x6x1
20x25x14
17x17x13
6x19x27
10x15x20
8x23x20
7x29x21
18x9x25
10x5x22
2x27x27
16x18x30
15x5x12
26x29x29
28x11x10
9x29x28
24x15x23
26x9x10
5x1x25
22x27x16
7x29x3
1x3x5
8x7x29
19x21x11
28x13x30
17x16x20
5x10x25
9x14x15
15x14x23
16x4x17
21x8x2
9x9x8
22x22x4
10x2x27
12x19x10
15x29x4
22x14x7
29x18x5
1x7x27
24x1x15
23x23x26
12x17x23
26x10x24
8x22x2
8x1x10
22x19x12
2x23x13
11x27x25
26x15x27
27x7x21
18x9x6
22x21x22
7x12x26
23x21x13
14x3x8
5x9x28
29x29x15
27x25x23
12x2x24
8x2x20
29x19x4
12x24x29
2x27x28
14x20x9
28x6x25
18x29x8
19x11x30
15x11x23
18x7x7
14x20x14
26x18x22
27x25x13
12x10x30
30x2x7
28x10x1
18x10x30
22x11x5
22x16x3
25x15x9
5x10x24
4x28x8
19x24x18
3x4x25
14x4x30
11x26x3
12x12x12
26x7x24
3x2x14
1x27x7
2x2x13
3x26x26
12x4x11
12x17x20
4x19x30
5x18x10
17x6x18
19x30x20
11x2x17
30x13x19
22x23x7
17x28x2
5x17x30
7x11x4
21x26x18
15x28x4
5x6x27
12x6x16
9x17x12
27x20x5
14x5x20
27x14x6
2x14x21
4x28x30
24x5x1
19x29x29
11x23x1
8x16x21
3x17x19
10x13x5
20x21x16
23x3x6
27x26x11
3x2x22
14x3x5
10x9x8"""

# ╔═╡ 0c139b96-8597-428f-b877-8ebe19f56b9b
presents = split(puzzle_input, "\n") |>
	Base.Fix1(map, line -> split(line, "x")) |>
	Base.Fix1(map, entry -> parse.(Int, entry))

# ╔═╡ fdebde6b-c7a2-46b7-936a-11c19fb037dc
# removes the entry with biggest values and returns the rest
shortest_side(dims) = dims[1:end .!= findfirst(x -> x == maximum(dims), dims)]

# ╔═╡ 5e2901d0-3731-441e-8f87-2db9b774760a
@assert prod(shortest_side([2,3,4])) == 6

# ╔═╡ 218a0a0f-4fa7-4b1b-8c3c-2ae96ef13c72
@assert prod(shortest_side([2,10,10])) == 20

# ╔═╡ a6a7817f-b41a-458e-a807-6cc44a82dc1f
begin
	ex = :(
		paper_required(l,w,h) = 2l*w + 2w*h + 2h*l
							+ prod(shortest_side((l,w,h))))
	eval(ex)
	latexify(ex)
end

# ╔═╡ ca06b273-ac98-4cde-94ec-6cd9a65e0807
@assert paper_required((2,3,4)...) == 52 + 6 == 58

# ╔═╡ 73275b5b-b4cd-47a7-ac0c-8b0ce25c6578
@assert paper_required((1,1,10)...) == 42 + 1 == 43

# ╔═╡ 77691a13-3967-4874-8d08-4f35ddb21369
part1 = sum(map(x->paper_required(x...), presents))

# ╔═╡ 4a9e893a-61d6-4c44-b032-50ce16aa245b
@test part1 == 1586300

# ╔═╡ b0d9c03f-6b5c-4d5f-a88d-bc1d641bd86a
md"Your puzzle answer was `1586300`."

# ╔═╡ c086f72d-ba34-42e5-bc38-4f3d40a6a0a2
md"""
# Part Two
The elves are also running low on ribbon. Ribbon is all the same width, so they only have to worry about the length they need to order, which they would again like to be exact.

The ribbon required to wrap a present is the shortest distance around its sides, or the smallest perimeter of any one face. Each present also requires a bow made out of ribbon as well; the feet of ribbon required for the perfect bow is equal to the cubic feet of volume of the present. Don't ask how they tie the bow, though; they'll never tell.

For example:

-   A present with dimensions `2x3x4` requires `2+2+3+3 = 10` feet of ribbon to wrap the present plus `2*3*4 = 24` feet of ribbon for the bow, for a total of `34` feet.
-   A present with dimensions `1x1x10` requires `1+1+1+1 = 4` feet of ribbon to wrap the present plus `1*1*10 = 10` feet of ribbon for the bow, for a total of `14` feet.

**How many total feet of ribbon should they order?**
"""

# ╔═╡ d23cc5bf-0748-489e-9da7-2054051482c9
begin
	ex2 = :(
		ribbon_required(l,w,h) = l*w*h + 2*sum(shortest_side((l,w,h)))
		)
	eval(ex2)
	latexify(ex2)
end

# ╔═╡ c1eef318-9168-4400-aec1-defd790d28db
@assert ribbon_required(2,3,4) == 10 + 24 == 34

# ╔═╡ 9bc722bb-f1ae-4cd9-9664-5ab7c282c57c
@assert ribbon_required(1,1,10) == 4 + 10 == 14

# ╔═╡ 630aad38-6f6f-4902-be97-8836fd7b8ed7
part2 = sum(map(x->ribbon_required(x...), presents))

# ╔═╡ 0bdf8ff7-fa34-429b-b7e4-899fd52c4f22
@test part2 == 3737498

# ╔═╡ 819beeb8-7b2a-4a88-aa50-032ef21dbae8
md"Your puzzle answer was `3737498`."

# ╔═╡ Cell order:
# ╠═b0602873-cda1-4d30-8ec8-1b9f7fdcdb07
# ╟─eac55bba-cc75-4632-832a-dbb06983be45
# ╟─0568b194-e70a-42b5-a9e8-89e710e61648
# ╠═0c139b96-8597-428f-b877-8ebe19f56b9b
# ╠═fdebde6b-c7a2-46b7-936a-11c19fb037dc
# ╠═5e2901d0-3731-441e-8f87-2db9b774760a
# ╠═218a0a0f-4fa7-4b1b-8c3c-2ae96ef13c72
# ╠═a6a7817f-b41a-458e-a807-6cc44a82dc1f
# ╠═ca06b273-ac98-4cde-94ec-6cd9a65e0807
# ╠═73275b5b-b4cd-47a7-ac0c-8b0ce25c6578
# ╠═77691a13-3967-4874-8d08-4f35ddb21369
# ╠═4a9e893a-61d6-4c44-b032-50ce16aa245b
# ╟─b0d9c03f-6b5c-4d5f-a88d-bc1d641bd86a
# ╟─c086f72d-ba34-42e5-bc38-4f3d40a6a0a2
# ╠═d23cc5bf-0748-489e-9da7-2054051482c9
# ╠═c1eef318-9168-4400-aec1-defd790d28db
# ╠═9bc722bb-f1ae-4cd9-9664-5ab7c282c57c
# ╠═630aad38-6f6f-4902-be97-8836fd7b8ed7
# ╠═0bdf8ff7-fa34-429b-b7e4-899fd52c4f22
# ╟─819beeb8-7b2a-4a88-aa50-032ef21dbae8


text:
using StatsBase

tests = ["weights",
         "moments",
         "scalarstats",
         "deviation",
         "cov",
         "counts",
         "ranking",
         "empirical",
         "hist",
         "rankcorr",
         "signalcorr",
         "misc",
         "sampling",
         "wsampling",
         "statmodels",
         "statquiz"]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
SIZE = 3
SEPARATOR = repeat("-",(4 * SIZE + 1))

function empty()
    # return an empty grid
    grid = Any[]

    return grid = [repeat(Any[" "],SIZE)  for i in range(1, stop=SIZE)] ##
end

function render(grid)

    #println("\n\n\n\t\t",SEPARATOR)

    for i in range(1,stop=SIZE)
        print("\n\t\t\t\t|")

        for j in range(1,stop=SIZE)

            if grid[i][j] == " "
                print(" ",(SIZE^2)-(i*SIZE-j)," |")

            elseif grid[i][j] == "X"
                print(" ")
                printstyled(grid[i][j],color=:green)
               print(" |")

            else
                print(" ")
                printstyled(grid[i][j],color=:red)
                print(" |")
            end

        end
println()
       # println("\n\t\t",SEPARATOR)
    end

end

function cellnum()
    return (SIZE^2)-(i*SIZE-j)
end

function win_player(grid, char)
    # check if a player wins the game
    # check rows

    for i = 1:3
        if all(grid[i][j] == char for j = 1:3)
            return true
        end
    end

    # check columns
    for j = 1:3
        if all(grid[i][j] == char for i = 1:3)
            return true
        end
    end

    # check diagonals
    if all(grid[i][i] == char for i = 1:3)
        return true
    end

    if all(grid[i][4 - i] == char for i = 1:3)
        return true
    end

    return false
end

function coordinates(choice)
    # return the row and coloumn corresponding to user input
    # get row
    row = ceil(Int,(3-((choice + 2) // 3)))
    #get col
    col = (choice -1) % 3

    return row+1,col+1

end

function terminal(grid)

    # check if the game is at a terminal state
    # player wins
    if win_player(grid, "X")
        return true
    end

    # computer wins
    if win_player(grid, "O")
        return true
    end

    # tie
    if all(grid[i][j] != " " for i = 1:3 for j = 1:3) ##
        return true
    end

    # otherwise the game isn't over yet
    return false

end

function utility(grid)
    # return the score corresponding to the terminal state
    if win_player(grid, "X")
        return -1
    end

    if win_player(grid, "O")
        return 1
        ###

    end
return 0
end

function actions(grid)
    # return possible actions a player can take at each state
    result = Any[]

    for i in range(1,stop =SIZE)
        for j in range(1,stop =SIZE)
            if grid[i][j] == " "
                push!(result, (i, j))

            end
        end

    end
    #shuffle!(result)
    return result

end

function minimax(grid, computer, alpha, beta, depth)

    # return the maximum value a player can obtain at each step
    if terminal(grid)
        return utility(grid), depth
    end

    if computer
        m = -Inf
        char = "O"
    else
        m = Inf
        char = "X"
    end

    for action in actions(grid)
        i , j = action
        grid[i][j] = char

        value , depth = minimax(grid, !computer, alpha, beta, depth + 1)

        if computer
         m = Int(max(m, value))
        else
         m = Int(min(m, value))
        end

        # undo the move
        grid[i][j] = " "

        # alpha-beta pruning
        if computer
            alpha = max(alpha, m)
        else
            beta = min(beta, m)
        end

        if beta <= alpha
            break
        end
    end
    return m, depth

end

function best_move(grid)
     # find all empty cells and compute the minimax for each one
    m = alpha = -Inf
    d = beta = Inf

result = 0;

    for action in actions(grid)
        i , j = action
        grid[i][j] = "O"
        value , depth = minimax(grid, false, alpha, beta, 0)

        if value > m || (value == m && depth < d)
            result = i ,j
            m = value
            d = depth
        end

        # undo the move
        grid[i][j] = " "

    end

    return result

end

function get_user_input(grid)
    # ask for player input : a number between 1 and 9

    while true
      #  render(grid)

        print("\n\nEnter choice (between 1-9): ")
        choice = readline()

        # validate input
        if !(choice in [string(i) for i in 1:9])##
            printstyled("\n\n\t\t\tinvalid choice",color=:red)
            sleep(1)
            continue
        end

         choice = parse(Int64,choice)
             i,j = coordinates(choice)

        if grid[i][j] != " "
            printstyled("\n\n\t\tThe chosen cell is already filled.",color=:red)
            sleep(1)
            render(grid)
            continue
        end

        # if cell is not full : fill with 'X'
        grid[i][j] = "X"
        break

    end

end

function game_loop(grid)

    #=bgn = true
    if rand(1:10) > 5=#

    while true

        printstyled("\n\n\t\t\t\tIt's  your turn\n\n",color=:green)

        # player turn
        get_user_input(grid)

        render(grid)

        # check if the player wins
        if win_player(grid, "X")
            print("\n\t\t\t\tYou win!\n",color=:cyan)
            break
        end

        # check if it's a tie
        if terminal(grid)
            println("\n\t\t\t\t    Tie!")
            break
        end

        # computer turn
        # find the best move to play
        i , j = best_move(grid)
        grid[i][j] = "O"

        printstyled("\n\n\t\t\t   It's the computer's turn\n\n",color=:yellow)
        sleep(2)
        # display the grid
        render(grid)

        # check if the computer wins with this choice
        if win_player(grid, "O")
            printstyled("\n\n\t\t\t\tYou lose!\n",color=:red)
            break
        end
    end

end

function play()
    while true
        grid = empty()
        render(grid)
        println("\n\n\t\twelcome to tic tac toe, good luck defeating the master")
        sleep(1)
        game_loop(grid)

         print("\n\n\t\t\tPlay again ? [y/n] ")
         again = readline()

        if uppercase(again) != "Y"
            break
        end
    end
end

play()

 function shuffle!(arr::AbstractArray)

    ind = rand(1:length(arr))
    checked = Any[]

      for i = 1:ind

         done = false
         chkd = true

         while !done

	        if length(checked) >= length(arr)
                return arr
            end

           tmp2 = rand(1:length(arr))
           tmp = rand(1:length(arr))

            while chkd
                   if !tmp in checked

                       if !tmp2 in checked

                           chkd = false

                       else
                           tmp2 = rand(1:length(arr))
                       end

                   else
                       tmp = rand(1:length(arr))
                   end
               end

                if tmp != tmp2

                   arr[tmp],arr[tmp2] = arr[tmp2], arr[tmp]
                   push!(checked,tmp,tmp2)
                   done = true

               end
            end
          end

     return arr
 end


text:
using PolyChaosODE, BenchmarkTools, DifferentialEquations, Distributions, PolyChaos, Plots

### This package ###

x0 = 2.0
μ, σ = -0.5, 0.05
tend = 10.0
L = 6

function exp_ode!(du, u, p, x, t)
    du[1] = (p[1] + p[2] * x[1]) * u[1]
end

prob_dim = 1
vars = (Normal(0.0, sqrt(2) / 2),)
num_polys = (L + 1,)
stoch_galerkin_ode = StochGalerkinODE(exp_ode!, prob_dim, vars; num_polys = num_polys)

x0_ = [x0]
tspan = (0, tend)
# sqrt(2) factor in the second term because we are using the Physicist's Hermite polynomials instead of the Probabilist's
p = [-0.5, 0.05 * sqrt(2)]
sol = stoch_galerkin_ode(x0_, tspan, p; alg = VCABM())

t_step = 0.01
interval_t = 0.0:t_step:tend
std_this =
    (
        x -> sqrt.(x)
    ).(
        compute_expectation_and_diag_variance(
            stoch_galerkin_ode,
            interval_t,
            [sol(t) for t in interval_t],
        )[2],
    )

@btime stoch_galerkin_ode($x0_, $tspan, $p; alg = VCABM())
@btime (
    x -> sqrt.(x)
).(
    compute_expectation_and_diag_variance(
        stoch_galerkin_ode,
        interval_t,
        [sol(t) for t in interval_t],
    )[2],
)

### PolyChaos (based on https://timueh.github.io/PolyChaos.jl/stable/random_ode/)###

Nrec = 40
opq = GaussOrthoPoly(L; Nrec = Nrec, addQuadrature = true)

a = [convert2affinePCE(μ, σ, opq); zeros(Float64, L - 1)]
xinit = [x0; zeros(Float64, L)]

const t2_ = Tensor(2, opq);
const t3_ = Tensor(3, opq);

function ODEgalerkin(du, u, p, t)
    du[:] = [
        sum(p[j+1] * u[k+1] * t3_.get([j, k, m]) / t2_.get([m, m]) for j = 0:L for k = 0:L) for m = 0:L
    ]
end

probgalerkin = ODEProblem(ODEgalerkin, xinit, tspan, a)
solgalerkin = solve(probgalerkin; alg = VCABM())

t, x = solgalerkin.t, solgalerkin.u;
mean_poly_chaos = [mean(x_, opq) for x_ in x]
std_poly_chaos = [std(x_, opq) for x_ in x]

@btime solve(probgalerkin; alg = VCABM())
@btime [mean(x_, opq) for x_ in x]
@btime [std(x_, opq) for x_ in x]

plot(t, mean_pce .- std_pce)
plot!(t, mean_pce)
plot!(t, mean_pce .+ std_pce)

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
using DeformableRegistration: Transformation, ImageProcessing, Interpolation, Distance
 using Base.Test
 using Interpolations: BSpline, Linear, Cubic, Line, Free
#  using Logging
#  Logging.configure(level=INFO)

@testset "Interpolation" begin
##
@testset "Interpolation: interpolation" begin
    ## load reference image

    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    img = loadImage(testimage)
    centeredGrid = getCellCenteredGrid(img)
    affineParameters = [0.5,0.4,50,0,0.5,100]
    deformationField = zeros(prod(size(img.data))*2)
    deformationField[prod(size(img.data))+1:end] = 10*sin.(0.01*centeredGrid.data[1:prod(size(img.data))])
    deformationField = scaledArray(deformationField, size(img.data), [1, 1.0], [0.0, 0])
    transformedGrid = transformGridAffine(centeredGrid,affineParameters) + deformationField
    transformedImageOnly,dx,dy = interpolateImage(img,transformedGrid,doDerivative=true)
    transformedImage,dY_transformedImage,dX_transformedImage = interpolateImage(img,transformedGrid,doDerivative=true, interpolationScheme=InterpLinearFast)
    transformedImageRef,dY_transformedImageRef,dX_transformedImageRef = interpolateImage(img,transformedGrid,doDerivative=true,interpolationScheme=BSpline(Linear()))

    ## using PyPlot
    #dY_transformedImage = reshape(dY_transformedImage, size(img.data))
    #dY_transformedImageRef = reshape(dY_transformedImageRef, size(img.data))
    #imshow(dY_transformedImage-dY_transformedImageRef)
    @test norm(transformedImageOnly - transformedImageRef) < 1.0
    @test transformedImageOnly[1:10] ≈ transformedImageRef[1:10] atol=1e-3

    @test norm(transformedImage - transformedImageRef) < 1.0
    @test transformedImage[1:10] ≈ transformedImageRef[1:10] atol=1e-3

    @test norm(dY_transformedImage - dY_transformedImageRef) ≈ 0 atol=0.5 #20
    @test dY_transformedImage[1:10] ≈ dY_transformedImageRef[1:10]

    @test norm(dX_transformedImage - dX_transformedImageRef) ≈ 0 atol=0.5 #20
    @test dX_transformedImage[1:10] ≈ dX_transformedImageRef[1:10]
end

@testset "Interpolation: interpolation with target grid" begin
    ## load reference image

    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    img = loadImage(testimage)
    img = restrictResolutionToLevel(img,1)
    centeredGrid = getCellCenteredGrid(img)
    affineParameters = [0.5,0.4,50,0,0.5,100]
    deformationField = zeros(prod(size(img.data))*2)
    deformationField[prod(size(img.data))+1:end] = 10*sin.(0.01*centeredGrid.data[1:prod(size(img.data))])
    deformationField = scaledArray(deformationField, size(img.data), img.voxelsize, img.shift)
    transformedGrid = transformGridAffine(centeredGrid,affineParameters) + deformationField

    targetGrid = getCellCenteredGrid(deformationField)

    transformedImage,dx,dy = interpolateImage(img,transformedGrid,doDerivative=true, interpolationScheme=BSpline(Linear()))
    tic()
    transformedImageTG,dxTG,dyTG = interpolateImage(img,transformedGrid,targetGrid, doDerivative=true, interpolationScheme=BSpline(Linear()))
    timing = toc()

    # @test timing < 0.3 not working on CI server :(

    @test norm(transformedImage - transformedImageTG) < 1e-10
    @test norm(dx - dxTG) < 1e-10
    @test norm(dy - dyTG) < 1e-10

    newSize = size(img.data)./2
    newSize = (Int(ceil(newSize[1])), Int(ceil(newSize[2])))
    targetGrid = getCellCenteredGrid(img.voxelsize*2, img.shift, newSize)
    @time transformedImageTG,dxTG,dyTG = interpolateImage(img,transformedGrid,targetGrid, doDerivative=true)

    transformedImageCoarse = transformedImage[1:2:end,1:2:end]
    @test norm(transformedImageCoarse[:] - transformedImageTG[:]) < 4
end

@testset "Interpolation: interpolation of part of image with target grid" begin
    ## load reference image

    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    img = loadImage(testimage)
    img = restrictResolutionToLevel(img,1)
    centeredGrid = getCellCenteredGrid(img.voxelsize*3, img.shift, (200,200))
    affineParameters = [0.5,0.4,25,0,0.5,50]
    transformedGrid = transformGridAffine(centeredGrid,affineParameters)

    # using PyPlot

    @time transformedImage,dx,dy = interpolateImage(img,transformedGrid,doDerivative=true, interpolationScheme=BSpline(Linear()))
    # subplot(221)
    # imshow(transformedImage)
    # subplot(222)
    # imshow(reshape(dx[1:200*200],transformedGrid.dimensions))

    targetGrid = getCellCenteredGrid(img.voxelsize*6, img.shift, (100,100))
    @time transformedImageTG,dxTG,dyTG = interpolateImage(img,transformedGrid,targetGrid, doDerivative=true)
    # subplot(223)
    # imshow(transformedImageTG)
    # subplot(224)
    # imshow(reshape(dxTG[1:100*100],targetGrid.dimensions))
    @test norm(transformedImage[1:2:200,1:2:200] - transformedImageTG) < 1

    targetGrid = getCellCenteredGrid(img.voxelsize*3, img.shift, (100,100))
    @time transformedImageTG,dxTG,dyTG = interpolateImage(img,transformedGrid,targetGrid, doDerivative=true)
    @test norm(transformedImage[1:100,1:100] - transformedImageTG) < 1

end

@testset "Interpolation: check derivatives with target grid 1" begin
    for interpolationScheme=[BSpline(Linear()), BSpline(Cubic(Line()))]

        include("helpers/checkDerivative.jl")
        testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
        img = loadImage(testimage)

        targetGrid = getCellCenteredGrid(img)
        # check image derivative
        transformedGrid = getCellCenteredGrid(img)
        transformedGrid.data[:] = transformedGrid.data[:] + 0.1
        Ifunc(x) = interpolateImage(img, x, targetGrid, interpolationScheme=interpolationScheme)[1][:]
        transformedImage, dX_transformedImage, dY_transformedImage =
            interpolateImage(img,transformedGrid,targetGrid,doDerivative=true,
                interpolationScheme=interpolationScheme)
        dTransformedImage = spdiagm((dX_transformedImage[:],dY_transformedImage[:]),[0, prod(size(img.data))])
        @time errlin,errquad = checkDerivative(Ifunc,dTransformedImage,transformedGrid)
        @test checkErrorDecay(errquad)
    end
end

##
@testset "Interpolation: check derivatives (all interpolation types)" for
    interpolationScheme=[InterpLinearFast, BSpline(Linear()), BSpline(Cubic(Line()))]
    include("helpers/checkDerivative.jl")
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    img = loadImage(testimage)

    # check image derivative
    centeredGrid = getCellCenteredGrid(img).data + 0.1
    Ifunc(x) = interpolateImage(img,x, interpolationScheme=interpolationScheme)[1][:]
    @time transformedImage, dX_transformedImage, dY_transformedImage =
        interpolateImage(img,centeredGrid,doDerivative=true,
            interpolationScheme=interpolationScheme)
    dTransformedImage = spdiagm((dX_transformedImage[:],dY_transformedImage[:]),[0, prod(size(img.data))])
    @time errlin,errquad = checkDerivative(Ifunc,dTransformedImage,centeredGrid)
    @test checkErrorDecay(errquad)
end

##
@testset "compare different interpolation types on real data" begin
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    refImg = loadImage(testimage)
    # define a cell centered grid, transform it and create a template image
    centeredGrid = getCellCenteredGrid(refImg)
    affineParametersInitial = [1.2,0.0,-50,0.0,1.2,-0]
    transformedGrid = transformGridAffine(centeredGrid,affineParametersInitial)
    temImg = interpolateImage(refImg,transformedGrid,interpolationScheme=InterpLinearFast)[1]
    temImg = createImage(temImg)
    ##
    targetGrid = getCellCenteredGrid(refImg)
    targetGrid.data[:] = targetGrid.data[:] - 0
    temImg2 = interpolateImage(refImg,transformedGrid,targetGrid,interpolationScheme=BSpline(Cubic(Free())))[1]
    temImg2 = createImage(temImg2)

    # if Logging.LogLevel == Logging.DEBUG
    # using PyPlot
    # figure(); PyPlot.imshow(Array(temImg2.data)-temImg.data)
    # end

    # border pixels differe in both interpolation schemes, most likely due
    # to InterpLinearFast padding the image too far...
    # allowing 2 pixels around the border of the image to be different
    @test norm(Array(temImg2.data)[:]-temImg.data[:])  < 32
    nDifferentPixels = sum(abs.(Array(temImg2.data)[:]-temImg.data[:]) .> 0.01)
    @test nDifferentPixels < maximum(size(refImg.data))*8
end

@testset "compare different interpolation types on real data" begin
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    refImg = loadImage(testimage)
    # define a cell centered grid, transform it and create a template image
    centeredGrid = getCellCenteredGrid(refImg)
    affineParametersInitial = [1.2,0.0,-50,0.0,1.2,-0]
    transformedGrid = transformGridAffine(centeredGrid,affineParametersInitial)
    temImg = interpolateImage(refImg,transformedGrid,interpolationScheme=InterpLinearFast)[1]
    temImg = createImage(temImg)

    temImg2 = interpolateImage(refImg,transformedGrid,interpolationScheme=BSpline(Linear()))[1]
    temImg2 = createImage(temImg2)

    # if Logging.LogLevel == Logging.DEBUG
    # #using PyPlot
    # #figure(); PyPlot.imshow(Array(temImg2.data)-temImg.data)
    # end

    # border pixels differe in both interpolation schemes, most likely due
    # to InterpLinearFast padding the image too far...
    # allowing 2 pixels around the border of the image to be different
    @test norm(Array(temImg2.data)[:]-temImg.data[:])  < 32
    nDifferentPixels = sum(abs.(Array(temImg2.data)[:]-temImg.data[:]) .> 0.01)
    @test nDifferentPixels < maximum(size(refImg.data))*4
end

##
@testset "interpolation with identity grid results in orgiginal data: $scheme" for scheme in (InterpLinearFast, BSpline(Linear()), BSpline(Cubic(Free())))
    data = [1.0 2; 3 4]
    testimage = createImage(data)
    g = getCellCenteredGrid(testimage)
    dataInterpolated = interpolateImage(testimage, g, interpolationScheme=scheme)[1]
    @test data ≈ dataInterpolated
end

@testset "smoothing" begin
    for k=[3,5,7,11]
        k_half = Int(ceil(k/2))
        for x = [ones(10,11)*2.0, rand(63,26)]
            x[1:k_half,:] = 0
            x[:,1:k_half] = 0
            x[end-k_half:end,:] = 0
            x[:,end-k_half:end] = 0

            #x = ones(10,11)*2.0;
            xnew = smoothArray(x,k, Float64(k))
            @test sum(x[:]) ≈ sum(xnew[:])
        end
    end

end

end


text:
using Test
import LessOLS: Process, make_sample, sampler_normal, ols, desc, show, yhat, equation

# # use profiler
# # https://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/

gen = sampler_normal(a=0, b=20, β_0=5, β=[1, 2.5], sd_e=0.5)
sam = gen(100)
lm = ols(sam, intercept=true)
show(lm)

@assert equation([-1/3, -2/7, 1/29],true,4) == "Y = -0.3333 - 0.2857*X1 + 0.0345*X2"


text:
module Slack

using HTTP
using CodecZlib
using JSON

include("send.jl")
export
    sendtoslack,
    sendattachmenttoslack

include("channel.jl")
export
    getchannels

# include("utils.jl")
# export
#     convert_HTTP_Response_To_JSON

end # module


text:
## CLOSURE RELATIONS
#============================================#

# Momentum shape for density: H
H_Whitfield(H, M_e) = (H - 0.290M_e^2)/(1 + 0.113M_e^2)

## Laminar closure for walls - Falkner-Skan profiles
#============================================#

# Laminar kin-ergy shape: H*
lam_H_star(H_k) = 1.515 + ifelse(H_k < 4., 0.076, 0.040) * (H_k - 4)^2 / H_k

# Laminar skin-friction: Re_θ CF/2
lam_CF(H_k) = -0.067 + ifelse(H_k < 7.4, 0.01977(7.4 - H_k)^2 / (H_k - 1), 0.022(1. - 1.4 / (H_k - 6.))^2)

# Laminar dissipation: Re_θ 2CD / H_k*
lam_CD(H_k) = 0.207 + ifelse(H_k < 4., 0.00205(4 - H_k)^5.5, -0.003(H_k - 4.)^2 / (1. + 0.02(H_k - 4.)^2))

## Laminar closure for wakes
#============================================#

# Laminar wake kin-ergy shape: H*
lam_wake_H_star(H_k) = 1.50 + ifelse(H_k < 3.5, 0.025(3.5 - H_k)^3 + 0.001(3.5 - H_k)^5, 0.015(H_k - 3.5)^2 / H_k)

# Laminar wake dissipation: Re_θ 2CD / H*
lam_wake_CD(H_k) = 1.52(H_k - 1)^2 / (3 + H_k^3)

## Turbulent closure relations
#============================================#

# Skin-friction
turb_CF(H_k, Re_θ, M_e) = let F_c = √(1 + 0.2M_e^2); 0.3exp(-1.33H_k) * log(10, Re_θ / F_c)^(-1.74 - 0.31H_k) / F_c + 1.1e-4(tanh(4 - H_k/0.875) - 1) end

# Density thickness (negligible at low speeds): H**
turb_H_star_star(H_k, M_e) = (0.064 / (H_k - 0.8) + 0.251)M_e^2

# Y+
y_plus(u_t, μ_e, ρ_e, η) = ρ_e * u_t * η / μ_e

# Velocity profile
function turb_ue(Cf, θ, y_p, η, a, b)
    uτ_ue = √(Cf / 2)
    s     = Cf / abs(Cf)

    uτ_ue * s / 0.09 * atan(0.09y_p) + (1 - uτ_ue * s * π / 0.18) * √(tanh(a * (η, θ)^b))
end

# H*?
function turb_H_star(H_k, Re_θ)
    H_0 = ifelse(Re_θ < 400, 4, 3 + 400/Re_θ)

    1.505 + 4/Re_θ + ifelse(H_k < H_0, (0.165 - 1.6/√Re_θ) * (H_0 - H_k)^1.6 / H_k, (H_k - H_0)^2 * (0.04 / H_k + 0.007log(10, Re_θ)/(H_k - H_0 + 4/log(10, Re_θ))^2))
end

# CD
turb_CD(cf, u_s, c_τ) = cf * u_s / 2 + c_τ * (1 - u_s)

# U_s
wall_slip(H_star, H_k, H) = H_star / 2 * (1 - 4/3(H_k - 1) / H)

## Boundary layer thickness δ from Green's correlation
δ_thick(H_k, θ, δ_star) = θ * (3.15 + 1.72 / (H_k - 1)) + δ_star

# Equilibirum shear stress coefficient C_τ_EQ
c_τ_equilibrium(H, H_k, H_star, u_s) = 0.015H_star / (1 - u_s) * (H_k - 1)^3 / (H_k^2 * H)

## Transition
#============================================#

# Rate slope dñ/dReθ
dñdReθ(H_k) = 0.01 * √((2.4H_k - 3.7 + 2.5tanh(1.5H_k - 4.65))^2 + 0.25)

# Critical Reynolds number: log₁₀ Re_θ
log10Reθ(H_k) = (1.415/(H_k - 1) - 0.489)tanh(20/(H_k - 1) - 12.9) + 3.295/(H_k - 1) + 0.44

# G-β locus?
G(β) = 6.7 * √(1 + 0.75β)

# Weird amplification ratio DE
l(H_k) = (6.54H_k - 14.07) / H_k^2
m(H_k) = (0.058(H_k - 4)^2 / (H_k - 1) - 0.068) / l(H_k)
# New correlation: -0.05 + 2.7*H_k - 5.5H_k^2 + 3.0H_k^3

## Generic conversions
#============================================#

# Cf = 2lam_CF(H) / Re
cf(U_e, θ, H, ν = 1.5e-5) = 2 * lam_CF(H) / reynolds_number(U_e, θ, ν)

# CD = lam_H_star(H) * lam_CD(H) / 2Re
cD(U_e, θ, H, ν) = lam_H_star(H) * lam_CD(H) / (2 * reynolds_number(U_e, θ, ν) )

# Squire-Young formula for drag coefficient
cd_squire_young(θ_wake, u_e_wake, H_wake, V_inf) = 2θ_wake * (u_e_wake / V_inf)^((5 + H) / 2)


text:
# structs and methods for basic neural accumulator / neural arithmetic
# logic unit.

# - Neural Accumulator - #
"""
A linear layer whose weights are soft constrained to be near one of
{-1, 0, 1}. The weights are calculated by tanh.(W) .* σ.(M).
"""
struct NeuralAccumulator{R <: AbstractMatrix}
    W::R
    M::R
    function NeuralAccumulator{R}(w, m) where {R <: AbstractMatrix}
        size(w) == size(m) || error("sizes of weight matrices must match")
        new{R}(w, m)
    end
end
NeuralAccumulator(W::R, M::R) where {R <: AbstractMatrix} = NeuralAccumulator{R}(W, M)

function NeuralAccumulator(::Type{T}, in::Integer, out::Integer, init_fn) where {T <: Number}
    return NeuralAccumulator(init_fn(T, out, in), init_fn(T, out, in))
end

NeuralAccumulator(in::Integer, out::Integer, init_fn=randn) = NeuralAccumulator(Float64, in, out, init_fn)

operator(nac::NeuralAccumulator) = (2 .* σ.(nac.W) .- 1) .* σ.(nac.M)

(nac::NeuralAccumulator)(x) = operator(nac) * x

# - Neural Arithmetic Logic Unit - #
struct NeuralALU{R <: AbstractMatrix}
    nac::NeuralAccumulator{R}
    G::R
    b
end

NALU(in::Integer, out::Integer; initW = Flux.initn, initb = zeros) =
    NALU(NAC(in, out, initW=initW), param(initW(out, in)), param(initb(out)))

Flux.@treelike(NALU)

function (nalu::NALU)(x)
    # gate
    g = σ_stable.(nalu.G*x .+ nalu.b)

    # addition
    a = nalu.nac(x)

    # multiplication
    m = exp.(nalu.nac(log.(abs.(x) .+ eps())))

    # nalu
    return g .* a + (1 .- g) .* m
end

max_stars_repo_name: dillondaudert/neural-alu

max_issues_repo_name: dillondaudert/neural-alu

max_forks_repo_name: dillondaudert/neural-alu


text:
# Parallel network simulation

using DifferentialEquations
using DynamicalNetworks
using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using Distributed
using Logging
using LoggingExtras
using JLD2
using Dates

# Includes
include("getclargs.jl")
include("setlogger.jl")

# Settings
clargs = getclargs()
mode = clargs["mode"]
numexp = clargs["ntrials"]
minsnr = clargs["minsnr"]
stepsnr = clargs["stepsnr"]
maxsnr = clargs["maxsnr"]
nbits = clargs["nbits"]
tbit = clargs["tbit"]
ncores = clargs["ncores"]
high = clargs["high"]
low = clargs["low"]
ϵ = clargs["coupling-strength"]
srcnode = clargs["srcnode"]
dstnode = clargs["dstnode"]
simdir = clargs["simdir"]
simprefix = clargs["simprefix"]
loglevel = clargs["loglevel"]
savenoise = clargs["savenoise"]
P = diagm([1.0, 0.0, 0.0])
H = [
    1 0 0 0
    -1 0 0 0
    0 1 0 0
    0 -1 0 0
    0 0 1 0
    0 0 -1 0
    0 0 0 1
    0 0 0 -1
]

# Construct simulation path
simpath = joinpath(simdir, simprefix * string(Dates.format(now(), "yy-mm-dd-HH-MM-SS")))
isdir(simpath) || mkpath(simpath)

# Set the logger
logger = setlogger(simpath, loglevel)

# Start simulation
tinit = time()
@info "Started simulation with settings" clargs

# Initializations
@info "Initializing simulation..."
ti = 0.0
tf = nbits * tbit
params = (
    mode = mode,
    nbits = nbits,
    tbit = tbit,
    bits = rand(Bool, nbits),
    high = high,
    low = low,
    ϵ = ϵ,
    tswitch = ti:tbit:tf-tbit,
    srcnode = srcnode,
    dstnode = dstnode,
    ctrlidx = findall(!iszero, H),
    snrrange = minsnr:stepsnr:maxsnr,
    noiserange = vcat(fill((minsnr:stepsnr:maxsnr)', numexp)...),
    ncores = ncores,
    simdir = simdir,
    simprefix = simprefix
)
@info "Done."

# Network
@info "Constructing network..."
n = 8
graph = cycle_graph(n) |> SimpleWeightedGraph
add_edge!(graph, 2, 7)
add_edge!(graph, 3, 6)
net = SDENetwork(Lorenz, ϵ, graph, H, P)
@info "Done."

# Loads
if params.mode == "distributed"
    @info "Loading processess..."
    using Distributed
    na = length(Sys.cpu_info()) - 1 - nprocs()
    params.ncores ≤ na ? addprocs(params.ncores) : addprocs(na)
    @everywhere begin
        using Pkg
        dev_env_path = joinpath(Pkg.envdir(), "dev-env")
        dirname(Pkg.project().path) == dev_env_path || Pkg.activate(dev_env_path)
        include("load.jl")
    end
    @info "Done."
else
    include("load.jl")
end

# Get base problem
@info "Constructing base problem..."
prob = netproblem(net, (ti, tf), params, noise_rate_prototype = get_noise_rate_prototype(net),
    callback = DiscreteCallback(condition, action!), tstops = params.tswitch)
@info "Done."

# Get ensemble problem
@info "Getting ensemble problem..."
ensembleprob = EnsembleProblem(prob, prob_func = prob_func)
@info "Done."

# Simulate ensemle problem
@info "Solving ensemble problems..."
if params.mode == "distributed"
    ensemblealg = EnsembleDistributed()
elseif params.mode == "threads"
    ensemblealg = EnsembleThreads()
else
    ensemblealg = EnsembleSerial()
end
sim = solve(ensembleprob, SRA(), ensemblealg, trajectories = length(params.noiserange))
@info "Done."

# Save simulation
@info "Writing simulation..."
cartidx = CartesianIndices(params.noiserange)
for (i, sol) in enumerate(sim)
    trial, snr = cartidx[i].I
    snrpath = joinpath(simpath, "$snr-dB")
    isdir(snrpath) || mkpath(snrpath)
    trialpath = joinpath(snrpath, "Traial-$trial")
    isdir(trialpath) || mkpath(trialpath)
    @save joinpath(trialpath, "sol.jld2") sol
end
@save joinpath(simpath, "bits.jld2") params.bits
@info "Done."